text:
module NL

import MathOptInterface
const MOI = MathOptInterface

include("NLExpr.jl")

### ============================================================================
### Nonlinear constraints
### ============================================================================

struct _NLConstraint
    lower::Float64
    upper::Float64
    opcode::Int
    expr::_NLExpr
end

"""
    _NLConstraint(expr::Expr, bound::MOI.NLPBoundsPair)

Convert a constraint in the form of a `expr` into a `_NLConstraint` object.

See `MOI.constraint_expr` for details on the format.

As a validation step, the right-hand side of each constraint must be a constant
term that is given by the `bound`. (If the constraint is an interval constraint,
both the left-hand and right-hand sides must be constants.)

The six NL constraint types are:

    l <= g(x) <= u : 0
         g(x) >= l : 1
         g(x) <= u : 2
         g(x)      : 3  # We don't support this
         g(x) == c : 4
     x ⟂ g(x)      : 5  # TODO(odow): Complementarity constraints
"""
function _NLConstraint(expr::Expr, bound::MOI.NLPBoundsPair)
    if expr.head == :comparison
        @assert length(expr.args) == 5
        if !(expr.args[1] ≈ bound.lower && bound.upper ≈ expr.args[5])
            _warn_invalid_bound(expr, bound)
        end
        return _NLConstraint(
            expr.args[1],
            expr.args[5],
            0,
            _NLExpr(expr.args[3]),
        )
    else
        @assert expr.head == :call
        @assert length(expr.args) == 3
        if expr.args[1] == :(<=)
            if !(-Inf ≈ bound.lower && bound.upper ≈ expr.args[3])
                _warn_invalid_bound(expr, bound)
            end
            return _NLConstraint(-Inf, expr.args[3], 1, _NLExpr(expr.args[2]))
        elseif expr.args[1] == :(>=)
            if !(expr.args[3] ≈ bound.lower && bound.upper ≈ Inf)
                _warn_invalid_bound(expr, bound)
            end
            return _NLConstraint(expr.args[3], Inf, 2, _NLExpr(expr.args[2]))
        else
            @assert expr.args[1] == :(==)
            if !(expr.args[3] ≈ bound.lower ≈ bound.upper)
                _warn_invalid_bound(expr, bound)
            end
            return _NLConstraint(
                expr.args[3],
                expr.args[3],
                4,
                _NLExpr(expr.args[2]),
            )
        end
    end
end

function _warn_invalid_bound(expr::Expr, bound::MOI.NLPBoundsPair)
    return @warn(
        "Invalid bounds detected in nonlinear constraint. Expected " *
        "`$(bound.lower) <= g(x) <= $(bound.upper)`, but got the constraint " *
        "$(expr)",
    )
end

### ============================================================================
### Nonlinear models
### ============================================================================

@enum(_VariableType, _BINARY, _INTEGER, _CONTINUOUS)

mutable struct _VariableInfo
    # Variable lower bound.
    lower::Float64
    # Variable upper bound.
    upper::Float64
    # Whether variable is binary or integer.
    type::_VariableType
    # Primal start of the variable.
    start::Union{Float64,Nothing}
    # Number of constraints that the variable appears in.
    jacobian_count::Int
    # If the variable appears in the objective.
    in_nonlinear_objective::Bool
    # If the variable appears in a nonlinear constraint.
    in_nonlinear_constraint::Bool
    # The 0-indexed column of the variable. Computed right at the end.
    order::Int
    function _VariableInfo()
        return new(-Inf, Inf, _CONTINUOUS, nothing, 0, false, false, 0)
    end
end

"""
    Model()

Create a new Optimizer object.
"""
mutable struct Model <: MOI.ModelLike
    # Store MOI.Name().
    name::String
    # The objective expression.
    f::_NLExpr
    sense::MOI.OptimizationSense
    # Number of nonlinear constraints in NLPBlock
    nlpblock_dim::Int
    # A vector of nonlinear constraints
    g::Vector{_NLConstraint}
    # A vector of linear constraints
    h::Vector{_NLConstraint}
    # A dictionary of info for the variables.
    x::Dict{MOI.VariableIndex,_VariableInfo}
    # A struct to help sort the mess that is variable ordering in NL files.
    types::Vector{Vector{MOI.VariableIndex}}
    # A vector of the final ordering of the variables.
    order::Vector{MOI.VariableIndex}

    function Model()
        return new(
            "",
            _NLExpr(false, _NLTerm[], Dict{MOI.VariableIndex,Float64}(), 0.0),
            MOI.FEASIBILITY_SENSE,
            0,
            _NLConstraint[],
            _NLConstraint[],
            Dict{MOI.VariableIndex,_VariableInfo}(),
            [MOI.VariableIndex[] for _ in 1:9],
            MOI.VariableIndex[],
        )
    end
end

Base.show(io::IO, ::Model) = print(io, "An AMPL (.nl) model")

MOI.get(model::Model, ::MOI.SolverName) = "AmplNLWriter"

MOI.supports(::Model, ::MOI.NLPBlock) = true

MOI.supports(::Model, ::MOI.Name) = true
MOI.get(model::Model, ::MOI.Name) = model.name
MOI.set(model::Model, ::MOI.Name, name::String) = (model.name = name)

function MOI.empty!(model::Model)
    model.f = _NLExpr(false, _NLTerm[], Dict{MOI.VariableIndex,Float64}(), 0.0)
    empty!(model.g)
    model.nlpblock_dim = 0
    empty!(model.h)
    empty!(model.x)
    for i in 1:9
        empty!(model.types[i])
    end
    empty!(model.order)
    return
end

function MOI.is_empty(model::Model)
    return isempty(model.g) && isempty(model.h) && isempty(model.x)
end

const _SCALAR_FUNCTIONS = Union{
    MOI.VariableIndex,
    MOI.ScalarAffineFunction{Float64},
    MOI.ScalarQuadraticFunction{Float64},
}

const _SCALAR_SETS = Union{
    MOI.LessThan{Float64},
    MOI.GreaterThan{Float64},
    MOI.EqualTo{Float64},
    MOI.Interval{Float64},
}

function MOI.supports_constraint(
    ::Model,
    ::Type{<:_SCALAR_FUNCTIONS},
    ::Type{<:_SCALAR_SETS},
)
    return true
end

function MOI.supports_constraint(
    ::Model,
    ::Type{MOI.VariableIndex},
    ::Type{<:Union{MOI.ZeroOne,MOI.Integer}},
)
    return true
end

MOI.supports(::Model, ::MOI.ObjectiveSense) = true
MOI.supports(::Model, ::MOI.ObjectiveFunction{<:_SCALAR_FUNCTIONS}) = true

# ==============================================================================

function MOI.supports(
    ::Model,
    ::MOI.VariablePrimalStart,
    ::Type{MOI.VariableIndex},
)
    return true
end

function MOI.set(
    model::Model,
    ::MOI.VariablePrimalStart,
    x::MOI.VariableIndex,
    v::Union{Nothing,Real},
)
    model.x[x].start = v === nothing ? nothing : convert(Float64, v)::Float64
    return
end

function MOI.get(model::Model, ::MOI.VariablePrimalStart, x::MOI.VariableIndex)
    return model.x[x].start
end

# ==============================================================================

struct _LinearNLPEvaluator <: MOI.AbstractNLPEvaluator end
MOI.initialize(::_LinearNLPEvaluator, ::Vector{Symbol}) = nothing

function MOI.copy_to(dest::Model, model::MOI.ModelLike)
    if !MOI.is_empty(dest)
        MOI.empty!(dest)
    end
    mapping = MOI.Utilities.IndexMap()
    nlp_block =
        MOI.NLPBlockData(MOI.NLPBoundsPair[], _LinearNLPEvaluator(), false)
    for attr in MOI.get(model, MOI.ListOfModelAttributesSet())
        if attr == MOI.NLPBlock()
            nlp_block = MOI.get(model, MOI.NLPBlock())
            if !(:ExprGraph in MOI.features_available(nlp_block.evaluator))
                error(
                    "Unable to use AmplNLWriter because the nonlinear " *
                    "evaluator does not supply expression graphs.",
                )
            end
        elseif attr == MOI.ObjectiveSense()
            dest.sense = MOI.get(model, MOI.ObjectiveSense())
        elseif attr isa MOI.ObjectiveFunction
            dest.f = _NLExpr(MOI.get(model, attr))
        else
            throw(MOI.UnsupportedAttribute(attr))
        end
    end
    MOI.initialize(nlp_block.evaluator, [:ExprGraph])
    if nlp_block.has_objective  # Nonlinear objective takes precedence.
        dest.f = _NLExpr(MOI.objective_expr(nlp_block.evaluator))
    end
    for (i, bound) in enumerate(nlp_block.constraint_bounds)
        push!(
            dest.g,
            _NLConstraint(MOI.constraint_expr(nlp_block.evaluator, i), bound),
        )
    end
    dest.nlpblock_dim = length(dest.g)
    x_src = MOI.get(model, MOI.ListOfVariableIndices())
    for x in x_src
        dest.x[x] = _VariableInfo()
        mapping[x] = x
    end
    MOI.Utilities.pass_attributes(dest, model, mapping, x_src)
    resize!(dest.order, length(dest.x))
    # Now deal with the normal MOI constraints.
    for (F, S) in MOI.get(model, MOI.ListOfConstraintTypesPresent())
        if !MOI.supports_constraint(dest, F, S)
            throw(MOI.UnsupportedConstraint{F,S}())
        end
        _process_constraint(dest, model, F, S, mapping)
    end
    # Correct bounds of binary variables. Mainly because AMPL doesn't have the
    # concept of binary nonlinear variables, but it does have binary linear
    # variables! How annoying.
    for (_, v) in dest.x
        if v.type == _BINARY
            v.lower = max(0.0, v.lower)
            v.upper = min(1.0, v.upper)
        end
    end
    # Jacobian counts. The zero terms for nonlinear constraints should have
    # been added when the expression was constructed.
    for g in dest.g, v in keys(g.expr.linear_terms)
        dest.x[v].jacobian_count += 1
    end
    for h in dest.h, v in keys(h.expr.linear_terms)
        dest.x[v].jacobian_count += 1
    end
    # Now comes the confusing part.
    #
    # AMPL, in all its wisdom, orders variables in a _very_ specific way.
    # The only hint in "Writing NL files" is the line "Variables are ordered as
    # described in Tables 3 and 4 of [5]".
    #
    # Reading these
    #
    # https://cfwebprod.sandia.gov/cfdocs/CompResearch/docs/nlwrite20051130.pdf
    # https://ampl.com/REFS/hooking2.pdf
    #
    # leads us to the following order
    #
    # 1) Continuous variables that appear in a
    #       nonlinear objective AND a nonlinear constraint
    # 2) Discrete variables that appear in a
    #       nonlinear objective AND a nonlinear constraint
    # 3) Continuous variables that appear in a
    #       nonlinear constraint, but NOT a nonlinear objective
    # 4) Discrete variables that appear in a
    #       nonlinear constraint, but NOT a nonlinear objective
    # 5) Continuous variables that appear in a
    #       nonlinear objective, but NOT a nonlinear constraint
    # 6) Discrete variables that appear in a
    #       nonlinear objective, but NOT a nonlinear constraint
    # 7) Continuous variables that DO NOT appear in a
    #       nonlinear objective or a nonlinear constraint
    # 8) Binary variables that DO NOT appear in a
    #       nonlinear objective or a nonlinear constraint
    # 9) Integer variables that DO NOT appear in a
    #       nonlinear objective or a nonlinear constraint
    #
    # Yes, nonlinear variables are broken into continuous/discrete, but linear
    # variables are partitioned into continuous, binary, and integer. (See also,
    # the need to modify bounds for binary variables.)
    if !dest.f.is_linear
        for x in keys(dest.f.linear_terms)
            dest.x[x].in_nonlinear_objective = true
        end
        for x in dest.f.nonlinear_terms
            if x isa MOI.VariableIndex
                dest.x[x].in_nonlinear_objective = true
            end
        end
    end
    for con in dest.g
        for x in keys(con.expr.linear_terms)
            dest.x[x].in_nonlinear_constraint = true
        end
        for x in con.expr.nonlinear_terms
            if x isa MOI.VariableIndex
                dest.x[x].in_nonlinear_constraint = true
            end
        end
    end
    types = dest.types
    for (x, v) in dest.x
        if v.in_nonlinear_constraint && v.in_nonlinear_objective
            push!(v.type == _CONTINUOUS ? types[1] : types[2], x)
        elseif v.in_nonlinear_constraint
            push!(v.type == _CONTINUOUS ? types[3] : types[4], x)
        elseif v.in_nonlinear_objective
            push!(v.type == _CONTINUOUS ? types[5] : types[6], x)
        elseif v.type == _CONTINUOUS
            push!(types[7], x)
        elseif v.type == _BINARY
            push!(types[8], x)
        else
            @assert v.type == _INTEGER
            push!(types[9], x)
        end
    end
    # However! Don't let Tables 3 and 4 fool you, because the ordering actually
    # depends on whether the number of nonlinear variables in the objective only
    # is _strictly_ greater than the number of nonlinear variables in the
    # constraints only. Quoting:
    #
    #   For all versions, the first nlvc variables appear nonlinearly in at
    #   least one constraint. If nlvo > nlvc, the first nlvc variables may or
    #   may not appear nonlinearly in an objective, but the next nlvo – nlvc
    #   variables do appear nonlinearly in at least one objective. Otherwise
    #   all of the first nlvo variables appear nonlinearly in an objective.
    #
    # However, even this is slightly incorrect, because I think it should read
    # "For all versions, the first nlvb variables appear nonlinearly." The "nlvo
    # - nlvc" part is also clearly incorrect, and should probably read "nlvo -
    # nlvb."
    #
    # It's a bit confusing, so here is the relevant code from Couenne:
    # https://github.com/coin-or/Couenne/blob/683c5b305d78a009d59268a4bca01e0ad75ebf02/src/readnl/readnl.cpp#L76-L87
    #
    # They interpret this paragraph to mean the switch on nlvo > nlvc determines
    # whether the next block of variables are the ones that appear in the
    # objective only, or the constraints only.
    #
    # That makes sense as a design choice, because you can read them in two
    # contiguous blocks.
    #
    # Essentially, what all this means is if !(nlvo > nlvc), then swap 3-4 for
    # 5-6 in the variable order.
    nlvc = length(types[3]) + length(types[4])
    nlvo = length(types[5]) + length(types[6])
    order_i = if nlvo > nlvc
        [1, 2, 3, 4, 5, 6, 7, 8, 9]
    else
        [1, 2, 5, 6, 3, 4, 7, 8, 9]
    end
    # Now we can order the variables.
    n = 0
    for i in order_i
        # Since variables come from a dictionary, there may be differences in
        # the order depending on platform and Julia version. Sort by creation
        # time for consistency.
        for x in sort!(types[i]; by = y -> y.value)
            dest.x[x].order = n
            dest.order[n+1] = x
            n += 1
        end
    end
    return mapping
end

_set_to_bounds(set::MOI.Interval) = (0, set.lower, set.upper)
_set_to_bounds(set::MOI.LessThan) = (1, -Inf, set.upper)
_set_to_bounds(set::MOI.GreaterThan) = (2, set.lower, Inf)
_set_to_bounds(set::MOI.EqualTo) = (4, set.value, set.value)

function _process_constraint(
    dest::Model,
    model,
    ::Type{F},
    ::Type{S},
    mapping,
) where {F,S}
    ci_src = MOI.get(model, MOI.ListOfConstraintIndices{F,S}())
    for ci in ci_src
        f = MOI.get(model, MOI.ConstraintFunction(), ci)
        s = MOI.get(model, MOI.ConstraintSet(), ci)
        op, l, u = _set_to_bounds(s)
        con = _NLConstraint(l, u, op, _NLExpr(f))
        if isempty(con.expr.linear_terms) && isempty(con.expr.nonlinear_terms)
            if !(l <= con.expr.constant <= u)
                error(
                    "Malformed constraint. There are no variables and the " *
                    "function constant $(con.expr.constant) is not in [$l, $u]",
                )
            end
            # Just use a placeholder for the constraint index. It's not going to
            # be used.
            mapping[ci] = MOI.ConstraintIndex{F,S}(-abs(ci.value))
        elseif con.expr.is_linear
            push!(dest.h, con)
            mapping[ci] = MOI.ConstraintIndex{F,S}(length(dest.h))
        else
            push!(dest.g, con)
            mapping[ci] = MOI.ConstraintIndex{F,S}(length(dest.g))
        end
    end
    MOI.Utilities.pass_attributes(dest, model, mapping, ci_src)
    return
end

function _process_constraint(
    dest::Model,
    model,
    F::Type{MOI.VariableIndex},
    S::Type{<:_SCALAR_SETS},
    mapping,
)
    ci_src = MOI.get(model, MOI.ListOfConstraintIndices{F,S}())
    for ci in ci_src
        mapping[ci] = ci
        f = MOI.get(model, MOI.ConstraintFunction(), ci)
        s = MOI.get(model, MOI.ConstraintSet(), ci)
        _, l, u = _set_to_bounds(s)
        if l > -Inf
            dest.x[f].lower = l
        end
        if u < Inf
            dest.x[f].upper = u
        end
    end
    MOI.Utilities.pass_attributes(dest, model, mapping, ci_src)
    return
end

function _process_constraint(
    dest::Model,
    model,
    F::Type{MOI.VariableIndex},
    S::Type{<:Union{MOI.ZeroOne,MOI.Integer}},
    mapping,
)
    ci_src = MOI.get(model, MOI.ListOfConstraintIndices{F,S}())
    for ci in ci_src
        mapping[ci] = ci
        f = MOI.get(model, MOI.ConstraintFunction(), ci)
        dest.x[f].type = S == MOI.ZeroOne ? _BINARY : _INTEGER
    end
    MOI.Utilities.pass_attributes(dest, model, mapping, ci_src)
    return
end

function _str(x::Float64)
    if isinteger(x) && (typemin(Int) <= x <= typemax(Int))
        return string(round(Int, x))
    end
    return string(x)
end

_write_term(io, x::Float64, ::Any) = println(io, "n", _str(x))
_write_term(io, x::Int, ::Any) = println(io, "o", x)
function _write_term(io, x::MOI.VariableIndex, nlmodel)
    return println(io, "v", nlmodel.x[x].order)
end

_is_nary(x::Int) = x in _NARY_OPCODES
_is_nary(x) = false

function _write_nlexpr(io::IO, expr::_NLExpr, nlmodel::Model)
    if expr.is_linear || length(expr.nonlinear_terms) == 0
        # If the expression is linear, just write out the constant term.
        _write_term(io, expr.constant, nlmodel)
        return
    end
    if !iszero(expr.constant)
        # If the constant term is non-zero, we need to write it out.
        _write_term(io, OPPLUS, nlmodel)
        _write_term(io, expr.constant, nlmodel)
    end
    last_nary = false
    for term in expr.nonlinear_terms
        if last_nary
            println(io, term::Int)
            last_nary = false
        else
            _write_term(io, term, nlmodel)
            last_nary = _is_nary(term)
        end
    end
    return
end

function _write_linear_block(io::IO, expr::_NLExpr, nlmodel::Model)
    elements = [(c, nlmodel.x[v].order) for (v, c) in expr.linear_terms]
    for (c, x) in sort!(elements; by = i -> i[2])
        println(io, x, " ", _str(c))
    end
    return
end

function Base.write(io::IO, nlmodel::Model)
    # ==========================================================================
    # Header
    # Line 1: Always the same
    # Notes:
    #  * I think these are magic bytes used by AMPL internally for stuff.
    println(io, "g3 1 1 0")

    # Line 2: vars, constraints, objectives, ranges, eqns, logical constraints
    # Notes:
    #  * We assume there is always one objective, even if it is just `min 0`.
    n_con, n_ranges, n_eqns = 0, 0, 0
    for cons in (nlmodel.g, nlmodel.h), c in cons
        n_con += 1
        if c.opcode == 0
            n_ranges += 1
        elseif c.opcode == 4
            n_eqns += 1
        end
    end
    println(io, " $(length(nlmodel.x)) $(n_con) 1 $(n_ranges) $(n_eqns) 0")

    # Line 3: nonlinear constraints, objectives
    # Notes:
    #  * We assume there is always one objective, even if it is just `min 0`.
    n_nlcon = length(nlmodel.g)
    println(io, " ", n_nlcon, " ", 1)

    # Line 4: network constraints: nonlinear, linear
    # Notes:
    #  * We don't support network constraints. I don't know how they are
    #    represented.
    println(io, " 0 0")

    # Line 5: nonlinear vars in constraints, objectives, both
    # Notes:
    #  * This order is confusingly different to the standard "b, c, o" order.
    nlvb = length(nlmodel.types[1]) + length(nlmodel.types[2])
    nlvc = nlvb + length(nlmodel.types[3]) + length(nlmodel.types[4])
    nlvo = nlvb + length(nlmodel.types[5]) + length(nlmodel.types[6])
    println(io, " ", nlvc, " ", nlvo, " ", nlvb)

    # Line 6: linear network variables; functions; arith, flags
    # Notes:
    #  * I don't know what this line means. It is what it is. Apparently `flags`
    #    is set to 1 to get suffixes in .sol file.
    println(io, " 0 0 0 1")

    # Line 7: discrete variables: binary, integer, nonlinear (b,c,o)
    # Notes:
    #  * The order is
    #    - binary variables in linear only
    #    - integer variables in linear only
    #    - binary or integer variables in nonlinear objective and constraint
    #    - binary or integer variables in nonlinear constraint
    #    - binary or integer variables in nonlinear objective
    nbv = length(nlmodel.types[8])
    niv = length(nlmodel.types[9])
    nl_both = length(nlmodel.types[2])
    nl_cons = length(nlmodel.types[4])
    nl_obj = length(nlmodel.types[6])
    println(io, " ", nbv, " ", niv, " ", nl_both, " ", nl_cons, " ", nl_obj)

    # Line 8: nonzeros in Jacobian, gradients
    # Notes:
    #  * Make sure to include a 0 element for every variable that appears in an
    #    objective or constraint, even if the linear coefficient is 0.
    nnz_jacobian = 0
    for g in nlmodel.g
        nnz_jacobian += length(g.expr.linear_terms)
    end
    for h in nlmodel.h
        nnz_jacobian += length(h.expr.linear_terms)
    end
    nnz_gradient = length(nlmodel.f.linear_terms)
    println(io, " ", nnz_jacobian, " ", nnz_gradient)

    # Line 9: max name lengths: constraints, variables
    # Notes:
    #  * We don't add names, so this is just 0, 0.
    println(io, " 0 0")

    # Line 10: common exprs: b,c,o,c1,o1
    # Notes:
    #  * We don't add common subexpressions (i.e., V blocks).
    #  * I assume the notation means
    #     - b = in nonlinear objective and constraint
    #     - c = in nonlinear constraint
    #     - o = in nonlinear objective
    #     - c1 = in linear constraint
    #     - o1 = in linear objective
    println(io, " 0 0 0 0 0")
    # ==========================================================================
    # Constraints
    # Notes:
    #  * Nonlinear constraints first, then linear.
    #  * For linear constraints, write out the constant term here.
    for (i, g) in enumerate(nlmodel.g)
        println(io, "C", i - 1)
        _write_nlexpr(io, g.expr, nlmodel)
    end
    for (i, h) in enumerate(nlmodel.h)
        println(io, "C", i - 1 + n_nlcon)
        _write_nlexpr(io, h.expr, nlmodel)
    end
    # ==========================================================================
    # Objective
    # Notes:
    #  * NL files support multiple objectives, but we're just going to write 1,
    #    so it's always `O0`.
    #  * For linear objectives, write out the constant term here.
    println(io, "O0 ", nlmodel.sense == MOI.MAX_SENSE ? "1" : "0")
    _write_nlexpr(io, nlmodel.f, nlmodel)
    # ==========================================================================
    # VariablePrimalStart
    # Notes:
    #  * Make sure to write out the variables in order.
    println(io, "x", length(nlmodel.x))
    for (i, x) in enumerate(nlmodel.order)
        start = MOI.get(nlmodel, MOI.VariablePrimalStart(), x)
        println(io, i - 1, " ", start === nothing ? 0 : _str(start))
    end
    # ==========================================================================
    # Constraint bounds
    # Notes:
    #  * Nonlinear constraints go first, then linear.
    #  * The constant term for linear constraints gets written out in the
    #    "C" block.
    if n_con > 0
        println(io, "r")
        # Nonlinear constraints
        for g in nlmodel.g
            print(io, g.opcode)
            if g.opcode == 0
                println(io, " ", _str(g.lower), " ", _str(g.upper))
            elseif g.opcode == 1
                println(io, " ", _str(g.upper))
            elseif g.opcode == 2
                println(io, " ", _str(g.lower))
            else
                @assert g.opcode == 4
                println(io, " ", _str(g.lower))
            end
        end
        # Linear constraints
        for h in nlmodel.h
            print(io, h.opcode)
            if h.opcode == 0
                println(io, " ", _str(h.lower), " ", _str(h.upper))
            elseif h.opcode == 1
                println(io, " ", _str(h.upper))
            elseif h.opcode == 2
                println(io, " ", _str(h.lower))
            else
                @assert h.opcode == 4
                println(io, " ", _str(h.lower))
            end
        end
    end
    # ==========================================================================
    # Variable bounds
    # Notes:
    #  * Not much to note, other than to make sure you iterate the variables in
    #    the correct order.
    println(io, "b")
    for x in nlmodel.order
        v = nlmodel.x[x]
        if v.lower == v.upper
            println(io, "4 ", _str(v.lower))
        elseif -Inf < v.lower && v.upper < Inf
            println(io, "0 ", _str(v.lower), " ", _str(v.upper))
        elseif -Inf == v.lower && v.upper < Inf
            println(io, "1 ", _str(v.upper))
        elseif -Inf < v.lower && v.upper == Inf
            println(io, "2 ", _str(v.lower))
        else
            println(io, "3")
        end
    end
    # ==========================================================================
    # Jacobian block
    # Notes:
    #  * If a variable appears in a constraint, it needs to have a corresponding
    #    entry in the Jacobian block, even if the linear coefficient is zero.
    #    AMPL uses this to determine the Jacobian sparsity.
    #  * As before, nonlinear constraints go first, then linear.
    #  * Don't write out the `k` entry for the last variable, because it can be
    #    inferred from the total number of elements in the Jacobian as given in
    #    the header.
    if n_con > 0
        println(io, "k", length(nlmodel.x) - 1)
        total = 0
        for i in 1:length(nlmodel.order)-1
            total += nlmodel.x[nlmodel.order[i]].jacobian_count
            println(io, total)
        end
        for (i, g) in enumerate(nlmodel.g)
            println(io, "J", i - 1, " ", length(g.expr.linear_terms))
            _write_linear_block(io, g.expr, nlmodel)
        end
        for (i, h) in enumerate(nlmodel.h)
            println(io, "J", i - 1 + n_nlcon, " ", length(h.expr.linear_terms))
            _write_linear_block(io, h.expr, nlmodel)
        end
    end
    # ==========================================================================
    # Gradient block
    # Notes:
    #  * You only need to write this out if there are linear terms in the
    #    objective.
    if nnz_gradient > 0
        println(io, "G0 ", nnz_gradient)
        _write_linear_block(io, nlmodel.f, nlmodel)
    end
    return nlmodel
end

end


text:
"""
        BoundaryFunction{B, X1, X2}(func, parameters=nothing)

A wrapper for the user-defined boundary condition function `func`, on the
boundary specified by symbol `B` and at location `(X1, X2)`, and with `parameters`.

Example
=======

julia> using Oceananigans, Oceananigans.BoundaryConditions, Oceananigans.Fields

julia> top_tracer_flux = BoundaryFunction{:z, Cell, Cell}((x, y, t) -> cos(2π*x) * cos(t))
(::BoundaryFunction{:z,Cell,Cell,var"#7#8",Nothing}) (generic function with 1 method)

julia> top_tracer_bc = BoundaryCondition(Flux, top_tracer_flux);

julia> flux_func(x, y, t, p) = cos(p.k * x) * cos(p.ω * t); # function with parameters

julia> parameterized_u_velocity_flux = BoundaryFunction{:z, Face, Cell}(flux_func, (k=4π, ω=3.0))
(::BoundaryFunction{:z,Face,Cell,typeof(flux_func),NamedTuple{(:k, :ω),Tuple{Float64,Float64}}}) (generic function with 1 method)

julia> top_u_bc = BoundaryCondition(Flux, parameterized_u_velocity_flux);
"""
struct BoundaryFunction{B, X1, X2, F, P} <: Function
          func :: F
    parameters :: P

    function BoundaryFunction{B, X1, X2}(func; parameters=nothing) where {B, X1, X2}

        B ∈ (:x, :y, :z) || throw(ArgumentError("The boundary B at which the BoundaryFunction is
                                                to be applied must be either :x, :y, or :z."))

        return new{B, X1, X2, typeof(func), typeof(parameters)}(func, parameters)
    end
end

@inline call_boundary_function(func, ξ, η, t, ::Nothing) = func(ξ, η, t)
@inline call_boundary_function(func, ξ, η, t, parameters) = func(ξ, η, t, parameters)

@inline (bc::BoundaryFunction{:x, Y, Z})(j, k, grid, clock, state) where {Y, Z} =
    call_boundary_function(bc.func, ynode(Y, j, grid), znode(Z, k, grid), clock.time, bc.parameters)

@inline (bc::BoundaryFunction{:y, X, Z})(i, k, grid, clock, state) where {X, Z} =
    call_boundary_function(bc.func, xnode(X, i, grid), znode(Z, k, grid), clock.time, bc.parameters)

@inline (bc::BoundaryFunction{:z, X, Y})(i, j, grid, clock, state) where {X, Y} =
    call_boundary_function(bc.func, xnode(X, i, grid), ynode(Y, j, grid), clock.time, bc.parameters)

ext: jl


text:
using Test
using DataFrames
using ..ChemistryFeaturization: elements
using ..ChemistryFeaturization.FeatureDescriptor: OrbitalFeatureDescriptor
using SparseArrays

@testset "OrbitalFeatureDescriptor" begin
    ofd = OrbitalFeatureDescriptor()

    @testset "get_value" begin
        # individual elements
        @test get_value(ofd, "He") == "1s2"
        @test get_value(ofd, "C") == "2s2.2p2"
        @test get_value(ofd, "La") == "5d1.6s2"
        @test get_value(ofd, "Po") == "4f14.5d10.6s2.6p4"

        # and on an AtomGraph
        He_mol = AtomGraph(Float32.([1 0; 0 1]), ["He", "He"])
        @test get_value(ofd, He_mol) == ["1s2", "1s2"]
    end

    @testset "Encode" begin
        silica = AtomGraph(Float32.([0 1 0; 1 0 1; 0 1 0]), ["O", "Si", "O"])
        encoded_silica = encode(ofd, silica)
        I, J, V = findnz(encoded_silica)

        @test I == [2, 3, 4, 5, 2, 3]
        @test J == [1, 1, 2, 2, 3, 3]
        @test V == [2, 4, 2, 2, 2, 4]

        # non existent element - doesn't really make sense
        nonexistent = AtomGraph(Float32.([1 1; 1 1]), ["Non", "Existent"])
        @test_throws AssertionError ofd(nonexistent)
    end

    @testset "Decode" begin
        alumina = AtomGraph(
            Float32.([0 1 0 0 0; 1 0 1 0 0; 0 1 0 1 0; 0 0 1 0 1; 0 0 0 1 0]),
            ["O", "Al", "O", "Al", "O"],
        )
        @test all(decode(ofd, encode(ofd, alumina)) .== elements(alumina))
    end
end


text:
# RSATools.jl
#=
General notation
p   : large prime number
q   : another large prime number
m   : p x q
k   : (p-1) x (q-1)
e   : a large number relativley prime to k
W   : converted message, relativley prime to m
      and less than m

Base equation
eq  : ex - ky = 1

Public shared
e, m
=#

__precompile__()
module RSATools
    using Primes

    function EncryptionParams()

        function gen_primes(lo::Int, hi::Int)
            return primes(lo,hi)
        end

        function prime_pair(generated_primes)
            prime_one = rand(generated_primes)
            prime_two = rand(generated_primes)
            while prime_two == prime_one
                prime_two = rand(generated_primes)
            end
            return (prime_one,prime_two)
        end

        function gen_e(k)
            e = rand(100000:200000)
            while gcd(e,k) != 1
                e = rand(100000:200000)
            end
            return e
        end

        secret_primes = prime_pair(gen_primes(10000,20000))
        p = secret_primes[1]
        q = secret_primes[2]
        m = p*q
        k = (p-1)*(q-1)
        e = gen_e(k)

        return p, q, m, k, e
    end

    p, q, m, k, e = EncryptionParams()

    function encrypt(Msg::String)
        W = Vector{Int}()
        C = Vector{Int}()
        for chr in Msg
            w = Int(chr)
            append!(W, w)
            append!(C, mod(w^e,m))
        end

        return C
    end

    function decrypt()
        Nothing
    end

    export encrypt, decrypt, m, e
end


text:
# Mutation operators
# ==================

# Genetic mutations
# =================

# Binary mutations
# ---------------------

"""
    flip(recombinant)

Returns an in-place mutated binary `recombinant` with a bit flips at random positions.
"""
function flip(
    recombinant::T;
    rng::AbstractRNG = Random.default_rng(),
) where {T<:AbstractVector{Bool}}
    s = length(recombinant)
    pos = rand(rng, 1:s)
    recombinant[pos] = !recombinant[pos]
    return recombinant
end
flip(recombinant::T; rng::AbstractRNG = Random.default_rng(),) where {T<:AbstractIndividual} =
    T(flip(variables(recombinant), rng=rng))

"""
    bitinversion(recombinant)

Returns an in-place mutated binary `recombinant` with its bits inverted.
"""
bitinversion(recombinant::T) where {T<:AbstractVector{Bool}} =
    map!(!, recombinant, recombinant)
bitinversion(recombinant::T) where {T<:AbstractIndividual} =
    T(bitinversion(variables(recombinant)))

# Real-valued mutations
# ---------------------

"""
    uniform(r = 1.0)

Returns an in-place real valued mutation function that performs the uniform distributed mutation [^1].

The mutation operator randomly chooses a number ``z`` in from the uniform distribution on the interval ``[-r,r]``, the mutation range.
The mutated individual is given by

- ``x_i^\\prime = x_i + z_i``

"""
function uniform(r::Real = 1.0)
    function mutation(
        recombinant::T;
        rng::AbstractRNG = Random.default_rng(),
    ) where {T<:AbstractVector}
        d = length(recombinant)
        recombinant .+= 2r .* rand(rng, d) .- r
        return recombinant
    end
    function mutation(
        recombinant::T;
        rng::AbstractRNG = Random.default_rng(),
    ) where {T<:AbstractIndividual}
        T(mutation(variables(recombinant), rng=rng))
    end
    return mutation
end

"""
    gaussian(σ = 1.0)

Returns an in-place real valued mutation function that performs the normal distributed mutation [^1].

The mutation operator randomly chooses a number ``z`` in from the normal distribution ``\\mathcal{N}(0,\\sigma)`` with standard deviation ``\\sigma``.
The mutated individual is given by

- ``x_i^\\prime = x_i + z_i``

"""
function gaussian(σ::Real = 1.0)
    function mutation(
        recombinant::T;
        rng::AbstractRNG = Random.default_rng(),
    ) where {T<:AbstractVector}
        d = length(recombinant)
        recombinant .+= σ .* randn(rng, d)
        return recombinant
    end
    function mutation(
        recombinant::T;
        rng::AbstractRNG = Random.default_rng(),
    ) where {T<:AbstractIndividual}
        T(mutation(variables(recombinant), rng=rng))
    end
    return mutation
end

"""
    BGA(valrange, m = 20)

Returns an in-place real valued mutation function that performs the BGA mutation scheme with the mutation range `valrange` and the mutation probability `1/m` [^1].
"""
function BGA(valrange::Vector, m::Int = 20)
    prob = 1.0 / m
    function mutation(
        recombinant::T;
        rng::AbstractRNG = Random.default_rng(),
    ) where {T<:AbstractVector}
        d = length(recombinant)
        @assert length(valrange) == d "Range matrix must have $(d) columns"
        δ = zeros(m)
        for i = 1:length(recombinant)
            for j = 1:m
                δ[j] = (rand(rng) < prob) ? δ[j] = 2.0^(-j) : 0.0
            end
            if rand(rng, Bool)
                recombinant[i] += sum(δ) * valrange[i]
            else
                recombinant[i] -= sum(δ) * valrange[i]
            end
        end
        return recombinant
    end
    function mutation(
        recombinant::T;
        rng::AbstractRNG = Random.default_rng(),
    ) where {T<:AbstractIndividual}
        T(mutation(variables(recombinant), rng=rng))
    end
    return mutation
end

"""
    PM(lower, upper, p = 2)

Returns an in-place real valued mutation function that performs the Power Mutation (PM)
scheme within `lower` and `upper` bound, and an index of
mutation `p`[^3].

*Note:* The implementation is a degenerate case of Mixed Integer Power Mutation ([`MIPM`](@ref))
"""
function PM(lower::Vector, upper::Vector, p::Float64 = 5.0) # index of distribution p
    return mipmmutation(lower, upper, p)
end

"""
    MIPM(lower, upper, p_real = 10, p_int = 4)

Returns an in-place real valued mutation function that performs the Mixed Integer Power Mutation (MI-PM) scheme within `lower` and `upper` bound, and an index of mutation `p_real` for real value and `p_int` for integer values[^4].
"""
function MIPM(
    lowerBounds::Vector,
    upperBounds::Vector,
    p_real::Float64 = 10.0,
    p_int::Float64 = 4.0,
) # index of distribution p
    return mipmmutation(lowerBounds, upperBounds, p_real, p_int)
end

function mipmmutation(
    lowerBounds::Vector,
    upperBounds::Vector,
    p_real::Float64,
    p_int::Union{Nothing,Float64} = nothing,
)
    function pm_mutation(rng, x, l, u, s, d)
        x̄ = d < rand(rng) ? x - s * (x - l) : x + s * (u - x)
        if isa(x, Integer)
            if isinteger(x̄)
                Int(x̄)
            else
                floor(Int, x̄) + (rand(rng) > 0.5)
            end
        else
            x̄
        end
    end
    function mutation(
        recombinant::T;
        rng::AbstractRNG = Random.default_rng(),
    ) where {T<:Vector}
        d = length(recombinant)
        @assert length(lowerBounds) == d "Bounds vector must have $(d) columns"
        @assert length(upperBounds) == d "Bounds vector must have $(d) columns"
        @assert !(p_int === nothing && any(isa(x, Integer) for x in recombinant)) "Need to set p_int for integer variables"
        u = rand(rng)
        P = (isa(x, Integer) ? p_int : p_real for x in recombinant)
        S = u .^ (1 ./ P) # random var following power distribution
        D = (recombinant - lowerBounds) ./ (upperBounds - lowerBounds)
        broadcast!(
            (x, l, u, s, d) -> pm_mutation(rng, x, l, u, s, d),
            recombinant,
            recombinant,
            lowerBounds,
            upperBounds,
            S,
            D,
        )
        return recombinant
    end
    function mutation(
        recombinant::T;
        rng::AbstractRNG = Random.default_rng(),
    ) where {T<:AbstractIndividual}
        T(mutation(variables(recombinant), rng=rng))
    end
    return mutation
end

"""
    PLM(lower, upper, η = 2)

Returns an in-place real valued mutation function that performs the Polynomial Mutation (PLM) scheme
within `lower` and `upper` bounds, and a mutation distribution index `η`[^9].
"""
function PLM(Δ::Union{Real,Vector} = 1.0; η = 2, pm::Real = NaN) # index of distribution p
    function mutation(
        recombinant::T;
        rng::AbstractRNG = Random.default_rng(),
    ) where {T<:AbstractVector}
        d = length(recombinant)
        pm = isnan(pm) ? 1 / d : pm
        mask = rand(rng, d) .< pm
        u = rand(rng, d)
        mask_p = u .<= 0.5
        mask_n = u .> 0.5
        δpw = 1 / (η + 1)
        δ = similar(recombinant)
        δ[mask_p] .= (2 * u[mask_p]) .^ δpw .- 1
        δ[mask_n] .= 1 .- (2(1 .- u[mask_n])) .^ δpw
        recombinant[mask] .+= (Δ.*δ)[mask]
        return recombinant
    end
    function mutation(
        recombinant::T;
        rng::AbstractRNG = Random.default_rng(),
    ) where {T<:AbstractIndividual}
        T(mutation(variables(recombinant), rng=rng))
    end
    return mutation
end
PLM(lower::Vector, upper::Vector; η::Real = 2, pm::Real = NaN) =
    PLM(upper - lower; η = η, pm = pm)

# Combinatorial mutations (applicable to binary vectors)
# ------------------------------------------------------

"""
    inversion(recombinant)

Returns an in-place mutated individual with a random arbitrary length segment of the genome in the reverse order.
"""
function inversion(
    recombinant::T;
    rng::AbstractRNG = Random.default_rng(),
) where {T<:AbstractVector}
    l = length(recombinant)
    from, to = randseg(rng, l)
    l = round(Int, (to - from) / 2)
    if from + 1 == to
        swap!(recombinant, from, to)
    else
        for i = 0:(l-1)
            swap!(recombinant, from + i, to - i)
        end
    end
    return recombinant
end
inversion(recombinant::T; rng::AbstractRNG = Random.default_rng(),) where {T<:AbstractIndividual} =
    T(inversion(variables(recombinant), rng = rng))

"""
    insertion(recombinant)

Returns an in-place mutated individual with an arbitrary element of the genome moved in a random position.
"""
function insertion(
    recombinant::T;
    rng::AbstractRNG = Random.default_rng(),
) where {T<:AbstractVector}
    l = length(recombinant)
    from, to = randseg(rng, l)
    val = recombinant[from]
    deleteat!(recombinant, from)
    return insert!(recombinant, to, val)
end
insertion(recombinant::T; rng::AbstractRNG = Random.default_rng(),) where {T<:AbstractIndividual} =
    T(insertion(variables(recombinant), rng = rng))

"""
    swap2(recombinant)

Returns an in-place mutated individual with a two random elements of the genome are swapped.
"""
function swap2(
    recombinant::T;
    rng::AbstractRNG = Random.default_rng(),
) where {T<:AbstractVector}
    l = length(recombinant)
    from, to = randseg(rng, l)
    swap!(recombinant, from, to)
    return recombinant
end
swap2(recombinant::T; rng::AbstractRNG = Random.default_rng(),) where {T<:AbstractIndividual} =
    T(swap2(variables(recombinant), rng = rng))

"""
    scramble(recombinant)

Returns an in-place mutated individual with elements, on a random arbitrary length segment of the genome, been scrambled.
"""
function scramble(
    recombinant::T;
    rng::AbstractRNG = Random.default_rng(),
) where {T<:AbstractVector}
    l = length(recombinant)
    from, to = randseg(rng, l)
    diff = to - from + 1
    if diff > 1
        patch = recombinant[from:to]
        idx = randperm(rng, diff)
        for i = 1:diff
            recombinant[from+i-1] = patch[idx[i]]
        end
    end
    return recombinant
end
scramble(recombinant::T; rng::AbstractRNG = Random.default_rng(),) where {T<:AbstractIndividual} =
    T(scramble(variables(recombinant), rng = rng))

"""
    shifting(recombinant)

Returns an in-place mutated individual with a random arbitrary length segment of the genome been shifted to an arbitrary position.
"""
function shifting(
    recombinant::T;
    rng::AbstractRNG = Random.default_rng(),
) where {T<:AbstractVector}
    l = length(recombinant)
    from, to, where = sort(rand(rng, 1:l, 3))
    patch = recombinant[from:to]
    diff = where - to
    if diff > 0
        # move values after tail of patch to the patch head position
        for i = 1:diff
            recombinant[from+i-1] = recombinant[to+i]
        end
        # place patch values in order
        start = from + diff
        for i = 1:length(patch)
            recombinant[start+i-1] = patch[i]
        end
    end
    return recombinant
end
shifting(recombinant::T; rng::AbstractRNG = Random.default_rng(),) where {T<:AbstractIndividual} =
    T(shifting(variables(recombinant), rng = rng))

"""
    replace(pool,[minchange=1])(recombinant)

Replacement mutation operator changes an arbitrary number, no smaller then `minchange`,
of elements in the individual by replacing them with elements from the predefined `pool` that are not in the individual.
"""
function replace(pool::Vector{P}; minchange = 1) where {P}
    function rplc(
        recombinant::T;
        rng::AbstractRNG = Random.default_rng(),
    ) where {T<:AbstractVector}
        l = length(recombinant)
        p = length(pool)
        # how many values to change
        nchg = max(minchange, rand(rng, 1:min(l, p - l)))
        # select new elements
        idxs = randperm(rng, p)
        new_vals = P[]
        for i in idxs
            if pool[i] ∉ recombinant
                push!(new_vals, pool[i])
            end
            length(new_vals) == nchg && break
        end
        # update arbitrary positions with new values
        new_idxs = randperm(rng, l)[1:nchg]
        recombinant[new_idxs] .= new_vals
        return recombinant
    end
    function rplc(
        recombinant::T;
        rng::AbstractRNG = Random.default_rng(),
    ) where {T<:AbstractIndividual}
        T(rplc(variables(recombinant), rng = rng))
    end
    return rplc
end

# Differential Evolution
# ======================

"""
    differentiation(recombinant, mutators; F = 1.0)

Returns an in-place differently mutated individual ``x^\\prime`` from `recombinant` ``x``  by `mutators` ``\\{\\xi_1, \\ldots, \\xi_n \\}`` as follows

- ``x^\\prime = x + \\sum_{i=1}^{n/2} F (\\xi_{2i-1} - \\xi_{2i})``

"""
function differentiation!(
    recombinant::AbstractVector,
    mutators::AbstractMatrix;
    F::Real = 1.0,
)
    r, c = size(mutators)
    @assert r == 2 "mutators must be 2"
    recombinant .+= F .* (mutators[1,:] .- mutators[2,:])
    return recombinant
end
differentiation!(recombinant::T, mutators::Vector{T}; F::Real = 1.0,) where {T<:AbstractIndividual} =
    T(differentiation!(variables(recombinant), variables(mutators), F=F))


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

JLLWrappers.@generate_main_file_header("libsndfile")
JLLWrappers.@generate_main_file("libsndfile", UUID("5bf562c0-5a39-5b4f-b979-f64ac885830c"))
end  # module libsndfile_jll


text:
#=
algebra:
- Julia version: 1.4.2
- Author: akhil
- Date: 2020-05-28
=#
"""
Only Algebra Module is defined.

Please include the algo file in this Module.
"""
module Algebra
    include("TestAlgebra.jl")
end


text:
@printf("\nTesting deep operations\n")

x = ([2.0, 3.0], [4.0, 5.0, 6.0], [1.0 2.0 3.0; 4.0 5.0 6.0])

lengths_x = (2, 3, 6)
deeplength_x = 11
deepvecnorm_x = 13.45362404707371

@test length.(x) == lengths_x
@test ProximalOperators.deeplength(x) == deeplength_x

y = ProximalOperators.deepsimilar(x)

@test ProximalOperators.deeplength(y) == deeplength_x
@test length.(y) == lengths_x

ProximalOperators.deepcopy!(y, x)

@test y == x
@test ProximalOperators.deepvecnorm(x) ≈ deepvecnorm_x
@test ProximalOperators.deepvecdot(x, y) ≈ deepvecnorm_x^2
@test ProximalOperators.deepmaxabs(x .- y) == 0

t1 = ProximalOperators.deepzeros((Float32, Float64), ((3, ), (4, )) )

t1 = (randn(20), randn(20))
t2 = (randn(20), randn(20))
t3 = ProximalOperators.deepsimilar(t1)
ProximalOperators.deepaxpy!(t3, t1, 0.5, t2)
t4 = (t1[1]+0.5*t2[1],t1[2]+0.5*t2[2])
@test max(norm(t4[1]-t3[1]),norm(t4[2]-t3[2])) <= 1e-12


text:
using Gadfly, DataArrays, RDatasets

custom_theme = Theme(panel_fill=color("black"),
                     default_color=color("pink"))

plot(dataset("datasets", "iris"), x="SepalLength", y="SepalWidth", Geom.point,
     custom_theme)


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

@recipe function f(ray::Ray)
  seriestype --> :path
  linecolor --> :black
  arrow --> true
  [ray(0), ray(1)]
end

@recipe function f(segment::Segment)
  seriestype --> :path
  linecolor --> :black
  vertices(segment)
end

@recipe function f(sphere::Sphere, nsamples=100)
  seriestype --> :path
  linecolor --> :black
  samples = sample(sphere, RegularSampling(nsamples))
  points  = collect(samples)
  [points; first(points)]
end

@recipe function f(ball::Ball, nsamples=100)
  seriestype --> :shape
  linecolor --> :black
  samples = sample(boundary(ball), RegularSampling(nsamples))
  points  = collect(samples)
  [points; first(points)]
end

@recipe function f(polygon::Polygon)
  seriestype --> :path
  seriescolor --> :auto
  fill --> true
  points = vertices(polygon)
  [points; first(points)]
end

@recipe function f(chain::Chain)
  seriestype --> :path
  linecolor --> :black
  points = vertices(chain)
  isclosed(chain) ? [points; first(points)] : points
end

@recipe function f(polyarea::PolyArea)
  if hasholes(polyarea)
    mesh  = discretize(polyarea, FIST())
    geoms = collect(mesh)
  else
    geoms = [first(chains(polyarea))]
  end

  seriestype --> :path
  @series begin
    geoms
  end
  for chain in chains(polyarea)
    @series begin
      primary --> false
      chain
    end
  end
end

@recipe function f(multi::Multi)
  collect(multi)
end

@recipe function f(geometries::AbstractVector{G}) where {G<:Geometry}
  @series begin
    first(geometries)
  end
  for geometry in Iterators.drop(geometries, 1)
    @series begin
      primary --> false
      geometry
    end
  end
end


text:
export iauHd2pa
"""
Parallactic angle for a given hour angle and declination.

This function is part of the International Astronomical Union's
SOFA (Standards of Fundamental Astronomy) software collection.

Status:  support function.

Given:
   ha     double     hour angle
   dec    double     declination
   phi    double     site latitude

Returned (function value):
         double     parallactic angle

Notes:

   1.  All the arguments are angles in radians.

   2.  The parallactic angle at a point in the sky is the position
   angle of the vertical, i.e. the angle between the directions to
   the north celestial pole and to the zenith respectively.

   3.  The result is returned in the range -pi to +pi.

   4.  At the pole itself a zero result is returned.

   5.  The latitude phi is pi/2 minus the angle between the Earth's
   rotation axis and the adopted zenith.  In many applications it
   will be sufficient to use the published geodetic latitude of the
   site.  In very precise (sub-arcsecond) applications, phi can be
   corrected for polar motion.

   6.  Should the user wish to work with respect to the astronomical
   zenith rather than the geodetic zenith, phi will need to be
   adjusted for deflection of the vertical (often tens of
   arcseconds), and the zero point of the hour angle ha will also
   be affected.

Reference:
   Smart, W.M., "Spherical Astronomy", Cambridge University Press,
   6th edition (Green, 1977), p49.

Last revision:   2017 September 12

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# double iauHd2pa (double ha, double dec, double phi)

function iauHd2pa(ha::Real, dec::Real, phi::Real)
   return ccall((:iauHd2pa, libsofa_c), Cdouble,
         (Cdouble, Cdouble, Cdouble),
          convert(Float64, ha), convert(Float64, dec),
          convert(Float64, phi))
end


text:
using testing_pkg
using Test

@testset "testing_pkg.jl" begin
    @test add2(2) == 4
    @test add2(5) == 7

end

@testset "my_f Tests" begin
    @test my_f(4) == 14
    @test my_f(1) == 5
end

@testset "d_my_f Tests" begin
    @test d_my_f(1) == 3
    @test d_my_f(99) == 3
end


text:
using SparseArrays
#------------------------------------------------------------------------------
struct SparseData
    i::Vector{Int}
    j::Vector{Int}
    e::Vector{Int}
    k0::Vector{Float64}
    k::Vector{Float64}
end

#------------------------------------------------------------------------------
struct FEM
    NNode::Int
    Node::Array{Float64,2}
    NElem::Int
    Elem::Array{Int,2}
    ElemLen::Vector{Float64}
    ElemNDof::Vector{Int}
    CumSumElemNDof2::Vector{Int}
    ElemDof::Vector{Vector{Int}}
    Supp::Support
    Load::NodalLoad
    LocalIntMap::Array{Float64,2}
    YoungModulus::Float64
    Force::Vector{Float64}
    FreeDofs::Vector{Int}
    Sparse::SparseData
end

#------------------------------------------------------------------------------
function InitFEA(Node::Array{Float64,2}, Bars::Array{Int,2}, Supp::Support, Load::NodalLoad, E::Float64)::FEM

    NNode, NElem = size(Node, 1), size(Bars, 1)

    D = Node[Bars[:,2],:] - Node[Bars[:,1],:]
    L = vec(sqrt.(sum(D.^2, dims=2))); N = D./L
    LocalInternalMap = [-N N]

    ElemNDof = 4*ones(Int, NElem); CumSumElemNDof2 = cumsum(ElemNDof.^2)
    Ke = zeros(4, 4); Nval::Int = sum(ElemNDof.^2)
    i = zeros(Int, Nval); j = zeros(Int, Nval); e = zeros(Int, Nval); k0 = zeros(Nval)

    idx = 0; ElemDof = Vector{Int}[]; sizehint!(ElemDof, NElem)
    for el = 1:NElem
        LocalK!(Ke, LocalInternalMap[el,:], L[el], E)

        NDof = ElemNDof[el]; Idx = idx+1:idx+NDof^2
        EDof = permutedims([2Bars[el,:].-1 2Bars[el,:]])[:]
        push!(ElemDof, EDof)

        IdxI = repeat(EDof, outer=(1, NDof))
        IdxJ = permutedims(IdxI)
        i[Idx] = IdxI[:]; j[Idx] = IdxJ[:]; k0[Idx] = Ke[:]; e[Idx] .= el
        idx += NDof^2
    end
    k = zeros(Nval)

    F = zeros(2NNode)
    for iL = 1:length(Load.node)
        F[2Load.node[iL] - 1] += Load.val[iL, 1]
        F[2Load.node[iL]] += Load.val[iL, 2]
    end

    NSupp = length(Supp.node)
    FixedDofs, idx = zeros(Int, 2NSupp), 0
    for iS = 1:NSupp
        if Supp.type[iS,1] == 1; idx += 1; FixedDofs[idx] = 2Supp.node[iS]-1; end
        if Supp.type[iS,2] == 1; idx += 1; FixedDofs[idx] = 2Supp.node[iS];   end
    end
    FreeDofs = setdiff(1:2NNode, FixedDofs[1:idx])

    return FEM(NNode, Node, NElem, Bars, L, ElemNDof, CumSumElemNDof2, ElemDof,
    Supp, Load, LocalInternalMap, E, F, FreeDofs, SparseData(i, j, e, k0, k))
end

#------------------------------------------------------------------------------
function LocalK!(Ke::Array{Float64,2}, LocalIntMap::Vector{Float64}, L::Float64, E::Float64)

    NDof = length(LocalIntMap)
    for j = 1:NDof, i = 1:NDof
        Ke[i, j] = E / L * (LocalIntMap[i] * LocalIntMap[j])
    end
end

#------------------------------------------------------------------------------
function FEAnalysis!(U::Vector{Float64}, x::Vector{Float64}, FEA::FEM)

    NDof, Free = 2FEA.NNode, FEA.FreeDofs

    for idx in eachindex(FEA.Sparse.i)
        @fastmath FEA.Sparse.k[idx] = FEA.Sparse.k0[idx] * x[FEA.Sparse.e[idx]]
    end
    @fastmath K = sparse(FEA.Sparse.i, FEA.Sparse.j, FEA.Sparse.k, NDof, NDof)

    @fastmath U[Free] = K[Free, Free] \ FEA.Force[Free]
end

#------------------------------------------------------------------------------
function Stress!(σ::Vector{Float64}, U::Vector{Float64}, FEA::FEM)

    for e = 1:FEA.NElem
        Δ = 0.
        for j = 1:FEA.ElemNDof[e]
            Δ += FEA.LocalIntMap[e, j] * U[FEA.ElemDof[e]][j]
        end
        σ[e] = FEA.YoungModulus / FEA.ElemLen[e] * Δ
    end
end

max_stars_repo_name: ricardoaf/jutop

max_issues_repo_name: ricardoaf/jutop

max_forks_repo_name: ricardoaf/jutop


text:
```

ДАНО: Робот - рядом с горизонтальной перегородкой под ней, бесконечно продолжающейся в обе стороны, в которой имеется проход шириной в одну клетку.
РЕЗУЛЬТАТ: Робот - в клетке под проходом

```
using HorizonSideRobots
originRobot = Robot(animate = true, "8example.sit")
include("../AbstractType.jl")
include("../functions.jl")

function master8!(robot::Robot)
    steps = 0
    flag = true
    while flag
        for side in [Ost, West]
            moving_defsteps!(robot, side, steps)
            if !isborder(robot, Nord)
                move!(robot, Nord)
                flag = false
                break
            else
                steps += 1
            end
        end
    end
end

master8!(originRobot)

max_stars_repo_name: Igetsuga/AlgorithmicsHomework

max_issues_repo_name: Igetsuga/AlgorithmicsHomework

max_forks_repo_name: Igetsuga/AlgorithmicsHomework


text:
using BenchmarkTools

# Get the number of iterations from the command line.
n, = size(ARGS)
if n < 1
    println("Usage: fibonnaci.jl N")
    println("       ---> Please specify the number of iterations.")
    exit()
end

N = parse(Int, ARGS[1])

"""
    recursive_fibonacci(n)

Finds the nth Fibonacci number using recursion.
"""
function recursive_fibonacci(n)
    if n <= 2
        1.0
    else
        recursive_fibonacci(n-1) + recursive_fibonacci(n-2);
    end
end

"""
    iterative_fibonacci(n)

Finds the nth Fibonacci number using iteration.
"""
function iterative_fibonacci(n)
    x, y = (0, 1)
    for i = 1:n x, y = (y, x + y) end
    x
end

println("--------------------------")
println(@sprintf "Iterative - Fibonnaci %d" N)
println("--------------------------")

@btime iterative_fibonacci(N)

println("--------------------------")
println(@sprintf "Recursive - Fibonnaci %d" N)
println("--------------------------")

@btime recursive_fibonacci(N)


text:
@testset "576.out-of-boundary-paths.jl" begin
    @test find_paths(2, 2, 2, 1, 1) == 6

    @test find_paths(1, 3, 1, 1, 2) == 2
    @test find_paths(1, 3, 2, 1, 2) == 8
    @test find_paths(1, 3, 3, 1, 2) == 12
end


text:
InputFolder = './Images/NucleiCytoo/';
OutputFolder ='./Results/Images/NucleiCytoo/';
Fill = 1;
Lbl = 1;

@iA = '*C00*.tif';

@fxg_mGradWaterTiles [iA] > [L];
params.GaussianRadInt = 2;
params.ExtendedMinThr = 2;
/endf

@fxm_lTilesMerge [L, iA] > [L2];
params.GaussianRad = 2;
params.MinObjArea = 175;
params.MinSal = 0;
params.MaxValleyness = 1.1;
params.ConcavityThresh = 0.4;
/endf

/show iA > L2;
/keep L2 > tif;

max_stars_repo_name: SebastienTs/LOBSTER

max_issues_repo_name: SebastienTs/LOBSTER

max_forks_repo_name: SebastienTs/LOBSTER


text:
## These problems are quadratics and transformed quadratics.
# They correspond to Problems A-C in
# Hans De Sterck - Steepest descent preconditioning for nonlinear GMRES optimization

function quad(x::Vector, param)
    mat = param.mat
    xt = x-param.vec
    return 0.5*dot(xt, mat*xt)
end

function quad_gradient!(storage::Vector, x::Vector, param)
    mat = param.mat
    xt = x-param.vec
    storage .= mat*xt
end

function quad_fun_gradient!(storage::Vector, x::Vector, param)
    mat = param.mat
    xt = x-param.vec
    storage .= mat*xt
    return 0.5*dot(xt, mat*xt)
end

function quad_hessian!(storage::Matrix, x::Vector, param)
    storage .= param.mat
end

struct MatVecHolder{Tv <: AbstractVector,
                    Tm <: AbstractArray}
    mat::Tm
    vec::Tv
end

function _quadraticproblem(N::Int; mat::AbstractArray{T,2} = sparse(Diagonal(float(1:N))),
                           x0::AbstractVector{T} = ones(N),
                           initial_x::AbstractVector{T} = zeros(N),
                           name::AbstractString = "Quadratic Diagonal ($N)") where T <: Number
    # Note that _quadraticproblem is a special case of
    # _paraboloidproblem, where param.alpha = 0.0
    OptimizationProblem(name,
                        quad,
                        quad_gradient!,
                        quad_fun_gradient!,
                        quad_hessian!,
                        nothing, # Constraints
                        initial_x,
                        x0,
                        zero(T),
                        true,
                        true,
                        MatVecHolder(mat, x0))
end

examples["Quadratic Diagonal"] = _quadraticproblem(100)

#######################
# Paraboloid. Similar to Rosenbrock
#######################

struct ParaboloidStruct{T, Tm <: AbstractArray{T,2},
                           Tv <: AbstractArray{T}} <: Any where T<:Number
    mat::Tm
    vec::Tv
    xt::Tv
    alpha::T
end

function paraboloid(x::AbstractArray, param::ParaboloidStruct)
    mat = param.mat
    xt = param.xt
    @. xt = x - param.vec
    xt[2:end] .-= param.alpha*xt[1]^2

    return 0.5*dot(xt, mat*xt)
end

function paraboloid_gradient!(storage::AbstractArray, x::AbstractArray, param::ParaboloidStruct)
    mat = param.mat
    xt = param.xt

    @. xt = x - param.vec
    xt[2:end] .-= param.alpha*xt[1]^2

    storage .= mat*xt
    storage[1] -= 2.0*param.alpha*xt[1]*sum(storage[2:end])
end

function paraboloid_fun_gradient!(storage::AbstractArray, x::AbstractArray, param::ParaboloidStruct)
    mat = param.mat
    xt = param.xt

    @. xt = x - param.vec
    xt[2:end] .-= param.alpha*xt[1]^2

    storage .= mat*xt
    storage[1] -= 2.0*param.alpha*xt[1]*sum(storage[2:end])

    return 0.5*dot(xt, mat*xt)
end

function paraboloid_hessian!(storage,x,param)
    error("Hessian not implemented for Paraboloid")
end

function _paraboloidproblem(N::Int; mat::AbstractArray{T,2} = sparse(Diagonal(float(1:N))),
                            x0::AbstractVector{T} = ones(N),
                            initial_x::AbstractVector{T} = zeros(N),
                            alpha::T = 10.0,
                            name::AbstractString = "Paraboloid Diagonal ($N)") where T <: Number
    OptimizationProblem(name,
                        paraboloid,
                        paraboloid_gradient!,
                        paraboloid_fun_gradient!,
                        paraboloid_hessian!,
                        nothing, # Constraints
                        initial_x,
                        x0, # x0 means the solution ...
                        zero(T),
                        true,
                        false,
                        ParaboloidStruct(mat, x0, similar(x0), alpha))
end

examples["Paraboloid Diagonal"] = _paraboloidproblem(100)

function _randommatrix(N::Int, scaling::Bool=true)
    F = qr(randn(N,N))
    if scaling
        retval = F.Q'*sparse(Diagonal(float(1:N)))*F.Q
    else
        retval = F.Q'*F.Q
    end
    retval
end

guardseed(0) do
examples["Paraboloid Random Matrix"] = _paraboloidproblem(100;
                                                          name = "Paraboloid Random Matrix (100)",
                                                          mat = _randommatrix(100))
end


text:
using Printf
using NPZ
include("util.jl")
include("kmeans_filters.jl")
include("quantum_filters.jl")

log_file = open("run_filters.log", "a")

for name in ARGS
    target_label = parse(Int, split(name, "-")[3][end:end])
    reps = npzread("./Spectre/output/$(name)/label_$(target_label)_reps.npy")'
    n = size(reps)[2]
    eps = parse(Int, match(r"[0-9]+$", name).match)
    removed = round(Int, 1.5*eps)

    @printf("%s: Running PCA filter\n", name)
    reps_pca, U = pca(reps, 1)
    pca_poison_ind = k_lowest_ind(-abs.(mean(reps_pca[1, :]) .- reps_pca[1, :]), round(Int, 1.5*eps))
    poison_removed = sum(pca_poison_ind[end-eps+1:end])
    clean_removed = removed - poison_removed
    @show poison_removed, clean_removed
    @printf(log_file, "%s-pca: %d, %d\n", name, poison_removed, clean_removed)
    npzwrite("./Spectre/output/$(name)/mask-pca-target.npy", pca_poison_ind)

    @printf("%s: Running kmeans filter\n", name)
    kmeans_poison_ind = .! kmeans_filter2(reps, eps)
    poison_removed = sum(kmeans_poison_ind[end-eps+1:end])
    clean_removed = removed - poison_removed
    @show poison_removed, clean_removed
    @printf(log_file, "%s-kmeans: %d, %d\n", name, poison_removed, clean_removed)
    npzwrite("./Spectre/output/$(name)/mask-kmeans-target.npy", kmeans_poison_ind)

    @printf("%s: Running quantum filter\n", name)
    quantum_poison_ind = .! rcov_auto_quantum_filter(reps, eps)
    poison_removed = sum(quantum_poison_ind[end-eps+1:end])
    clean_removed = removed - poison_removed
    @show poison_removed, clean_removed
    @printf(log_file, "%s-quantum: %d, %d\n", name, poison_removed, clean_removed)
    npzwrite("./Spectre/output/$(name)/mask-rcov-target.npy", quantum_poison_ind)
end

max_stars_repo_name: TDteach/spectre-defense

max_issues_repo_name: TDteach/spectre-defense

max_forks_repo_name: TDteach/spectre-defense


text:
# Now, discard the corrupted lines. The remaining lines are incomplete.
#
# Incomplete lines don't have any incorrect characters - instead, they're
# missing some closing characters at the end of the line. To repair the
# navigation subsystem, you just need to figure out the sequence of closing
# characters that complete all open chunks in the line.
#
# You can only use closing characters (), ], }, or >), and you must add them in
# the correct order so that only legal pairs are formed and all chunks end up
# closed.
#
# In the example above, there are five incomplete lines:
#
# [({(<(())[]>[[{[]{<()<>> - Complete by adding }}]])})].
# [(()[<>])]({[<{<<[]>>( - Complete by adding )}>]}).
# (((({<>}<{<{<>}{[]{[]{} - Complete by adding }}>}>)))).
# {<[[]]>}<{[{[{[]{()[[[] - Complete by adding ]]}}]}]}>.
# <{([{{}}[<[[[<>{}]]]>[]] - Complete by adding ])}>.
#
# Did you know that autocomplete tools also have contests? It's true! The score
# is determined by considering the completion string character-by-character.
# Start with a total score of 0. Then, for each character, multiply the total
# score by 5 and then increase the total score by the point value given for the
# character in the following table:
#
# ): 1 point.
# ]: 2 points.
# }: 3 points.
# >: 4 points.
#
# So, the last completion string above - ])}> - would be scored as follows:
#
# - Start with a total score of 0.
# - Multiply the total score by 5 to get 0, then add the value of ] (2) to get
#   a new total score of 2.
# - Multiply the total score by 5 to get 10, then add the value of ) (1) to get
#   a new total score of 11.
# - Multiply the total score by 5 to get 55, then add the value of } (3) to get
#   a new total score of 58.
# - Multiply the total score by 5 to get 290, then add the value of > (4) to
#   get a new total score of 294.
#
# The five lines' completion strings have total scores as follows:
#
# }}]])})] - 288957 total points.
# )}>]}) - 5566 total points.
# }}>}>)))) - 1480781 total points.
# ]]}}]}]}> - 995444 total points.
# ])}> - 294 total points.
#
# Autocomplete tools are an odd bunch: the winner is found by sorting all of
# the scores and then taking the middle score. (There will always be an odd
# number of scores to consider.) In this example, the middle score is 288957
# because there are the same number of scores smaller and larger than it.
#
# Find the completion string for each incomplete line, score the completion
# strings, and sort the scores. What is the middle score?

const scores = Dict(')' => 1, ']' => 2, '}' => 3, '>' => 4)
const match = Dict('(' => ')', '[' => ']', '{' => '}', '<' => '>')
stack = Array{Char}(undef, 100)
cnts = map(eachline("input.txt")) do l
  p = 1
  for c in l
    if haskey(match, c)
      stack[p] = match[c]
      p += 1
    elseif p > 1 && stack[p-1] == c
      p -= 1
    else
      return 0
    end
  end
  return foldr((i,t) -> t * 5 + scores[stack[i]], 1:(p-1); init=0)
end

# find the middle score
sort!(cnts)
idx = findfirst(>(0), cnts)
idx += (length(cnts) - idx) >>> 1
println("Result: ", cnts[idx])


text:
#
# Bayes Net
#
Base.showable(::MIME"image/svg+xml", bn::BayesNet) = success(`lualatex -v`)

function plot(bn::BayesNet)
	if !isempty(names(bn))
		plot(bn.dag, AbstractString[string(s) for s in names(bn)]) # NOTE: sometimes the same var shows up twice
	else
		plot(DiGraph(1), ["Empty Graph"])
	end
end

function Base.show(f::IO, a::MIME"image/svg+xml", bn::BayesNet)
 	show(f, a, plot(bn))
end


text:
function AEMWriter()
    return AEMWriter(())
end

function AEMWriter(arg0::JString, arg1::JString, arg2::JString)
    return AEMWriter((JString, JString, JString), arg0, arg1, arg2)
end

function write(obj::AEMWriter, arg0::Appendable, arg1::AEMFile)
    return jcall(obj, "write", void, (Appendable, AEMFile), arg0, arg1)
end

function write(obj::AEMWriter, arg0::JString, arg1::AEMFile)
    return jcall(obj, "write", void, (JString, AEMFile), arg0, arg1)
end


text:
# Statistical_Analysis
module Statistical_Analysis

#Statistical_Analysis related code here
include(statistical_analysis/barchart.jl)
include(statistical_analysis/combined.jl)
include(statistical_analysis/correlation.jl)
include(statistical_analysis/fJBTest.jl)
include(statistical_analysis/fKurt.jl)
include(statistical_analysis/fSkew.jl)
include(statistical_analysis/fStd.jl)
include(statistical_analysis/fVAI.jl)
include(statistical_analysis/fgetXLData.jl)
include(statistical_analysis/fhistogram.jl)
include(statistical_analysis/jarquebera.jl)
include(statistical_analysis/moments.jl)
include(statistical_analysis/normalplot.jl)
include(statistical_analysis/regression.jl)
include(statistical_analysis/vami.jl)
end  #  module


text:
# test function kurt
x = [-1,0,2,-1]
@test kurt(x) == 18/4 - 3*(6/4)^2

x = [0]
@test kurt(x) == 0

x = [2,1,3,-2,-1,-3]
@test kurt(x) == mean(x.^4) - 3*mean(x.^2)^2

# test function grad_neg_sum_kurt()
@testset "Test grad_neg_sum_kurt" begin
    for idx = 1: 10
        Z = [randn(100000) for i = 1: 3]
        W = randn(3,3)
        @test rel_error(ForwardDiff.gradient(x -> neg_abs_sum_kurt(x'*Z), W), grad_neg_abs_sum_kurt(W, Z)) < 1e-8
    end
end

@testset "Test kurtosis-based icf $opt" for opt in ["orth", "sphe"]
    for idx = 1: 10
        # set up
        T = 100000
        x = [rand(T) + 1*ones(T), rand(T) - 1*ones(T)]
        A = randn(2,2)
        z = A*x;

        # apply icf
        Aest = icf(z; opt = opt)[1]
        @test err_pd(pinv(Aest)*A) < 5e-2
    end
end

# test function grad_ent_sum()
@testset "Test grad_ent_sum" begin
    for idx = 1: 10
        Z = [randn(100000) for i = 1: 3]
        W = randn(3,3)
        @test rel_error(ForwardDiff.gradient(x -> ent_sum(x'*Z), W), grad_ent_sum(W, Z)) < 1e-8
    end
end

@testset "Test entropy-based icf" begin
    for idx = 1: 10
        # set up
        T = 100000
        x = [rand(T) + 1*ones(T), rand(T) - 1*ones(T)]
        A = randn(2,2)
        z = A*x;

        # apply icf
        Aest = icf(z; obj = "ent")[1]
        @test err_pd(pinv(Aest)*A) < 5e-2
    end
end


text:
importall  Base

abstract FixedSizeArray{T, SZ, N}
abstract FixedSizeVector{T, N} <: FixedSizeArray{T, (N,), 1}

getindex{T <: FixedSizeVector}(c::T, i::Int) = c.(i)

# Ugly workaround for not having triangular dispatch:
eltype{T <: FixedSizeVector}(A::Type{T})                = A.types[1]
length{T <: FixedSizeVector}(A::Type{T})                = length(A.types)
ndims{T <: FixedSizeVector}(A::Type{T})                 = 1
size{T <: FixedSizeVector}(A::Type{T})                  = (length(A),)
size{T <: FixedSizeVector}(A::Type{T}, d::Integer)      = (length(A),) # should throw an error!?

abstract Func{Name, N}
stagedfunction map{T <: FixedSizeVector}(f::Func{1}, a::T)
	typename  = symbol("$(a.name.name)")
	arguments = [:(f(a.($i))) for i=1:length(T)]
    :($typename($(arguments...)))
end

stagedfunction map{T <: FixedSizeVector, FName}(f::Func{FName, 2}, a::T, b::T)
	typename  = symbol("$(a.name.name)")
	arguments = [:($FName(a.($i), b.($i))) for i=1:length(T)]
    :($typename($(arguments...)))
end

const binaryOps = (:.+, :.-,:.*, :./, :.\, :.^,:*,:/,
                   :.==, :.!=, :.<, :.<=, :.>, :.>=, :+, :-,
                   :min, :max,
                   :div, :fld, :rem, :mod, :mod1, :cmp,
                   :atan2, :besselj, :bessely, :hankelh1, :hankelh2,
                   :besseli, :besselk, :beta, :lbeta)
for op in binaryOps
    const unicsymb = gensym()

    eval(quote
        immutable $unicsymb <: Func{symbol(string($op)), 2} end
        $op{T <: FixedSizeVector}(x::T, y::T) = map($unicsymb(), x, y)
    end)
end

gen_fixedsizevector_type(name::DataType, T::Symbol, N::Int) = gen_fixedsizevector_type(symbol(string(name.name.name)), T, N)
function gen_fixedsizevector_type(name::Symbol, T::Symbol, N::Int)
	fields = [Expr(:(::), symbol("I_$i"), T) for i = 1:N]
	typename = symbol("FS$name")
	eval(quote
		immutable $(typename){$T} <: $(name){$T}
			$(fields...)
		end
	end)
	typename
end
abstract Dimension{T} <: Number

stagedfunction getindex{FT <: FixedSizeVector, T <: Dimension}(a::FT, key::Type{T})
    index = fieldindex(a, T)
    :(T(a[$index]))
end

macro accessors(typ, fields)
    result = Any[]
    for elem in fields.args
        push!(result, quote
            fieldindex{T <: $typ}(::Type{T}, ::Type{$(elem.args[1])}) = $(elem.args[2])
        end)
    end
    esc(Expr(:block, result...))
end

call{T <: FixedSizeVector, ET}(t::Type{T}, data::ET...) = t(ntuple(x->data[x], length(data)))
stagedfunction call{T <: FixedSizeVector, N, ET}(t::Type{T}, data::NTuple{N, ET})
	Tsuper, Nsuper = super(T).parameters
	@assert Nsuper == N "not the right dimension"
	typename = gen_fixedsizevector_type(T, Tsuper.name, N)
	:($typename(data...))
end

immutable RGB{T} <: FixedSizeVector{T, 3}
  r::T
  g::T
  b::T
end

immutable Red{T} <: Dimension{T}
    val::T
end
immutable Green{T} <: Dimension{T}
    val::T
end
immutable Blue{T} <: Dimension{T}
    val::T
end
@accessors RGB (Red => 1, Green => 2, Blue => 3)

@show a = RGB(0.1f0,0.1f0,0.3f0)
@show a = RGB(0.1f0,0.1f0,0.3f0)
@show a = RGB(0.1f0,0.1f0,0.3f0)
@show a[Green]
@show a[2]
@show a + a


text:
function AdamsFieldIntegrator(arg0::Field, arg1::JString, arg2::jint, arg3::jint, arg4::jdouble, arg5::jdouble, arg6::Vector{jdouble}, arg7::Vector{jdouble})
    return AdamsFieldIntegrator((Field, JString, jint, jint, jdouble, jdouble, Vector{jdouble}, Vector{jdouble}), arg0, arg1, arg2, arg3, arg4, arg5, arg6, arg7)
end

function AdamsFieldIntegrator(arg0::Field, arg1::JString, arg2::jint, arg3::jint, arg4::jdouble, arg5::jdouble, arg6::jdouble, arg7::jdouble)
    return AdamsFieldIntegrator((Field, JString, jint, jint, jdouble, jdouble, jdouble, jdouble), arg0, arg1, arg2, arg3, arg4, arg5, arg6, arg7)
end

function integrate(obj::AdamsFieldIntegrator, arg0::FieldExpandableODE, arg1::FieldODEState, arg2::RealFieldElement)
    return jcall(obj, "integrate", FieldODEStateAndDerivative, (FieldExpandableODE, FieldODEState, RealFieldElement), arg0, arg1, arg2)
end


text:
# using Revise, Plots
using Test, ForwardDiff, Parameters, Setfield, LinearAlgebra
using BifurcationKit
const BK = BifurcationKit

norminf = x -> norm(x, Inf)
####################################################################################################
function COm(u, p)
	@unpack q1,q2,q3,q4,q5,q6,k = p
	x, y, s = u
	z = 1-x-y-s
	out = similar(u)
	out[1] = 2 * q1 * z^2 - 2 * q5 * x^2 - q3 * x * y
	out[2] = q2 * z - q6 * y - q3 * x * y
	out[3] = q4 * z - k * q4 * s
	out
end
dCOm = (z, p) -> ForwardDiff.jacobian(x -> COm(x, p), z)
jet = BK.get3Jet(COm, dCOm)

par_com = (q1 = 2.5, q2 = 2.0, q3 = 10., q4 = 0.0675, q5 = 1., q6 = 0.1, k = 0.4)
z0 = [0.07,0.2,05]

opts_br = ContinuationPar(pMin = 0.6, pMax = 2.5, ds = 0.002, dsmax = 0.01, nInversion = 4, detectBifurcation = 3, maxBisectionSteps = 25, nev = 2, maxSteps = 20000)
	br, = @time continuation(jet[1], jet[2], z0, par_com, (@lens _.q2), opts_br;
	printSolution = (x, p) -> (x = x[1], y = x[2]),
	plot = false, verbosity = 0, normC = norminf,
	bothside = true)

# plot(br, plotfold=false, markersize=4, legend=:topleft)
####################################################################################################
hp, = newton(jet[1:2]..., br, 1; options = NewtonPar( opts_br.newtonOptions; maxIter = 10),startWithEigen=true, d2F = jet[3])

hpnf = computeNormalForm(jet..., br, 4)
@test hpnf.nf.b |> real ≈ 1.070259e+01 rtol = 1e-3

hpnf = computeNormalForm(jet..., br, 1)
@test hpnf.nf.b |> real ≈ 4.332247e+00 rtol = 1e-2
####################################################################################################
# different tests for the Fold point
snpt = computeNormalForm(jet..., br, 2)
@set! opts_br.newtonOptions.verbose = false
@set! opts_br.newtonOptions.maxIter = 10

sn = newton(jet[1:2]..., br, 2; options = opts_br.newtonOptions, bdlinsolver = MatrixBLS())
# printstyled(color=:red, "--> guess for SN, p = ", br.specialpoint[2].param, ", psn = ", sn[1].p)
	# plot(br);scatter!([sn.x.p], [sn.x.u[1]])
@test sn[3] && sn[5] == 6

sn = newton(jet[1:2]..., br, 2; options = opts_br.newtonOptions, bdlinsolver = MatrixBLS(), d2F = jet[3])
@test sn[3] && sn[5] == 8

sn = newton(jet[1:2]..., br, 2; options = opts_br.newtonOptions, bdlinsolver = MatrixBLS(), d2F = jet[3], startWithEigen = true)
@test sn[3] && sn[5] == 8

sn_br, = continuation(jet[1:2]..., br, 2, (@lens _.k), ContinuationPar(opts_br, pMax = 1., pMin = 0., detectBifurcation = 1, maxSteps = 50, saveSolEveryStep = 1, detectEvent = 2), bdlinsolver = MatrixBLS(), d2F = jet[3], startWithEigen = true, updateMinAugEveryStep = 1, verbosity = 0, plot=false)
@test sn_br.specialpoint[1].type == :bt
@test ~isnothing(sn_br.eig)

# we test the jacobian and problem update
par_sn = set(br.params, br.lens, sn_br.sol[end].x.p)
par_sn = set(par_sn, sn_br.lens, sn_br.sol[end].p)
_J = dCOm(sn_br.sol[end].x.u, par_sn)
_eigvals, eigvec, = eigen(_J)
ind = argmin(abs.(_eigvals))
@test _eigvals[ind] ≈ 0 atol=1e-10
ζ = eigvec[:, ind]
@test sn_br.functional.b ./ norm(sn_br.functional.b) ≈ -ζ

_eigvals, eigvec, = eigen(_J')
ind = argmin(abs.(_eigvals))
ζstar = eigvec[:, ind]
@test sn_br.functional.a ≈ ζstar
####################################################################################################
# different tests for the Hopf point
hppt = computeNormalForm(jet..., br, 1)

@set! opts_br.newtonOptions.verbose = false

hp = newtonHopf(jet[1:2]..., br, 1; options = opts_br.newtonOptions, startWithEigen = true)
# printstyled(color=:red, "--> guess for HP, p = ", br.specialpoint[1].param, ", php = ", hp[1].p)
# plot(br);scatter!([hp[1].p[1]], [hp[1].u[1]])
@test hp[3] && hp[5] == 12

hp = newtonHopf(jet[1:2]..., br, 1; options = opts_br.newtonOptions, startWithEigen = false, bdlinsolver = MatrixBLS())
@test hp[3] && hp[5] == 12

hp = newtonHopf(jet[1:2]..., br, 1; options = opts_br.newtonOptions, startWithEigen = true, bdlinsolver = MatrixBLS(), verbose = true)
@test hp[3] && hp[5] == 8

# we check that we truly have a bifurcation point.
pb = hp[end]
ω = hp[1].p[2]
par_hp = set(br.params, br.lens, hp[1].p[1])
_J = dCOm(hp[1].u, par_hp)
_eigvals, eigvec, = eigen(_J)
ind = argmin(abs.(_eigvals .- Complex(0, ω)))
@test real(_eigvals[ind]) ≈ 0 atol=1e-9
@test abs(imag(_eigvals[ind])) ≈ abs(hp[1].p[2]) atol=1e-9
ζ = eigvec[:, ind]
# reminder: pb.b should be a null vector of (J+iω)
@test pb.b ≈ ζ atol = 1e-3

hp, = newton(jet[1:2]..., br, 1;
	options = NewtonPar( opts_br.newtonOptions; maxIter = 10),
	startWithEigen = true,
	bdlinsolver = MatrixBLS(),
	d2F = jet[3])
# printstyled(color=:red, "--> guess for HP, p = ", br.specialpoint[1].param, ", php = ", hp.p)
# plot(br);scatter!([hp.p[1]], [hp.u[1]])

hp, = newton(jet[1:2]..., br, 1; options = NewtonPar( opts_br.newtonOptions; maxIter = 10),startWithEigen=true, d2F = jet[3])

hp_br, = continuation(jet[1:2]..., br, 1, (@lens _.k), ContinuationPar(opts_br, ds = -0.001, pMax = 1., pMin = 0., detectBifurcation = 1, maxSteps = 50, saveSolEveryStep = 1, detectEvent = 2), bdlinsolver = MatrixBLS(), d2F = jet[3], d3F = jet[4], startWithEigen = true, updateMinAugEveryStep = 1, verbosity=0, plot=false)
@test hp_br.specialpoint[1].type == :gh
@test ~isnothing(hp_br.eig)


text:
# Copyright (c) 2019 Arpit Bhatia, Chris Coey, and contributors                  #src
#                                                                                #src
# Permission is hereby granted, free of charge, to any person obtaining a copy   #src
# of this software and associated documentation files (the "Software"), to deal  #src
# in the Software without restriction, including without limitation the rights   #src
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell      #src
# copies of the Software, and to permit persons to whom the Software is          #src
# furnished to do so, subject to the following conditions:                       #src
#                                                                                #src
# The above copyright notice and this permission notice shall be included in all #src
# copies or substantial portions of the Software.                                #src
#                                                                                #src
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR     #src
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,       #src
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE    #src
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER         #src
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,  #src
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE  #src
# SOFTWARE.                                                                      #src

# # Experiment Design

# **Originally Contributed by**: Arpit Bhatia, Chris Coey

# This tutorial covers experiment design examples (D-optimal, A-optimal, and
# E-optimal) from section 7.5 of the book [Convex Optimization](https://doi.org/10.1017/CBO9780511804441)
# by Boyd and Vandenberghe.

# The tutorial uses the following packages
using JuMP
import SCS
import LinearAlgebra
import Random

# We set a seed so the random numbers are repeatable:
Random.seed!(1234)

# ## Relaxed Experiment Design Problem

# The basic experiment design problem is as follows.

# Given the menu of possible choices for experiments, $v_{1}, \ldots, v_{p}$,
# and the total number $m$ of experiments to be carried out, choose the numbers
# of each type of experiment, i.e., $m_{1}, \ldots, m_{p}$ to make the error
# covariance $E$ small (in  some sense).

# The variables $m_{1}, \ldots, m_{p}$ must, of course, be integers and sum to
# $m$ the given total number of experiments. This leads to the optimization
# problem:

# ```math
# \begin{aligned}
# \min\left(\mathrm{w.r.t.} \mathbf{S}_{+}^{n}\right) & E=\left(\sum_{j=1}^{p} m_{j} v_{j} v_{j}^{T}\right)^{-1} \\
# \text{subject to} & m_{i} \geq 0 \\
# & \sum\limits_{i=1}^p m_{i} = m \\
# & m_{i} \in \mathbb{Z},\quad i=1,\ldots,p
# \end{aligned}
# ```

# The basic experiment design problem can be a hard combinatorial problem when
# $m$, the total number of experiments, is comparable to $n$, since in this case
# the $m_{i}$ are all small integers.

# In the case when $m$ is large compared to $n$, however, a good approximate
# solution can be found by ignoring, or relaxing, the constraint that the
# $m_{i}$ are integers.

# Let $\lambda_{i}=m_{i} / m,$ which is the fraction of the total number of
# experiments for which $a_{j}=v_{i},$ or the relative frequency of experiment
# $i$. We can express the error covariance in terms of $\lambda_{i}$ as:

# ```math
# E=\frac{1}{m}\left(\sum_{i=1}^{p} \lambda_{i} v_{i} v_{i}^{T}\right)^{-1}
# ```

# The vector $\lambda \in \mathbf{R}^{p}$ satisfies
# $\lambda \succeq 0, \mathbf{1}^{T} \lambda=1,$ and also, each
# $\lambda_{i}$ is an integer multiple of $1 / m$. By ignoring this last
# constraint, we arrive at the problem:

# ```math
# \begin{aligned}
# \min\left(\mathrm{w.r.t.} \mathbf{S}_{+}^{n}\right) & E=(1 / m)\left(\sum_{i=1}^{p} \lambda_{i} v_{i} v_{i}^{T}\right)^{-1} \\
# \text{subject to:}\quad & \lambda \succeq 0 \\
# & \mathbf{1}^{T} \lambda=1
# \end{aligned}
# ```

# Several scalarizations have been proposed for the experiment design problem,
# which is a vector optimization problem over the positive semidefinite cone.

q = 4 # dimension of estimate space
p = 8 # number of experimental vectors
nmax = 3 # upper bound on lambda
n = 12

V = randn(q, p)

eye = Matrix{Float64}(LinearAlgebra.I, q, q);

# ## A-optimal design

# In A-optimal experiment design, we minimize tr $E$, the trace of the
# covariance matrix. This objective is simply the mean of the norm of the error
# squared:

# ```math
# \mathbf{E}\|e\|_{2}^{2}=\mathbf{E} \operatorname{tr}\left(e e^{T}\right)=\operatorname{tr} E
# ```

# The A-optimal experiment design problem in SDP form is

# ```math
# \begin{aligned}
# \min & \mathbf{1}^{T} u \\
# \text{subject to} & \left[\begin{aligned}{\sum_{i=1}^{p} \lambda_{i} v_{i} v_{i}^{T}} & {e_{k}} \\ {e_{k}^{T}} & {u_{k}}\end{aligned}\right] \succeq 0, \quad k=1, \ldots, n \\
# & \lambda \succeq 0 \\
# &  \mathbf{1}^{T} \lambda=1
# \end{aligned}
# ```

aOpt = Model(SCS.Optimizer)
set_silent(aOpt)
@variable(aOpt, np[1:p], lower_bound = 0, upper_bound = nmax)
@variable(aOpt, u[1:q], lower_bound = 0)
@constraint(aOpt, sum(np) <= n)
for i = 1:q
    matrix = [
        V * LinearAlgebra.diagm(0 => np ./ n) * V' eye[:, i];
        eye[i, :]' u[i]
    ]
    @SDconstraint(aOpt, matrix >= 0)
end
@objective(aOpt, Min, sum(u))
optimize!(aOpt)
objective_value(aOpt)

#-

value.(np)

# ## E-optimal design

# In $E$ -optimal design, we minimize the norm of the error covariance matrix,
# i.e. the maximum eigenvalue of $E$.

# Since the diameter (twice the longest semi-axis) of the confidence ellipsoid
# $\mathcal{E}$ is proportional to $\|E\|_{2}^{1 / 2}$, minimizing $\|E\|_{2}$
# can be interpreted geometrically as minimizing the diameter of the confidence
# ellipsoid.

# E-optimal design can also be interpreted as minimizing the maximum variance of
# $q^{T} e$, over all $q$ with $\|q\|_{2}=1$. The E-optimal experiment design
# problem in SDP form is:

# ```math
# \begin{aligned}
# \min & t \\
# \text{subject to} & \sum_{i=1}^{p} \lambda_{i} v_{i} v_{i}^{T} \succeq t I \\
# & \lambda \succeq 0 \\
# & \mathbf{1}^{T} \lambda=1
# \end{aligned}
# ```

eOpt = Model(SCS.Optimizer)
set_silent(eOpt)
@variable(eOpt, 0 <= np[1:p] <= nmax)
@variable(eOpt, t)
@SDconstraint(eOpt, V * LinearAlgebra.diagm(0 => np ./ n) * V' - (t .* eye) >= 0)
@constraint(eOpt, sum(np) <= n)
@objective(eOpt, Max, t)
optimize!(eOpt)
objective_value(eOpt)
#-
value.(np)

# ## D-optimal design

# The most widely used scalarization is called $D$ -optimal design,
# in which we minimize the determinant of the error covariance matrix $E$.
# This corresponds to designing the experiment to minimize the volume of the
# resulting confidence ellipsoid (for a fixed confidence level). Ignoring the
# constant factor $1 / m$ in $E$, and taking the logarithm of the objective, we
# can pose this problem as convex optimization problem:

# ```math
# \begin{aligned}
# \min & \log \operatorname{det}\left(\sum_{i=1}^{p} \lambda_{i} v_{i} v_{i}^{T}\right)^{-1} \\
# \text{subject to} & \lambda \succeq 0 \\
# & \mathbf{1}^{T} \lambda=1
# \end{aligned}
# ```

dOpt = Model(SCS.Optimizer)
set_silent(dOpt)
@variable(dOpt, np[1:p], lower_bound = 0, upper_bound = nmax)
@variable(dOpt, t)
@objective(dOpt, Max, t)
@constraint(dOpt, sum(np) <= n)
E = V * LinearAlgebra.diagm(0 => np ./ n) * V'
@constraint(
    dOpt,
    [t, 1, (E[i, j] for i in 1:q for j in 1:i)...] in MOI.LogDetConeTriangle(q)
)
optimize!(dOpt)
objective_value(dOpt)
#-
value.(np)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
#!/usr/bin/env julia

# 1.6: Bubble sort
function bubbleSort!(a)
  n = length(a)
  for i in 1:(n - 1)
    for j in 1:(n - i)
      if a[j] > a[j + 1]
        a[j], a[j + 1] = a[j + 1], a[j]
      end
    end
  end
  return a
end

data = [65, 51, 32, 12, 23, 84, 68, 1]
bubbleSort!(data)


text:
#=
This file is auto-generated. Do not edit.
=#

"""Data Structure Operational Cost Data in Three parts fixed, variable cost and start - stop costs."""
mutable struct ThreePartCost <: OperationalCost
    variable::VariableCost
    fixed::Float64
    startup::Float64
    shutdn::Float64
    internal::InfrastructureSystemsInternal
end

function ThreePartCost(variable, fixed, startup, shutdn, )
    ThreePartCost(variable, fixed, startup, shutdn, InfrastructureSystemsInternal())
end

function ThreePartCost(; variable, fixed, startup, shutdn, )
    ThreePartCost(variable, fixed, startup, shutdn, )
end

# Constructor for demo purposes; non-functional.

function ThreePartCost(::Nothing)
    ThreePartCost(;
        variable=VariableCost((0.0, 0.0)),
        fixed=0.0,
        startup=0.0,
        shutdn=0.0,
    )
end

"""Get ThreePartCost variable."""
get_variable(value::ThreePartCost) = value.variable
"""Get ThreePartCost fixed."""
get_fixed(value::ThreePartCost) = value.fixed
"""Get ThreePartCost startup."""
get_startup(value::ThreePartCost) = value.startup
"""Get ThreePartCost shutdn."""
get_shutdn(value::ThreePartCost) = value.shutdn
"""Get ThreePartCost internal."""
get_internal(value::ThreePartCost) = value.internal


text:
using Immerse, Colors
using Base.Test

include("figs.jl")
include("properties.jl")
include("gadfly_panzoom.jl")
include("hittesting.jl")
include("faces.jl")


text:
using Documenter
using BioGraph

makedocs(
    sitename="BioGraph",
    format=Documenter.HTML(),
    modules=[BioGraph],
    pages=[
        "Home" => "index.md",
        "Function" => "function.md"
    ],
    authors="Nguyet Dang, Tuan Do, Francois Sabot and other contributors."
)

deploydocs(
    repo="github.com/nguyetdang/BioGraph.jl.git",
    devbranch = "main",
    devurl = "dev"
)


text:
# RANDOM FEASIBLE POPULATION #####################################################

# generate a totally random population of feasible solutions for the given momkp
function INIT_random_feasible_population(size_population::Int, momkp::_MOMKP; seeding_sols = Vector{Sol}(undef, 0))

    population = Vector{Sol}(undef, size_population)

    K = length(seeding_sols)

    # first insert all the seeding solutions that fit inside the population
    if K > 0
        perm = randperm(K) # shuffle the seeding solutions
        for i in 1:min(K, size_population)
            population[i] = seeding_sols[perm[i]]
        end
    end

    # complete the population by adding random individuals
    for i in K+1:size_population
        population[i] = INIT_random_feasible_solution(momkp)
    end

    # sort!(population, by = sol -> sol.z[1]) # O(nlog(n))

    debug && DEBUG_feasible_solutions(population)

    return population
end

# RANDOM FEASIBLE SOLUTION #####################################################

# generate a random feasible solution for the given momkp
function INIT_random_feasible_solution(momkp::_MOMKP)

    perm_var = randperm(momkp.n) # random order for the variables
    can_be_added = true # true if the current variable fit in the bag
    dim = 1
    sol = SOLUTION_empty_solution(momkp)
    for var in 1:momkp.n
        # make sure the object won't exceed one of the dimensions
        can_be_added = true
        dim = 1
        while can_be_added && dim <= momkp.m
            if sol.w[dim] + momkp.W[dim,perm_var[var]] > momkp.ω[dim]
                can_be_added = false
            end
            dim += 1
        end
        can_be_added && SOLUTION_add_item!(sol, perm_var[var])
    end

    debug && DEBUG_feasible_solution(sol)

    return sol
end


text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

#TODO this is redundant. Many of these systems are already defined in Examples.otherexamples.jl. We should only have one copy

function cooketriplet(::Type{T} = Float64) where {T<:Real}
    return AxisymmetricOpticalSystem{T}(
        DataFrame(
            SurfaceType = ["Object", "Standard", "Standard", "Standard", "Stop", "Standard", "Standard", "Image"],
            Radius = [Inf, 26.777, 66.604, -35.571, 35.571, 35.571, -26.777, Inf],
            Thickness = [Inf, 4.0, 2.0, 4.0, 2.0, 4.0, 44.748, missing],

            Material = [Air, OpticSim.Examples.Examples_N_SK16, Air, OpticSim.Examples.Examples_N_SF2, Air, OpticSim.Examples.Examples_N_SK16, Air, missing],
            SemiDiameter = [Inf, 8.580, 7.513, 7.054, 6.033, 7.003, 7.506, 15.0]
        )
    )
end

function cooketripletfirstelement(::Type{T} = Float64) where {T<:Real}
    return AxisymmetricOpticalSystem{T}(
        DataFrame(
            SurfaceType = ["Object", "Standard", "Standard", "Image"],
            Radius = [Inf, -35.571, 35.571, Inf],
            Thickness = [Inf, 4.0, 44.748, missing],
            Material = [Air, OpticSim.Examples.Examples_N_SK16, Air, missing],
            SemiDiameter = [Inf, 7.054, 6.033, 15.0]
        )
    )
end

function convexplano(::Type{T} = Float64) where {T<:Real}
    return AxisymmetricOpticalSystem{T}(
        DataFrame(
            SurfaceType = ["Object", "Standard", "Standard", "Image"],
            Radius = [Inf, 60.0, Inf, Inf],
            Thickness = [Inf, 10.0, 57.8, missing],
            Material = [Air, OpticSim.Examples.Examples_N_BK7, Air, missing],
            SemiDiameter = [Inf, 9.0, 9.0, 15.0]
        )
    )
end

function doubleconvex(
    frontradius::T, rearradius::T;
    temperature::Unitful.Temperature = GlassCat.TEMP_REF_UNITFUL, pressure::T = T(PRESSURE_REF)
) where {T<:Real}
    return AxisymmetricOpticalSystem{T}(
        DataFrame(
            SurfaceType = ["Object", "Standard", "Standard", "Image"],
            Radius = [T(Inf64), frontradius, rearradius, T(Inf64)],
            Thickness = [T(Inf64), T(10.0), T(57.8), missing],
            Material = [Air, OpticSim.Examples.Examples_N_BK7, Air, missing],
            SemiDiameter = [T(Inf64), T(9.0), T(9.0), T(15.0)]
        );
        temperature,
        pressure
    )
end

function doubleconvex(
    ::Type{T} = Float64; temperature::Unitful.Temperature = GlassCat.TEMP_REF_UNITFUL, pressure::T = T(PRESSURE_REF)
) where {T<:Real}
    return AxisymmetricOpticalSystem{T}(
        DataFrame(
            SurfaceType = ["Object", "Standard", "Standard", "Image"],
            Radius = [Inf64, 60, -60, Inf64],
            Thickness = [Inf64, 10.0, 57.8, missing],
            Material = [Air, OpticSim.Examples.Examples_N_BK7, Air, missing],
            SemiDiameter = [Inf64, 9.0, 9.0, 15.0]
        );
        temperature,
        pressure
    )
end

function doubleconcave(::Type{T} = Float64) where {T<:Real}
    return AxisymmetricOpticalSystem{T}(
        DataFrame(SurfaceType = ["Object", "Standard", "Standard", "Image"],
            Radius = [Inf64, -41.0, 41.0, Inf64],
            Thickness = [Inf64, 10.0, 57.8, missing],
            Material = [Air, OpticSim.Examples.Examples_N_BK7, Air, missing],
            SemiDiameter = [Inf64, 9.0, 9.0, 15.0]
        )
    )
end

function planoconcaverefl(::Type{T} = Float64) where {T<:Real}
    return AxisymmetricOpticalSystem{T}(
        DataFrame(SurfaceType = ["Object", "Standard", "Standard", "Image"],
            Radius = [Inf64, Inf64, -41.0, Inf64],
            Thickness = [Inf64, 10.0, -57.8, missing],
            Material = [Air, OpticSim.Examples.Examples_N_BK7, Air, missing],
            SemiDiameter = [Inf64, 9.0, 9.0, 25.0],
            Reflectance = [missing, missing, 1.0, missing]
        )
    )
end

function concaveplano(::Type{T} = Float64) where {T<:Real}
    return AxisymmetricOpticalSystem{T}(
        DataFrame(SurfaceType = ["Object", "Standard", "Standard", "Image"],
            Radius = [Inf64, -41.0, Inf64, Inf64],
            Thickness = [Inf64, 10.0, 57.8, missing],
            Material = [Air, OpticSim.Examples.Examples_N_BK7, Air, missing],
            SemiDiameter = [Inf64, 9.0, 9.0, 15.0]
        )
    )
end

function planoplano(::Type{T} = Float64) where {T<:Real}
    return AxisymmetricOpticalSystem{T}(
        DataFrame(SurfaceType = ["Object", "Standard", "Standard", "Image"],
            Radius = [Inf64, Inf64, Inf64, Inf64],
            Thickness = [Inf64, 10.0, 57.8, missing],
            Material = [Air, OpticSim.Examples.Examples_N_BK7, Air, missing],
            SemiDiameter = [Inf64, 9.0, 9.0, 15.0]
        )
    )
end

"""
    planarshapes()

Returns a CSGOpticalSystem containing a vertical stack (along the z-axis) of alternating planar shapes and a detector
placed at the end of the stack. This is for the purpose of benchmarking allocations when planar shapes are stored in a
shared Vector within LensAssembly.
"""
function planarshapes()
    interface = FresnelInterface{Float64}(OpticSim.Examples.Examples_N_BK7, Air; interfacemode=Transmit)
    s = 10.0
    surfacenormal = SVector(0., 0., 1.)

    elements = []
    for i = 1:9
        centrepoint = SVector(0., 0., -i)
        if i % 3 == 1
            push!(elements, Ellipse(s, s, surfacenormal, centrepoint; interface))
        elseif i % 3 == 2
            push!(elements, Hexagon(s, surfacenormal, centrepoint; interface))
        elseif i % 3 == 0
            push!(elements, Rectangle(s, s, surfacenormal, centrepoint; interface))
        end
    end

    lensassembly = LensAssembly(elements...)
    detector = Rectangle(s, s, surfacenormal, SVector(0., 0., -10.); interface = opaqueinterface())

    return CSGOpticalSystem(lensassembly, detector)
end

zernikesystem() = CSGOpticalSystem(LensAssembly(zernikelens()), Rectangle(40.0, 40.0, [0.0, 0.0, 1.0], [0.0, 0.0, -67.8], interface = opaqueinterface()))
asphericsystem() = CSGOpticalSystem(LensAssembly(asphericlens()), Rectangle(25.0, 25.0, [0.0, 0.0, 1.0], [0.0, 0.0, -67.8], interface = opaqueinterface()))
evenasphericsystem()=CSGOpticalSystem(LensAssembly(evenasphericlens()), Rectangle(25.0, 25.0, [0.0, 0.0, 1.0], [0.0, 0.0, -67.8], interface = opaqueinterface()))
oddasphericsystem()=CSGOpticalSystem(LensAssembly(oddasphericlens()), Rectangle(25.0, 25.0, [0.0, 0.0, 1.0], [0.0, 0.0, -67.8], interface = opaqueinterface()))
conicsystemA() = CSGOpticalSystem(LensAssembly(coniclensA()), Rectangle(25.0, 25.0, [0.0, 0.0, 1.0], [0.0, 0.0, -67.8], interface = opaqueinterface()))
conicsystemQ() = CSGOpticalSystem(LensAssembly(coniclensQ()), Rectangle(25.0, 25.0, [0.0, 0.0, 1.0], [0.0, 0.0, -67.8], interface = opaqueinterface()))
qtypesystem() = CSGOpticalSystem(LensAssembly(qtypelens()), Rectangle(25.0, 25.0, [0.0, 0.0, 1.0], [0.0, 0.0, -67.8], interface = opaqueinterface()))
chebyshevsystem() = CSGOpticalSystem(LensAssembly(chebyshevlens()), Rectangle(40.0, 40.0, [0.0, 0.0, 1.0], [0.0, 0.0, -67.8], interface = opaqueinterface()))
gridsagsystem() = CSGOpticalSystem(LensAssembly(gridsaglens()), Rectangle(40.0, 40.0, [0.0, 0.0, 1.0], [0.0, 0.0, -67.8], interface = opaqueinterface()))


text:
# julia1_4.1.jl  Basic Stream I/O
#
# well depth example with file i/o
# constants
using Printf
g= 9.807 # acceleration due to gravity m/s^2
vs= 343.21 # velocity of sound in air m/s 20˚C
# (depth,tv)=well(t)  judging well depth using a rock
function well(t)
  tdn = (-vs + sqrt(vs^2 + 2.0*g*vs*t))/g
  tup = t - tdn
  depth = vs*tup; termvel=g*tdn
  return (depth,termvel)
end
t=5.0       # seconds from drop to ear
(d,tv)=well(t)
println(t," ",d,"  ",tv)
#rounding results
rd=(round(d*10))/10.0
rtv=(round(tv*10))/10.0
println("depth $rd meters")
println("terminal velocity $rtv meters/second")
# could also use @printf() macro for C printf
@printf("depth %6.1f meters",d)
@printf("%6.1f,%6.1f,%6.1f\n",t,d,tv)
#
###
#
win=open("wellin.dat","r")
wout=open("wellout.csv","w")
@printf(wout,"   Time  Depth Speed  Place\n")
@printf(wout,"   sec     M    M/s   CA/USA\n")
while true
  s=readline(win)
  if(s==""); break; end
  st=split(s)
  t=parse(Float64,st[1])
  (d,tv)=well(t)
  @printf(wout,"%6.1f,%6.1f,%6.1f %-22s\n",t,d,tv,st[2])
end
println("done")
close(win)
close(wout)
# other printf formats designators:
#
# many system-like directory functions
# examples:
pwd() # path to current working directory
readdir() # list of files/directories in working directory

stat("wellin.dat") # yields string, list of file parameters
# mode=6 Owner R&W, =4 Group R only, =4 Others R only
# uperm(file) user permission, gperm(file) group permission
# operm(file) other permission

chmod("wellout.csv",0o664)
stat("wellout.csv")


text:
import Dates
import MLJMultivariateStatsInterface: _replace!
import MultivariateStats
import Random

using LinearAlgebra
using MLJBase
using MLJMultivariateStatsInterface
using StableRNGs
using Test

const MS = MultivariateStats

include("testutils.jl")
println("\nutils"); include("utils.jl")
println("\ncomponent_analysis"); include("models/decomposition_models.jl")
println("\ndiscriminant_analysis"); include("models/discriminant_analysis.jl")
println("\nlinear_models"); include("models/linear_models.jl")


text:
module LintExample3

include(joinpath(Pkg.dir(), "no_such_file"))

end