text:
################################################################################
#                                BASIC HELPERS
################################################################################
# Macro error function
# inspired from https://github.com/jump-dev/JuMP.jl/blob/709d41b78e56efb4f2c54414266b30932010bd5a/src/macros.jl#L923-L928
function _macro_error(macroname, args, source, str...)
    error("At $(source.file):$(source.line): `@$macroname($(join(args, ", ")))`: ",
          str...)
end

# Escape when needed
# taken from https://github.com/jump-dev/JuMP.jl/blob/709d41b78e56efb4f2c54414266b30932010bd5a/src/macros.jl#L895-L897
_esc_non_constant(x::Number) = x
_esc_non_constant(x::Expr) = isexpr(x,:quote) ? x : esc(x)
_esc_non_constant(x) = esc(x)

# Extract the name from a macro expression
# Inspired from https://github.com/jump-dev/JuMP.jl/blob/45ce630b51fb1d72f1ff8fed35a887d84ef3edf7/src/Containers/macro.jl#L8-L17
_get_name(c::Symbol) = c
_get_name(c::Nothing) = ()
_get_name(c::AbstractString) = c
function _get_name(c::Expr)
    if isexpr(c, :string)
        return c
    else
        return c.args[1]
    end
end

# Given a base_name and idxvars, returns an expression that constructs the name
# of the object.
# Inspired from https://github.com/jump-dev/JuMP.jl/blob/709d41b78e56efb4f2c54414266b30932010bd5a/src/macros.jl#L930-L946
function _name_call(base_name, idxvars)
    if isempty(idxvars) || base_name == ""
        return base_name
    end
    ex = Expr(:call, :string, base_name, "[")
    for i in eachindex(idxvars)
        # Converting the arguments to strings before concatenating is faster:
        # https://github.com/JuliaLang/julia/issues/29550.
        esc_idxvar = esc(idxvars[i])
        push!(ex.args, :(string($esc_idxvar)))
        i < length(idxvars) && push!(ex.args, ",")
    end
    push!(ex.args, "]")
    return ex
end

# Process macro arugments
function _extract_kwargs(args)
    arg_list = collect(args)
    if !isempty(args) && isexpr(args[1], :parameters)
        p = popfirst!(arg_list)
        append!(arg_list, p.args)
    end
    extra_kwargs = filter(x -> isexpr(x, :(=)) && x.args[1] != :container &&
                          x.args[1] != :base_name, arg_list)
    container_type = :Auto
    base_name = nothing
    for kw in arg_list
        if isexpr(kw, :(=)) && kw.args[1] == :container
            container_type = kw.args[2]
        elseif isexpr(kw, :(=)) && kw.args[1] == :base_name
            base_name = esc(kw.args[2])
        end
    end
    pos_args = filter!(x -> !isexpr(x, :(=)), arg_list)
    return pos_args, extra_kwargs, container_type, base_name
end

# Add on keyword arguments to a function call expression and escape the expressions
# Adapted from https://github.com/jump-dev/JuMP.jl/blob/d9cd5fb16c2d0a7e1c06aa9941923492fc9a28b5/src/macros.jl#L11-L36
function _add_kwargs(call, kwargs)::Nothing
    for kw in kwargs
        push!(call.args, esc(Expr(:kw, kw.args...)))
    end
    return
end

# Determine if an expression contains any index variable symbols
function _has_idxvars(expr, idxvars)::Bool
    expr in idxvars && return true
    if expr isa Expr
        return any(_has_idxvars(a, idxvars) for a in expr.args)
    end
    return false
end

# Ensure a model argument is valid
# Inspired from https://github.com/jump-dev/JuMP.jl/blob/d9cd5fb16c2d0a7e1c06aa9941923492fc9a28b5/src/macros.jl#L38-L44
_valid_model(_error::Function, model::InfiniteModel, name) = nothing
function _valid_model(_error::Function, model, name)
    _error("Expected $name to be an `InfiniteModel`, but it has type ",
           typeof(model))
end

# Check if a macro julia variable can be registered
# Adapted from https://github.com/jump-dev/JuMP.jl/blob/d9cd5fb16c2d0a7e1c06aa9941923492fc9a28b5/src/macros.jl#L66-L86
function _error_if_cannot_register(
    _error::Function,
    model::InfiniteModel,
    name::Symbol
    )::Nothing
    if haskey(JuMP.object_dictionary(model), name)
        _error("An object of name $name is already attached to this model. If ",
               "this is intended, consider using the anonymous construction ",
               "syntax, e.g., `x = @macro_name(model, ...)` where the name ",
               "of the object does not appear inside the macro. Alternatively, ",
               "use `unregister(model, :$(name))` to first unregister the ",
               "existing name from the model. Note that this will not delete ",
               "the object; it will just remove the reference at ",
               "`model[:$(name)]`")
    end
    return
end
function _error_if_cannot_register(_error::Function, model::InfiniteModel, name)
    return _error("Invalid name $name.")
end

# Update the creation code to register and assign the object to the name
# Inspired from https://github.com/jump-dev/JuMP.jl/blob/d9cd5fb16c2d0a7e1c06aa9941923492fc9a28b5/src/macros.jl#L88-L120
function _macro_assign_and_return(_error, code, name, model)
    return quote
        _error_if_cannot_register($_error, $model, $(quot(name)))
        $(esc(name)) = $code
        $model[$(quot(name))] = $(esc(name))
    end
end

# Wrap the macro generated code for better stacttraces (assumes model is escaped)
# Inspired from https://github.com/jump-dev/JuMP.jl/blob/d9cd5fb16c2d0a7e1c06aa9941923492fc9a28b5/src/macros.jl#L46-L64
function _finalize_macro(_error, model, code, source::LineNumberNode)
    return Expr(:block, source,
                :(_valid_model($_error, $model, $(quot(model.args[1])))), code)
end

################################################################################
#                          INFINITE PARAMETER MACRO
################################################################################
## Process a distribution input from a macro into a distribution domain
# Univariate Distribution
function _distribution_or_error(
    _error::Function,
    dist::D
    )::UniDistributionDomain{D} where {D <: Distributions.UnivariateDistribution}
    return UniDistributionDomain(dist)
end

# Multivariate Distribution
function _distribution_or_error(
    _error::Function,
    dist::D
    )::MultiDistributionDomain{D} where {D <: NonUnivariateDistribution}
    return MultiDistributionDomain(dist)
end

# Fallback
function _distribution_or_error(_error::Function, dist)
    _error("Expected distribution from `Distributions.jl`, but got input of ",
           "type `$(typeof(dist))`.")
end

## Process an infinite domain macro input and make sure it is correct
# AbstractInfiniteDomain
function _domain_or_error(
    _error::Function,
    domain::D
    )::D where {D <: AbstractInfiniteDomain}
    return domain
end

# Vector{Real}
function _domain_or_error(
    _error::Function,
    vect::Vector{<:Real}
    )::IntervalDomain
    if length(vect) == 2
        return IntervalDomain(vect...)
    else
        _error("Expected interval domain of format `[lb, ub]`, but got `$vect`.")
    end
end

# Distribution
function _domain_or_error(_error::Function, dist::Distributions.Distribution)
    _error("Distribution `$dist` was specified as a domain, but should " *
          "be given as a distribution via the syntax: `@infinite_parameter( " *
          "model, param_expr ~ distribution, kwargs...)`.")
end

# Fallback
function _domain_or_error(_error::Function, domain)
    _error("Expected an infinite domain, but got input of ",
           "type `$(typeof(domain))`.")
end

# Process an expr as a distribution
function _make_distribution_call(_error, expr)
    return :( _distribution_or_error($_error, $(_esc_non_constant(expr))) )
end

# Process an expr as an infinite domain
function _make_domain_call(_error, expr)
    return :( _domain_or_error($_error, $(_esc_non_constant(expr))) )
end

# Parse a distribution input
function _parse_parameter(_error::Function, ::Val{:~}, param, dist)
    return param, _make_distribution_call(_error, dist)
end

# Parse a domain input
function _parse_parameter(
    _error::Function,
    ::Union{Val{:in}, Val{:∈}},
    param,
    domain
    )
    return param, _make_domain_call(_error, domain)
end

# Fallback
function _parse_parameter(_error::Function, ::Val{S}, param, rhs) where S
    _error("Unexpected operator $S.")
end

# Define the keyword arg types that can vary for individual dependent parameters
const _ArrayKwargNames = (:supports, :derivative_method)

"""
    @infinite_parameter(model::InfiniteModel, kwargs...)

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

```julia
@infinite_parameter(model::InfiniteModel, expr, kwargs...)
```

Add an infinite parameter to the model `model` described by the expression `expr`,
and the keyword arguments `kw_args`. (Note that in the following the symbol `in`
can be used instead of `∈`) The expression `expr` can be of the form:
- `paramexpr` creating parameters described by `paramexpr`.
- `paramexpr ∈ [lb, ub]` creating parameters described by `paramexpr` characterized
   by a continuous interval domain with lower bound `lb` and upper bound `ub`.
- `paramexpr ~ dist` creating parameters described by `paramexpr` characterized
   by the `Distributions.jl` distribution object `dist`.
- `paramexpr ∈ domain` creating parameters described by `paramexpr` characterized
   by the `AbstractInfiniteDomain` object `domain`.

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 `kwargs` 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`.
- `domain`: The `InfiniteDomain` characterizing the parameters see subtypes of
            [`AbstractInfiniteDomain`](@ref).
- `distribution`: Sets the `Distributions.jl` distribution object that
   characterizes the parameters (specified instead of a domain).
- `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 0.
- `derivative_method`: Specify the numerical method to evaluate derivatives that
                       are taken with respect to the parameter.
- `sig_digits`: Specifies the number of significant digits that should be used
                in automatic support generation. Defaults to `DefaultSigDigits`.
- `independent`: Specifies if the each parameter is independent from each other
  or not. Defaults to false.
- `container`: Specify the container type. Defaults to `Auto`

**Examples**
```julia-repl
julia> @infinite_parameter(m, x in [0, 1])
x

julia> @infinite_parameter(m, y[i = 1:2] ~ MvNormal(ones(2)), num_supports = 10)
2-element Array{GeneralVariableRef,1}:
 y[1]
 y[2]

julia> z = @infinite_parameter(m, [["a", "b"]], distribution = Normal(),
                               independent = true)
1-dimensional DenseAxisArray{GeneralVariableRef,1,...} with index sets:
    Dimension 1, ["a", "b"]
And data, a 2-element Array{GeneralVariableRef,1}:
 noname[a]
 noname[b]
```
"""
macro infinite_parameter(model, args...)
    # define error message function
    _error(str...) = _macro_error(:infinite_parameter, (model, args...),
                                  __source__, str...)
    # parse the arguments
    pos_args, extra_kwargs, container_type, base_name = _extract_kwargs(args)

    # process the keyword arguments
    domain_kwarg = filter(kw -> kw.args[1] == :domain, extra_kwargs)
    dist_kwarg = filter(kw -> kw.args[1] == :distribution, extra_kwargs)
    indep_kwarg = filter(kw -> kw.args[1] == :independent, extra_kwargs)
    filter!(kw -> !(kw.args[1] in (:domain, :independent, :distribution)),
            extra_kwargs)

    # process the positional arguments
    if isempty(pos_args)
        p = gensym()
        is_anon_single = true
    elseif length(pos_args) == 1
        p = first(pos_args)
        is_anon_single = false
    else
        _error("Too many positional arguments. See docs for correct syntax.")
    end

    # extract the domain if given symbolically
    # There are 9 cases to consider:
    # p                                  | type of p | p.head
    # -----------------------------------+-----------+------------
    # name                               | Symbol    | NA
    # name[1:2]                          | Expr      | :ref
    # name[i = 1:2, j = 1:2; i + j >= 3] | Expr      | :typed_vcat
    # [1:2]                              | Expr      | :vect
    # [i = 1:2, j = 1:2; i + j >= 3]     | Expr      | :vcat
    # name in domain                     | Expr      | :call
    # name[...] in domain                | Expr      | :call
    # name ~ dist                        | Expr      | :call
    # name[...] ~ dist                   | Expr      | :call
    # In the 4 last cases, we call _parse_parameter
    if isexpr(p, :call)
        param, domain = _parse_parameter(_error, Val(p.args[1]), p.args[2:end]...)
        if !isempty(domain_kwarg) || !isempty(dist_kwarg)
            _error("Cannot double specify the infinite domain and/or ",
                   "distribution.")
        end
    elseif !isempty(domain_kwarg) && isempty(dist_kwarg)
        param = p
        domain = _make_domain_call(_error, domain_kwarg[1].args[2])
    elseif isempty(domain_kwarg) && !isempty(dist_kwarg)
        param = p
        domain = _make_distribution_call(_error, dist_kwarg[1].args[2])
    elseif !isempty(domain_kwarg) && !isempty(dist_kwarg)
        _error("Cannot specify both a domain and a distribution.")
    else
        _error("Must specify the infinite domain and/or the given syntax is ",
               "unrecognized. See docs for accepted forms.")
    end

    # determine if it is a single parameter
    is_single = param isa Symbol

    # determine if independent
    is_independent = false
    if is_single
        is_independent = true
    elseif !isempty(indep_kwarg)
        if !isa(indep_kwarg[1].args[2], Bool)
            _error("Can only specify boolean literals (`false` or `true`) ",
                   "with keyword `independent`.")
        end
        is_independent = indep_kwarg[1].args[2]
    end

    # determine if anonymous
    is_anon = isexpr(param, (:vect, :vcat)) || is_anon_single

    # process the naming
    name = _get_name(param)
    if !isa(name, Symbol) && !is_anon
        _error("Expression $name should not be used as a parameter name. Use ",
               "the \"anonymous\" syntax ",
               "`$name = @infinite_parameter(model, kwargs...)`` instead.")
    end
    if isnothing(base_name)
        base_name = is_anon ? "" : string(name)
    end

    # prepare the name code if it is multi-dimensional
    if !is_single
        idxvars, inds = JuMPC.build_ref_sets(_error, param)
        name_code = _name_call(base_name, idxvars)
        if model in idxvars
            _error("Index $(model) is the same symbol as the model. Use a ",
                   "different name for the index.")
        end
    end

    # make the build call
    if is_independent
        # we only have a single parameter or an array of independent parameters
        build_call = :( build_parameter($_error, $domain) )
    else
        # we have multi-dimensional dependent parameters
        # let's first process the kwargs accordingly (array ones are made into
        # containers)
        array_kwargs = filter(kw -> kw.args[1] in _ArrayKwargNames, extra_kwargs)
        filter!(kw -> !(kw.args[1] in _ArrayKwargNames), extra_kwargs)
        # check that the non array keywords do not contain index variables used
        # for building containers
        for kw in extra_kwargs
            if _has_idxvars(kw.args[2], idxvars)
                _error("Cannot use index variables with `$(kw.args[1])` keyword.")
            end
        end
        # now let's make the build call
        inds_var = gensym() # used as a placeholder variable for the ContainerIndices used to vectorize
        domain_var = gensym() # used as a placeholder for the domain container
        vect_domain_call = :( Collections.vectorize($domain_var, $inds_var) )
        build_call = :( _build_parameters($_error, $vect_domain_call, $inds_var) )
        for kw in array_kwargs
            code = JuMPC.container_code(idxvars, inds,
                                        _esc_non_constant(kw.args[2]),
                                        container_type)
            vect_code = :( Collections.vectorize($code, $inds_var) )
            push!(build_call.args, Expr(:kw, kw.args[1], vect_code))
        end
    end
    _add_kwargs(build_call, extra_kwargs)

    # make the creation code
    esc_model = esc(model)
    if is_single
        # only single independent parameter
        creation_code = :( add_parameter($esc_model, $build_call, $base_name) )
    elseif is_independent
        # multi-dimensional independent parameters
        parameter_call = :( add_parameter($esc_model, $build_call, $name_code) )
        creation_code = JuMPC.container_code(idxvars, inds, parameter_call,
                                             container_type)
    else
        # multi-dimensional dependent parameters
        domain_call = JuMPC.container_code(idxvars, inds, domain, container_type)
        name_code = JuMPC.container_code(idxvars, inds, name_code, container_type)
        name_code = :( Collections.vectorize($name_code, $inds_var) )
        vect_add_call = :( add_parameters($esc_model, $build_call, $name_code) )
        creation_code = quote
            $domain_var = $domain_call
            $inds_var = Collections.indices($domain_var)
            Collections.unvectorize($vect_add_call, $inds_var)
        end
    end

    # finalize the code and register the parameter name if needed
    if is_anon
        macro_code = creation_code
    else
        macro_code = _macro_assign_and_return(_error, creation_code, name,
                                              esc_model)
    end
    return _finalize_macro(_error, esc_model, macro_code, __source__)
end

################################################################################
#                           FINITE PARAMETER MACRO
################################################################################

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

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, kwargs...)
```
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 `kwargs` 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, defaults to `Auto`.

**Examples**
```julia-repl
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...)
    # process the inputs
    esc_model = esc(model)
    pos_args, kwargs, container_type, base_name = _extract_kwargs(args)

    # make an error function
    _error(str...) = _macro_error(:finite_parameter, (model, args...),
                                  __source__, str...)

    # process the positional arguments
    if length(pos_args) == 1
        expr = popfirst!(pos_args)
    else
        _error("Incorrect number of arguments. Must be of form ",
               "@finite_parameter(model, name_expr == value).")
    end
    if isexpr(expr, :call)
        if expr.args[1] !== :(==)
            _error("Unrecognized operator. Must be of form ",
                   "@finite_parameter(model, name_expr == value).")
        end
        param = expr.args[2]
        value = _esc_non_constant(expr.args[3])
        is_anon = isexpr(param, (:vect, :vcat))
    else
        param = gensym()
        value = _esc_non_constant(expr)
        is_anon = true
    end

    # process the name
    name = _get_name(param)
    if isnothing(base_name)
        base_name = is_anon ? "" : string(name)
    end
    if !isa(name, Symbol) && !is_anon
        _error("Expression $name should not be used as a parameter name. Use " *
               "the \"anonymous\" syntax $name = @finite_parameter(model, " *
               "...) instead.")
    end

    # make the build call
    build_call = :( build_parameter($_error, $value) )
    _add_kwargs(build_call, kwargs)

    # make the creation code
    if isa(param, Symbol)
        # easy case with single parameter
        creation_code = :( add_parameter($esc_model, $build_call, $base_name) )
    else
        # we have a container of parameters
        idxvars, inds = JuMPC.build_ref_sets(_error, param)
        if model in idxvars
            _error("Index $(model) is the same symbol as the model. Use a ",
                   "different name for the index.")
        end
        name_code = _name_call(base_name, idxvars)
        parameter_call = :( add_parameter($esc_model, $build_call, $name_code) )
        creation_code = JuMPC.container_code(idxvars, inds, parameter_call,
                                             container_type)
    end

    # finalize the macro
    if is_anon
        macro_code = creation_code
    else
        macro_code = _macro_assign_and_return(_error, creation_code, name,
                                              esc_model)
    end
    return _finalize_macro(_error, esc_model, macro_code, __source__)
end

################################################################################
#                         PARAMETER FUNCTION MACRO
################################################################################
const _BadOperators = (:in, :<=, :>=, :(==), :≥, :≤, ∈)

## Make methods to find and replace a part of an expression
# Expr
function _expr_replace!(ex::Expr, old, new)
    for (i, a) in enumerate(ex.args)
        if a == old
            ex.args[i] = new
        elseif a isa Expr
            _expr_replace!(a, old, new)
        end
    end
    return ex
end

# Symbol
function _expr_replace!(ex::Symbol, old, new)
    return ex == old ? new : ex
end

## Safely extract the parameter references (copy as needed)
# Symbol
function _extract_parameters(ex::Symbol)
    return esc(ex)
end

# Expr
function _extract_parameters(ex::Expr)
    return esc(copy(ex))
end

# Helper method to process parameter function expressions
function _process_func_expr(_error::Function, raw_expr)
    if isexpr(raw_expr, :call)
        # check that the call is not some operator
        if raw_expr.args[1] in _BadOperators
            _error("Invalid input syntax.")
        end
        func_expr = esc(raw_expr.args[1])
        is_anon = false
        # check for keywords
        if isexpr(raw_expr.args[2], :parameters) ||
            any(isexpr(a, :kw) for a in raw_expr.args[2:end])
            _error("Cannot specify keyword arguements directly, try using an ",
                   "anonymous function.")
        end
        # extract the parameter inputs
        pref_expr = esc(Expr(:tuple, raw_expr.args[2:end]...))
    elseif isexpr(raw_expr, :(->))
        # extract the parameter inputs
        pref_expr = _extract_parameters(raw_expr.args[1]) # this will create a copy if needed
        # fix the function if the parameter arguments were given as references
        if isexpr(raw_expr.args[1], :ref)
            raw_expr = _expr_replace!(raw_expr, raw_expr.args[1], gensym())
        elseif isexpr(raw_expr.args[1], :tuple)
            for a in raw_expr.args[1].args
                if isexpr(a, :ref)
                    raw_expr = _expr_replace!(raw_expr, a, gensym())
                end
            end
        end
        # extract the function expression
        func_expr = esc(raw_expr)
        is_anon = true
    else
        _error("Unrecognized syntax.")
    end
    return func_expr, pref_expr, is_anon
end

"""
    @parameter_function(model::InfiniteModel, func_expr, kwargs...)

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

```julia
@parameter_function(model::InfiniteModel, var_expr == func_expr, kwargs...)
```

Add a parameter function to `model` described by the expression `var_expr`, the
function expression `func_expr`, and the keyword arguments `kwargs`. The
expression `var_expr` is used to define the parameter function references of the
form `varname[...]` where the indexing matches the container syntax of other
macros.

The expression `func_expr` determines the concrete Julia function that defines the
behavior of parameter function and also specifies the infinite parameters it
depends on. The accepted forms are:
- `func(params...)`: where `func` is the function that takes supports of the
  infinite parameters `params` as input and outputs a scalar value.
- `(params...) -> my_func_expr`: where `params` are the infinite parameters and
  `my_func_expr` is the source code of the anonymous function.

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

**Examples**
```julia-repl
julia> @parameter_function(model, sin(t))
sin(t)

julia> func_vect = [sin, cos];

julia> @parameter_function(model, [i = 1:2] == func_vect[i](t))
2-element Array{GeneralVariableRef,1}:
 sin(t)
 cos(t)

julia> f(t_val, x_vals) = t_val + sum(x_vals)
f (generic function with 1 method)

julia> @parameter_function(model, pf == f(t, x))
pf(t, x)

julia> g(t_val, a; b = 0) = t_val + a + b
g (generic function with 1 method)

julia> @parameter_function(model, pf2[i = 1:2] == t -> g(t, i, b = 2 * i ))
2-element Array{GeneralVariableRef,1}:
 pf2[1](t)
 pf2[2](t)
```
"""
macro parameter_function(model, args...)
    # prepare the model
    esc_model = esc(model)

    # define error message function
    _error(str...) = _macro_error(:parameter_function, (model, args...),
                                  __source__, str...)

    # parse the arguments
    pos_args, kwargs, container_type, base_name = _extract_kwargs(args)

    # process the positional arguements
    if length(pos_args) == 1
        expr = popfirst!(pos_args)
    else
        _error("Incorrect amount of positional arguments given.")
    end
    if isexpr(expr, :call) && expr.args[1] === :(==)
        var = expr.args[2]
        func, prefs, is_anon_func = _process_func_expr(_error, expr.args[3])
        is_anon = isexpr(var, (:vect, :vcat))
    elseif isexpr(expr, (:call, :(->)))
        var = gensym()
        func, prefs, is_anon_func = _process_func_expr(_error, expr)
        is_anon = true
    else
        _error("Unrecognized syntax.")
    end

    # parse the name appropriately
    name = _get_name(var)
    if !isa(name, Symbol) && !is_anon
        _error("Expression $name should not be used as a parameter function name. Use " *
               "the \"anonymous\" syntax $name = @parameter_function(model, " *
               "...) instead.")
    end
    idxvars, inds = JuMPC.build_ref_sets(_error, var)
    if model in idxvars
        _error("Index $(model) is the same symbol as the model. Use a ",
               "different name for the index.")
    end
    if isnothing(base_name) && is_anon && !is_anon_func
        name_code = :( string(nameof($func)) )
    elseif isnothing(base_name) && is_anon
        name_code = _name_call("", idxvars)
    elseif isnothing(base_name)
        name_code = _name_call(string(name), idxvars)
    else
        name_code = _name_call(base_name, idxvars)
    end

    # make the build call
    build_call = :( build_parameter_function($_error, $func, $prefs) )
    _add_kwargs(build_call, kwargs)

    # make the creation code
    if isa(var, Symbol)
        # create a single parameter function
        creation_code = :( add_parameter_function($esc_model, $build_call, $name_code) )
    else
        # create code for a container
        pfunction_call = :( add_parameter_function($esc_model, $build_call,
                                                   $name_code) )
        creation_code = JuMPC.container_code(idxvars, inds, pfunction_call,
                                             container_type)
    end

    # finalize the macro
    if is_anon
        # anonymous creation
        macro_code = creation_code
    else
        # want to register the object/container reference
        macro_code = _macro_assign_and_return(_error, creation_code, name,
                                              esc_model)
    end
    return _finalize_macro(_error, esc_model, macro_code, __source__)
end


text:
mutable struct OnlineSVGP{
    T<:Real,
    TLikelihood<:AbstractLikelihood,
    TInference<:AbstractInference,
    TData<:AbstractDataContainer,
    N,
} <: AbstractGP{T,TLikelihood,TInference,N}
    data::TData
    f::NTuple{N,OnlineVarLatent{T}}
    likelihood::TLikelihood
    inference::TInference
    verbose::Int64
    atfrequency::Int64
    trained::Bool
end

"""
    OnlineSVGP(args...; kwargs...)

Online Sparse Variational Gaussian Process

## Arguments
- `kernel::Kernel` : Covariance function, can be any kernel from KernelFunctions.jl
- `likelihood` : Likelihood of the model. For compatibilities, see [`Likelihood Types`](@ref likelihood_user)
- `inference` : Inference for the model, see the [`Compatibility Table`](@ref compat_table))
- `Zalg` : Algorithm selecting how inducing points are selected

## Keywords arguments
- `verbose` : How much does the model print (0:nothing, 1:very basic, 2:medium, 3:everything)
- `optimiser` : Optimiser used for the kernel parameters. Should be an Optimiser object from the [Flux.jl](https://github.com/FluxML/Flux.jl) library, see list here [Optimisers](https://fluxml.ai/Flux.jl/stable/training/optimisers/) and on [this list](https://github.com/theogf/AugmentedGaussianProcesses.jl/tree/master/src/inference/optimisers.jl). Default is `ADAM(0.001)`
- `atfrequency::Int=1` : Choose how many variational parameters iterations are between hyperparameters optimization
- `mean=ZeroMean()` : PriorMean object, check the documentation on it [`MeanPrior`](@ref meanprior)
- `Zoptimiser` : Optimiser for inducing points locations
- `T::DataType=Float64` : Hint for what the type of the data is going to be.
"""
function OnlineSVGP(
    kernel::Kernel,
    likelihood::AbstractLikelihood,
    inference::AbstractInference,
    Zalg::InducingPoints.OnIPSA=OIPS(0.9);
    verbose::Integer=0,
    optimiser=ADAM(0.01),
    atfrequency::Integer=1,
    mean::Union{<:Real,AbstractVector{<:Real},PriorMean}=ZeroMean(),
    Zoptimiser=nothing,
    T::DataType=Float64,
)
    data = OnlineDataContainer()
    inference isa AnalyticVI || error("The inference object should be of type `AnalyticVI`")
    implemented(likelihood, inference) ||
        error("The $likelihood is not compatible or implemented with the $inference")

    if isa(optimiser, Bool)
        optimiser = optimiser ? ADAM(0.01) : nothing
    end

    if typeof(mean) <: Real
        mean = ConstantMean(mean)
    elseif typeof(mean) <: AbstractVector{<:Real}
        mean = EmpiricalMean(mean)
    end

    nLatent = num_latent(likelihood)
    latentf = ntuple(_ -> OnlineVarLatent(T, 0, 0, [], Zalg, kernel, mean, optimiser, Zoptimiser), nLatent)
    inference = tuple_inference(inference, nLatent, 0, 0, 0, [], [])
    inference.nIter = 1
    return OnlineSVGP{T,typeof(likelihood),typeof(inference),typeof(data),nLatent}(
        data, latentf, likelihood, inference, verbose, atfrequency, false
    )
end

function Base.show(io::IO, model::OnlineSVGP) where {T}
    return print(
        io,
        "Online Variational Gaussian Process with a $(likelihood(model)) infered by $(inference(model)) ",
    )
end

@traitimpl IsSparse{OnlineSVGP}

Zviews(m::OnlineSVGP) = Zview.(m.f)
objective(m::OnlineSVGP) = ELBO(m)
MBIndices(m::OnlineSVGP) = 1:nSamples(m)
xview(m::OnlineSVGP) = input(m)
yview(m::OnlineSVGP) = output(m)
nFeatures(m::OnlineSVGP) = collect(dim.(getf(m)))

## Accessors to InducingPoints methods
# function InducingPoints.init(Z::OptimIP, m::OnlineSVGP, gp::OnlineVarLatent)
#     return InducingPoints.init(Z.Z, m, gp)
# end
# InducingPoints.init(Z::OIPS, m, gp) = InducingPoints.init(Z, input(m), kernel(gp))

# function InducingPoints.add_point!(Z::OptimIP, m::OnlineSVGP, gp::OnlineVarLatent)
#     return InducingPoints.add_point!(Z.Z, m, gp)
# end
# function InducingPoints.add_point!(Z::OIPS, m::OnlineSVGP, gp::OnlineVarLatent)
#     return InducingPoints.add_point!(Z, input(m), kernel(gp))
# end

# function InducingPoints.remove_point!(Z::OptimIP, m::OnlineSVGP, gp::OnlineVarLatent)
#     return InducingPoints.remove_point!(Z.Z, m, gp)
# end
# function InducingPoints.remove_point!(Z::OIPS, m::OnlineSVGP, gp::OnlineVarLatent)
#     return InducingPoints.remove_point!(Z, pr_cov(gp), kernel(gp))
# end

ext: jl


text:
################################################################################
# Drone Example
################################################################################
# using Algames
using StaticArrays
using LinearAlgebra

vis = Visualizer()
open(vis)

T = Float64

# Define the dynamics of the system
p = 2 # Number of players
d = 3 # Number of dimensions
model = DoubleIntegratorGame(p=p, d=d) # game with 3 players with double integrator dynamics in 3D
n = model.n
m = model.m

# Define the horizon of the problem
N = 20 # N time steps
dt = 0.1 # each step lasts 0.1 second
probsize = ProblemSize(N,model) # Structure holding the relevant sizes of the problem

# Define the objective of each player
# We use a LQR cost
Q = [Diagonal(10*ones(SVector{model.ni[i],T})) for i=1:p] # Quadratic state cost
R = [Diagonal(0.1*ones(SVector{model.mi[i],T})) for i=1:p] # Quadratic control cost
# Desrired state
xf = [SVector{model.ni[1],T}([2,+0.7,+0.3,0,0,0]),
      SVector{model.ni[2],T}([2,-0.7,-0.3,0,0,0]),
      ]
# Desired control
uf = [zeros(SVector{model.mi[i],T}) for i=1:p]
# Objectives of the game
game_obj = GameObjective(Q,R,xf,uf,N,model)
radius = 1.0*ones(p)
μ = 5.0*ones(p)
add_collision_cost!(game_obj, radius, μ)

# Define the constraints that each player must respect
game_con = GameConstraintValues(probsize)
# Add collision avoidance
radius = 0.08
add_spherical_collision_avoidance!(game_con, radius)
# Add control bounds
u_max =  5*ones(SVector{m,T})
u_min = -5*ones(SVector{m,T})
add_control_bound!(game_con, u_max, u_min)
# # Add state bounds for player 1
x_max =  5*ones(SVector{n,T})
x_min = -5*ones(SVector{n,T})
add_state_bound!(game_con, 1, x_max, x_min)
# Add wall constraint
# walls = [Wall3D([1.,0,0], [2.,0,0], [2.,1,0], [0.,0,1])]
room_walls = [
    # Wall3D([1.,0.5,0], [2.,0,1], [2.5,1.0,1], [0.,0,1]),
    Wall3D([-3.00, -1.00, -1.00], [3.00, -1.00, -1.00], [3.00,  1.00, -1.00], [0.00,  0.00, -1.00]),
    Wall3D([-3.00, -1.00,  1.00], [3.00, -1.00,  1.00], [3.00,  1.00,  1.00], [0.00,  0.00,  1.00]),
    Wall3D([-3.00, -1.00, -1.00], [3.00, -1.00, -1.00], [3.00, -1.00,  1.00], [0.00, -1.00,  0.00]),
    Wall3D([-3.00,  1.00, -1.00], [3.00,  1.00, -1.00], [3.00,  1.00,  1.00], [0.00,  1.00,  0.00]),
    ]
door_walls = [
    Wall3D([-0.75, -1.00,  1.00], [0.00, -1.00,  0.25], [0.00,  1.00,  0.25], [ sqrt(2),  0.00,  sqrt(2)]),
    Wall3D([ 0.75, -1.00,  1.00], [0.00, -1.00,  0.25], [0.00,  1.00,  0.25], [-sqrt(2),  0.00,  sqrt(2)]),

    Wall3D([-0.75, -1.00, -1.00], [0.00, -1.00, -0.25], [0.00,  1.00, -0.25], [ sqrt(2),  0.00, -sqrt(2)]),
    Wall3D([ 0.75, -1.00, -1.00], [0.00, -1.00, -0.25], [0.00,  1.00, -0.25], [-sqrt(2),  0.00, -sqrt(2)]),
    #
    Wall3D([-0.75, -1.00, -1.00], [0.00, -0.25, -1.00], [0.00, -0.25,  1.00], [ sqrt(2), -sqrt(2),  0.00]),
    Wall3D([ 0.75, -1.00, -1.00], [0.00, -0.25, -1.00], [0.00, -0.25,  1.00], [-sqrt(2), -sqrt(2),  0.00]),
    #
    Wall3D([-0.75,  1.00, -1.00], [0.00,  0.25, -1.00], [0.00,  0.25,  1.00], [ sqrt(2),  sqrt(2),  0.00]),
    Wall3D([ 0.75,  1.00, -1.00], [0.00,  0.25, -1.00], [0.00,  0.25,  1.00], [-sqrt(2),  sqrt(2),  0.00]),
    ]
add_wall_constraint!(game_con, room_walls)
add_wall_constraint!(game_con, door_walls)

build_wall!(vis, room_walls, α=0.1, name=:Room)
build_wall!(vis, door_walls, α=0.1, name=:Door1)

# Define the initial state of the system
x0 = SVector{model.n,T}([
   -1.2,-1.0,
   -0.4, 0.4,
    0.6, 0.6,
    0.0, 0.0,
    0.0, 0.0,
    0.0, 0.0,
    ])

# Define the Options of the solver
opts = Options()
# Define the game problem
prob = GameProblem(N,dt,x0,model,opts,game_obj,game_con)

# Solve the problem
@time newton_solve!(prob)

# Visualize the Results
# const Algames = Main
using Plots
plot(prob.model, prob.pdtraj.pr)
plot(prob.stats)

build_traj!(vis, model, prob.pdtraj.pr, α=1.0, name=:Traj)

# build_drone!(vis, model, name=:Robot, α=1.0)
# set_drone(vis, model, Algames.state(prob.pdtraj.pr[1]), name=:Robot)
# animate_drone(vis, model, α=1.4, name=:Robot)
anim = visualize_drone!(vis, model, prob.pdtraj.pr)


text:
import MittagLefflerFunctions.MLF: mlf1

α = 1/2
β = 1.0
z = Complex(-0.5, 0.4)
@assert ( α < 1 ) && ( abs(z) ≤ 1 )

val = mlf1(α, β, z, 10)
ref = erfcx(-z)

@test abs(val-ref) ≤ 1e-10

α = 1.5
β = 1.0
z = Complex(5.5, 0.0)
@assert ( 1 ≤ α < 2 ) && ( abs(z) ≤ floor(20/(2.1-α)^(5.5-2α)) )

val = mlf1(α, β, z, 10)
Eq = MLQuad(α, β, 10, :hyperbola)
ref = Eq(real(z))

@test abs(val-ref) ≤ 1e-10 * abs(ref)

α = 4.0
β = 1.0
z = Complex(15.0, 17.3)
@assert ( α ≥ 2 ) && ( abs(z) ≤ 50 )

val = mlf1(α, β, z, 10)
ref = ( cos(z^(1/4)) + cosh(z^(1/4)) ) / 2

@test abs(val-ref) ≤ 1e-10 * abs(ref)


text:
include("testBICGSTB.jl")
include("testBlockBiCGSTB.jl")
include("testCG.jl")
include("testBlockCG.jl")
include("testCGLS.jl")
include("testGS.jl")
include("testGMRES.jl")
include("testLANCZOS.jl")
include("testSSOR.jl")
include("testLSQR.jl")
include("testMINRES.jl")
include("testMatrixDepot.jl")


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:3] <= 1)

                     @variable(m, -24.652335523355834 <= q <= 21.32472335033785)

                     add_NL_constraint(m, :(softplus(-0.6982139593401175 + -0.19504968807074974*softplus(0.6793369023557885 + 0.6523875537646115*softplus(0.1880507757579264 + -0.8703531991881865*$(x[1]) + 0.8423062039785143*$(x[2]) + 0.27626604928751597*$(x[3])) + -0.5142038874540145*softplus(0.8004597365811503 + 0.4093664270553261*$(x[1]) + -0.7983161406747903*$(x[2]) + -0.755209029087681*$(x[3])) + -0.0456439541419269*softplus(-0.8428376848805481 + -0.4971126206674761*$(x[1]) + -0.0245959153233013*$(x[2]) + 0.4386350201089879*$(x[3])) + 0.6096854328317232*softplus(-0.26516744777719925 + 0.12012704290515952*$(x[1]) + -0.4833932117044406*$(x[2]) + -0.25804110713494666*$(x[3]))) + -0.41803298572193937*softplus(-0.7048442383777926 + 0.47781322839334983*softplus(0.1880507757579264 + -0.8703531991881865*$(x[1]) + 0.8423062039785143*$(x[2]) + 0.27626604928751597*$(x[3])) + -0.713369260368808*softplus(0.8004597365811503 + 0.4093664270553261*$(x[1]) + -0.7983161406747903*$(x[2]) + -0.755209029087681*$(x[3])) + -0.44964139924772173*softplus(-0.8428376848805481 + -0.4971126206674761*$(x[1]) + -0.0245959153233013*$(x[2]) + 0.4386350201089879*$(x[3])) + -0.9876898123617992*softplus(-0.26516744777719925 + 0.12012704290515952*$(x[1]) + -0.4833932117044406*$(x[2]) + -0.25804110713494666*$(x[3]))) + 0.9542284451941447*softplus(0.9969387928417852 + -0.011395717836928387*softplus(0.1880507757579264 + -0.8703531991881865*$(x[1]) + 0.8423062039785143*$(x[2]) + 0.27626604928751597*$(x[3])) + 0.8103499615759371*softplus(0.8004597365811503 + 0.4093664270553261*$(x[1]) + -0.7983161406747903*$(x[2]) + -0.755209029087681*$(x[3])) + -0.12667586987008894*softplus(-0.8428376848805481 + -0.4971126206674761*$(x[1]) + -0.0245959153233013*$(x[2]) + 0.4386350201089879*$(x[3])) + 0.6252342222679741*softplus(-0.26516744777719925 + 0.12012704290515952*$(x[1]) + -0.4833932117044406*$(x[2]) + -0.25804110713494666*$(x[3]))) + 0.7916680056214345*softplus(-0.5156629799588863 + -0.22745979325095078*softplus(0.1880507757579264 + -0.8703531991881865*$(x[1]) + 0.8423062039785143*$(x[2]) + 0.27626604928751597*$(x[3])) + -0.8222624364306839*softplus(0.8004597365811503 + 0.4093664270553261*$(x[1]) + -0.7983161406747903*$(x[2]) + -0.755209029087681*$(x[3])) + -0.4937419854366709*softplus(-0.8428376848805481 + -0.4971126206674761*$(x[1]) + -0.0245959153233013*$(x[2]) + 0.4386350201089879*$(x[3])) + -0.1046356011041838*softplus(-0.26516744777719925 + 0.12012704290515952*$(x[1]) + -0.4833932117044406*$(x[2]) + -0.25804110713494666*$(x[3])))) + softplus(-0.4539405071565561 + 0.5320780591355936*softplus(0.6793369023557885 + 0.6523875537646115*softplus(0.1880507757579264 + -0.8703531991881865*$(x[1]) + 0.8423062039785143*$(x[2]) + 0.27626604928751597*$(x[3])) + -0.5142038874540145*softplus(0.8004597365811503 + 0.4093664270553261*$(x[1]) + -0.7983161406747903*$(x[2]) + -0.755209029087681*$(x[3])) + -0.0456439541419269*softplus(-0.8428376848805481 + -0.4971126206674761*$(x[1]) + -0.0245959153233013*$(x[2]) + 0.4386350201089879*$(x[3])) + 0.6096854328317232*softplus(-0.26516744777719925 + 0.12012704290515952*$(x[1]) + -0.4833932117044406*$(x[2]) + -0.25804110713494666*$(x[3]))) + -0.8489788623944037*softplus(-0.7048442383777926 + 0.47781322839334983*softplus(0.1880507757579264 + -0.8703531991881865*$(x[1]) + 0.8423062039785143*$(x[2]) + 0.27626604928751597*$(x[3])) + -0.713369260368808*softplus(0.8004597365811503 + 0.4093664270553261*$(x[1]) + -0.7983161406747903*$(x[2]) + -0.755209029087681*$(x[3])) + -0.44964139924772173*softplus(-0.8428376848805481 + -0.4971126206674761*$(x[1]) + -0.0245959153233013*$(x[2]) + 0.4386350201089879*$(x[3])) + -0.9876898123617992*softplus(-0.26516744777719925 + 0.12012704290515952*$(x[1]) + -0.4833932117044406*$(x[2]) + -0.25804110713494666*$(x[3]))) + -0.25620536200506905*softplus(0.9969387928417852 + -0.011395717836928387*softplus(0.1880507757579264 + -0.8703531991881865*$(x[1]) + 0.8423062039785143*$(x[2]) + 0.27626604928751597*$(x[3])) + 0.8103499615759371*softplus(0.8004597365811503 + 0.4093664270553261*$(x[1]) + -0.7983161406747903*$(x[2]) + -0.755209029087681*$(x[3])) + -0.12667586987008894*softplus(-0.8428376848805481 + -0.4971126206674761*$(x[1]) + -0.0245959153233013*$(x[2]) + 0.4386350201089879*$(x[3])) + 0.6252342222679741*softplus(-0.26516744777719925 + 0.12012704290515952*$(x[1]) + -0.4833932117044406*$(x[2]) + -0.25804110713494666*$(x[3]))) + 0.9264345634068367*softplus(-0.5156629799588863 + -0.22745979325095078*softplus(0.1880507757579264 + -0.8703531991881865*$(x[1]) + 0.8423062039785143*$(x[2]) + 0.27626604928751597*$(x[3])) + -0.8222624364306839*softplus(0.8004597365811503 + 0.4093664270553261*$(x[1]) + -0.7983161406747903*$(x[2]) + -0.755209029087681*$(x[3])) + -0.4937419854366709*softplus(-0.8428376848805481 + -0.4971126206674761*$(x[1]) + -0.0245959153233013*$(x[2]) + 0.4386350201089879*$(x[3])) + -0.1046356011041838*softplus(-0.26516744777719925 + 0.12012704290515952*$(x[1]) + -0.4833932117044406*$(x[2]) + -0.25804110713494666*$(x[3])))) + softplus(-0.4017849565754932 + 0.7092320156791541*softplus(0.6793369023557885 + 0.6523875537646115*softplus(0.1880507757579264 + -0.8703531991881865*$(x[1]) + 0.8423062039785143*$(x[2]) + 0.27626604928751597*$(x[3])) + -0.5142038874540145*softplus(0.8004597365811503 + 0.4093664270553261*$(x[1]) + -0.7983161406747903*$(x[2]) + -0.755209029087681*$(x[3])) + -0.0456439541419269*softplus(-0.8428376848805481 + -0.4971126206674761*$(x[1]) + -0.0245959153233013*$(x[2]) + 0.4386350201089879*$(x[3])) + 0.6096854328317232*softplus(-0.26516744777719925 + 0.12012704290515952*$(x[1]) + -0.4833932117044406*$(x[2]) + -0.25804110713494666*$(x[3]))) + -0.20657569919094643*softplus(-0.7048442383777926 + 0.47781322839334983*softplus(0.1880507757579264 + -0.8703531991881865*$(x[1]) + 0.8423062039785143*$(x[2]) + 0.27626604928751597*$(x[3])) + -0.713369260368808*softplus(0.8004597365811503 + 0.4093664270553261*$(x[1]) + -0.7983161406747903*$(x[2]) + -0.755209029087681*$(x[3])) + -0.44964139924772173*softplus(-0.8428376848805481 + -0.4971126206674761*$(x[1]) + -0.0245959153233013*$(x[2]) + 0.4386350201089879*$(x[3])) + -0.9876898123617992*softplus(-0.26516744777719925 + 0.12012704290515952*$(x[1]) + -0.4833932117044406*$(x[2]) + -0.25804110713494666*$(x[3]))) + 0.6260976974386958*softplus(0.9969387928417852 + -0.011395717836928387*softplus(0.1880507757579264 + -0.8703531991881865*$(x[1]) + 0.8423062039785143*$(x[2]) + 0.27626604928751597*$(x[3])) + 0.8103499615759371*softplus(0.8004597365811503 + 0.4093664270553261*$(x[1]) + -0.7983161406747903*$(x[2]) + -0.755209029087681*$(x[3])) + -0.12667586987008894*softplus(-0.8428376848805481 + -0.4971126206674761*$(x[1]) + -0.0245959153233013*$(x[2]) + 0.4386350201089879*$(x[3])) + 0.6252342222679741*softplus(-0.26516744777719925 + 0.12012704290515952*$(x[1]) + -0.4833932117044406*$(x[2]) + -0.25804110713494666*$(x[3]))) + 0.026773003462286393*softplus(-0.5156629799588863 + -0.22745979325095078*softplus(0.1880507757579264 + -0.8703531991881865*$(x[1]) + 0.8423062039785143*$(x[2]) + 0.27626604928751597*$(x[3])) + -0.8222624364306839*softplus(0.8004597365811503 + 0.4093664270553261*$(x[1]) + -0.7983161406747903*$(x[2]) + -0.755209029087681*$(x[3])) + -0.4937419854366709*softplus(-0.8428376848805481 + -0.4971126206674761*$(x[1]) + -0.0245959153233013*$(x[2]) + 0.4386350201089879*$(x[3])) + -0.1046356011041838*softplus(-0.26516744777719925 + 0.12012704290515952*$(x[1]) + -0.4833932117044406*$(x[2]) + -0.25804110713494666*$(x[3])))) + softplus(0.8541923459688499 + -0.7979821016654509*softplus(0.6793369023557885 + 0.6523875537646115*softplus(0.1880507757579264 + -0.8703531991881865*$(x[1]) + 0.8423062039785143*$(x[2]) + 0.27626604928751597*$(x[3])) + -0.5142038874540145*softplus(0.8004597365811503 + 0.4093664270553261*$(x[1]) + -0.7983161406747903*$(x[2]) + -0.755209029087681*$(x[3])) + -0.0456439541419269*softplus(-0.8428376848805481 + -0.4971126206674761*$(x[1]) + -0.0245959153233013*$(x[2]) + 0.4386350201089879*$(x[3])) + 0.6096854328317232*softplus(-0.26516744777719925 + 0.12012704290515952*$(x[1]) + -0.4833932117044406*$(x[2]) + -0.25804110713494666*$(x[3]))) + 0.2172669330681649*softplus(-0.7048442383777926 + 0.47781322839334983*softplus(0.1880507757579264 + -0.8703531991881865*$(x[1]) + 0.8423062039785143*$(x[2]) + 0.27626604928751597*$(x[3])) + -0.713369260368808*softplus(0.8004597365811503 + 0.4093664270553261*$(x[1]) + -0.7983161406747903*$(x[2]) + -0.755209029087681*$(x[3])) + -0.44964139924772173*softplus(-0.8428376848805481 + -0.4971126206674761*$(x[1]) + -0.0245959153233013*$(x[2]) + 0.4386350201089879*$(x[3])) + -0.9876898123617992*softplus(-0.26516744777719925 + 0.12012704290515952*$(x[1]) + -0.4833932117044406*$(x[2]) + -0.25804110713494666*$(x[3]))) + -0.6626207135044284*softplus(0.9969387928417852 + -0.011395717836928387*softplus(0.1880507757579264 + -0.8703531991881865*$(x[1]) + 0.8423062039785143*$(x[2]) + 0.27626604928751597*$(x[3])) + 0.8103499615759371*softplus(0.8004597365811503 + 0.4093664270553261*$(x[1]) + -0.7983161406747903*$(x[2]) + -0.755209029087681*$(x[3])) + -0.12667586987008894*softplus(-0.8428376848805481 + -0.4971126206674761*$(x[1]) + -0.0245959153233013*$(x[2]) + 0.4386350201089879*$(x[3])) + 0.6252342222679741*softplus(-0.26516744777719925 + 0.12012704290515952*$(x[1]) + -0.4833932117044406*$(x[2]) + -0.25804110713494666*$(x[3]))) + -0.5735396155482761*softplus(-0.5156629799588863 + -0.22745979325095078*softplus(0.1880507757579264 + -0.8703531991881865*$(x[1]) + 0.8423062039785143*$(x[2]) + 0.27626604928751597*$(x[3])) + -0.8222624364306839*softplus(0.8004597365811503 + 0.4093664270553261*$(x[1]) + -0.7983161406747903*$(x[2]) + -0.755209029087681*$(x[3])) + -0.4937419854366709*softplus(-0.8428376848805481 + -0.4971126206674761*$(x[1]) + -0.0245959153233013*$(x[2]) + 0.4386350201089879*$(x[3])) + -0.1046356011041838*softplus(-0.26516744777719925 + 0.12012704290515952*$(x[1]) + -0.4833932117044406*$(x[2]) + -0.25804110713494666*$(x[3])))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
using DifferentialEquations, ModelingToolkit
using HoJMTK.Electrical
using Plots

# Based on https://mtk.sciml.ai/dev/tutorials/acausal_components/
R = 1.0
C = 1.0
V = 1.0
@named resistor = Resistor(R=R)
@named capacitor = Capacitor(C=C)
@named source = DCVoltageSource(V=V)
@named ground = Ground()

rc_eqs = [
          connect(source.p, resistor.p)
          connect(resistor.n, capacitor.p)
          connect(capacitor.n, source.n, ground.g)
         ]

@variables t
@named _rc_model = ODESystem(rc_eqs, t)
@named rc_model = compose(_rc_model, [resistor, capacitor, source, ground])
sys = structural_simplify(rc_model)
u0 = [
      capacitor.v => 0.0
     ]
prob = ODAEProblem(sys, u0, (0, 10.0))
sol = solve(prob, Tsit5())

plot(sol)


text:
#----------------------------------------------------------------------------------------------#
#                               Fundamental Phys-Chem Constants                                #
#----------------------------------------------------------------------------------------------#

"""
`const _NA = _Amt(measurement("6.0221415(10)e+23") / u"mol")`\n
The Avogadro constant, \$N_A\$, [Lide, D. R., 2006], as a `_Amt{Float64,MM}`.
"""
const _NA = _Amt(measurement("6.0221415(10)e+23") / u"mol")

"""
`NA(𝖯::Type{𝗽}=Float64, 𝖷::Type{𝘅}=EX) where {𝗽<:PREC,𝘅<:EXAC}`\n
Returns the Avogadro constant as a `_Amt{𝖯,𝖷}`.\n
Arguments `𝖯` and `𝖷` can be ommitted and/or be supplied in any order.
"""
NA(𝖯::Type{𝗽}=Float64, 𝖷::Type{𝘅}=EX) where {𝗽<:PREC,𝘅<:EXAC} = _Amt{𝖯,𝖷}(_NA)
NA(𝖷::Type{𝘅}, 𝖯::Type{𝗽}=Float64) where {𝗽<:PREC,𝘅<:EXAC} = _Amt{𝖯,𝖷}(_NA)

"""
`const _mu = m(measurement("1.66053886(28)e-27") ,SY)`\n
The atomic mass constant, \$m_u = (1/12)m(¹²C)\$, [Lide, D. R., 2006], as a
`mAmt{Float64,MM,SY}`.
"""
const _mu = m(measurement("1.66053886(28)e-27") ,SY)

"""
`mu(𝖯::Type{𝗽}=Float64, 𝖷::Type{𝘅}=EX) where {𝗽<:PREC,𝘅<:EXAC}`\n
Returns the atomic mass constant as a `mAmt{𝖯:𝖷}`.\n
Arguments `𝖯` and `𝖷` can be ommitted and/or be supplied in any order.
"""
mu(𝖯::Type{𝗽}=Float64, 𝖷::Type{𝘅}=EX) where {𝗽<:PREC,𝘅<:EXAC} = mAmt{𝖯,𝖷}(_mu)
mu(𝖷::Type{𝘅}, 𝖯::Type{𝗽}=Float64) where {𝗽<:PREC,𝘅<:EXAC} = mAmt{𝖯,𝖷}(_mu)

"""
`const _R̄ = R(measurement("8.314472(15)"), MO)`\n
The molar gas constant, \$R̄\$, [Lide, D. R., 2006], as a `RAmt{Float64,MM,MO}`.
"""
const _R̄ = R(measurement("8.314472(15)"), MO)

"""
`R(𝖯::Type{𝗽}=Float64, 𝖷::Type{𝘅}=EX) where {𝗽<:PREC,𝘅<:EXAC}`\n
Returns the molar gas constant as a `RAmt{𝖯:𝖷,MO}`.\n
Arguments `𝖯` and `𝖷` can be ommitted and/or be supplied in any order.
"""
R(𝖯::Type{𝗽}=Float64, 𝖷::Type{𝘅}=EX) where {𝗽<:PREC,𝘅<:EXAC} = RAmt{𝖯,𝖷}(_R̄)
R(𝖷::Type{𝘅}, 𝖯::Type{𝗽}=Float64) where {𝗽<:PREC,𝘅<:EXAC} = RAmt{𝖯,𝖷}(_R̄)

"""
`const _kB = _Amt(measurement("1.3806505(24)e-23") * u"J/K")`\n
The Boltzmann constant, \$k_B = R̄/N_A\$, [Lide, D. R., 2006], as a `_Amt{Float64,MM}`.
"""
const _kB = _Amt(measurement("1.3806505(24)e-23") * u"J/K")

"""
`kB(𝖯::Type{𝗽}=Float64, 𝖷::Type{𝘅}=EX) where {𝗽<:PREC,𝘅<:EXAC}`\n
Returns the Boltzmann constant as a `_Amt{𝖯,𝖷}`.\n
Arguments `𝖯` and `𝖷` can be ommitted and/or be supplied in any order.
"""
kB(𝖯::Type{𝗽}=Float64, 𝖷::Type{𝘅}=EX) where {𝗽<:PREC,𝘅<:EXAC} = _Amt{𝖯,𝖷}(_kB)
kB(𝖷::Type{𝘅}, 𝖯::Type{𝗽}=Float64) where {𝗽<:PREC,𝘅<:EXAC} = _Amt{𝖯,𝖷}(_kB)

export NA, mu, kB   # as R is already exported on "amounts.jl"

#----------------------------------------------------------------------------------------------#
#                                     Reference Constants                                      #
#----------------------------------------------------------------------------------------------#

"""
`const _stdT = sysT{Float64,MM}(T(25u"°C"))`\n
The `sysT{Float64,MM}` representation of the exact standard temperature, \$T_0 ≡ 25°C\$, [Lide,
D. R., 2006].
"""
const _stdT = sysT{Float64,MM}(T(25u"°C"))

"""
`T(𝖯::Type{𝗽}=Float64, 𝖷::Type{𝘅}=EX) where {𝗽<:PREC,𝘅<:EXAC}`\n
Returns the standard temperature as a `sysT{𝖯,𝖷}`.\n
Arguments `𝖯` and `𝖷` can be ommitted and/or be supplied in any order.
"""
T(𝖯::Type{𝗽}=Float64, 𝖷::Type{𝘅}=EX) where {𝗽<:PREC,𝘅<:EXAC} = sysT{𝖯,𝖷}(_stdT)
T(𝖷::Type{𝘅}, 𝖯::Type{𝗽}=Float64) where {𝗽<:PREC,𝘅<:EXAC} = sysT{𝖯,𝖷}(_stdT)

"""
`const _stdP = sysP{Float64,MM}(P(101350u"Pa"))`\n
The `sysP{Float64,MM}` representation of the exact standard atmosphere, \$P_0 ≡ 101350Pa\$,
[Lide, D. R., 2006].
"""
const _stdP = sysP{Float64,MM}(P(101350u"Pa"))

"""
`P(𝖯::Type{𝗽}=Float64, 𝖷::Type{𝘅}=EX) where {𝗽<:PREC,𝘅<:EXAC}`\n
Returns the standard atmosphere as a `sysP{𝖯,𝖷}`.\n
Arguments `𝖯` and `𝖷` can be ommitted and/or be supplied in any order.
"""
P(𝖯::Type{𝗽}=Float64, 𝖷::Type{𝘅}=EX) where {𝗽<:PREC,𝘅<:EXAC} = sysP{𝖯,𝖷}(_stdP)
P(𝖷::Type{𝘅}, 𝖯::Type{𝗽}=Float64) where {𝗽<:PREC,𝘅<:EXAC} = sysP{𝖯,𝖷}(_stdP)

"""
`const _gn = grav{Float64,MM}(grav(9_806_650u"μm/s^2"))`\n
The `grav{Float64,MM}` representation of the exact standard gravity, \$g_n ≡ 9.80665 m/s^2\$,
[Lide, D. R., 2006].
"""
const _gn = grav{Float64,MM}(grav(9_806_650u"μm/s^2"))

"""
`grav(𝖯::Type{𝗽}=Float64, 𝖷::Type{𝘅}=E𝖷) where {𝗽<:PREC,𝘅<:EXAC}`\n
Returns the standard gravity as a `grav{𝖯,𝖷}`.\n
Arguments `𝖯` and `𝖷` can be ommitted and/or be supplied in any order.
"""
grav(𝖯::Type{𝗽}=Float64, 𝖷::Type{𝘅}=EX) where {𝗽<:PREC,𝘅<:EXAC} = grav{𝖯,𝖷}(_gn)
grav(𝖷::Type{𝘅}, 𝖯::Type{𝗽}=Float64) where {𝗽<:PREC,𝘅<:EXAC} = grav{𝖯,𝖷}(_gn)


text:
export order, roots, any_root, issquare, next_prime

######################################################################
# tools.jl: useful tools for EllipticCurves
######################################################################

######################################################################
# Rational fractions
######################################################################

function (P::MPoly)(args...)
	return evaluate(P, [args...])
end

function (P::Frac)(args...)
	n = P.num
	d = P.den
	return n(args...) // d(args...)
end

######################################################################
# Useful functions for finite fields
######################################################################

#Multiplicative order in (small) finite fields

function order(x::FinFieldElem)
	K = parent(x)
	if x == Nemo.zero(K)
		throw(ArgumentError("Argument must be non-zero"))
	else
		k = 1
		z = x
		while z != Nemo.one(K)
			z *= x
			k += 1
		end
		return k
	end
end

#Computing roots of polynomials over finite fields

function roots(f::PolyElem{T}) where T<:FinFieldElem
	fac = Nemo.factor(f)
	res = Dict{T, Int}()
	for index in fac
		P, exp = index
		if degree(P) == 1
			res[- coeff(P, 0)] = exp
		end
	end
	res = collect(res)
	return res
end

function any_root(f::PolyElem{T}) where T<:FinFieldElem
	r = roots(f)
	if any(_->true, r) #tests whether r contains any element
		y0, _ = r[1]
		return (true, y0)
	else
		return (false, Nemo.zero(base_ring(f)))
	end
end

function issquare(x::FinFieldElem)
	K = parent(x)
	A, y = PolynomialRing(K, "y")
	return any_root(y^2 - x)
end

#Conversions between finite fields

"""
Convert an element of a finite field to another one.

No control is made on the input.
"""
function convert(x::FinFieldElem, K::FinField)
	K1 = parent(x)
	p1 = characteristic(K1)
	p = characteristic(K)
	(p1 == p) || throw(ArgumentError("Fields must have the same characteristic"))
	y = deepcopy(x)
	y.parent = K
	return y
end

"""
When used with polynomials, convert each coefficient.
"""
function convert(P::PolyElem{T}, K::FinField) where T<:FinFieldElem
	A, Y = PolynomialRing(K, "Y")
	poly = Nemo.zero(A)
	for i = 0:degree(P)
		Nemo.setcoeff!(poly, i, convert(coeff(P, i), K))
	end
	return poly
end

function sqrt(P::PolyElem)
	return gcd(P, derivative(P))
end


text:
module TestCompare
    using Benchmark

    f1() = svd(zeros(100, 100))
    f2() = svd(ones(100, 100))
    f3() = svd(eye(100, 100))

    benchmark(f1, "Linear Algebra", "svd(zeros(100, 100))", 10)
    benchmark(f1, "svd(zeros(100, 100))", 10)

    benchmark(f2, "Linear Algebra", "svd(ones(100, 100))", 10)
    benchmark(f2, "svd(ones(100, 100))", 10)

    compare([f1, f2, f3], 100)
end


text:
using SIMDPirates
using LinearAlgebra: Transpose, Adjoint
using Base.Cartesian: @nexprs, @nif
using ..LoopInfo
using LoopVectorization: @avx
using SIMDPirates: vfmadd231
using VectorizationBase: VectorizationBase
using TimerOutputs: TimerOutputs, @timeit_debug, TimerOutput
using UnPack: @unpack

# General direction: we want to avoid pointer arithmetic as much as possible,
# also, don't strict type the container type

# I'll try to implement everything without defining a type, and all variables
# should stay local if possible. We can always clean up the code later when we
# find a good direction/structure.

###
### Prelude
###

const N_THREADS_BUFFERS = map(_->Vector{UInt8}(undef, 0), 1:Threads.nthreads())
const PAGE_SIZE = ccall(:jl_getpagesize, Int, ())

const BLAS_TIMER = TimerOutput()
reset_timer!() = TimerOutputs.reset_timer!(BLAS_TIMER)
get_timer() = BLAS_TIMER
enable_timer() = TimerOutputs.enable_debug_timings(@__MODULE__)
disable_timer() = TimerOutputs.disable_debug_timings(@__MODULE__)

# alias scopes from https://github.com/JuliaLang/julia/pull/31018
struct Const{T<:Array}
    a::T
end

@eval Base.getindex(A::Const, i1::Int) = Core.const_arrayref($(Expr(:boundscheck)), A.a, i1)
@eval Base.getindex(A::Const, i1::Int, i2::Int, I::Int...) = (Base.@_inline_meta; Core.const_arrayref($(Expr(:boundscheck)), A.a, i1, i2, I...))
constarray(a::Array) = Const(a)
constarray(a) = a

macro aliasscope(body)
    sym = gensym()
    esc(quote
        $(Expr(:aliasscope))
        $sym = $body
        $(Expr(:popaliasscope))
        $sym
    end)
end

###
### User-level API
###

function mul!(C, A, B, α=true, β=false; cache_params=(cache_m=72, cache_k=256, cache_n=4080), kernel_params=(Val(8), Val(6)), packing=nothing)
    m, k, n = checkmulsize(C, A, B)
    α = convert(eltype(C), α)
    β = convert(eltype(C), β)
    iszeroα = iszero(α)
    if iszero(β) && iszeroα
        return fill!(C, zero(eltype(C)))
    elseif isone(β) && iszeroα
        return C
    elseif iszeroα
        # TODO: optimize this
        for ii in eachindex(C)
            Cii = @inbounds C[ii]
            Cii = β * Cii
            @inbounds C[ii] = Cii
        end
        return C
    end
    # TODO: packing strategy
    packing === nothing && (packing = (false, false))
    packa, packb = packing
    if packa isa Bool && packb isa Bool
        if packa # manual union split
            if packb
                _mul!(C, A, B, α, β, (Val(true),Val(true)), cache_params, kernel_params)
            else
                _mul!(C, A, B, α, β, (Val(true),Val(false)), cache_params, kernel_params)
            end
        else
            if packb
                _mul!(C, A, B, α, β, (Val(false),Val(true)), cache_params, kernel_params)
            else
                _mul!(C, A, B, α, β, (Val(false),Val(false)), cache_params, kernel_params)
            end
        end
    else
        _mul!(C, A, B, α, β, packing, cache_params, kernel_params)
    end
    return C
end

partition_k(k, cache_k) = cld(k, cld(k, cache_k))

###
### Lower-level `_mul!`
###

# use lowercase here because these are not actually types
struct GemmContext{packa,packb,micro_m,micro_n,C1,C2,A1,A2,B1,B2}
    cache_m::Int
    cache_k::Int
    cache_n::Int
    cs1::C1
    cs2::C2
    as1::A1
    as2::A2
    bs1::B1
    bs2::B2
end
function GemmContext{packa,packb,micro_m,micro_n}(cache_m, cache_k, cache_n, cs1, cs2, as1, as2, bs1, bs2) where {packa,packb,micro_m,micro_n}
    return GemmContext{packa,packb,micro_m,micro_n,
                       typeof(cs1),typeof(cs2),
                       typeof(as1),typeof(as2),
                       typeof(bs1),typeof(bs2)}(
                                                cache_m, cache_k, cache_n,
                                                cs1, cs2,
                                                as1, as2,
                                                bs1, bs2
                                               )
end

function _mul!(C, A, B, α, β, packing::Tuple{Val{packa},Val{packb}}, (cache_m, cache_k, cache_n), kernel_params::Tuple{Val{micro_m},Val{micro_n}}) where {packa,packb,micro_m,micro_n}
    cs1, cs2 = _strides(C)
    as1, as2 = _strides(A)
    bs1, bs2 = _strides(B)
    cs1 == 1 || throw(ArgumentError("Packing kernel doesn't support nonunit leading stride C matrix. Got stride(C, 1) = $cs1."))
    if !packa
        as1 == 1 || throw(ArgumentError("Kernel with packed A doesn't support nonunit leading stride C and A matrices. Got stride(C, 1) = $cs1, stride(A, 1) = $as1, and stride(B, 1) = $bs1."))
    end
    ctx = GemmContext{packa,packb,micro_m,micro_n}(cache_m, cache_k, cache_n,
                                                   cs1, cs2, as1, as2, bs1, bs2)
    m, k, n = checkmulsize(C, A, B)
    A′, B′ = canonicalize(A), canonicalize(B)
    if !packa && !packb
        macrokernel!(C, A′, B′, α, β, (1, m), (1, k), (1, n), ctx)
        return C
    end

    cache_k = partition_k(k, cache_k) # More evenly divide K

    # get buffer
    Abuffersize = cache_m * cache_k
    Bbuffersize = cache_k * cache_n
    T = eltype(C)
    buffer = N_THREADS_BUFFERS[Threads.threadid()]
    resize!(buffer, (Abuffersize + Bbuffersize) * sizeof(T) + 2PAGE_SIZE)
    ptrbuffer = align(convert(Ptr{T}, pointer(buffer)), PAGE_SIZE)
    Abuffer = packa ? unsafe_wrap(Array, ptrbuffer, Abuffersize) : A′
    ptrbuffer = align(ptrbuffer + Abuffersize * sizeof(T), PAGE_SIZE)
    Bbuffer = packb ? unsafe_wrap(Array, ptrbuffer, Bbuffersize) : B′
    ptrbuffer = align(ptrbuffer + Bbuffersize * sizeof(T), PAGE_SIZE)

    for cachej₁ in 1:cache_n:n; cachej₂ = min(cachej₁ + cache_n - 1, n)
        for cachep₁ in 1:cache_k:k; cachep₂ = min(cachep₁ + cache_k - 1, k) # TODO: add adaptive packing?
            β′ = cachep₁ == 1 ? β : one(β)
            packb && @timeit_debug BLAS_TIMER "Pack B" packBbuffer!(Bbuffer, B, cachep₁, cachep₂, cachej₁, cachej₂, Val(micro_n))
            for cachei₁ in 1:cache_m:m; cachei₂ = min(cachei₁ + cache_m - 1, m)
                packa && @timeit_debug BLAS_TIMER "Pack A" packAbuffer!(Abuffer, A, cachei₁, cachei₂, cachep₁, cachep₂, Val(micro_m))
                # macrokernel
                macrokernel!(C, Abuffer, Bbuffer, α, β′, (cachei₁, cachei₂), (cachep₁, cachep₂), (cachej₁, cachej₂), ctx)
            end # cachei
        end # cachep
    end # cachej
    return C
end

###
### Macro kernel
###

@generated function macrokernel!(C, A, B, α, β, cacheis, cacheps, cachejs, ctx::GemmContext{packa,packb,micro_m,micro_n})::Nothing where {packa,packb,micro_m,micro_n}
    T = eltype(C)
    N = VectorizationBase.pick_vector_width(T)
    micro_m % N == 0 || error("`micro_m` must be an integer multiple of vector register width for efficient microkernel generation, got width = $N, micro_m = $micro_m.")
    mregister = max(1, micro_m ÷ N)

    mloopexpr = quote
        ii = cachei₁
        m′ = $micro_m
        while (m = cachei₂ - ii + 1) >= $micro_m
            Coffset = (ii - 1) * cs1 + (jj - 1) * cs2
            Aoffset′ = Aoffset + (packa ? (ii - cachei₁) * ps : (ii - 1) * as1)
            Boffset′ = Boffset + (packb ? (jjfull - cachej₁) * ps + micro_n_offset : (jj - 1) * bs2)
            microkernel!(C, A, B, α, β, Coffset, Aoffset′, Boffset′, ps, (Val(m′), Val(n′)), Val(4), nothing, ctx) # no mask
            ii += m′
        end
        m = cachei₂ - ii + 1
        mask = VectorizationBase.mask($T, m)
        Coffset = (ii - 1) * cs1 + (jj - 1) * cs2
        Aoffset′ = Aoffset + (packa ? (ii - cachei₁) * ps : (ii - 1) * as1)
        Boffset′ = Boffset + (packb ? (jjfull - cachej₁) * ps + micro_n_offset : (jj - 1) * bs2)
        @nif $(mregister+1) midx -> begin
            m′ = $micro_m - (midx - 1) * $N
            m > m′ - $N
        end midx -> begin
            microkernel!(C, A, B, α, β, Coffset, Aoffset′, Boffset′, ps, (Val(m′), Val(n′)), Val(1), mask, ctx) # mask
        end midx -> nothing
        jj += n′
    end

    nloopexpr = quote
        cachei₁, cachei₂ = cacheis
        cachep₁, cachep₂ = cacheps
        cachej₁, cachej₂ = cachejs
        ps = cachep₂ - cachep₁ + 1
        @unpack cs1, cs2, as1, as2, bs1, bs2 = ctx
        jj = jjfull = cachej₁
        micro_n_offset = 0
        Aoffset = packa ? 0 : (cachep₁ - 1) * as2
        Boffset = packb ? 0 : (cachep₁ - 1) * bs1
        n′ = $micro_n
        while (n = cachej₂ - jj + 1) >= n′
            $mloopexpr
            jjfull = jj
        end # full n
        jjfull = jj # set `jjfull` after completing a full `micro_n`

        @nif $micro_n nidx -> begin
            n′ = $micro_n - nidx
            n = cachej₂ - jj + 1
            n >= n′
        end nidx -> begin
            $mloopexpr
            micro_n_offset += n′
        end nidx -> nothing
        return nothing
    end
    # macroexpand(@__MODULE__, nloopexpr)
    return nloopexpr
end

###
### Packing
###

function packAbuffer!(Abuffer, A, cachei₁, cachei₂, cachep₁, cachep₂, ::Val{micro_m}) where micro_m
    # terminology:
    # a `micro_m × 1` vector is a panel
    # a `micro_m × cache_k` matrix is a block
    l = 0 # write position in packing buffer
    @inbounds @aliasscope for i₁ in cachei₁:micro_m:cachei₂ # iterate over blocks
        i₂ = i₁ + micro_m - 1
        if i₂ <= cachei₂ # full panel
            # copy a panel to the buffer contiguously
            @simd for p in cachep₁:cachep₂
                @unroll for i in i₁:i₂ # has constant loop count
                    Abuffer[l += 1] = constarray(A)[i, p]
                end
            end
        else # a panel is not full
            for p in cachep₁:cachep₂ # iterate through panel columns
                for i in i₁:cachei₂  # iterate through live panel rows
                    Abuffer[l += 1] = constarray(A)[i, p]
                end
                for i in (cachei₂ + 1):i₂ # pad the rest of the panel with zero
                    Abuffer[l += 1] = zero(eltype(Abuffer))
                end
            end
        end
    end
    return nothing
end

function packBbuffer!(Bbuffer, B, cachep₁, cachep₂, cachej₁, cachej₂, ::Val{micro_n}) where micro_n
    # terminology:
    # a `1 × micro_n` vector is a panel
    # a `cache_k × micro_n` matrix is a block
    l = 0 # write position in packing buffer
    @inbounds @aliasscope for j₁ in cachej₁:micro_n:cachej₂ # iterate over blocks
        j₂ = j₁ + micro_n - 1
        if j₂ <= cachej₂ # full panel
            # copy a panel to the buffer contiguously
            @simd for p in cachep₁:cachep₂
                @unroll for j in j₁:j₂ # has constant loop count
                    Bbuffer[l += 1] = constarray(B)[p, j]
                end
            end
        else # a panel is not full
            for p in cachep₁:cachep₂ # iterate through panel columns
                for j in j₁:cachej₂  # iterate through live panel rows
                    Bbuffer[l += 1] = constarray(B)[p, j]
                end
                for j in (cachej₂ + 1):j₂ # pad the rest of the panel with zero
                    Bbuffer[l += 1] = zero(eltype(Bbuffer))
                end
            end
        end
    end
    return nothing
end

###
### Micro kernel
###

"""
    microkernel!(C, Abuffer, Bbuffer, α, β, Coffset, Aoffset, Boffset, m, ps,
        kernel_params::Tuple{Val{micro_m},Val{micro_n}}) where {micro_m,micro_n}

A kernel that computes
```julia
Ĉ[1:m, 1:n] = α * Â[1:m, 1:ps] * B̂[1:ps, 1:n] + β * Ĉ[1:m, 1:n]
```
where ``{⋅̂}`` denotes the matrix with offset.
"""
@generated function microkernel!(C, A, B, α, β, Coffset, Aoffset, Boffset, ps, ::Tuple{Val{micro_m},Val{micro_n}},::Val{unroll}, mask,
                                 ctx::GemmContext{packa,packb,fullmicro_m,fullmicro_n})::Nothing where {micro_m,micro_n,unroll,packa,packb,fullmicro_m,fullmicro_n}
    T = eltype(C)
    TA = eltype(A)
    TB = eltype(B)
    T === TA === TB || throw(ArgumentError("eltype(C) === eltype(A) === eltype(B) must hold in the microkernel! Got eltype(C) = $T, eltype(A) = $TA, eltype(B) = $TB."))
    N = VectorizationBase.pick_vector_width(T)
    mregister = max(1, micro_m ÷ N)
    V = SVec{N,T}
    dounroll = unroll > 1
    matA = !packa
    matB = !packb
    usemask = mask != Nothing

    nounroll_loopexpr = quote
        prefetcht0(ptrÂ + 64 * $fullmicro_m + 2 * $fullmicro_n)
        @nexprs $mregister m̂ -> begin
            # assumption: A has unit leading stride
            A_m̂ = vload($V, ptrÂ + (m̂ - 1) * $N * st, mask_m̂)
        end
        @nexprs $micro_n n̂ -> begin
            if $matB
                B_n̂ = $V(vload(ptrB̂ + (n̂ - 1) * bs2))
            else
                B_n̂ = $V(vload(ptrB̂ + (n̂ - 1) * st))
            end
            @nexprs $mregister m̂ -> begin
                AB_m̂_n̂ = vfmadd231(A_m̂, B_n̂, AB_m̂_n̂)
            end
        end
        ptrÂ += $matA ? as2 : $fullmicro_m * st
        ptrB̂ += $matB ? bs1 : $fullmicro_n * st
    end

    unroll_loopexpr = quote
        @nexprs $unroll u -> begin
            #TODO
            if !$matA
                u == 1 && prefetcht0(ptrÂ + 64 * $micro_m + 2 * $fullmicro_n)
                u == 3 && prefetcht0(ptrÂ + 76 * $micro_m + 2 * $fullmicro_n)
            end
            ## unroll variable: u
            @nexprs $mregister m̂ -> begin
                # assumption: A has unit leading stride, as1 == 1
                ptrA′ = ptrÂ + (u - 1) * ($matA ? as2 : $fullmicro_m * st)
                A_m̂ = vload($V, ptrA′ + (m̂ - 1) * $N * st, mask_m̂)
            end
            @nexprs $micro_n n̂ -> begin
                ptrB′ = ptrB̂ + (n̂ - 1) * ($matB ? bs2 : st) + (u - 1) * ($matB ? bs1 : $fullmicro_n * st)
                B_n̂ = $V(unsafe_load(ptrB′))
                @nexprs $mregister m̂ -> begin
                    AB_m̂_n̂ = vfmadd231(A_m̂, B_n̂, AB_m̂_n̂)
                end
            end
        end
        ptrÂ += $matA ? $unroll * as2 : $unroll * $fullmicro_m * st
        ptrB̂ += $matB ? $unroll * bs1 : $unroll * $fullmicro_n * st
    end

    kernel_code = quote
        # tell the compiler that the iteration is nonempty
        SIMDPirates.assume(ps > 0)

        st = sizeof($T)
        cs1 = 1 * st # ctx.cs1 * st
        cs2 = ctx.cs2 * st
        if $matA # A is not packed
            as1 = 1 * st # ctx.as1 * st
            as2 = ctx.as2 * st
        end
        if $matB # B is not packed
            bs1 = ctx.bs1 * st
            bs2 = ctx.bs2 * st
        end
        @nexprs $mregister m̂ -> begin
            mask_m̂ = ($usemask && m̂ == $mregister) ? mask : # nontrival mask
                                                     VectorizationBase.max_mask($T) # trival mask that should be optimized away
        end

        ptrĈ = pointer(C) + Coffset * st
        ptrÂ = pointer(A) + Aoffset * st
        ptrB̂ = pointer(B) + Boffset * st

        # prefetch C matrix
        @nexprs $micro_n n̂ -> begin
            prefetcht0(ptrĈ + (n̂ - 1) * cs2 + 7 * 8)
        end

        # intializing AB registers
        @nexprs $micro_n n̂ -> @nexprs $mregister m̂ -> AB_m̂_n̂ = zero($V)

        # rank-1 updates
        if $dounroll
            punroll, pleft = divrem(ps, $unroll)
            # tell the compiler that the iteration is nonempty
            SIMDPirates.assume(punroll > 0 || pleft > 0)
            for _ in 1:punroll
                $unroll_loopexpr
            end

            for _ in 1:pleft
                $nounroll_loopexpr
            end
        else
            for _ in 1:ps
                $nounroll_loopexpr
            end
        end

        _α = $V(α)
        @nexprs $micro_n n̂ -> @nexprs $mregister m̂ -> begin
            C_m̂_n̂ = _α * AB_m̂_n̂
        end

        if iszero(β)
            @nexprs $micro_n n̂ -> @nexprs $mregister m̂ -> begin
                addr = ptrĈ + (m̂ - 1) * $N * cs1 + (n̂ - 1) * cs2
                vstore!(addr, C_m̂_n̂, mask_m̂)
            end
        else
            _β = $V(β)
            @nexprs $micro_n n̂ -> @nexprs $mregister m̂ -> begin
                addr = ptrĈ + (m̂ - 1) * $N * cs1 + (n̂ - 1) * cs2
                C_m̂_n̂ = vfmadd231(_β, vload($V, addr, mask_m̂), C_m̂_n̂)
                vstore!(addr, C_m̂_n̂, mask_m̂)
            end
        end
    end

    expr = quote
        begin
            $(Expr(:meta, usemask ? :noinline : :inline))
            $kernel_code
            return nothing
        end
    end
    return expr
end

###
### Fast stride handling
###

struct One
end
@inline Base.:*(::One, x::Int) = x
@inline Base.:*(x::Int, ::One) = x
@inline Base.:+(::One, x::Int) = x + 1
@inline Base.:+(x::Int, ::One) = x + 1
@inline Base.:(==)(x::Int, ::One) = isone(x)
@inline Base.:(==)(::One, x::Int) = isone(x)
@inline _strides(A) = strides(A)
@inline _strides(A::DenseArray) = (One(), stride(A, 2))
@inline _strides(A::SubArray{<:Any,2,<:DenseArray,<:Tuple{Union{Base.Slice, AbstractUnitRange}, Vararg{Any}}}) = (One(), stride(A, 2))
@inline _strides(A::Union{Adjoint,Transpose}) = reverse(_strides(parent(A)))
@inline function _strides(A::PermutedDimsArray{T,N,perm}) where {T,N,perm}
    s = _strides(parent(A))
    return ntuple(d->s[perm[d]], Val(N))
end
@inline canonicalize(A) = A
@inline canonicalize(A::Union{Adjoint,Transpose,PermutedDimsArray}) = canonicalize(parent(A))

###
### Utilities
###

function checkmulsize(C, A, B)
    cm, cn = size(C)
    am, ak = size(A)
    bk, bn = size(B)
    (cm == am && ak == bk && cn == bn) || throw(DimensionMismatch("C has dimensions ($cm, $cn), A has dimensions ($am, $ak), but B has dimensions ($bk, $bn)"))
    return cm, ak, bn
end

prefetcht0(ptr::Ptr) = __prefetch(ptr, Val(:read), Val(3), Val(:data))
# From https://github.com/vchuravy/GPUifyLoops.jl/pull/5
@generated function __prefetch(ptr::T, ::Val{RW}, ::Val{Locality}, ::Val{Cache}) where {T, RW, Locality, Cache}
    decls = """
    declare void @llvm.prefetch(i8*, i32, i32, i32)
    """

    if RW == :read
        f_rw = 0
    elseif RW == :write
        f_rw = 1
    end

    f_locality = Locality

    if Cache == :data
        f_cache = 1
    elseif Cache == :instruction
        f_cache = 0
    end

    ir = """
        %ptr = inttoptr i64 %0 to i8*
        call void @llvm.prefetch(i8* %ptr, i32 $f_rw, i32 $f_locality, i32 $f_cache)
        ret void
    """

    quote
        Base.@_inline_meta
        Base.llvmcall(($decls, $ir), Nothing, Tuple{T}, ptr)
    end
end

@inline align(x::Ptr{T}, n) where {T} = reinterpret(Ptr{T}, align(reinterpret(UInt, x), n))
@inline align(x, n) = (nm1 = n - 1; (x + nm1) & -n)


text:
using Random

"""
    roll_dice(n_dice::Int=5, sides=(1,2,3,4,5,6))

Returns array of n_dice random integers corresponding to the sides of dice
"""
function roll_dice(n_dice::Int=5, sides=(1,2,3,4,5,6))
    return [rand(sides) for i in 1:n_dice]
end

"""
    yahtzee()

Play Yahtzee with 5 6-sided dice and 3 throws.
Collect as many of the same dice side as possible.
Returns a tuple with the collected side (e.g. 4's) and
how many of that side (between 1 and 5).
"""
function yahtzee()
    sides = (1,2,3,4,5,6)
    n_dice = 5
    # we first throw all dice
    first_throw = roll_dice(n_dice, sides)
    # count how many times each side comes up
    side_counts = [count(x->x==i,first_throw) for i in sides]
    # collected_side is the dice side we will start collecting
    n_collected, collected_side = findmax(side_counts)
    if n_collected == 5
        return collected_side, n_collected
    end

    # now we throw n_dice-n_collected dice and hope to get more collected_side
    second_throw = roll_dice(n_dice-n_collected, sides)
    n_new_matches = count(x->x==collected_side,second_throw)
    n_collected += n_new_matches
    if n_collected == 5
        return collected_side, n_collected
    end

    # final throw...
    third_throw = roll_dice(n_dice-n_collected, sides)
    n_new_matches = count(x->x==collected_side,third_throw)
    n_collected += n_new_matches

    return collected_side, n_collected
end

for i in 1:10
    collected_side, n_collected = yahtzee()
    println("We got $n_collected $collected_side's in this round")
    if n_collected == 5
        println("Yay, it's a Yahtzee!")
    end
end


text:
@testset "restrict" begin
    img = repeat(range(0, stop=1, length=100), 1, 100)
    tfm = recenter(RotMatrix(pi/8), center(img))

    for wv in [
        WarpedView(img, tfm),
        InvWarpedView(img, tfm)
    ]
        out = restrict(wv)
        ref = restrict(OffsetArray(collect(wv), axes(wv)))
        @test nearlysame(out, ref)
        @test out isa OffsetArray

        for dims in [1, (2, )]
            out = restrict(wv, dims)
            ref = restrict(OffsetArray(collect(wv), axes(wv)), dims)
            @test nearlysame(out, ref)
            @test out isa OffsetArray
        end
    end
end

@testset "Image resize" begin
    testcolor = (RGB,Gray)
    testtype = (Float32, Float64, N0f8, N0f16)

    @testset "Interface" begin
        function test_imresize_interface(img, outsz, args...; kargs...)
            img2 = @test_broken @inferred imresize(img, args...; kargs...) # FIXME: @inferred failed
            img2 = @test_nowarn imresize(img, args...; kargs...)
            @test size(img2) == outsz
            @test eltype(img2) == eltype(img)
        end
        for C in testcolor, T in testtype
            img = rand(C{T},10,10)

            test_imresize_interface(img, (5,5), (5,5))
            test_imresize_interface(img, (5,5), (1:5,1:5)) # FIXME: @inferred failed
            test_imresize_interface(img, (5,5), 5,5)
            test_imresize_interface(img, (5,5), 1:5,1:5) # FIXME: @inferred failed
            test_imresize_interface(img, (5,5), ratio = 0.5)
            test_imresize_interface(img, (20,20), ratio = 2)
            test_imresize_interface(img, (20,20), ratio = (2, 2))
            test_imresize_interface(img, (20,10), ratio = (2, 1))
            test_imresize_interface(img, (10,10), ())
            test_imresize_interface(img, (5,10), 5)
            test_imresize_interface(img, (5,10), (5,))
            test_imresize_interface(img, (5,10), 1:5) # FIXME: @inferred failed
            test_imresize_interface(img, (5,10), (1:5,)) # FIXME: @inferred failed

            @test_throws MethodError imresize(img,5.0,5.0)
            @test_throws MethodError imresize(img,(5.0,5.0))
            @test_throws MethodError imresize(img,(5, 5.0))
            @test_throws MethodError imresize(img,[5,5])
            @test_throws UndefKeywordError imresize(img)
            @test_throws DimensionMismatch imresize(img,(5,5,5))
            @test_throws ArgumentError imresize(img, ratio = -0.5)
            @test_throws ArgumentError imresize(img, ratio = (-0.5, 1))
            @test_throws DimensionMismatch imresize(img, ratio=(5,5,5))
            @test_throws DimensionMismatch imresize(img, (5,5,1))
        end
    end

    @testset "Numerical" begin
        # RGB is skipped
        for C in (Gray,), T in testtype
            A = [1 0; 0 1] .|> C{T}
            @test imresize(A, (1, 1)) ≈ reshape([0.5], 1, 1) .|> C{T}
            # if only 1 of 4 pixels is nonzero, check that you get 1/4 the
            # value (no matter where it is)
            for i = 1:4
                fill!(A, 0)
                A[i] = 1
                A = C{T}.(A)
                @test imresize(A, (1, 1)) ≈ reshape([0.25], 1, 1) .|> C{T}
            end
            A = [1 0; 0 0] .|> C{T}
            @test imresize(A, (2, 1)) ≈ reshape([0.5, 0], (2, 1)) .|> C{T}
            @test imresize(A, (1, 2)) ≈ reshape([0.5, 0], (1, 2)) .|> C{T}
            A = [1 0 0; 0 0 0; 0 0 0] .|> C{T}
            @test imresize(A, (2, 2)) ≈ [0.75^2 0; 0 0] .|> C{T}

            A = ones(5,5) .|> C{T}
            for l2 = 3:7, l1 = 3:7
                R = imresize(A, (l1, l2))
                @test all(x->x==1, R)
            end
            for l2 = 3:7, l1 = 3:7
                R = imresize(A, (0:l1-1, 0:l2-1))
                @test all(x->x==1, R)
                @test axes(R) == (0:l1-1, 0:l2-1)
            end
            for l1 = 3:7
                R = imresize(A, (0:l1-1,))
                @test all(x->x==1, R)
                @test axes(R) == (0:l1-1, 1:5)
            end

            R = imresize(A, (2:6, 0:4))
            @test all(x->x==1, R)
            @test axes(R) == (2:6, 0:4)
            R = imresize(A, ())
            @test all(x->x==1, R)
            @test axes(R) == axes(A)
            @test !(R === A)
            Ao = OffsetArray(A, -2:2, 0:4)
            R = imresize(Ao, (5,5))
            @test axes(R) === axes(A)
            R = imresize(Ao, axes(Ao))
            @test axes(R) === axes(Ao)
            @test !(R === A)
            img = reshape([0.5])
            R = imresize(img, ())
            @test ndims(R) == 0
            @test !(R === A)
        end

    @testset "Interpolation" begin
        img = rand(16,16)
        out = imresize(img, (128,128), method=Lanczos4OpenCV())
        @test size(out) == (128,128)

        out = imresize(img, (128,128), method=Linear())
        @test size(out) == (128,128)

        out = imresize(img, (16,16), method=Lanczos4OpenCV())
        @test size(out) == (16,16)

        #test default behavior
        @test imresize(img, (128,128)) == imresize(img, (128,128), method=Linear())

        # test that out-of-place and in-place produce the same result.
        resized = zeros(eltype(img), 128, 128)
        @test imresize(img, (128, 128)) == imresize!(resized, img)
        @test imresize(img, (128, 128), method=Linear()) == imresize!(resized, img, method=Linear())
        @test imresize(img, (128, 128), method=BSpline(Linear())) == imresize!(resized, img, method=BSpline(Linear()))
        @test imresize(img, (128, 128), method=Lanczos4OpenCV()) == imresize!(resized, img, method=Lanczos4OpenCV())

        #check error handling
        @test_throws ArgumentError imresize(img, BSpline(Linear()))
        @test_throws ArgumentError imresize(img, Linear())

        #consisency checks
        @test imresize(img, (128, 128), method=Linear()) == imresize(OffsetArray(img, -1, -1), (128, 128), method=Linear())
        @test imresize(img, (128, 128), method=BSpline(Linear())) == imresize(OffsetArray(img, -1, -1), (128, 128), method=BSpline(Linear()))
        @test imresize(img, (128, 128), method=Lanczos4OpenCV()) == imresize(OffsetArray(img, -1, -1), (128, 128), method=Lanczos4OpenCV())

        out = imresize(img, (0:127, 0:127), method=Linear())
        @test axes(out) == (0:127, 0:127)
        @test OffsetArrays.no_offset_view(out) == imresize(img, (128, 128), method=Linear())

        out = imresize(img, (0:127, 0:127), method=BSpline(Linear()))
        @test axes(out) == (0:127, 0:127)
        @test OffsetArrays.no_offset_view(out) == imresize(img, (128, 128), method=BSpline(Linear()))

        out = imresize(img, (0:127, 0:127), method=Lanczos4OpenCV())
        @test axes(out) == (0:127, 0:127)
        @test OffsetArrays.no_offset_view(out) == imresize(img, (128, 128), method=Lanczos4OpenCV())
        end
    end

    @testset "special RGB/Gray types (#97)" begin
        ori = repeat(distinguishable_colors(10), inner=(1, 10))
        for T in (
            RGB, BGR, RGBX, XRGB,
            ARGB, RGBA,
            RGB24, ARGB32,
        )
            img = T.(ori)
            out = @inferred imresize(img, ratio=2)
            @test eltype(out) <: T
            ref = imresize(ori, ratio=2)
            @test ref ≈ RGB.(out)
        end
        for T in (Gray, AGray, GrayA, Gray24)
            img = T.(ori)
            out = @inferred imresize(img, ratio=2)
            @test eltype(out) <: T
            ref = imresize(Gray.(ori), ratio=2)
            @test ref ≈ Gray.(out)
        end
    end
end


text:
using SafeTestsets

@safetestset "Normalization Test Sets" begin include("normalizationtest.jl") end


text:
struct IntegralPotential{T<:AbstractKernel,S}
    kernel::T
    surface::S
end

function IntegralPotential(op::AbstractPDE,surf,ktype::Symbol)
    if ktype == :singlelayer
        k = SingleLayerKernel(op)
    elseif ktype == :doublelayer
        k = DoubleLayerKernel(op)
    else
        error("unknown kernel type: valid options are :singlelayer or :doublelayer")
    end
    return IntegralPotential(k,surf)
end

kernel_type(pot::IntegralPotential) = kernel_type(pot.kernel)

(pot::IntegralPotential)(σ::AbstractVector,x) = pot(kernel_type(pot),σ,x)
function (pot::IntegralPotential)(::SingleLayer,σ::AbstractVector,x)
    f = pot.kernel
    Γ = pot.surface
    iter = zip(qnodes(Γ),qweights(Γ),σ)
    out = mapreduce(+,iter) do (y,w,σ)
        f(x,y)*σ*w
    end
    return out
end
function (pot::IntegralPotential)(::DoubleLayer,σ::AbstractVector,x)
    f = pot.kernel
    Y = pot.surface
    iter = zip(qnodes(Y),qweights(Y),qnormals(Y),σ)
    out  = mapreduce(+,iter) do (y,w,ny,σ)
        f(x,y,ny)*σ*w
    end
    return out
end

Base.getindex(pot::IntegralPotential,σ::AbstractVector) = (x) -> pot(σ,x)

SingleLayerPotential(op::AbstractPDE,surf) = IntegralPotential(SingleLayerKernel(op),surf)
DoubleLayerPotential(op::AbstractPDE,surf) = IntegralPotential(DoubleLayerKernel(op),surf)

max_stars_repo_name: WaveProp/WaveProp

max_issues_repo_name: WaveProp/WaveProp

max_forks_repo_name: WaveProp/WaveProp


text:
@testset "getindex" begin
    headers = Headers("Content-Type" => "application/json")
    @test headers["Content-Type"] == "application/json"
    @test headers["content-type"] == "application/json"
    @test headers["CoNtEnT-tYpE"] == "application/json"

    headers = Headers("CoNtEnT-tYpE" => "text/csv")

    @test headers["Content-Type"] == "text/csv"
    @test headers["content-type"] == "text/csv"
    @test headers["CoNtEnT-tYpE"] == "text/csv"

    @test_throws BoundsError headers["Bad-Header"]
end

@testset "get" begin
    headers = Headers("Content-Type" => "application/json")

    @test get(headers, "Content-Type", "default") == "application/json"
    @test get(headers, "content-type", "default") == "application/json"
    @test get(headers, "Fake-Header", "default") == "default"
end

@testset "haskey" begin
    headers = Headers("Content-Type" => "application/json")

    @test haskey(headers, "Content-Type") == true
    @test haskey(headers, "content-TYPE") == true
    @test haskey(headers, "Fake-Header") == false
end

@testset "set" begin
    headers = Headers("Content-Type" => "application/json")

    modified = set(headers, "X-Powered-By", "Spirit")

    @test haskey(modified, "Content-Type") == true
    @test haskey(modified, "X-Powered-By") == true
    @test haskey(headers, "X-Powered-By") == false
    @test modified["X-Powered-By"] == "Spirit"
    @test modified["Content-Type"] == "application/json"

    modified = set(headers, "content-type", "application/yaml")

    @test headers["content-type"] == "application/json"
    @test modified["Content-Type"] == "application/yaml"
end

@testset "append" begin
    headers = Headers()

    headers = append(headers, "Accept", "application/json; q=0.6")
    @test headers["accept"] == "application/json; q=0.6"

    headers = append(headers, "aCcEpT", "application/yaml; q=0.8")

    @test headers["accept"] == "application/json; q=0.6, application/yaml; q=0.8"
end

@testset "delete" begin
    headers = Headers("Content-Type" => "application/json")

    modified = delete(headers, "content-type")

    @test haskey(headers, "content-type") == true
    @test haskey(modified, "content-type") == false
end


text:
include("../utils/Utils.jl")
using .Utils

function align(read::String, cigar::String, md::String)

    ref_seq = ""
    newread = ""

    # # if empty md tag
    # if md == ""
    #     md = string(length(read))
    # end

    MDlist = collect(m.match for m in eachmatch(r"(\d+|\D+)", md))

    MDcounter = 1
    alignment = ""

    for e in MDlist

        if is_integer(e)
            i = tryparse(Int, e)
            alignment *= "."^i
            ref_seq *= read[MDcounter:MDcounter+i-1]
            newread *= read[MDcounter:MDcounter+i-1]
            MDcounter += i

        elseif occursin("^", e)
            ef = lstrip(e, '^')
            alignment *= ef
            continue

        elseif is_alpha(e)
            alignment *= e
            ref_seq *= e
            newread *= read[MDcounter]
            MDcounter += length(e)
        end

    end

    if occursin("I", cigar) | occursin("S", cigar)

        # find insertions and clips in cigar
        insertions = Int[]
        softclips = Int[]
        cigarcount = 1

        cigarcomponents = collect(c.captures for c in eachmatch(r"([0-9]*)([A-Z])", cigar))

        for p in cigarcomponents

            cigaraddition = tryparse(Int, p[1])

            if occursin("I", p[2])
                for c = cigarcount:cigarcount+cigaraddition-1
                    push!(insertions, c)
                end

            elseif occursin("S", p[2])
                for c = cigarcount:cigarcount+cigaraddition-1
                    push!(softclips, c)
                end
            end

            cigarcount += cigaraddition

        end
        # end cigar parsing

        # redo the read and ref using indel and clip info
        ref_seq = ""
        # newread = ""
        alignmentcounter = 1

        for (x, r) in zip(1:length(read), read)

            if x in insertions
                ref_seq *= "-"
                # newread *= read[x]

            elseif x in softclips
                ref_seq *= "-"
                # newread *= read[x]

            else
                # newread *= read[x]
                refbasealn = alignment[alignmentcounter]
                if refbasealn == '.'
                    ref_seq *= read[x]
                else
                    ref_seq *= refbasealn
                end
                alignmentcounter += 1

            end
        end

    end

    return read, ref_seq

end

function align(read::LongDNASeq, cigar::String, md::String)
    return align(convert(String, read), cigar, md)
end

max_stars_repo_name: ChristianMichelsen/reference-free

max_issues_repo_name: ChristianMichelsen/reference-free

max_forks_repo_name: ChristianMichelsen/reference-free


text:
using Documenter, DemoCards
using DIPTutorials

templates, theme = cardtheme("list")
tutorials, tutorials_cb = makedemos("tutorials", templates)
dip_examples, dip_examples_cb = makedemos("dip", templates)

format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true",
                         assets = [theme])

makedocs(format = format,
         modules = [DIPTutorials],
         sitename = "Julia 数字图像处理",
         pages = [
            "Home" => "index.md",
            "JuliaImages 入门" => tutorials,
            "DIP 示例代码" => dip_examples,
            "Reference" => "reference.md"
         ]
)

tutorials_cb()
dip_examples_cb()

deploydocs(repo="github.com/johnnychen94/Image-Processing-in-Julia.git")


text:
module TraceTest

using Cassette
using Test
using SemanticModels.Dubstep
using LinearAlgebra

struct ODEProblem{T,U,V,W}
    ode::T
    init::U
    tspan::V
    parms::W
end

# Define == for ODEProblems to mean they are semantically equivalent. The
# normal definition for structs is that equality of fields is done with ===
# which means structs that contain arrays, cannot be == if those arrays are
# not === perhaps we should use the AutoHashEquals.jl package here.
function Base.:(==)(p::ODEProblem, q::ODEProblem)
    return p.ode == q.ode && p.init == q.init && p.tspan == q.tspan && p.parms == q.parms
end

"""   sir_ode2(du,u,p,t)

        computes the du/dt array for the SIR system. parameters p is b,g = beta,gamma.
        """
sir_ode2(du,u,p,t) = begin
    S,I,R = u
    b,g = p
    du[1] = -b*S*I
    du[2] = b*S*I-g*I
    du[3] = g*I
end

# Cassette.@context TraceCtx

# function Cassette.overdub(ctx::TraceCtx, args...)
#     subtrace = Any[]
#     push!(ctx.metadata, args => subtrace)
#     if Cassette.canrecurse(ctx, args...)
#         newctx = Cassette.similarcontext(ctx, metadata = subtrace)
#         return Cassette.recurse(newctx, args...)
#     else
#         return Cassette.fallback(ctx, args...)
#     end
# end

trace = Any[]
x, y, z = rand(3)
f(x, y, z) = x*y + y*z
Cassette.recurse(Dubstep.TraceCtx(metadata = trace), () -> f(x, y, z))
@testset "Cassette" begin
        @test trace == Any[
            (f,x,y,z) => Any[
                (*,x,y) => Any[(Base.mul_float,x,y)=>Any[]]
                (*,y,z) => Any[(Base.mul_float,y,z)=>Any[]]
                (+,x*y,y*z) => Any[(Base.add_float,x*y,y*z)=>Any[]]
            ]
        ]
end

function g()
    parms = [0.1,0.05]
    init = [0.99,0.01,0.0]
    tspan = (0.0,200.0)
    sir_prob2 = construct(ODEProblem,sir_ode2,init,tspan,parms)
    # sir_sol = solve(sir_prob2,saveat = 0.1)
    return sir_prob2
end

function h()
    parms = [0.1,0.05]
    init = [0.99,0.01,0.0]
    tspan = (0.0,200.0)
    sir_prob2 = ODEProblem(sir_ode2,init,tspan,parms)
    # sir_sol = solve(sir_prob2,saveat = 0.1)
    return sir_prob2
end

# these are equal, because we wrote a custom == function for ODEProblem.
@testset "Dubstep" begin
    @testset "SIR" begin

    @info "Tracing implementation with construct call"
    trace1 = Dubstep.trace(g)
    @info "Tracing implementation without construct call"
    trace2 = Dubstep.trace(h)
    @test trace2.value == trace1.value

    for i in 1:3
        @test trace1.trace[i] == trace2.trace[i]
    end

    @info "Trace step for construct call"
    @show trace1.trace[4]
    @test trace1.trace[4][1] == construct
    @info "Trace step for direct apply type"
    @show trace2.trace[4]
    @test last(trace2.trace[4])[1][1] == Core.apply_type

end #SIR

end #Dubstep

subg(x,y) = norm([x x x]/6 - [y y y]/2, 2)
function g()
    a = 5+7
    b = 3+4
    c = subg(a,b)
    return c
end
function h()
    a = 5+7
    c = subg(a,(3+4))
    return c
end

@testset "LP" begin
@test 2.5980 < g() < 2.599
ctx = Dubstep.LPCtx(metadata=Dict(1=>2, 2=>1, Inf=>1))
@test Cassette.recurse(ctx, g) == 4.5
end #LP

end #module


text:
using MLDatasets
using Flux
using Flux: Data.DataLoader
using Flux: onehotbatch, onecold, @epochs
using Statistics: mean

x_train, y_train = MNIST.traindata()
x_test, y_test = MNIST.testdata()
# reshape image to 1D array
x_train = reshape(x_train, (28 * 28, :))
x_test = reshape(y_train, (28 * 28, :))
# one hot batch
y_train = onehotbatch(y_train, 0:9)
y_test = onehotbatch(y_test, 0:9)
data = DataLoader((x_train, y_train), batchsize=40, shuffle=true)

model = Chain(
  Dense(28 * 28, 40, relu),
  Dense(40, 10),
  softmax
)
model_with_batch_normalization = Chain(
  Dense(28 * 28, 40, relu),
  BatchNorm(40, relu),
  Dense(40, 10),
  BatchNorm(10),
  softmax
)

loss(x, y) = Flux.crossentropy(model(x), y)
acc(x, y) = mean(Flux.onecold(model(x)) .== Flux.onecold(y))
# opt = Descent()
# opt = ADADelta()
# opt = Momentum()
opt = ADAM()

@epochs 10 begin
  Flux.train!(loss, params(model), data, opt)
  @show(loss(x_test, y_test))
  @show(acc(x_test, y_test))
end


text:
module ApplyTests

using Test
using InplaceArrays.Arrays
using InplaceArrays.Fields
using FillArrays
using InplaceArrays.TensorValues

using InplaceArrays.Fields: MockField, MockBasis

k = bcast(-)

v = 3.0
d = 2
f = MockField{d}(v)
g = apply_kernel_to_field(k,f)

p = 4
p1 = Point(1,2)
p2 = Point(2,1)
p3 = Point(4,3)
p4 = Point(6,1)
x = [p1,p2,p3,p4]

fx = evaluate(f,x)
∇fx = evaluate(∇(f),x)
gx = apply_kernel(k,fx)
∇gx = apply_kernel(k,∇fx)
test_field(g,x,gx,grad=∇gx)

fi = 3.0
gi = 4.5
d = 2
f = MockField{d}(fi)
g = MockField{d}(gi)
h = apply_kernel_to_field(k,f,g)
hx = apply_kernel(k,evaluate(f,x),evaluate(g,x))
∇hx = apply_kernel(k,evaluate(∇(f),x),evaluate(∇(g),x))
test_field(h,x,hx,grad=∇hx)

fi = 3.0
gi = VectorValue(4,5)
d = 2
f = MockField{d}(fi)
g = gi
h = apply_kernel_to_field(k,f,g)
hx = apply_kernel(k,evaluate(f,x),evaluate_field(g,x))
∇hx = apply_kernel(k,evaluate(∇(f),x),evaluate_field(field_gradient(g),x))
test_field(h,x,hx,grad=∇hx)

fi = 3.0
gi = [1,2,3,4,5]
d = 2
f = MockField{d}(fi)
g = gi
h = apply_kernel_to_field(k,f,g)
hx = apply_kernel(k,evaluate(f,x),evaluate_field(g,x))
∇hx = apply_kernel(k,evaluate(∇(f),x),evaluate_field(field_gradient(g),x))
test_field(h,x,hx,grad=∇hx)

end # module


text:
# module classical

export Particles, CavityMode, MultimodeSystem,
        ClassicalState, x, v, αn, splitstate,
        orderparameter, orderparameters

using ArrayViews

# Define system

type Particles
    N::Int
    mass::Float64
end

type CavityMode
    k::Float64
    delta::Float64
    eta::Float64
    kappa::Float64
    U0::Float64
end

type MultimodeSystem
    particles::Particles
    modes::Vector{CavityMode}
end

# Define state of system

type ClassicalState
    particlenumber::Int
    cavitymodenumber::Int
    data::Vector{Float64}
    function ClassicalState(particlenumber::Int, cavitymodenumber::Int, data::Vector{Float64})
        @assert particlenumber > 0
        @assert cavitymodenumber > 0
        @assert length(data)==2*(particlenumber+cavitymodenumber)
        x = new()
        x.particlenumber = particlenumber
        x.cavitymodenumber = cavitymodenumber
        x.data = data
        return x
    end
end

ClassicalState(particlenumber::Int, cavitymodenumber::Int) = ClassicalState(particlenumber::Int, cavitymodenumber::Int, zeros(Float64, 2*(particlenumber + cavitymodenumber)))
ClassicalState(s::MultimodeSystem, data::Vector{Float64}) = ClassicalState(s.particles.N, length(s.modes), data)
ClassicalState(s::MultimodeSystem) = ClassicalState(s.particles.N, length(s.modes))

x(state::ClassicalState) = view(state.data, 1:state.particlenumber)
v(state::ClassicalState) = view(state.data, state.particlenumber+1:2*state.particlenumber)
αn(state::ClassicalState, n::Int) = view(state.data, 2*state.particlenumber+2*(n-1)+1:2*state.particlenumber+2*n)

splitstate(N::Int, M::Int, data::Vector{Float64}) = view(data, 1:N), view(data, N+1:2*N), view(data, 2*N+1:2*(N+M))

orderparameter(k::Float64, state::ClassicalState) = 1./state.particlenumber*sum(sin(k*x(state)))
orderparameters(K::Vector{Float64}, state::ClassicalState) = Float64[orderparameter(k, state) for k=K]
orderparameters(s::MultimodeSystem, state::ClassicalState) = orderparameters(Float64[mode.k for mode=s.modes], state)

# end # module