text:
################################################################################
#                          PARAMETER FUNCTION METHODS
################################################################################
# Extend dispatch_variable_ref
function dispatch_variable_ref(
    model::InfiniteModel,
    index::ParameterFunctionIndex
    )::ParameterFunctionRef
    return ParameterFunctionRef(model, index)
end

# Extend _add_data_object
function _add_data_object(
    model::InfiniteModel,
    object::ParameterFunctionData
    )::ParameterFunctionIndex
    return MOIUC.add_item(model.param_functions, object)
end

# Extend _data_dictionary (reference based)
function _data_dictionary(
    fref::ParameterFunctionRef
    )::MOIUC.CleverDict{ParameterFunctionIndex, ParameterFunctionData{<:ParameterFunction}}
    return JuMP.owner_model(fref).param_functions
end

# Extend _data_object
function _data_object(fref::ParameterFunctionRef)
  object = get(_data_dictionary(fref), JuMP.index(fref), nothing)
  object === nothing && error("Invalid parameter function reference, cannot find " *
                              "corresponding object in the model. This is likely " *
                              "caused by using the reference of a deleted function.")
  return object
end

# Extend _core_variable_object
function _core_variable_object(fref::ParameterFunctionRef)
    return _data_object(fref).func
end

# Extend _object_numbers
function _object_numbers(fref::ParameterFunctionRef)::Vector{Int}
    return _core_variable_object(fref).object_nums
end

# Extend _parameter_numbers
function _parameter_numbers(fref::ParameterFunctionRef)::Vector{Int}
    return _core_variable_object(fref).parameter_nums
end

"""
    build_parameter_function(
        _error::Function,
        func::Function,
        parameter_refs::Union{GeneralVariableRef, AbstractArray{<:GeneralVariableRef}, Tuple}
        )::ParameterFunction

Build an [`ParameterFunction`](@ref) object that employs a parameter function
`func` that takes instances of the infinite parameter(s) as input. This can
ultimately by incorporated into expressions to enable nonlinear infinite parameter
behavior and/or incorporate data over infinite domains.

Here `func` should be of the form:
```
func(paramvals...)::Float64
```
where the formatting of `paramvals` is analagous to point variables (and will be
based on the tuple of infinite parameter references given in `parameter_refs`).

Errors if the infinite parameter tuple is formatted incorrectly. The allowed
format follows that of infinite variables. Also errors if the function doesn't
accept a support realization of the `parameter_refs` as input.

**Example**
```julia-repl
julia> f = build_parameter_function(error, sin, t);
```
"""
function build_parameter_function(
    _error::Function,
    func::Function,
    parameter_refs::Union{GeneralVariableRef, AbstractArray{<:GeneralVariableRef}, Tuple};
    extra_kwargs...
    )
    # check for unneeded keywords
    for (kwarg, _) in extra_kwargs
        _error("Keyword argument $kwarg is not for use with parameter functions.")
    end
    # process the parameter reference inputs and check
    prefs = Collections.VectorTuple(parameter_refs)
    _check_parameter_tuple(_error, prefs)
    # check the function
    supp_type = typeof(Tuple(Vector{Float64}(undef, length(prefs)), prefs))
    if !hasmethod(func, supp_type)
        _error("Parameter function method `$(func)` is not defined for `$(func)(",
               join(Tuple(prefs), ", ") * ")`. Note that the infinite parameter ",
               "arguments `" * join(Tuple(prefs), ", ") * "` are checked via a ",
               "numeric support (each parameter is a `Float64`) of the same format.")
    end
    # get the parameter object numbers
    object_nums = Int[]
    for pref in prefs
        union!(object_nums, _object_number(pref))
    end
    # get the parameter numbers
    param_nums = [_parameter_number(pref) for pref in prefs]
    # make the variable and return
    return ParameterFunction(func, prefs, object_nums, param_nums)
end

# Fallback for weird macro inputs
function build_parameter_function(_error::Function, func, prefs; kwargs...)
    _error("Invalid syntax causing unexpected parsing of the function and " *
           "infinite parameter arguments.")
end

# Used to update parameter-parameter function mappings
function _update_param_var_mapping(
    fref::ParameterFunctionRef,
    prefs::Collections.VectorTuple
    )::Nothing
    for pref in prefs
        dependency_list = _parameter_function_dependencies(pref)
        if !(JuMP.index(fref) in dependency_list)
            push!(dependency_list, JuMP.index(fref))
        end
    end
    return
end

"""
    add_parameter_function(model::InfiniteModel, pfunc::ParameterFunction,
                           [name::String])::GeneralVariableRef

Add an [`ParameterFunction`](@ref) `pfunc` to the `model` using `name` for
printing and return a `GeneralVariableRef` such that it can be embedded in
expressions. Errors if the parameter function `pfunc` points to do not belong to
`model`. Note that `pfunc` should be created using [`build_parameter_function`](@ref).

**Example**
```julia-repl
julia> f = build_parameter_function(error, sin, t);

julia> fref = add_parameter_function(model, f)
sin(t)
```
"""
function add_parameter_function(
    model::InfiniteModel,
    pfunc::ParameterFunction,
    name::String = String(nameof(pfunc.func))
    )::GeneralVariableRef
    _check_parameters_valid(model, pfunc.parameter_refs)
    data_object = ParameterFunctionData(pfunc, name)
    findex = _add_data_object(model, data_object)
    fref = ParameterFunctionRef(model, findex)
    _update_param_var_mapping(fref, pfunc.parameter_refs)
    return _make_variable_ref(model, findex)
end

"""
    JuMP.name(fref::ParameterFunctionRef)::String

Extend `JuMP.name` to return the base name of `fref`.

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

"""
    JuMP.set_name(fref::ParameterFunctionRef, name::String)::Nothing

Extend `JuMP.set_name` to set the name of a parameter function.

**Example**
```julia-repl
julia> set_name(fref, "func_name")

julia> name(fref)
"func_name"
```
"""
function JuMP.set_name(fref::ParameterFunctionRef, name::String)::Nothing
    _data_object(fref).name = name
    return
end

"""
    parameter_function(func::Function,
                       pref_inputs::Union{GeneralVariableRef, AbstractArray{GeneralVariableRef}, Tuple};
                       [name::String = [the name of `func`]]
                       )::GeneralVariableRef

Make a parameter function and return a `GeneralVariableRef` that can be
embedded in InfiniteOpt expressions. This serves as a convenient wrapper for
[`build_parameter_function`](@ref) and [`add_parameter_function`](@ref). For an
even more convenient definition method see [`@parameter_function`](@ref).

Here `func` denotes the function that will take a support of infinite parameters as
input (formatted like `pref_inputs`) and will return a scalar value. Specifically,
`func` should be of the form:
```
func(paramvals...)::Float64
```
where the formatting of `paramvals` is analagous to point variables (and will be
based on the tuple of infinite parameter references given in `parameter_refs`).
Moreover, `func` must be a function that returns a scalar numeric value.

Errors if `func` will not take a support formatted like `pref_inputs` in
combination with the `fargs` and `fkwargs` specified. Also errors if `pref_inputs`
follow an invalid input format.

**Example**
```julia-repl
julia> p_func = parameter_function(sin, t)
sin(t)

julia> p_func3 = parameter_function((t_supp) -> 2 * sin(2 * t_supp), t, name = "mysin")
mysin(t)

julia> p_func4 = parameter_function(t, name = "mysin") do t_supp
                    if t_supp <= 5
                        return sin(t_supp)
                    else
                        return 2 * sin(2 * t_supp)
                    end
                 end

mysin(t)
```
"""
function parameter_function(
    func::Function,
    pref_inputs;
    name::String = String(nameof(func))
    )::GeneralVariableRef
    f = build_parameter_function(error, func, pref_inputs)
    model = JuMP.owner_model(first(f.parameter_refs))
    return add_parameter_function(model, f, name)
end

"""
    raw_parameter_refs(fref::ParameterFunctionRef)::VectorTuple

Return the raw [`VectorTuple`](@ref InfiniteOpt.Collections.VectorTuple) of the
parameter references that `fref` depends on. This is primarily an internal method
where [`parameter_refs`](@ref parameter_refs(::ParameterFunctionRef))
is intended as the preferred user function.
"""
function raw_parameter_refs(fref::ParameterFunctionRef)
    return _core_variable_object(fref).parameter_refs
end

"""
    parameter_refs(fref::ParameterFunctionRef)::Tuple

Return the parameter references associated with `fref`. This
is formatted as a Tuple of containing the parameter references as they inputted
to define `fref`.

**Example**
```julia-repl
julia> parameter_refs(p_func)
(t,)
```
"""
function parameter_refs(fref::ParameterFunctionRef)
    return Tuple(raw_parameter_refs(fref))
end

"""
    parameter_list(fref::ParameterFunctionRef)::Vector{GeneralVariableRef}

Return a vector of the parameter references that `fref` depends on. This is
primarily an internal method where [`parameter_refs`](@ref parameter_refs(::ParameterFunctionRef))
is intended as the preferred user function.
"""
function parameter_list(fref::ParameterFunctionRef)::Vector{GeneralVariableRef}
    return raw_parameter_refs(fref).values
end

"""
    raw_function(fref::ParameterFunctionRef)::Function

Returns the raw function behind `fref` that takes a particular support of `fref`'s
infinite parameters as input.
"""
function raw_function(fref::ParameterFunctionRef)
    return _core_variable_object(fref).func
end

"""
    call_function(fref::ParameterFunctionRef, support...)::Float64

Safely evaluates the [`raw_function`](@ref) of `fref` at a particular support `support`
point that matches the format of the infinite parameter tuple given when the `fref`
was defined. This is essentially equivalent to `raw_function(fref)(supps...)`.
"""
function call_function(fref::ParameterFunctionRef, supps...)::Float64
    pfunc = _core_variable_object(fref)
    return pfunc.func(supps...)
end

# Extend _semi_infinite_variable_dependencies
function _semi_infinite_variable_dependencies(
    fref::ParameterFunctionRef
     )::Vector{SemiInfiniteVariableIndex}
    return _data_object(fref).semi_infinite_var_indices
end

# Extend _derivative_dependencies
function _derivative_dependencies(
    fref::ParameterFunctionRef
     )::Vector{DerivativeIndex}
    return _data_object(fref).derivative_indices
end

# Extend _measure_dependencies
function _measure_dependencies(fref::ParameterFunctionRef)::Vector{MeasureIndex}
    return _data_object(fref).measure_indices
end

# Extend _constraint_dependencies
function _constraint_dependencies(
    fref::ParameterFunctionRef
    )::Vector{InfOptConstraintIndex}
    return _data_object(fref).constraint_indices
end

"""
    used_by_semi_infinite_variable(fref::ParameterFunctionRef)::Bool

Return a `Bool` indicating if `fref` is used by a semi-infinite infinite variable.

**Example**
```julia-repl
julia> used_by_semi_infinite_variable(fref)
false
```
"""
function used_by_semi_infinite_variable(fref::ParameterFunctionRef)::Bool
    return !isempty(_semi_infinite_variable_dependencies(fref))
end

"""
    used_by_derivative(fref::ParameterFunctionRef)::Bool

Return a `Bool` indicating if `fref` is used by a derivative.

**Example**
```julia-repl
julia> used_by_derivative(vref)
true
```
"""
function used_by_derivative(fref::ParameterFunctionRef)::Bool
    return !isempty(_derivative_dependencies(fref))
end

"""
    used_by_measure(fref::ParameterFunctionRef)::Bool

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

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

"""
    used_by_constraint(fref::ParameterFunctionRef)::Bool

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

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

"""
    is_used(fref::ParameterFunctionRef)::Bool

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

**Example**
```julia-repl
julia> is_used(fref)
true
```
"""
function is_used(fref::ParameterFunctionRef)::Bool
    return used_by_measure(fref) || used_by_constraint(fref) ||
           used_by_semi_infinite_variable(fref) || used_by_derivative(fref)
end

"""
    JuMP.delete(model::InfiniteModel, fref::ParameterFunctionRef)::Nothing

Extend `JuMP.delete` to delete parameter functions and their dependencies. Errors
if `fref` is invalid, meaning it has already been deleted or it belongs to
another model.
"""
function JuMP.delete(model::InfiniteModel, fref::ParameterFunctionRef)::Nothing
    @assert JuMP.is_valid(model, fref) "Parameter function is invalid."
    # update the optimizer model status
    if is_used(fref)
        set_optimizer_model_ready(model, false)
    end
    # update parameter mapping
    all_prefs = parameter_list(fref)
    for pref in all_prefs
        filter!(e -> e != JuMP.index(fref), _parameter_function_dependencies(pref))
    end
    gvref = _make_variable_ref(model, JuMP.index(fref))
    # delete associated semi-infinite variables and mapping
    for index in _semi_infinite_variable_dependencies(fref)
        JuMP.delete(model, dispatch_variable_ref(model, index))
    end
    # delete associated derivative variables and mapping
    for index in _derivative_dependencies(fref)
        JuMP.delete(model, dispatch_variable_ref(model, index))
    end
    # remove from measures if used
    for mindex in _measure_dependencies(fref)
        mref = dispatch_variable_ref(model, mindex)
        func = measure_function(mref)
        data = measure_data(mref)
        if func isa GeneralVariableRef
            new_func = zero(JuMP.GenericAffExpr{Float64, GeneralVariableRef})
            new_meas = Measure(new_func, data, Int[], Int[], true)
        else
            _remove_variable(func, gvref)
            new_meas = build_measure(func, data)
        end
        _set_core_variable_object(mref, new_meas)
    end
    # remove from constraints if used
    for cindex in _constraint_dependencies(fref)
        cref = _make_constraint_ref(model, cindex)
        func = JuMP.jump_function(JuMP.constraint_object(cref))
        if func isa GeneralVariableRef
            set = JuMP.moi_set(JuMP.constraint_object(cref))
            new_func = zero(JuMP.GenericAffExpr{Float64, GeneralVariableRef})
            new_constr = JuMP.ScalarConstraint(new_func, set)
            _set_core_constraint_object(cref, new_constr)
            empty!(_object_numbers(cref))
        elseif func isa AbstractArray{GeneralVariableRef}
            JuMP.delete(model, cref)
        else
            _remove_variable(func, gvref)
            _data_object(cref).object_nums = sort(_object_numbers(func))
        end
    end
    # delete the data object
    _delete_data_object(fref)
    return
end

################################################################################
#                              COMPARISON METHODS
################################################################################
## Extend for better comparisons than default
# GenericAffExpr
function Base.:(==)(aff1::JuMP.GenericAffExpr{C, V},
    aff2::JuMP.GenericAffExpr{C, V}
    )::Bool where {C, V <: GeneralVariableRef}
    return aff1.constant == aff2.constant && aff1.terms == aff2.terms
end

# GenericQuadExpr
function Base.:(==)(quad1::JuMP.GenericQuadExpr{C, V},
    quad2::JuMP.GenericQuadExpr{C, V}
    )::Bool where {C, V <: GeneralVariableRef}
    pairs1 = collect(pairs(quad1.terms))
    pairs2 = collect(pairs(quad2.terms))
    if length(pairs1) != length(pairs2)
        return false
    end
    for i in eachindex(pairs1)
        if pairs1[i][1].a != pairs2[i][1].a || pairs1[i][1].b != pairs2[i][1].b ||
           pairs1[i][2] != pairs2[i][2]
            return false
        end
    end
    return quad1.aff == quad2.aff
end

################################################################################
#                            VARIABLE LIST MAKING
################################################################################
## Determine which variables are present in a function
# GeneralVariableRef
function _all_function_variables(
    f::GeneralVariableRef
    )::Vector{GeneralVariableRef}
    return [f]
end

# GenericAffExpr
function _all_function_variables(
    f::JuMP.GenericAffExpr{C, V}
    )::Vector{V} where {C, V <: GeneralVariableRef}
    return collect(keys(f.terms))
end

# GenericQuadExpr
function _all_function_variables(
    f::JuMP.GenericQuadExpr{C, V}
    )::Vector{V} where {C, V <: GeneralVariableRef}
    vref_set = Set(keys(f.aff.terms))
    for pair in keys(f.terms)
        push!(vref_set, pair.a)
        push!(vref_set, pair.b)
    end
    return collect(vref_set)
end

# Array of expressions
function _all_function_variables(
    arr::AbstractArray
    )::Vector{GeneralVariableRef}
    vref_set = Set{GeneralVariableRef}()
    for f in arr
        union!(vref_set, _all_function_variables(f))
    end
    return collect(vref_set)
end

# Fallback
function _all_function_variables(f)
    error("`_all_function_variables` not defined for expression of type $(typeof(f)).")
end

################################################################################
#                            OBJECT NUMBER METHODS
################################################################################
## Return the unique set of object numbers in an expression
# Dispatch fallback (--> should be defined for each non-empty variable type)
_object_numbers(expr::DispatchVariableRef)::Vector{Int} = Int[]

# GeneralVariableRef
function _object_numbers(expr::GeneralVariableRef)::Vector{Int}
    return _object_numbers(dispatch_variable_ref(expr))
end

# GenericAffExpr
function _object_numbers(expr::JuMP.GenericAffExpr)::Vector{Int}
    obj_nums = Set{Int}()
    for vref in keys(expr.terms)
        union!(obj_nums, _object_numbers(vref))
    end
    return collect(obj_nums)
end

# GenericQuadExpr
function _object_numbers(expr::JuMP.GenericQuadExpr)::Vector{Int}
    obj_nums = Set{Int}()
    for vref in keys(expr.aff.terms)
        union!(obj_nums, _object_numbers(vref))
    end
    for pair in keys(expr.terms)
        union!(obj_nums, _object_numbers(pair.a))
        union!(obj_nums, _object_numbers(pair.b))
    end
    return collect(obj_nums)
end

# Variable list
function _object_numbers(vrefs::Vector{GeneralVariableRef})::Vector{Int}
    obj_nums = Set{Int}()
    for vref in vrefs
        union!(obj_nums, _object_numbers(vref))
    end
    return collect(obj_nums)
end

################################################################################
#                             PARAMETER NUMBER METHODS
################################################################################
## Return the unique set of parameter numbers in an expression
# Dispatch fallback (--> should be defined for each non-empty variable type)
_parameter_numbers(expr::DispatchVariableRef)::Vector{Int} = Int[]

# GeneralVariableRef
function _parameter_numbers(expr::GeneralVariableRef)::Vector{Int}
    return _parameter_numbers(dispatch_variable_ref(expr))
end

# GenericAffExpr
function _parameter_numbers(expr::JuMP.GenericAffExpr)::Vector{Int}
    par_nums = Set{Int}()
    for vref in keys(expr.terms)
        union!(par_nums, _parameter_numbers(vref))
    end
    return collect(par_nums)
end

# GenericQuadExpr
function _parameter_numbers(expr::JuMP.GenericQuadExpr)::Vector{Int}
    par_nums = Set{Int}()
    for vref in keys(expr.aff.terms)
        union!(par_nums, _parameter_numbers(vref))
    end
    for pair in keys(expr.terms)
        union!(par_nums, _parameter_numbers(pair.a))
        union!(par_nums, _parameter_numbers(pair.b))
    end
    return collect(par_nums)
end

################################################################################
#                             MODEL EXTRACTION METHODS
################################################################################
## Get the model from an expression
# GeneralVariableRef
function _model_from_expr(expr::GeneralVariableRef)::InfiniteModel
    return JuMP.owner_model(expr)
end

# AffExpr
function _model_from_expr(
    expr::JuMP.GenericAffExpr
    )::Union{InfiniteModel, Nothing}
    if isempty(expr.terms)
        return
    else
        return JuMP.owner_model(first(keys(expr.terms)))
    end
end

# QuadExpr
function _model_from_expr(
    expr::JuMP.GenericQuadExpr
    )::Union{InfiniteModel, Nothing}
    result = _model_from_expr(expr.aff)
    if result !== nothing
        return result
    elseif isempty(expr.terms)
        return
    else
        return JuMP.owner_model(first(keys(expr.terms)).a)
    end
end

# Vector{GeneralVariableRef}
function _model_from_expr(
    vrefs::Vector{GeneralVariableRef}
    )::Union{InfiniteModel, Nothing}
    if isempty(vrefs)
        return
    else
        return JuMP.owner_model(first(vrefs))
    end
end

# Fallback
function _model_from_expr(expr)
    error("`_model_from_expr` not defined for expr of type $(typeof(expr)).")
end

################################################################################
#                            VARIABLE REMOVAL BOUNDS
################################################################################
## Delete variables from an expression
# GenericAffExpr
function _remove_variable(
    f::JuMP.GenericAffExpr,
    vref::GeneralVariableRef
    )::Nothing
    if haskey(f.terms, vref)
        delete!(f.terms, vref)
    end
    return
end

# GenericQuadExpr
function _remove_variable(
    f::JuMP.GenericQuadExpr,
    vref::GeneralVariableRef
    )::Nothing
    _remove_variable(f.aff, vref)
    for (pair, coef) in f.terms
        if pair.a == vref || pair.b == vref
            delete!(f.terms, pair)
        end
    end
    return
end

# AbstractArray
function _remove_variable(
    arr::AbstractArray,
    vref::GeneralVariableRef
    )::Nothing
    for f in arr
        _remove_variable(f, vref)
    end
    return
end

################################################################################
#                             COEFFICIENT METHODS
################################################################################
## Modify linear coefficient of variable in expression
# GeneralVariableRef
function _set_variable_coefficient!(expr::GeneralVariableRef,
    var::GeneralVariableRef,
    coeff::Real
    )::JuMP.GenericAffExpr{Float64, GeneralVariableRef}
    # Determine if variable is that of the expression and change accordingly
    if expr == var
        return Float64(coeff) * var
    else
        return expr + Float64(coeff) * var
    end
end

# GenericAffExpr
function _set_variable_coefficient!(expr::JuMP.GenericAffExpr{C, V},
    var::V,
    coeff::Real
    )::JuMP.GenericAffExpr{C, V} where {C, V <: GeneralVariableRef}
    # Determine if variable is in the expression and change accordingly
    if haskey(expr.terms, var)
        expr.terms[var] = coeff
        return expr
    else
        return expr + coeff * var
    end
end

# GenericQuadExpr
function _set_variable_coefficient!(expr::JuMP.GenericQuadExpr{C, V},
    var::V,
    coeff::Real
    )::JuMP.GenericQuadExpr{C, V} where {C, V <: GeneralVariableRef}
    # Determine if variable is in the expression and change accordingly
    if haskey(expr.aff.terms, var)
        expr.aff.terms[var] = coeff
        return expr
    else
        return expr + coeff * var
    end
end

# Fallback
function _set_variable_coefficient!(expr, var::GeneralVariableRef, coeff::Real)
    error("Unsupported function type for coefficient modification.")
end

## Query the coefficient of a variable
# GeneralVariableRef
function _affine_coefficient(
    func::GeneralVariableRef,
    var::GeneralVariableRef
    )::Float64
    return func == var ? 1.0 : 0.0
end

# GenericAffExpr
function _affine_coefficient(
    func::GenericAffExpr,
    var::GeneralVariableRef
    )::Float64
    return get(func.terms, var, 0.0)
end

# GenericQuadExpr
function _affine_coefficient(
    func::GenericQuadExpr,
    var::GeneralVariableRef
    )::Float64
    return get(func.aff.terms, var, 0.0)
end

# Fallback
function _affine_coefficient(func, var::GeneralVariableRef)
    error("Unsupported function type for coefficient queries.")
end

################################################################################
#                         PARAMETER REFERENCE METHODS
################################################################################
## Return an element of a parameter reference tuple given the model and index
# IndependentParameterIndex
function _make_param_tuple_element(model::InfiniteModel,
    idx::IndependentParameterIndex,
    )::GeneralVariableRef
    return _make_parameter_ref(model, idx)
end

# DependentParametersIndex
function _make_param_tuple_element(model::InfiniteModel,
    idx::DependentParametersIndex,
    )::Vector{GeneralVariableRef}
    dpref = DependentParameterRef(model, DependentParameterIndex(idx, 1))
    num_params = _num_parameters(dpref)
    return [GeneralVariableRef(model, idx.value, DependentParameterIndex, i)
            for i in 1:num_params]
end

"""
    parameter_refs(expr)::Tuple

Return the tuple of parameter references that determine the infinite
dependencies of `expr`.

**Example**
```julia-repl
julia> parameter_refs(my_expr)
(t,)
```
"""
function parameter_refs(
    expr::Union{JuMP.GenericAffExpr, JuMP.GenericQuadExpr}
    )::Tuple
    model = _model_from_expr(expr)
    if model === nothing
        return ()
    else
        obj_nums = sort!(_object_numbers(expr))
        obj_indices = _param_object_indices(model)[obj_nums]
        return Tuple(_make_param_tuple_element(model, idx) for idx in obj_indices)
    end
end

################################################################################
#                             EXPRESSION SEARCHING
################################################################################
# Check expression for a particular variable type via a recursive search
#=
function _has_variable(vrefs::Vector{GeneralVariableRef},
                       vref::GeneralVariableRef; prior=GeneralVariableRef[]
                       )::Bool
    if vrefs[1] == vref
        return true
    elseif _index_type(vrefs[1]) == MeasureIndex
        dvref = dispatch_variable_ref(vrefs[1])
        if length(vrefs) > 1
            return _has_variable(_all_function_variables(measure_function(dvref)),
                                 vref, prior = append!(prior, vrefs[2:end]))
        else
            return _has_variable(_all_function_variables(measure_function(dvref)),
                                 vref, prior = prior)
        end
    elseif length(vrefs) > 1
        return _has_variable(vrefs[2:end], vref, prior = prior)
    elseif length(prior) > 0
        return _has_variable(prior, vref)
    else
        return false
    end
end
=#


text:
struct SingleSpecies <: DiscreteSpeciesPool
    species::Vector{Symbol}

end

"""
    DiscreteUnipartiteSpeciesPool
"""
struct DiscreteUnipartiteSpeciesPool <: DiscreteSpeciesPool
    species::Vector{Symbol}
    metaweb::Matrix
    function DiscreteUnipartiteSpeciesPool(
        s::Vector{Symbol},
        adjmat::T,
    ) where {U<:Number,T<:AbstractArray{U,2}}
        new(s, adjmat)
    end
end

function DiscreteUnipartiteSpeciesPool(
    en::T,
) where {T<:EcologicalNetworks.UnipartiteNetwork}
    DiscreteUnipartiteSpeciesPool(en.S, Matrix(en.edges))
end

function DiscreteUnipartiteSpeciesPool(; numspecies = 30, connectance = 0.1)
    nm = EcologicalNetworks.nichemodel(numspecies, connectance)
    return DiscreteUnipartiteSpeciesPool(nm)
end
function DiscreteUnipartiteSpeciesPool(s::Vector{ST}) where {ST<:String}
    DiscreteUnipartiteSpeciesPool(Symbol.(s), zeros(Int32, length(s), length(s)))
end
function DiscreteUnipartiteSpeciesPool(s::Vector{ST}) where {ST<:Symbol}
    DiscreteUnipartiteSpeciesPool(s, zeros(Int32, length(s), length(s)))
end

function DiscreteUnipartiteSpeciesPool(
    s::Vector{String},
    adjmat::T,
) where {U<:Number,T<:AbstractArray{U,2}}
    DiscreteUnipartiteSpeciesPool(Symbol.(s), adjmat)
end

struct DiscreteKpartiteSpeciesPool <: DiscreteSpeciesPool
    species::Vector{Symbol}
    metaweb::Matrix
    partitions::Int
end

struct ContinuousUnipartiteSpeciesPool <: ContinuousSpeciesPool
    traitdists::Any
end

struct ContinuousKpartiteSpeciesPool <: ContinuousSpeciesPool
    traitdists::Any
    partitions::Int
end

numspecies(sp::T) where {T<:DiscreteSpeciesPool} = length(species(sp))
species(sp::T) where {T<:DiscreteSpeciesPool} = sp.species
metaweb(sp::T) where {T<:DiscreteSpeciesPool} = sp.metaweb

# some interfaces to EN generators
"""
    layers(speciespool::ST; element=Biomass, dims=(100,100))

    Converts a species pool to a set of layers of type `element`.
    Should be able to apply generic function. Mostly used to generate initial conditions.
"""
function layers(
    speciespool::ST;
    element = Biomass,
    dims = (100, 100),
) where {ST<:DiscreteSpeciesPool}
    specieslist = species(speciespool)
    layerlist = []
    for sp in specieslist
        thislayer = zeros(element, dims...)
        push!(layerlist, thislayer)
    end
    return layerlist
end
layers(
    speciespool::SPT,
    element::ET;
    dims = (100, 100),
) where {SPT<:DiscreteSpeciesPool,ET<:Measurement} =
    layers(speciespool, element = element, dims = dims)

function layernames(sp::SPT) where {SPT<:DiscreteSpeciesPool}
    return Symbol.(sp.species)
end

"""
    speciespool(net::UnipartiteNetwork)
"""
function speciespool(net::UnipartiteNetwork)
    DiscreteUnipartiteSpeciesPool([Symbol(i) for i in net.S], Matrix(net.edges))
end

"""
    interesting idea but getting names is better for filtering
    because returning a sub-metaweb losess information about interactions

"""
function Base.filter(f::Function, sp::T) where {T<:DiscreteSpeciesPool}

    newspecieslist = []
    indecies = []

    for (i, spec) in enumerate(species(sp))
        if f(String(spec)) == true
            push!(newspecieslist, spec)
            push!(indecies, i)
        end
    end

    newmetaweb = sp.metaweb[indecies, indecies]

    return T(newspecieslist, newmetaweb)
end

ext: jl


text:
mutable struct Skin
    inverseBindMatrices::Union{Nothing,Int}
    skeleton::Union{Nothing,Int}
    joints::Vector{Int}
    name::Union{Nothing,String}
    extensions::Union{Nothing,Dict}
    extras::Union{Nothing,Dict}
    function Skin()
        obj = new()
        obj.inverseBindMatrices = nothing
        obj.skeleton = nothing
        obj.name = nothing
        obj.extensions = nothing
        obj.extras = nothing
        obj
    end
end

function Base.setproperty!(obj::Skin, sym::Symbol, x)
    if sym === :inverseBindMatrices && x !== nothing
        x ≥ 0 || throw(ArgumentError("the index of the accessor containing the inverse-bind matrices should be ≥ 0"))
    elseif sym === :skeleton && x !== nothing
        skeleton ≥ 0 || throw(ArgumentError("the index of the node used as a skeleton root should be ≥ 0"))
    elseif sym === :joints
        isempty(x) && throw(ArgumentError("joints should not be empty."))
        allunique(x) || throw(ArgumentError("joint number must be unique."))
    end
    setfield!(obj, sym, x)
end

JSON3.StructType(::Type{Skin}) = JSON3.Mutable()
JSON3.omitempties(::Type{Skin}) = (:inverseBindMatrices, :skeleton, :name, :extensions, :extras)


text:
abstract type MonteCarloSampler <: Sampler end
abstract type FullSpaceSampler <: Sampler end

parallel_type(s::Sampler) = true

abstract type SamplerCache{T}  end

"""
    caches(sampler, problem, net, [parallel_type=NotParallel()])

Creates a `SamplerCache` for this problem and sampler.
"""
cache(s::Sampler, prob::AbstractProblem, net, par=NotParallel()) =
    cache(s, state(prob, net), net, par)

"""
    caches(sampler, v, net, [parallel_type=NotParallel()])

Creates a `SamplerCache` for this sampler and state.
"""
cache(s::Sampler, v::State, net, par=NotParallel()) =
    _sampler_cache(s, v, net, par)

"""
    init_sampler!(sampler, net, σ, [c=cache(sampler, σ, net)]) -> c

Initializes the sampler `sampler` and state `σ`. If no `SamplerCache` is
provided, one will be initialized and returned. The state σ is the first in the
list of sampled states.
"""
init_sampler!(s::Sampler, net, σ) = init_sampler!(s, net, σ, cache(s, σ, net))


text:
# ----------------------------------------------------------------------------------------------------------------------
# Question 5
# Fomulation of Matrix A
mat_a1_2 = zeros(numNodes^2, numNodes^2)
for mn = 1: numNodes^2
    mat_a1_2[mn, mn] = -1
end
mat_a3_2 = zeros(numNodes, numNodes^2)
seq_zeroTenNinety = collect(0:10:90)
for m = 1: numNodes
    mat_a3_2[m, (10 * m - 9): (10 * m)] = repeat([-1], 10)  # sent out
    mat_a3_2[m, (seq_zeroTenNinety + repeat([m],10))] = repeat([1], 10)  # sent in
    mat_a3_2[:, (11 * m - 10)] = repeat([0], 10)  # self-sending doesn't count
end
mat_a_2 = vcat(mat_a1_2, mat_a3_2)
# Fomulation of Matrix B
mat_b1_2 = zeros(numNodes^2, numNodes^2)
for mn = 1: numNodes^2
    mat_b1_2[mn, mn] = 170
end
mat_b3_2 = zeros(numNodes, numNodes^2)
mat_b_2 = vcat(mat_b1_2, mat_b3_2)
# Fomulation of Vector b
vec_b1_2 = zeros(numNodes^2)
vec_b3_2 = zeros(numNodes)
for n = 1: numNodes
    vec_b3_2[n] = vec_netInject[n] - vec_netEject[n]
end
vec_b_2 = vcat(vec_b1_2, vec_b3_2)
vec_b_2 = hcat(vec_b_2)
#
vec_max_y_2 = repeat([1], numNodes^2)
vec_max_y_2 = hcat(vec_max_y_2)
#
vec_c_2 = zeros(numNodes^2)
for m = 1: numNodes
    for n = 1: numNodes
        vec_c_2[10 * (m - 1) + n] = mat_distance[m, n]
    end
end
vec_c_2 = hcat(vec_c_2)
#
vec_f_2 = zeros(numNodes^2)
for m = 1: numNodes
    for n = 1: numNodes
        vec_f_2[10 * (m - 1) + n] = mat_fixedCost[m, n]
    end
end
vec_f_2 = hcat(vec_f_2)
#
Benders.milp(n_x = numNodes^2,
            n_y = numNodes^2,
            vec_min_y = repeat([0], numNodes^2),
            vec_max_y = vec_max_y_2,
            vec_c = vec_c_2,
            vec_f = vec_f_2,
            vec_b = vec_b_2,
            mat_a = mat_a_2,
            mat_b = mat_b_2,
            epsilon = 0.0001,
            timesIterationMax = 1000)
# obj - sum(mat_arcTwoNodes .* mat_fixedCost)
# ----------------------------------------------------------------------------------------------------------------------
# using LightGraphs
# using GraphPlot
# using Compose
# using Cairo
# using Fontconfig
# mat_result_y = [0 1 0 0 0 0 0 0 0 0;
#                 0 0 0 0 0 1 0 0 0 0;
#                 0 0 0 0 0 0 0 0 0 1;
#                 1 0 0 0 1 0 0 0 0 0;
#                 0 0 0 0 0 1 0 0 0 0;
#                 1 1 0 1 1 0 0 1 0 0;
#                 0 1 1 0 0 0 0 0 0 0;
#                 0 0 0 0 0 0 0 0 1 0;
#                 0 0 0 0 1 0 0 0 0 0;
#                 0 0 1 0 0 1 1 1 1 0]
# graph_initial = DiGraph(mat_arcTwoNodes)
# plot_graph_initial = gplot(graph_initial, nodelabel = vec_nameNodes)
# graph_result = graph_initial
# mat_arcTwoNodes_new = mat_result_y - mat_arcTwoNodes
# for m = 1: numNodes
#     for n = 1: numNodes
#         if mat_arcTwoNodes_new[m,n] == 1
#             add_edge!(graph_result, m, n)
#         end
#     end
# end
# gplot(graph_result, nodelabel = vec_nameNodes)
# draw(PNG("graph_initial.png", 16cm, 16cm), gplot(graph_initial))
# draw(PNG("graph_result.png", 16cm, 16cm), gplot(result))

# Final Corrected
vec_x= [0 0 0 0 0 0 0 0 0 0
        0 0 0 0 0 0 0 0 0 0
        0 0 0 0 0 0 0 0 0 0
        0 0 0 0 0 0 0 0 0 0
        0 0 0 0 0 0 0 0 0 0
        10 10 0 10 40 0 0 20 0 0
        0 20 20 0 0 0 0 0 0 0
        0 0 0 0 0 0 0 0 0 0
        0 0 0 0 0 0 0 0 0 0
        0 0 0 0 0 0 65 0 15 0]
vec_y= [0 1 0 0 0 0 0 0 0 0
        0 0 0 0 0 1 0 0 0 0
        0 0 0 0 0 0 0 0 0 1
        1 0 0 0 1 0 0 0 0 0
        0 0 0 0 0 1 0 0 0 0
        1 1 0 1 1 0 0 1 0 0
        0 1 1 0 0 0 0 0 0 0
        0 0 0 0 0 0 0 0 1 0
        0 0 0 0 1 0 0 0 0 0
        0 0 0 0 0 1 1 0 1 0]

vec_y: [0 1 0 0 0 0 0 0 0 0
        0 0 0 0 0 1 0 0 0 0
        0 0 0 0 0 0 0 0 0 1
        1 0 0 0 1 0 0 0 0 0
        0 0 0 0 0 1 0 0 0 0
        1 1 0 1 1 0 0 1 0 0
        0 1 1 0 0 0 0 0 0 0
        0 0 0 0 0 0 0 0 1 0
        0 0 0 0 1 0 0 0 0 0
        0 0 1 0 0 1 1 1 1 0]
result_q: 600.0

vec_x= [0 0 0 0 0 0 0 0 0 0
        0 0 0 0 0 0 0 0 0 0
        0 0 0 0 0 0 0 0 0 0
        0 0 0 0 0 0 0 0 0 0
        0 0 0 0 0 0 0 0 0 0
        10 30 0 10 40 0 0 0 0 0
        0 0 0 0 0 0 0 0 0 0
        0 0 0 0 0 0 0 0 15 0
        0 0 0 0 0 0 0 0 0 0
        0 0 20 0 0 0 25 35 0 0]

# vec_result_x = reshape(transpose(vec_x), numNodes^2)
vec_result_node = zeros(numNodes)
for m = 1: numNodes
    vec_result_node[m] = sum(vec_x[(10 * m - 9): (10 * m)]) - sum(vec_x[(seq_zeroTenNinety + repeat([m],10))])
end
vec_y= [0 1 0 0 0 0 0 0 0 0
        0 0 0 0 0 1 0 0 0 0
        0 0 0 0 0 0 0 0 0 1
        1 0 0 0 1 0 0 0 0 0
        0 0 0 0 0 1 0 0 0 0
        1 1 0 1 1 0 0 1 0 0
        0 1 1 0 0 0 0 0 0 0
        0 0 0 0 0 0 0 0 1 0
        0 0 0 0 1 0 0 0 0 0
        0 0 1 0 0 1 1 1 1 0]


text:
#
# IID case
#
function min(Xs::IIDRandomSequence{Continuous,Exponential})
  N,λ = length(Xs),rate(Xs.d)
  return Exponential(1./(N*λ))
end

function spacing(Xs::IIDRandomSequence,space::Integer)
  N,λ = length(Xs),rate(Xs.d)
  return Exponential(1./((N-space+1)*λs))
end

#
# INID case
#
function min(Xs::INIDRandomSequence{Continuous,Exponential})
  N,λs = length(Xs),[rate(X) for X in Xs.distributions]
  return Exponential(1./sum(λs))
end


text:
using Documenter, Impute

makedocs(
    modules=[Impute],
    format=Documenter.HTML(assets=["assets/invenia.css"]),
    pages=[
        "Home" => "index.md",
        "Walkthroughs" => [
            "Spatiotemporal" => "walkthroughs/spatiotemporal.md",
            "SVD" => "walkthroughs/svd.md",
        ],
        "API" => [
            "Declare Missings" => "api/declaremissings.md",
            "Validation" => "api/validation.md",
            "Filtering" => "api/filtering.md",
            "Imputation" => "api/imputation.md",
            "Chaining" => "api/chaining.md",
        ],
    ],
    repo="https://github.com/invenia/Impute.jl/blob/{commit}{path}#L{line}",
    sitename="Impute.jl",
    authors="Invenia Technical Computing Corporation",
    strict=true,
    checkdocs=:exports,
)

deploydocs(
    repo = "github.com/invenia/Impute.jl",
    devbranch = "master",
    push_preview = true,
)


text:
# --- testingChannelApplication.jl
# ---
# Description
#    Testing how wo apply time varying channel model
# ---
# Syntax
#    include("testingChannelApplication.jl");
# ---
# v 1.0 - Robin Gerzaguet.

module example_multipathChannel

# ----------------------------------------------------
# --- Modules
# ----------------------------------------------------
# --- External
using Plots
using Printf
using DSP
using FFTW
# --- Custom
using DigitalComm

# ----------------------------------------------------
## --- Functions
# ----------------------------------------------------
# --- Tx scheme
function minimalPhy(nbSymb,mcs,ofdm)
	# Minimal transmitter scheme
	nbSubcarriers = length(ofdm.allocatedSubcarriers);		  # Subcarrier length
	nbBits		  = Int(nbSymb * nbSubcarriers * log2(mcs));		  # nbBits
	bitSeq	      = genBitSequence(nbBits);					  # Binary sequence
	qamSeq		  = bitMappingQAM(mcs,bitSeq);				  # QAM generation
	qamMat		  = reshape(qamSeq,nbSubcarriers,nbSymb);	  # QAM matrix
	sigId		  = genSig(qamMat,ofdm);					  # OFDM signal
	return sigId;
end
function zeroPad(x,n)
	if n > length(x)
		out = [x;zeros(n-length(x))];
	else
		out = x;
	end
	return out;
end
# --- Receiver equalizer
function minimalEqualizer(qamRx,ofdm,strucChan)
	# Simple ZF equalisation scheme:
	qamEq	= zeros(Complex{Float64},size(qamRx));
	firAll	= zeros(Complex{Float64},ofdm.nFFT,size(qamRx,2));
	for iN = 1 : 1 : size(qamRx,2)
		# --- Getting all channel samples in current symbol
		# We get a nFFT x nbTaps matrix (all channel taps in current symbol)
		# We will consider that channel is constant, and that channel of interest
		# is at the middle of the OFDM symbol
		if Bool(strucChan.timeVarying)
			currChan  = strucChan.cir[(iN-1)*(ofdm.nFFT+ofdm.nCP) + ofdm.nCP + ofdm.nFFT÷2,:];
		else
			# COnstant FIR
			currChan = strucChan.cir;
		end
		# Switch to frequency domain
		firFreq   = fft(zeroPad(currChan,ofdm.nFFT));
		qamEq[:,iN] = qamRx[:,iN] ./ firFreq[ofdm.allocatedSubcarriers];
		firAll[:,iN] = firFreq;
	end
	return qamEq,firAll;
end

function plotFIR(fir,ofdm,samplingFreq)
	# --- Calculate freq axis
	freqAx	= (((0:1:ofdm.nFFT-1) ./ ofdm.nFFT) .- 0.5) .* samplingFreq;
	firC	= circshift(10*log10.(abs2.(fir)), ofdm.nFFT ÷2);
	pObj	= plot(freqAx,firC[:,1],label="First Channel");
	plot!(pObj,freqAx,firC[:,end],label="Last Channel");
	xlabel!("Frequency [Hz]");
	ylabel!("Channel magnitude");
	return pObj;
end

function main()
	# ----------------------------------------------------
	# --- Parameters
	# ----------------------------------------------------
	sizeFFT	  		= 1024;
	carrierFreq		= 5.2e9;
	speed	  		= 1300;
	samplingFreq	= 15.36e6;
	allocatedSubcarriers= getLTEAlloc(sizeFFT);

	# ----------------------------------------------------
	# --- Minimal Phy
	# ----------------------------------------------------
	# --- Init OFDM structure
	ofdm  = initOFDM(
								 sizeFFT,						    # --- nFFT                 : FFT size
								 72,						        # --- nCP                  : CP size
								 allocatedSubcarriers,		        # --- allocatedSubcarriers : Subcarrier allocation
								 );
	# --- Init parameters
	sigId	  = minimalPhy(14,4,ofdm);
	# ----------------------------------------------------
	# --- Multipath Channel
	# ----------------------------------------------------
	# --- profile
	profile         = "etu";
	# --- Channel object
	channelObj      = initChannel(profile,carrierFreq,samplingFreq,speed);
	# --- Channel instance
	strucChan0       = getChannel(length(sigId),channelObj,-1);
	# ----------------------------------------------------
	# --- Channel application
	# ----------------------------------------------------
	sigChan	  = applyChannel(sigId,strucChan0);
	qamRx	  = ofdmSigDecode(sigChan,ofdm);
	scatter(qamRx[:])
	# ----------------------------------------------------
	# --- Setting a minimal equalizer
	# ----------------------------------------------------
	# --- Assuming a perfect CSI
	qamEq,firFreq	= minimalEqualizer(qamRx,ofdm,strucChan0)
	#
	pObj = scatter(qamRx[:],xlims=(-1,1),ylims=(-1,1),label="before eq.");
	scatter!(pObj,qamEq[:],xlims=(-1,1),ylims=(-1,1),label="after eq.");
	display(pObj);
	#
	pObj2 = plotFIR(firFreq,ofdm,samplingFreq);
	display(pObj2);

end

end


text:
module Plot

import PyPlot; plt = PyPlot
plt.svg(true)
using PyCall

export plt, make_bar_plot

"""
    make_bar_plot(lookup_value, groups, keys; ...)

Create a bar plot.

# Examples

```julia
fig, ax = make_bar_plot(
    ["A", "B", "C", "D", "E"], ["1", "2", "3", "4"],
) do group, key
    rand()
end
ax.set_title("...")
ax.set_xlabel("...")
ax.set_ylabel("...")
```
```julia
fig, ax = make_bar_plot(
    ["A"=>1, "B"=>2, "C"=>3, "D"=>4, "E"=>5], [1, 2, 3, 4, 5],
    label_rotation=45, label_anchor="right",
) do group, key
    group * key
end
ax.set_title("...")
ax.set_xlabel("...")
ax.set_ylabel("...")
```
"""
function make_bar_plot(lookup_value::Function, groups, keys;
                       colors=nothing, patterns=nothing, ax=nothing,
                       width::Real=0.7, gap::Real=0.2, label_rotation::Real=0.0,
                       label_anchor=nothing)
    is_latex = ax isa PyObject && startswith(ax.__module__, "latextools")

    n_groups = length(groups)
    n_bars = length(keys)
    width /= n_bars
    gap /= n_bars
    step = width + gap
    x_arr = collect(1:n_groups)
    offset = -(step * (n_bars+1)) / 2

    if ax === nothing
        fig, ax = plt.subplots()
    else
        fig = nothing
    end

    group_strs = [g isa Pair ? g[1] : string(g) for g in groups]
    groups = [g isa Pair ? g[2] : g for g in groups]
    for (i, key) in enumerate(keys)
        key_str, key = (key isa Pair ? (key[1], key[2])
                        : (string(key), key))
        vals = [lookup_value(g isa Pair ? g[2] : g, key) for g in groups]
        c = colors !== nothing && i <= length(colors) ? colors[i] : nothing
        p = (patterns !== nothing && i <= length(patterns)
             ? patterns[i] : nothing)
        if is_latex
            ax.bar(group_strs, vals, color=c, pattern=p, legend=key_str,
                   extra_config="pattern color=$c,")
        else
            ax.bar(x_arr .+ (offset + i*step), vals, width, color=c,
                   label=key_str)
        end
    end

    if !is_latex
        ax.set_xticks(collect(1:n_groups))
        ax.set_xticklabels(group_strs, rotation=label_rotation, ha=label_anchor)
        ax.legend()
    end
    fig === nothing ? nothing : (fig, ax)
end

end


text:
module TestBayes
using Test

using Random
using IterativeSolvers

using Catlab, Catlab.Doctrines, Catlab.WiringDiagrams, Catlab.Programs
using Catlab.LinearAlgebra.Bayes

import LinearAlgebra: norm, svd

import LinearMaps: BlockDiagonalMap, UniformScaling
using LinearOperators

using BayesNets
import BayesNets:
  BayesNet
# Evaluation
############

# BayesMap instance
#--------------------

#Create Equations
a_equation = randn(100)
b_equation = randn(100) .+ 2*a_equation .+ 3
c_equation = randn(100) .+ 2*b_equation .+ 3

#Build Rain, Sprinkler, GrassWet Example BayesNet
data = DataFrame(rain=a_equation, sprinkler=b_equation, grasswet=c_equation)
cpdA = fit(StaticCPD{Normal}, data, :rain)
cpdB = fit(LinearGaussianCPD, data, :sprinkler, [:rain])
cpdC = fit(LinearGaussianCPD, data, :grasswet, [:rain,:sprinkler])

bn1 = BayesNet([cpdA, cpdB, cpdC])
bn1_dom = [1]
bn1_codom = [2,3]

openbn1 = OpenBayesNet(bn1, bn1_dom, bn1_codom)

#Create Equations
a_equation = randn(100)
b_equation = randn(100) .+ 2*a_equation .+ 3

#Build Rain, Sprinkler, GrassWet Example BayesNet
data = DataFrame(mosquito=a_equation, bird=b_equation)
cpdA = fit(StaticCPD{Normal}, data, :mosquito)
cpdB = fit(LinearGaussianCPD, data, :bird, [:mosquito])

bn2 = BayesNet([cpdA, cpdB])
bn2_dom = [1]
bn2_codom = [2]

openbn2 = OpenBayesNet(bn2, bn2_dom, bn2_codom)

openbn_otimes = Bayes.otimes(openbn1, openbn2)
@test dom(openbn_otimes).N == dom(openbn1).N + dom(openbn2).N
@test codom(openbn_otimes).N == codom(openbn1).N + codom(openbn2).N

# @test (f⋅h)*x == h*(f*x)
# @test (f⊕g)*[x;y] == [f*x; g*y]
# @test braid(dom(f),dom(g)) * [x;y] == [y;x]

# @test mcopy(dom(f))*x == [x;x]
# @test delete(dom(f))*x == []
# @test plus(dom(f))*[x;x] == 2x
# @test zero(dom(f))*Float64[] == zero(x)

# @test (h+k)*x == h*x + k*x
# @test scalar(dom(f),3)*x == 3*x
# @test antipode(dom(f))*x == -1*x
# @test adjoint(g)*y == g'*y
end


text:
using PyPlot
nu = 3
mu = linspace(0,1,10000)
ntrial = 100

function iofmu(u,mu)
  intensity = ones(mu)
  for i=1:length(u)
    intensity .-= u[i]*(1 .-mu).^i
  end
  return intensity
end

function diofmudmu(u,mu)
  didmu = zeros(mu)
  for i=1:length(u)
    didmu .+= i*u[i]*(1 .-mu).^(i-1)
  end
  return didmu
end

coeff = zeros(nu,ntrial)
u = zeros(nu)
clf()
for i=1:ntrial
  bad = true
  while bad
    u = -1.0+2.0*rand(nu)
    test = minimum(iofmu(u,mu)) < 0.0 || minimum(diofmudmu(u,mu)) <= 0.0
    if test
      bad = true
#      plot(mu,iofmu(u,mu),linestyle=":")
    else
      bad = false
      plot(mu,iofmu(u,mu))
      read(STDIN,Char)
      clf()
    end
  end
end


text:
module Vasilek

include(joinpath("MaxwellSolver", "FDTD1D.jl"))
include(joinpath("MaxwellSolver", "PoissonFourier1D.jl"))

include(joinpath("VlasovSolver", "LaxWendroff.jl"))
include(joinpath("VlasovSolver", "Upwind.jl"))
include(joinpath("VlasovSolver", "Godunov.jl"))
include(joinpath("VlasovSolver", "SemiLagrangian.jl"))
include(joinpath("VlasovSolver", "PFC.jl"))
include(joinpath("VlasovSolver", "PFCNonUniform.jl"))

include(joinpath("BoltzmannSolver", "BGK.jl"))
include(joinpath("BoltzmannSolver", "Landau1P.jl"))

end  # module Vasilek


text:
@testset "variable.jl" begin
    @testset ":tallied" begin
        global tallied = Variable{Int64}(data=:tally, history=true)

        update_list = []
        final_value = nothing

        @process test_process(value_durations) begin
            for (value, duration) in value_durations
                tallied.value = final_value = value
                push!(update_list, tallied.prev_update)
                work(duration)
            end
        end

        @simulation begin
            value_durations = [(1, 2.0), (2, 1.0), (3, 2.0), (4, 3.0)]

            @schedule now test_process(value_durations)
            @schedule at 1.0 test_process(value_durations)
            @schedule in 3.0 test_process(value_durations)

            start_simulation()
        end

        want = [0.0, 1.0, 2.0, 3.0, 3.0, 3.0, 4.0, 5.0, 5.0, 6.0, 6.0, 8.0]
        @test update_list == want
        @test typeof(tallied.history) == SimLynx.TalliedHistory{Int64}
        @test typeof(tallied.stats) == SimLynx.TalliedStats{Int64}

    end

    @testset ":accumulate" begin
        global accumulated = Variable{Int64}(data=:accumulate, history=true)

        update_list = []
        final_value = nothing

        @process test_process(value_durations) begin
            for (value, duration) in value_durations
                accumulated.value = final_value = value
                push!(update_list, accumulated.prev_update)
                work(duration)
            end
        end

        @simulation begin
            value_durations = [(1, 2.0), (2, 1.0), (3, 2.0), (4, 3.0)]

            @schedule now test_process(value_durations)
            @schedule at 1.0 test_process(value_durations)
            @schedule in 3.0 test_process(value_durations)

            start_simulation()
        end

        want = [0.0, 1.0, 2.0, 3.0, 3.0, 3.0, 4.0, 5.0, 5.0, 6.0, 6.0, 8.0]
        @test update_list == want
        @test typeof(accumulated.history) == SimLynx.AccumulatedHistory{Int64}
        @test typeof(accumulated.stats) == SimLynx.AccumulatedStats{Int64}
    end
end


text:
using TimeTrees
using Base.Test

# Test newick parsing:
t = TimeTree("((A:1,B:1):1,C:2):0;")
r = t.root

# Test getNodes and getLeaves
nodes = getNodes(t)
leaves = getLeaves(t)
l = leaves[1]

# Test basic methods
@test isRoot(r) == true
@test isRoot(l) == false
@test isLeaf(r) == false
@test isLeaf(l) == true

# Test tree sorting
@test getDecendentCount(r) == 5
s = getSorted(t)
@test s.root.children[1].label == "C"


text:
function DSSTPropagatorBuilder(arg0::Orbit, arg1::ODEIntegratorBuilder, arg2::jdouble, arg3::PropagationType, arg4::PropagationType)
    return DSSTPropagatorBuilder((Orbit, ODEIntegratorBuilder, jdouble, PropagationType, PropagationType), arg0, arg1, arg2, arg3, arg4)
end

function DSSTPropagatorBuilder(arg0::Orbit, arg1::ODEIntegratorBuilder, arg2::jdouble, arg3::PropagationType, arg4::PropagationType, arg5::AttitudeProvider)
    return DSSTPropagatorBuilder((Orbit, ODEIntegratorBuilder, jdouble, PropagationType, PropagationType, AttitudeProvider), arg0, arg1, arg2, arg3, arg4, arg5)
end

function add_force_model(obj::DSSTPropagatorBuilder, arg0::DSSTForceModel)
    return jcall(obj, "addForceModel", void, (DSSTForceModel,), arg0)
end

function build_kalman_model(obj::DSSTPropagatorBuilder, arg0::List, arg1::List, arg2::ParameterDriversList)
    return jcall(obj, "buildKalmanModel", DSSTKalmanModel, (List, List, ParameterDriversList), arg0, arg1, arg2)
end

function build_ls_model(obj::DSSTPropagatorBuilder, arg0::Vector{IntegratedPropagatorBuilder}, arg1::List, arg2::ParameterDriversList, arg3::ModelObserver)
    return jcall(obj, "buildLSModel", DSSTBatchLSModel, (Vector{IntegratedPropagatorBuilder}, List, ParameterDriversList, ModelObserver), arg0, arg1, arg2, arg3)
end

function build_propagator(obj::DSSTPropagatorBuilder, arg0::Vector{jdouble})
    return jcall(obj, "buildPropagator", DSSTPropagator, (Vector{jdouble},), arg0)
end

function copy(obj::DSSTPropagatorBuilder)
    return jcall(obj, "copy", DSSTPropagatorBuilder, ())
end

function get_all_force_models(obj::DSSTPropagatorBuilder)
    return jcall(obj, "getAllForceModels", List, ())
end

function get_integrator_builder(obj::DSSTPropagatorBuilder)
    return jcall(obj, "getIntegratorBuilder", ODEIntegratorBuilder, ())
end

function get_mass(obj::DSSTPropagatorBuilder)
    return jcall(obj, "getMass", jdouble, ())
end

function set_mass(obj::DSSTPropagatorBuilder, arg0::jdouble)
    return jcall(obj, "setMass", void, (jdouble,), arg0)
end


text:
using Documenter
using BaseExtensions

DocMeta.setdocmeta!(BaseExtensions, :DocTestSetup, :(using BaseExtensions); recursive=true)

makedocs(;
    modules=[BaseExtensions],
    authors="Jishnu Bhattacharya",
    repo="https://github.com/jishnub/BaseExtensions.jl/blob/{commit}{path}#L{line}",
    sitename="BaseExtensions.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://jishnub.github.io/BaseExtensions.jl",
        assets=String[],
    ),
    pages=[
        "API" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/jishnub/BaseExtensions.jl",
    devbranch = "main",
)


text:
module CreateTableItems

import SearchLight.Migrations: create_table, column, primary_key, add_index, drop_table

function up()
  create_table(:items) do
    [
      primary_key()
      column(:a, :string, limit=100)
      column(:b, :int)
    ]
  end

  add_index(:items, :a)
end

function down()
  drop_table(:items)
end

end

max_stars_repo_name: jaantollander/genie-webapp-db

max_issues_repo_name: jaantollander/genie-webapp-db

max_forks_repo_name: jaantollander/genie-webapp-db


text:
"""
This module provides a collection of multivariate polynomial bases.

The exported names are:

$(EXPORTS)
"""
module Polynomials

using DocStringExtensions
using LinearAlgebra: mul!
using Gridap.Helpers
using Gridap.Arrays
using Gridap.TensorValues
using Gridap.Fields

import Gridap.Fields: evaluate!
import Gridap.Fields: return_cache
import Gridap.Arrays: return_type

export MonomialBasis
export QGradMonomialBasis
export QCurlGradMonomialBasis
export PCurlGradMonomialBasis
export get_exponents

export get_order
export get_orders
export num_terms

include("MonomialBases.jl")

include("QGradMonomialBases.jl")

include("QCurlGradMonomialBases.jl")

include("PCurlGradMonomialBases.jl")

end # module


text:
using CSV, StateSpaceModels, Plots, Dates

# Load the Nile annual flow dataset
nile = CSV.read("Nile.csv")

# Convert data to an array of Float64
flow = Float64.(nile.Flow)

# Plot the data
p1 = plot(nile.Year, flow, label = "Annual flow", legend = :topright, color = :black)

# Specify the state-space model
model = local_level(flow)

# Estimate the state-space model
ss = statespace(model)

# Innovations
p2 = plot(nile.Year, ss.filter.v, label = "Innovations (v_t)", legend = :topright, color = :black)

# Plot filtered state and 90% confidence interval
std_ptt = sqrt.(ss.filter.Ptt[1, 1, :])
att = vec(ss.filter.att)
lb = vec(att + 1.64 * std_ptt)
ub = vec(att - 1.64 * std_ptt)
p3 = plot!(p1, nile.Year, att, ribbon = [ub - att, lb + att],
           label = "Filtered state (att)", color = :red, fillalpha = 0.1)


text:
using EffectiveWaves
using JLD2
using FileIO
# using OffsetArrays
# using Plots; pyplot()
# Plots.scalefontsizes(1.5)

# physical parameters
medium = Medium(1.0,1.0+0.0im)
specie = Specie(ρ=0.1, r=0.4, c=0.1, volfrac=0.04)
specie = Specie(ρ=0.2, r=0.4, c=0.2, volfrac=0.2)
# specie = Specie(ρ=20.0, r=0.4, c=20.0, volfrac=0.10)
species = [specie]

k = 1.0
# k = 0.4;
θin = 0.4
# θin = 0.0
ω = real(k*medium.c)
radius_multiplier = 1.005

kws = ()
T = Float64
tol = 1e-7
ho = maximum_basis_order(ω, medium, [specie]; tol=tol) -1
basis_order = ho

# for X = 0.0:0.05:12. the relative convergence error for
# Specie(ρ=0.5, r=0.5, c=0.2, volfrac=0.1), k = 0.8, and θin = 0.0,
#  inbetween 6 and 8 is 0.05%
a12 = radius_multiplier*2.0*specie.r
num_wavenumbers = 90;
num_wavenumbers = 40;
# X = 0.0:(a12k/31):15.
# X = 0.0:0.04:12.

wave_effs = WaveModes(ω, medium, [specie];
    basis_order=ho, tol = 1e-7,  θin = θin, verbose = true,
    num_wavenumbers = num_wavenumbers,
    apply_meshing = true,
    radius_multiplier = radius_multiplier,
    mesh_refine = 0.18, mesh_points = 10, mesh_size=2.0,
    extinction_rescale = false
)

# k_effs = [w.wavenumber for w in wave_effs];
# mesh_ks = wavenumbers_mesh(ω, k_effs, medium, [specie];
#     basis_order=ho, tol = 1e-7,  θin = θin, verbose = true,
#     num_wavenumbers = num_wavenumbers,
#     radius_multiplier = radius_multiplier, mesh_refine = 0.2
#     );

using Plots; pyplot()
plot(wave_effs, color=:black)
plot!(wave_effs[1:num_wavenumbers], markercolor=:white)

collect(x_mesh_match(wave_effs[1:num_wavenumbers])[2])

# Calculate discretised average wave directly from governing equations
L, x = x_mesh_match(wave_effs[1:num_wavenumbers]; a12 = a12, max_size = 220, tol=1e-9);

x = x_mesh(wave_effs[1], wave_effs[2]; a12 = a12, max_size = 1600, tol=1e-8)

avg_wave_x = DiscretePlaneWaveMode(ω, medium, specie;
    basis_order=ho, θin=θin, x = x,
    wave_effs = wave_effs[1:num_wavenumbers],
    radius_multiplier = radius_multiplier, tol = 1e-8)
X = avg_wave_x.x.*k
# x = X./k

scatter(avg_wave_x)
# scatter!(avg_wave_x, hankel_indexes = 0:ho, apply=abs)

L, x = x_mesh_match(wave_effs[1:num_wavenumbers]; a12 = a12, max_size = 220, tol=1e-9);
x_short = x[1:Int(round(length(x)/3))]

match_w = MatchPlaneWaveMode(ω, medium, specie;
    x = x_short, L_match = 20,
    radius_multiplier = radius_multiplier,
    basis_order=ho,
    tol = 1e-7, θin = θin, wave_effs = wave_effs[1:num_wavenumbers]);

    ω = real(k*medium.c)
    radius_multiplier = 1.005

    kws = ()
    T = Float64
    tol = 1e-7
    ho = maximum_basis_order(ω, medium, [specie]; tol=tol) -1
    basis_order = ho

plot(match_w, hankel_indexes = 0:1)
plot(avg_wave_x, hankel_indexes = 0:1, seriestype=:scatter, xlims = (0.,5.))
plot!(match_w, hankel_indexes = 0:1, seriestype=:line, xlims = (0.,5.))

avg_wave_small = DiscretePlaneWaveMode(ω, medium, specie; basis_order=ho, θin=θin, x=match_w.discrete_wave.x,
    radius_multiplier = radius_multiplier, tol = 1e-6, max_size=500)
plot!(avg_wave_small, seriestype=:line, linestyle=:dash)

# X = avg_wave_x.x
X0 = 2. # match from X[L]
L = findmin(abs.(X .- X0))[2]
X_max = X[Int(round(0.65*length(X)))] # maximum X considered to have the correct numerical solution
i_max = findmin(abs.(X .- X_max))[2]

# # If we assume the average wave is a sum of plane waves everywhere, then we can estimate which of these effectives waves will have decayed too much already, and therefore discard them
avg_wave_effs = [DiscretePlaneWaveMode(X[L:L+1], wave) for wave in wave_effs]
ref_norm = norm(avg_wave_effs[1].amplitudes[end,:,:])/norm(wave_effs[1].eigenvectors)
eff_inds = find(
    norm(avg_wave_effs[i].amplitudes[end,:,:])/norm(wave_effs[i].eigenvectors) > ref_norm*tol
for i in eachindex(wave_effs))

k_effs = k_effs[eff_inds]
wave_effs = wave_effs[eff_inds];

J = L + 3*length(k_effs) # double number of points than waves to avoid overfitting
X_match = X[L:J]
avg_wave = DiscretePlaneWaveMode(ho, X[1:J], avg_wave_x.amplitudes[1:J,:,:])

# rescale to lessen roundoff errors
for i in eachindex(k_effs)[2:end]
    wave_effs[i].eigenvectors = wave_effs[i].eigenvectors * exp(-im*k_effs[i]*cos(wave_effs[i].θ_eff)*X_match[1]/k) / norm(wave_effs[i].eigenvectors)
end

include("match_and_extinction.jl")

Y = match_only_arrays(ω, wave_effs, L, X[1:J], medium, [specie]; θin = θin)
(w_vec, G_arr) = extinc_arrays(ω, wave_effs, L, X[1:J], medium, [specie]; θin = θin)
(LT_mat, E_mat, b_eff) = match_arrays(ω, wave_effs, L, X[1:J], medium, [specie]; θin = θin)

LA = LT_mat*avg_wave.amplitudes[:]
qA = sum(G_arr[ind] * avg_wave.amplitudes[ind] for ind in eachindex(avg_wave.amplitudes))

αs = LA + b_eff
scatter([abs.(exp.(im*k_effs./k)), sign.(real.(k_effs)).*abs.(αs)], lab = ["exp(im*k_eff)" "sign(Re k_eff) * abs αs"])

# check extinction
    norm((sum(w_vec .* αs) - (qA + im*k^2))/sum(w_vec .* αs)) < 1e-10

# For only one effective wave
    (w_vec, q_arr) = extinc_arrays(ω, wave_effs[1:1], L, X[1:J], medium, [specie]; θin = θin);
    qA = sum(q_arr[ind] * avg_wave.amplitudes[ind] for ind in eachindex(avg_wave.amplitudes));

    α1 = (qA + im*k^2)/w_vec[1]

# check matching
    αs_match = Y*avg_wave.amplitudes[:]
    # Calculate the discretised wave from these effective wave
    avg_wave_effs = [DiscretePlaneWaveMode(X[L:i_max], wave) for wave in wave_effs]

    amplitudes_eff = sum(αs[i] * avg_wave_effs[i].amplitudes[:,:,:] for i in eachindex(avg_wave_effs))
    avg_wave_eff = DiscretePlaneWaveMode(ho, avg_wave_effs[1].x, amplitudes_eff)

    amplitudes_eff = sum(αs_match[i] * avg_wave_effs[i].amplitudes[:,:,:] for i in eachindex(avg_wave_effs))
    avg_wave_eff_match = DiscretePlaneWaveMode(ho, avg_wave_effs[1].x, amplitudes_eff)

    ps = map(0:ho) do n
        max_w = maximum(real.(avg_wave.amplitudes[L:end,n+ho+1]))
        min_w = minimum(real.(avg_wave.amplitudes[L:end,n+ho+1]))
        max_w = (max_w > 0) ? 1.2*max_w : 0.0
        min_w = (min_w < 0) ? 1.2*min_w : 0.0
        plot(x -> x,y->max_w,X[L],X[J],line=0,fill=(0,:orange), fillalpha=0.4, lab="match region")
        plot!(x -> x,y->min_w,X[L],X[J],line=0,fill=(0,:orange), fillalpha=0.4, lab="")
        plot!(xlab = "X", ylab = "Re amplitude", title="Hankel order $n")
        plot!(avg_wave_eff.x, real.(avg_wave_eff.amplitudes[:,n+ho+1]), lab = "Wave eff", linewidth=2)
        # plot!(avg_wave_eff_match.x, real.(avg_wave_eff_match.amplitudes[:,n+ho+1]), lab = "Wave eff match", linewidth=2)
        scatter!(X[L:i_max], real.(avg_wave_x.amplitudes[L:i_max,n+ho+1]), lab = "Avg wave", linewidth =2 )
        plot!(avg_wave_effs[1].x, real.(α1*avg_wave_effs[1].amplitudes[:,n+ho+1]), lab = "Wave eff 1", linewidth =2, linestyle=:dash)
    end
    plot(ps...)

    plot(xlab = "X", ylab = "Abs amplitude")
    for n = -ho:ho
        scatter!(avg_wave_eff.x, abs.(avg_wave_eff.amplitudes[:,n+ho+1]), lab = "Wave eff hankel $n")
        plot!(X[L:end], abs.(avg_wave_x.amplitudes[L:end,n+ho+1]), lab = "Avg wave hankel $n" )
    end
    gui()

max_stars_repo_name: arturgower/EffectiveWaves


text:
"""
Pseudopotential correction energy. TODO discuss the need for this.
"""
struct PspCorrection end
(::PspCorrection)(basis) = TermPspCorrection(basis)

struct TermPspCorrection{T} <: Term
    energy::T  # precomputed energy
end
function TermPspCorrection(basis::PlaneWaveBasis)
    model = basis.model
    if model.n_dim != 3 && any(attype isa ElementPsp for attype in model.atoms)
        error("The use of pseudopotentials is only sensible for 3D systems.")
    end
    TermPspCorrection(energy_psp_correction(model))
end

function ene_ops(term::TermPspCorrection, basis::PlaneWaveBasis, ψ, occ; kwargs...)
    (E=term.energy, ops=[NoopOperator(basis, kpt) for kpt in basis.kpoints])
end

"""
Compute the correction term for properly modelling the interaction of the pseudopotential
core with the compensating background charge induced by the `Ewald` term.
"""
function energy_psp_correction(lattice::AbstractMatrix{T}, atoms, atom_groups) where T
    psp_groups = [group for group in atom_groups if atoms[first(group)] isa ElementPsp]
    isempty(psp_groups) && return zero(T)

    n_electrons::Int = sum(n_elec_valence, atoms)
    correction_per_cell = sum(
        length(group) * eval_psp_energy_correction(T, atoms[first(group)].psp, n_electrons)
        for group in psp_groups
    )

    correction_per_cell / compute_unit_cell_volume(lattice)
end
function energy_psp_correction(model::Model)
    energy_psp_correction(model.lattice, model.atoms, model.atom_groups)
end


text:
using LazyArrays, LinearAlgebra, FillArrays, Test
    import LazyArrays: MemoryLayout, DenseRowMajor, DenseColumnMajor, StridedLayout,
                            ConjLayout, RowMajor, ColumnMajor, UnknownLayout,
                            SymmetricLayout, HermitianLayout, UpperTriangularLayout,
                            UnitUpperTriangularLayout, LowerTriangularLayout,
                            UnitLowerTriangularLayout, ScalarLayout,
                            hermitiandata, symmetricdata, FillLayout, ZerosLayout,
                            VcatLayout, BroadcastLayout, Add

struct FooBar end
struct FooNumber <: Number end

@testset "MemoryLayout" begin
    @testset "Trivial" begin
        @test MemoryLayout(1.0) == MemoryLayout(1) == MemoryLayout(FooNumber()) == ScalarLayout()
        @test MemoryLayout(FooBar()) == UnknownLayout()

        A = randn(6)
        @test MemoryLayout(A) == MemoryLayout(Base.ReshapedArray(A,(2,3),())) == MemoryLayout(reinterpret(Float32,A)) ==
            DenseColumnMajor()

    end

    @testset "adjoint and transpose MemoryLayout" begin
        A = [1.0 2; 3 4]
        @test MemoryLayout(A') == DenseRowMajor()
        @test MemoryLayout(transpose(A)) == DenseRowMajor()
        B = [1.0+im 2; 3 4]
        @test MemoryLayout(B') == ConjLayout(DenseRowMajor())
        @test MemoryLayout(transpose(B)) == DenseRowMajor()
        VA = view(A, 1:1, 1:1)
        @test MemoryLayout(VA') == RowMajor()
        @test MemoryLayout(transpose(VA)) == RowMajor()
        VB = view(B, 1:1, 1:1)
        @test MemoryLayout(VB') == ConjLayout(RowMajor())
        @test MemoryLayout(transpose(VB)) == RowMajor()
        VA = view(A, 1:2:2, 1:2:2)
        @test MemoryLayout(VA') == StridedLayout()
        @test MemoryLayout(transpose(VA)) == StridedLayout()
        VB = view(B, 1:2:2, 1:2:2)
        @test MemoryLayout(VB') == ConjLayout(StridedLayout())
        @test MemoryLayout(transpose(VB)) == StridedLayout()
        VA2 = view(A, [1,2], :)
        @test MemoryLayout(VA2') == UnknownLayout()
        @test MemoryLayout(transpose(VA2)) == UnknownLayout()
        VB2 = view(B, [1,2], :)
        @test MemoryLayout(VB2') == UnknownLayout()
        @test MemoryLayout(transpose(VB2)) == UnknownLayout()
        VAc = view(A', 1:1, 1:1)
        @test MemoryLayout(VAc) == RowMajor()
        VAt = view(transpose(A), 1:1, 1:1)
        @test MemoryLayout(VAt) == RowMajor()
        VBc = view(B', 1:1, 1:1)
        @test MemoryLayout(VBc) == ConjLayout(RowMajor())
        VBt = view(transpose(B), 1:1, 1:1)
        @test MemoryLayout(VBt) == RowMajor()
    end

    @testset "Symmetric/Hermitian MemoryLayout" begin
         A = [1.0 2; 3 4]
         @test MemoryLayout(Symmetric(A)) == SymmetricLayout(DenseColumnMajor(),'U')
         @test MemoryLayout(Hermitian(A)) == SymmetricLayout(DenseColumnMajor(),'U')
         @test MemoryLayout(Transpose(Symmetric(A))) == SymmetricLayout(DenseColumnMajor(),'U')
         @test MemoryLayout(Transpose(Hermitian(A))) == SymmetricLayout(DenseColumnMajor(),'U')
         @test MemoryLayout(Adjoint(Symmetric(A))) == SymmetricLayout(DenseColumnMajor(),'U')
         @test MemoryLayout(Adjoint(Hermitian(A))) == SymmetricLayout(DenseColumnMajor(),'U')
         @test MemoryLayout(view(Symmetric(A),:,:)) == SymmetricLayout(DenseColumnMajor(),'U')
         @test MemoryLayout(view(Hermitian(A),:,:)) == SymmetricLayout(DenseColumnMajor(),'U')
         @test MemoryLayout(Symmetric(A')) == SymmetricLayout(DenseRowMajor(),'U')
         @test MemoryLayout(Hermitian(A')) == SymmetricLayout(DenseRowMajor(),'U')
         @test MemoryLayout(Symmetric(transpose(A))) == SymmetricLayout(DenseRowMajor(),'U')
         @test MemoryLayout(Hermitian(transpose(A))) == SymmetricLayout(DenseRowMajor(),'U')

         @test symmetricdata(Symmetric(A)) ≡ A
         @test symmetricdata(Hermitian(A)) ≡ A
         @test symmetricdata(Transpose(Symmetric(A))) ≡ A
         @test symmetricdata(Transpose(Hermitian(A))) ≡ A
         @test symmetricdata(Adjoint(Symmetric(A))) ≡ A
         @test symmetricdata(Adjoint(Hermitian(A))) ≡ A
         @test symmetricdata(view(Symmetric(A),:,:)) ≡ A
         @test symmetricdata(view(Hermitian(A),:,:)) ≡ A
         @test symmetricdata(Symmetric(A')) ≡ A'
         @test symmetricdata(Hermitian(A')) ≡ A'
         @test symmetricdata(Symmetric(transpose(A))) ≡ transpose(A)
         @test symmetricdata(Hermitian(transpose(A))) ≡ transpose(A)

         B = [1.0+im 2; 3 4]
         @test MemoryLayout(Symmetric(B)) == SymmetricLayout(DenseColumnMajor(),'U')
         @test MemoryLayout(Hermitian(B)) == HermitianLayout(DenseColumnMajor(),'U')
         @test MemoryLayout(Transpose(Symmetric(B))) == SymmetricLayout(DenseColumnMajor(),'U')
         @test MemoryLayout(Transpose(Hermitian(B))) == UnknownLayout()
         @test MemoryLayout(Adjoint(Symmetric(B))) == UnknownLayout()
         @test MemoryLayout(Adjoint(Hermitian(B))) == HermitianLayout(DenseColumnMajor(),'U')
         @test MemoryLayout(view(Symmetric(B),:,:)) == SymmetricLayout(DenseColumnMajor(),'U')
         @test MemoryLayout(view(Hermitian(B),:,:)) == HermitianLayout(DenseColumnMajor(),'U')
         @test MemoryLayout(Symmetric(B')) == UnknownLayout()
         @test MemoryLayout(Hermitian(B')) == UnknownLayout()
         @test MemoryLayout(Symmetric(transpose(B))) == SymmetricLayout(DenseRowMajor(),'U')
         @test MemoryLayout(Hermitian(transpose(B))) == HermitianLayout(DenseRowMajor(),'U')

         @test symmetricdata(Symmetric(B)) ≡ B
         @test hermitiandata(Hermitian(B)) ≡ B
         @test symmetricdata(Transpose(Symmetric(B))) ≡ B
         @test hermitiandata(Adjoint(Hermitian(B))) ≡ B
         @test symmetricdata(view(Symmetric(B),:,:)) ≡ B
         @test hermitiandata(view(Hermitian(B),:,:)) ≡ B
         @test symmetricdata(Symmetric(B')) ≡ B'
         @test hermitiandata(Hermitian(B')) ≡ B'
         @test symmetricdata(Symmetric(transpose(B))) ≡ transpose(B)
         @test hermitiandata(Hermitian(transpose(B))) ≡ transpose(B)
     end

    @testset "triangular MemoryLayout" begin
        A = [1.0 2; 3 4]
        B = [1.0+im 2; 3 4]
        for (TriType, TriLayout, TriLayoutTrans) in ((UpperTriangular, UpperTriangularLayout, LowerTriangularLayout),
                              (UnitUpperTriangular, UnitUpperTriangularLayout, UnitLowerTriangularLayout),
                              (LowerTriangular, LowerTriangularLayout, UpperTriangularLayout),
                              (UnitLowerTriangular, UnitLowerTriangularLayout, UnitUpperTriangularLayout))
            @test MemoryLayout(TriType(A)) == TriLayout(DenseColumnMajor())
            @test MemoryLayout(TriType(transpose(A))) == TriLayout(DenseRowMajor())
            @test MemoryLayout(TriType(A')) == TriLayout(DenseRowMajor())
            @test MemoryLayout(transpose(TriType(A))) == TriLayoutTrans(DenseRowMajor())
            @test MemoryLayout(TriType(A)') == TriLayoutTrans(DenseRowMajor())

            @test MemoryLayout(TriType(B)) == TriLayout(DenseColumnMajor())
            @test MemoryLayout(TriType(transpose(B))) == TriLayout(DenseRowMajor())
            @test MemoryLayout(TriType(B')) == TriLayout(ConjLayout(DenseRowMajor()))
            @test MemoryLayout(transpose(TriType(B))) == TriLayoutTrans(DenseRowMajor())
            @test MemoryLayout(TriType(B)') == TriLayoutTrans(ConjLayout(DenseRowMajor()))
        end

        @test MemoryLayout(UpperTriangular(B)') == MemoryLayout(LowerTriangular(B'))
    end

    @testset "Fill and Vcat" begin
        @test MemoryLayout(Fill(1,10)) == FillLayout()
        @test MemoryLayout(Ones(10)) == FillLayout()
        @test MemoryLayout(Zeros(10)) == ZerosLayout()
        @test MemoryLayout(Vcat(Ones(10),Zeros(10))) == VcatLayout((FillLayout(), ZerosLayout()))
        @test MemoryLayout(Vcat([1.],Zeros(10))) == VcatLayout((DenseColumnMajor(), ZerosLayout()))
    end

    @testset "BroadcastArray" begin
        A = [1.0 2; 3 4]
        @test MemoryLayout(Add(A, Fill(0, (2, 2)), Zeros(2, 2))) ==
            BroadcastLayout(+, (DenseColumnMajor(), FillLayout(), ZerosLayout()))
    end
end


text:
module Car2DExamples

using OneVision
using OneVision: ℝ, 𝕋, ℕ, @kwdef, °, @set, @_, @asserteq
using ..Examples
using OneVision.NumericalIntegration
using OneVision.SymbolMaps
using OneVision.Examples
using Random
using StaticArrays
using AxisArrays
using Parameters
using AbstractPlotting
using Makie
using AbstractPlotting.MakieLayout
using Printf: @sprintf
using LinearAlgebra: norm
import ColorSchemes
import Dates
using Statistics: mean

include("common_2D.jl")
include("tracking_example.jl")
include("formation_example.jl")
include("interactive_example.jl")

end  # module Car2DExamples


text:
using Documenter
using Controlz

makedocs(
    root = joinpath(dirname(pathof(Controlz)), "..", "docs"),
    modules = [Controlz],
 #     source = "src",
    sitename = "Controlz.jl",
    clean = true,
    pages = ["Controlz" => "index.md",
             "Transfer Functions" => "tfs.md",
             "Simulation" => "sim.md",
             "Visualization" => "viz.md",
             "Conrol systems" => "controls.md",
             "FAQ" => "faq.md"],
    format = Documenter.HTML(assets = ["assets/flux.css"])
)

deploydocs(
    repo = "github.com/SimonEnsemble/Controlz.jl.git",
    # This is a link to the main repo and the master branch
    # target = "build",
    deps = Deps.pip("mkdocs", "mkdocs-material", "pymdown-extensions") # These are dependencies for the site, not the package
)


text:
export @ein_str, @ein, ein
"""
    ein"ij,jk -> ik"(A,B)

String macro interface which understands `numpy.einsum`'s notation.
Translates strings into `EinCode`-structs that can be called to evaluate
an `einsum`.
To control evaluation order, use parentheses - instead of an `EinCode`,
a `NestedEinsum` is returned which evaluates the expression
according to parens.
The valid character ranges for index-labels are `a-z` and `α-ω`.

# example

```jldoctest; setup = :(using OMEinsum)
julia> a, b, c = rand(10,10), rand(10,10), rand(10,1);

julia> ein"ij,jk,kl -> il"(a,b,c) ≈ ein"(ij,jk),kl -> il"(a,b,c) ≈ a * b * c
true
```
"""
macro ein_str(s::AbstractString)
    ein(s)
end

function ein(s::AbstractString)
    s = replace(replace(s, "\n" => ""), " "=>"")
    m = match(r"([\(\)a-z,α-ω]*)->([a-zα-ω]*)", s)
    m === nothing && throw(ArgumentError("invalid einsum specification $s"))
    sixs, siy = m.captures
    if '(' in sixs
        return parse_nested(sixs, collect(siy))
    else
        iy  = Tuple(siy)
        ixs = Tuple(Tuple(ix) for ix in split(sixs,','))
        return EinCode(ixs, iy)
    end
end

function (code::EinCode)(xs...; size_info=nothing)
    size_dict = get_size_dict(getixs(code), xs, size_info)
    einsum(code, xs, size_dict)
end

@doc raw"
    get_size_dict(ixs, xs, size_info=nothing)

return a dictionary that is used to get the size of an index-label
in the einsum-specification with input-indices `ixs` and tensors `xs` after
consistency within `ixs` and between `ixs` and `xs` has been verified.
"
function get_size_dict(ixs, xs, size_info=nothing)
    # check size of input tuples
    length(xs)<1 && error("empty input tensors")
    d = size_info === nothing ? Dict{promote_type(eltype.(ixs)...),Int}() : size_info
    length(ixs) != length(xs) && throw(ArgumentError("$(length(xs)) tensors labelled by $(length(ixs)) indices"))

    # check tensor orders
    foreach(ixs, xs) do ix, x
        length(ix) == ndims(x) || throw(
            ArgumentError("indices $ix invalid for tensor with ndims = $(ndims(x))"))
        for j = 1:length(ix)
            k = ix[j]
            if haskey(d, k)
                @assert size(x, j) == d[k] || throw(DimensionMismatch("$k = $(d[k]) or $(size(x,j))?"))
            else
                d[k] = size(x, j)
            end
        end
    end
    return d
end

using MacroTools
"""
    @ein A[i,k] := B[i,j] * C[j,k]     # A = B * C

Macro interface similar to that of other packages.

You may use numbers in place of letters for dummy indices, as in `@tensor`,
and need not name the output array. Thus `A = @ein [1,2] := B[1,ξ] * C[ξ,2]`
is equivalent to the above. This can also be written `A = ein"ij,jk -> ik"(B,C)`
using the numpy-style string macro.

# example

```jldoctest; setup = :(using OMEinsum)
julia> a, b = rand(2,2), rand(2,2);

julia> @ein c[i,k] := a[i,j] * b[j,k];

julia> c ≈ a * b
true

julia> c ≈ ein"ij,jk -> ik"(a,b)
true
```
"""
macro ein(exs...)
    _ein_macro(exs...)
end

primefix!(ind) = map!(i -> @capture(i, (j_)') ? Symbol(j, '′') : i, ind, ind)

function _ein_macro(ex; einsum=:einsum)
    @capture(ex, (left_ := right_)) || throw(ArgumentError("expected A[] := B[]... "))
    @capture(left, Z_[leftind__] | [leftind__] ) || throw(
        ArgumentError("can't understand LHS, expected A[i,j] etc."))
    if Z===nothing
        @gensym Z
    end
    primefix!(leftind)

    rightind, rightpairs = [], []
    @capture(right, *(factors__)) || (factors = Any[right])
    for fact in factors
        @capture(fact, A_[Aind__]) || return _nested_ein_macro(ex)
        primefix!(Aind)
        append!(rightind, Aind)
        push!(rightpairs, (A, Aind) )
    end
    unique!(rightind)
    isempty(setdiff(leftind, rightind)) || throw(
        ArgumentError("some indices appear only on the left"))

    lefttuple = Tuple(indexin(leftind, rightind))
    righttuples = [ Tuple(indexin(ind, rightind)) for (A, ind) in rightpairs ]
    rightnames = [ esc(A) for (A, ind) in rightpairs ]

    return :( $(esc(Z)) = $einsum( EinCode{($(righttuples...),), $lefttuple}(), ($(rightnames...),)) )
end

@doc raw"
    einsum(::EinCode{ixs, iy}, xs, size_dict) where {ixs, iy}

return the tensor that results from contracting the tensors `xs` according
to their indices `ixs`, where all indices that do not appear in the output `iy` are
summed over.
The result is permuted according to `out`.

- `ixs` - tuple of tuples of index-labels of the input-tensors `xs`

- `iy` - tuple of index-labels of the output-tensor

- `xs` - tuple of tensors

- `size_dict` - a dictionary that maps index-labels to their sizes

# example

```jldoctest; setup = :(using OMEinsum)
julia> a, b = rand(2,2), rand(2,2);

julia> einsum(EinCode((('i','j'),('j','k')),('i','k')), (a, b)) ≈ a * b
true

julia> einsum(EinCode((('i','j'),('j','k')),('k','i')), (a, b)) ≈ permutedims(a * b, (2,1))
true
```
"
@generated function einsum(code::EinCode{ixs, iy}, xs, size_dict::Dict{LT}) where {ixs, iy, LT}
    rule = match_rule(ixs, iy)
    if length(ixs) == 1
        :(einsum($rule, ixs[1], iy, xs[1], size_dict))
    else
        :(einsum($rule, ixs, iy, xs, size_dict))
    end
end

function einsum(code::EinCode{ixs, iy}, xs) where {ixs, iy}
    einsum(code, xs, get_size_dict(ixs, xs))
end

function dynamic_einsum(ixs, xs, iy; size_info=nothing)
    size_dict = get_size_dict(ixs, xs, size_info)
    dynamic_einsum(ixs, xs, iy, size_dict)
end

function dynamic_einsum(ixs, xs, iy, size_dict)
    rule = match_rule(ixs, iy)
    if length(ixs) == 1
        einsum(rule, ixs[1], iy, xs[1], size_dict)
    else
        einsum(rule, ixs, iy, (xs...,), size_dict)
    end
end

dynamic_einsum(::EinCode{ixs, iy}, xs; kwargs...) where {ixs, iy} = dynamic_einsum(ixs, xs, iy; kwargs...)

# the fallback
function einsum(::DefaultRule, ixs, iy, xs, size_dict)
    @debug "DefaultRule loop_einsum" ixs => iy size.(xs)
    loop_einsum(EinCode{ixs, iy}(), (xs...,), size_dict)
end


text:
---
title: shutdown
tags: shellcommand,advance
---

Shutdown computer in a duration given by the user and message to be displayed stating the reason for shutdown.

Use variable `Sys.KERNEL` to handle differences between operating systems.
Use `@static` macro at beginning of expression for mutually exclusive subsets of isunix (Optional).
Use `run` function which take an object of type `osCommand` defined by enclosing a command string in backticks and execute
the required command for the operating system based on the conditions.

```jl
function shutdown(time,message)

  #Operating system commands
  windowsCommand=`cmd /c shutdown/s /t $time /c $message`
  linuxCommand=`sudo shutdown $time $message`
  appleCommand=`sudo shutdown -h +$time $message`
  freebsdCommand=`shutdown -h +$time $message`

  #Condition before command is executed based on OS
  if Sys.iswindows()==true
    run(windowsCommand)
  elseif Sys.islinux()==true
    run(linuxCommand)
  elseif Sys.isapple()==true
    run(appleCommand)
  elseif Sys.isfreebsd()==true
    run(freebsdCommand)
  else
    println("Operating system not supported!")
   end
end
```

```jl
shutdown(400,"Fatal System Error! Blue Screen Loading..") # Dialogue box appear showing "Fatal System Error! Blue Screen Loading" and computer shutdown in 400sec.
```


text:
using MPI
using CLIMA
using CLIMA.Mesh.Topologies
using CLIMA.Mesh.Grids
using CLIMA.DGmethods
using CLIMA.DGmethods.NumericalFluxes
using CLIMA.MPIStateArrays
using CLIMA.LowStorageRungeKuttaMethod
using CLIMA.ODESolvers
using CLIMA.GenericCallbacks
using CLIMA.Atmos
using CLIMA.VariableTemplates
using CLIMA.MoistThermodynamics
using CLIMA.PlanetParameters
using LinearAlgebra
using StaticArrays
using Logging, Printf, Dates
using CLIMA.VTK

using CLIMA.Atmos: vars_state, vars_aux

using Random
const seed = MersenneTwister(0)

@static if haspkg("CuArrays")
  using CUDAdrv
  using CUDAnative
  using CuArrays
  CuArrays.allowscalar(false)
  const ArrayTypes = (CuArray,)
else
  const ArrayTypes = (Array,)
end

if !@isdefined integration_testing
  const integration_testing =
    parse(Bool, lowercase(get(ENV,"JULIA_CLIMA_INTEGRATION_TESTING","false")))
end

"""
  Initial Condition for DYCOMS_RF01 LES
@article{doi:10.1175/MWR2930.1,
author = {Stevens, Bjorn and Moeng, Chin-Hoh and Ackerman,
          Andrew S. and Bretherton, Christopher S. and Chlond,
          Andreas and de Roode, Stephan and Edwards, James and Golaz,
          Jean-Christophe and Jiang, Hongli and Khairoutdinov,
          Marat and Kirkpatrick, Michael P. and Lewellen, David C. and Lock, Adrian and
          Maeller, Frank and Stevens, David E. and Whelan, Eoin and Zhu, Ping},
title = {Evaluation of Large-Eddy Simulations via Observations of Nocturnal Marine Stratocumulus},
journal = {Monthly Weather Review},
volume = {133},
number = {6},
pages = {1443-1462},
year = {2005},
doi = {10.1175/MWR2930.1},
URL = {https://doi.org/10.1175/MWR2930.1},
eprint = {https://doi.org/10.1175/MWR2930.1}
}
"""
function Initialise_DYCOMS!(state::Vars, aux::Vars, (x,y,z), t)
  DT         = eltype(state)
  xvert::DT  = z

  epsdv::DT     = molmass_ratio
  q_tot_sfc::DT = 8.1e-3
  Rm_sfc::DT    = gas_constant_air(PhasePartition(q_tot_sfc))
  ρ_sfc::DT     = 1.22
  P_sfc::DT     = 1.0178e5
  T_BL::DT      = 285.0
  T_sfc::DT     = P_sfc/(ρ_sfc * Rm_sfc);

  q_liq::DT      = 0
  q_ice::DT      = 0
  zb::DT         = 600
  zi::DT         = 840
  dz_cloud       = zi - zb
  q_liq_peak::DT = 4.5e-4

  if xvert > zb && xvert <= zi
    q_liq = (xvert - zb)*q_liq_peak/dz_cloud
  end
  if ( xvert <= zi)
    θ_liq  = DT(289)
    q_tot  = DT(8.1e-3)
  else
    θ_liq = DT(297.5) + (xvert - zi)^(DT(1/3))
    q_tot = DT(1.5e-3)
  end

  q_pt = PhasePartition(q_tot, q_liq, DT(0))
  Rm    = gas_constant_air(q_pt)
  cpm   = cp_m(q_pt)
  #Pressure
  H = Rm_sfc * T_BL / grav;
  P = P_sfc * exp(-xvert/H);
  #Exner
  exner_dry = exner_given_pressure(P, PhasePartition(DT(0)))
  #Temperature
  T             = exner_dry*θ_liq + LH_v0*q_liq/(cpm*exner_dry);
  #Density
  ρ             = P/(Rm*T);
  #Potential Temperature
  θv     = virtual_pottemp(T, P, q_pt)
  # energy definitions
  u, v, w     = DT(7), DT(-5.5), DT(0)
  U           = ρ * u
  V           = ρ * v
  W           = ρ * w
  e_kin       = DT(1//2) * (u^2 + v^2 + w^2)
  e_pot       = grav * xvert
  E           = ρ * total_energy(e_kin, e_pot, T, q_pt)
  state.ρ     = ρ
  state.ρu    = SVector(U, V, W)
  state.ρe    = E
  state.moisture.ρq_tot = ρ * q_tot
end

function run(mpicomm, ArrayType, dim, topl, N, timeend, DT, dt, C_smag, LHF, SHF, C_drag, zmax, zsponge)

  grid = DiscontinuousSpectralElementGrid(topl,
                                          FloatType = DT,
                                          DeviceArray = ArrayType,
                                          polynomialorder = N,
                                         )
  model = AtmosModel(FlatOrientation(),
                     NoReferenceState(),
                     SmagorinskyLilly{DT}(C_smag),
                     EquilMoist(),
                     StevensRadiation{DT}(85, 1, 840, 1.22, 3.75e-6, 70, 22),
                     (Gravity(),
                      RayleighSponge{DT}(zmax, zsponge, 1),
                      Subsidence(),
                      GeostrophicForcing{DT}(7.62e-5, 7, -5.5)),
                     DYCOMS_BC{DT}(C_drag, LHF, SHF),
                     Initialise_DYCOMS!)

  dg = DGModel(model,
               grid,
               Rusanov(),
               CentralNumericalFluxDiffusive(),
               CentralGradPenalty())

  Q = init_ode_state(dg, DT(0))

  lsrk = LSRK54CarpenterKennedy(dg, Q; dt = dt, t0 = 0)

  eng0 = norm(Q)
  @info @sprintf """Starting
  norm(Q₀) = %.16e""" eng0

  # Set up the information callback
  starttime = Ref(now())
  cbinfo = GenericCallbacks.EveryXWallTimeSeconds(60, mpicomm) do (s=false)
    if s
      starttime[] = now()
    else
      energy = norm(Q)
      @info @sprintf("""Update
                     simtime = %.16e
                     runtime = %s
                     norm(Q) = %.16e""", ODESolvers.gettime(lsrk),
                     Dates.format(convert(Dates.DateTime,
                                          Dates.now()-starttime[]),
                                  Dates.dateformat"HH:MM:SS"),
                     energy)
    end
  end

  step = [0]
    cbvtk = GenericCallbacks.EveryXSimulationSteps(5000) do (init=false)
    mkpath("./vtk-dycoms/")
    outprefix = @sprintf("./vtk-dycoms/dycoms_%dD_mpirank%04d_step%04d", dim,
                           MPI.Comm_rank(mpicomm), step[1])
    @debug "doing VTK output" outprefix
    writevtk(outprefix, Q, dg, flattenednames(vars_state(model,DT)),
             dg.auxstate, flattenednames(vars_aux(model,DT)))

    step[1] += 1
    nothing
  end

  solve!(Q, lsrk; timeend=timeend, callbacks=(cbinfo, cbvtk))

  # Print some end of the simulation information
  engf = norm(Q)
  Qe = init_ode_state(dg, DT(timeend))

  engfe = norm(Qe)
  errf = euclidean_distance(Q, Qe)
  @info @sprintf """Finished
  norm(Q)                 = %.16e
  norm(Q) / norm(Q₀)      = %.16e
  norm(Q) - norm(Q₀)      = %.16e
  norm(Q - Qe)            = %.16e
  norm(Q - Qe) / norm(Qe) = %.16e
  """ engf engf/eng0 engf-eng0 errf errf / engfe
  engf/eng0
end

using Test
let
  MPI.Initialized() || MPI.Init()
  mpicomm = MPI.COMM_WORLD
  ll = uppercase(get(ENV, "JULIA_LOG_LEVEL", "INFO"))
  loglevel = ll == "DEBUG" ? Logging.Debug :
    ll == "WARN"  ? Logging.Warn  :
    ll == "ERROR" ? Logging.Error : Logging.Info
  logger_stream = MPI.Comm_rank(mpicomm) == 0 ? stderr : devnull
  global_logger(ConsoleLogger(logger_stream, loglevel))
  @static if haspkg("CUDAnative")
      device!(MPI.Comm_rank(mpicomm) % length(devices()))
  end
  @testset "$(@__FILE__)" for ArrayType in ArrayTypes
    # Problem type
    DT = Float64
    # DG polynomial order
    N = 4
    # SGS Filter constants
    C_smag = DT(0.15)
    LHF    = DT(115)
    SHF    = DT(15)
    C_drag = DT(0.0011)
    # User defined domain parameters
    brickrange = (grid1d(0, 2000, elemsize=DT(50)*N),
                  grid1d(0, 2000, elemsize=DT(50)*N),
                  grid1d(0, 1500, elemsize=DT(20)*N))
    zmax = brickrange[3][end]
    zsponge = DT(0.75 * zmax)

    topl = StackedBrickTopology(mpicomm, brickrange,
                                periodicity = (true, true, false),
                                boundary=((0,0),(0,0),(1,2)))
    dt = 0.02
    timeend = 100dt
    dim = 3
    @info (ArrayType, DT, dim)
    result = run(mpicomm, ArrayType, dim, topl,
                 N, timeend, DT, dt, C_smag, LHF, SHF, C_drag, zmax, zsponge)
    @test result ≈ DT(0.9999737848359238)
  end
end

#nothing

max_stars_repo_name: akshaysridhar/CLIMA

max_issues_repo_name: akshaysridhar/CLIMA

max_forks_repo_name: akshaysridhar/CLIMA


text:
# download 'atlas_data.csv' from "https://mega.nz/file/7sxx3ajJ#cgPEn1LGSutiJTCrNeeoVbuHGclxfYKZL2IhqhzpLHA"
# and put it in the same folder as this script

import MLJFlux
using DataFrames, DataFramesMeta, CSV, Alert, ProgressMeter, Plots, Flux, MLJ,  Random
using Flux: onehotbatch, onecold, @epochs, Data.DataLoader, Optimiser
using Chain: @chain
using MLDataUtils: splitobs, shuffleobs
using StatsBase: standardize, ZScoreTransform

function build_model(input, layers, output; activation = relu, use_softmax = true, use_last_activation = false)
	f = []
	in_layer = input

	for out_layer in layers
		append!(f, [Dense(in_layer, out_layer, activation)])
		in_layer = out_layer
	end

	if use_last_activation
		append!(f, [Dense(in_layer, output, activation)])
	else
		append!(f, [Dense(in_layer, output)])
	end

	if use_softmax
		append!(f, [softmax])
	end

	Chain(f...)
end

function score_accuracy(X_output, y_output, classes = [1, 0])
	X = onecold(X_output, classes)
	y = onecold(y_output, classes)
	comparison = X .== y
	return sum(a -> a > 0, comparison) / length(comparison)
end

function flux_run_model(df, dataset_frac, hidden_layers, input_loss, opt, n_epochs, batchsize, λ = 0.0001, tol = 1e-4, iter_tol = 10)
	# converts the 'Label' column to a matrix of 0 and 1
	# compatible with the neural newtork
	X, y = select(df, Not(:Label)), @chain df begin
		select(_, :Label)
		Flux.onehotbatch(_.Label, ["s", "b"])
	end

	N_input  = length(names(X))
	N_output = size(y, 1)

	X = transpose(standardize(ZScoreTransform, Matrix(X)))
	X_train, X_test = splitobs(X, at = dataset_frac)
	y_train, y_test = splitobs(y, at = dataset_frac)

	model = build_model(N_input, hidden_layers, N_output)
	loss(a, b) = input_loss(model(a), b)

	ps = Flux.params(model)

	loader = DataLoader(
		(X_train, y_train),
		batchsize = batchsize,
		shuffle = false
	)

	tol_counter = 0
	loss_train_values = []
	loss_test_values = []
	acc_train_values = []
	acc_test_values = []

	prev = loss(X_train, y_train)
	p = Progress(n_epochs, dt = 1)
	generate_showvalues(i, n) = () -> [(:current_epoch, i), (:tot_epochs, n)]

	for i in 1:n_epochs
		Flux.train!(loss, ps, loader, Optimiser(WeightDecay(λ), opt))

		l_train = loss(X_train, y_train)
		l_test = loss(X_test, y_test)
		append!(loss_train_values, l_train)
		append!(loss_test_values, l_test)
		append!(acc_train_values, score_accuracy(model(X_train), y_train))
		append!(acc_test_values, score_accuracy(model(X_test), y_test))

		if abs(l_train - prev) < tol
			tol_counter += 1
		else
			tol_counter = 0
		end

		if tol_counter == iter_tol
			break
		end

		prev = l_train

		ProgressMeter.next!(p; showvalues = generate_showvalues(i, n_epochs))
	end

	if tol_counter == iter_tol
		@warn "Terminated due to having reached tol = $tol for $iter_tol times in a row"
	end

	return loss_train_values, loss_test_values, acc_train_values, acc_test_values
end

mutable struct MyNetwork{F <: Function} <: MLJFlux.Builder
    layers :: Vector{Int64}
	activation :: F
	use_softmax :: Bool
	use_last_activation :: Bool
end

function MLJFlux.build(nn::MyNetwork, n_in, n_out)
	layers = nn.layers
	activation = nn.activation
	use_softmax = nn.use_softmax
	use_last_activation = nn.use_last_activation

    f = []
	in_layer = n_in

	for out_layer in layers
		append!(f, [Dense(in_layer, out_layer, activation)])
		in_layer = out_layer
	end

	if use_last_activation
		append!(f, [Dense(in_layer, n_out, activation)])
	else
		append!(f, [Dense(in_layer, n_out)])
	end

	if use_softmax
		append!(f, [softmax])
	end

	Chain(f...)
end

function mljflux_run_model(df, holdout_frac, hidden_layers, input_loss, opt, n_epochs, batchsize, λ = 0.0001, α = 0.0)
	y, X = unpack(df, ==(:Label), colname -> true)

	X = coerce(X, Count => Continuous)
	y = coerce(y, Multiclass)

	NeuralNetworkClassifier = @load NeuralNetworkClassifier

	clf = NeuralNetworkClassifier(
		builder = MyNetwork(
			hidden_layers,
			relu,
			false,
			false
		),
		finaliser = softmax,
		optimiser = opt,
		loss = input_loss,
		epochs = n_epochs,
		batch_size = batchsize,
		lambda = λ,
		alpha = α,
		optimiser_changes_trigger_retraining = false
	)

	mach = machine(clf, X, y)

	evaluate!(
		mach,
		resampling = Holdout(
			fraction_train = holdout_frac
		),
		operation = predict_mode,
		measure = accuracy,
		verbosity = 2
	)
end

# assuming 'atlas_data.csv' is in the same folder
filepath = joinpath(pwd(), "atlas_data.csv")

# after this, df contains 250k rows and 31 columns (labels included)
# this simply removes useless columns from the dataset
# and select data tagged with KaggleSet = "t"
df = @chain begin
	CSV.read(filepath, DataFrame)
	@where(_, :KaggleSet .== "t")
	select(_, Not([:Weight, :EventId, :KaggleSet, :KaggleWeight]))
end

# PARAMETERS ------------------------------------------------------------------------------------------------------------
hidden_layers = [20, 10, 5]

optimizer = "adam" # or "momentum"

η  = 1e-3 # learning rate
ρ  = 0.99 # for momentum opt
β₁ = 0.9  # for adam opt
β₂ = 0.999 # for adam opt
λ  = 0.0002 # L2 regularization
α  = 0.0 # used only in MLJ's NeuralNetworkClassifier

# for Flux model only
tol = 1e-4
iter_tol = 10

batchsize = 200
n_epochs = 10

# rows of dataset to use. Set it to `length(names(df))` to use the entire dataset
N_rows = 20000

# fraction of dataset for training (in Flux) or holdout fraction (in MLJFlux)
frac = 0.7

opt_dict = Dict([
	("momentum", Momentum(η, ρ)),
	("adam", ADAM(η, (β₁, β₂)))
])

#
# FLUX-----------------------------------------------------------------------------------------------------------------
#

loss_train_values, loss_test_values, acc_train_values, acc_test_values = @alert "Flux finished" @time flux_run_model(
	df[1:N_rows, :],
	frac,
	hidden_layers,
	Flux.Losses.crossentropy,
	opt_dict[optimizer],
	n_epochs,
	batchsize,
	λ,
	tol,
	iter_tol
)

p_loss = plot(
	loss_train_values,
	title = "Loss",
	xlabel = "Epoch",
	ylabel = "Loss",
	label = "Training",
)

plot!(
	p_loss,
	loss_test_values,
	label = "Testing"
)

p_acc = plot(
	acc_train_values,
	title = "Accuracy",
	xlabel = "Epoch",
	ylabel = "Accuracy",
	label = "Training",
	legend = :bottomright
)

plot!(
	p_acc,
	acc_test_values,
	label = "Testing"
)

#
# MLJFlux ---------------------------------------------------------------------------------------------------
#

ev = @alert "MLJFlux Finished" @time mljflux_run_model(
	df[1:N_rows, :],
	frac,
	hidden_layers,
	Flux.crossentropy,
	opt_dict[optimizer],
	n_epochs,
	batchsize,
	λ,
	α
)