text:
using LinearAlgebra, ForwardDiff

# test Mehrotra predictor-corrector algorithm
n = 10
me = 5
mi = n
p = 0

idx_ineq = collect(1:mi)

x0 = rand(n)
p0 = zeros(p)
a0 = zeros(mi)

# objective
c = rand(n)
function obj(x)
    transpose(c) * x
end

# constraints
A = rand(me, n)
b = A * x0
function con(x)
    A * x - b
end

# unpack
function unpack(w)
    x = view(w, 1:n)
    y = view(w, n .+ (1:me))
    z = view(w, n + me .+ (1:mi))

    return x, y, z
end

# barrier check
function check_variables(w)
    x, y, z = unpack(w)
    if any(x .<= 0.0)
        return true
    end
    if any(z .<= 0.0)
        return true
    end
    return false
end

function r(w, θ, a)
    x, y, z = unpack(w)

    lag = c + transpose(A) * y
    lag[idx_ineq] -= z

    c_eq = con(x)

    comp = x[idx_ineq] .* z
    comp -= a

    return [lag; c_eq; comp]
end

function Rz(w, θ)
    return ForwardDiff.jacobian(q -> r(q, θ, a0), w)
end

function Rθ(w, θ)
    return ForwardDiff.jacobian(q -> r(w, q, a0), θ)
end

function interior_point(x0, θ;
    x_init_tol = 0.1,
    r_tol = 1.0e-8,
    μ_tol = 1.0e-5,
    max_iter_inner = 100,
    max_iter_ls = 25)

    # initialization
    x = max.(x0, x_init_tol)
    y = zeros(me)
    z = ones(mi)
    w = [x; y; z]
    μ = 1.0
    e = ones(mi)

    flag = false
    total_iter = 0

    while μ > μ_tol
        for i = 1:max_iter_inner

            # compute residual, residual Jacobian
            res = r(w, θ, μ * e)
            res_norm = norm(res)

            # println("iteration: $i, residual: $res_norm")

            if res_norm < r_tol
                flag = true
                continue
            end

            jac = Rz(w, θ)

            # compute step
            Δ = -jac \ res

            # line search the step direction
            α = 1.0

            iter = 0
            while check_variables(w + α * Δ) # backtrack inequalities
                α = 0.5 * α
                # println("   α = $α")
                iter += 1
                if iter > max_iter_ls
                    @error "backtracking line search fail"
                    flag = false
                    return w, zeros(n + me + mi, p), false, total_iter
                end
            end

            while norm(r(w + α * Δ, θ, μ * e))^2.0 >= (1.0 - 0.001 * α) * res_norm^2.0
                α = 0.5 * α
                # println("   α = $α")
                iter += 1
                if iter > max_iter_ls
                    @error "line search fail"
                    flag = false
                    return w, zeros(n + me + mi, p), false, total_iter
                end
            end

            # update
            w .+= α * Δ

            total_iter += 1
        end

        μ = 0.1 * μ
        println("μ: $μ")
    end

    δw = -1.0 * Rz(w, θ) \ Rθ(w, θ)

    return w, δw, flag, total_iter
end

w_soli, δw_soli, flag, total_iter = interior_point(rand(n), p0)
@show total_iter

function mehrotra(x0, θ;
    x_init_tol = 0.1,
    r_tol = 1.0e-8,
    μ_tol = 1.0e-5,
    max_iter_inner = 100,
    max_iter_ls = 25)

    # initialization
    x = max.(x0, x_init_tol)
    y = zeros(me)
    z = ones(mi)
    w = [x; y; z]
    μ = 1.0
    e = ones(mi)

    flag = false
    total_iter = 0

    while μ > μ_tol
        for i = 1:max_iter_inner

            # unpack
            x, y, z = unpack(w)

            # compute residual, residual Jacobian
            res_aff = r(w, θ, a0)
            res_aff_norm = norm(res_aff, 1)

            # println("iteration: $i, residual: $res_norm")

            if res_aff_norm < r_tol
                flag = true
                continue
            end

            jac = Rz(w, θ)

            # compute affine step
            Δ_aff = -jac \ res_aff

            Δx_aff, Δy_aff, Δz_aff = unpack(Δ_aff)

            α_aff_pr = min(1.0, minimum([Δxi < 0.0 ? -x[i] / Δxi : Inf for (i, Δxi) in enumerate(Δx_aff)]))
            α_aff_du = min(1.0, minimum([Δzi < 0.0 ? -z[i] / Δzi : Inf for (i, Δzi) in enumerate(Δz_aff)]))
            μ_aff = transpose(x[idx_ineq] + α_aff_pr * Δx_aff[idx_ineq]) * (z + α_aff_du * Δz_aff) / mi

            # compute corrector step
            σ = (μ_aff / μ)^3.0

            res_corr = r(w, θ, -Δx_aff[idx_ineq] .* Δz_aff + σ * μ * e)

            Δ_corr = -jac \ res_corr

            # line search the step direction
            η = 1.0
            α_pr = α_aff_pr
            α_du = α_aff_du

            α = η .* [α_pr * ones(n); α_du * ones(me + mi)]

            iter = 0
            # while check_variables(w + α .* Δ_corr) # backtrack inequalities
            #     η = 0.99 * η
            #     α = η .* [α_pr * ones(n); α_du * ones(me + mi)]
            #
            #     # println("   α = $α")
            #     iter += 1
            #     if iter > max_iter_ls
            #         @error "backtracking line search fail"
            #         flag = false
            #         return w, zeros(n + me + mi, p), false, total_iter
            #     end
            # end

            # while norm(r(w + α .* Δ_corr, θ, a0))^2.0 >= (1.0 - 0.001 * η) * res_aff_norm^2.0
            #     η = 0.99 * η
            #     α = η .* [α_pr * ones(n); α_du * ones(me + mi)]
            #     # println("   α = $α")
            #     iter += 1
            #     if iter > max_iter_ls
            #         @error "line search fail"
            #         flag = false
            #         return w, zeros(n + me + mi, p), false, total_iter
            #     end
            # end

            # update
            w .+= α .* Δ_corr

            total_iter += 1
        end

        μ = 0.1 * μ
        println("μ: $μ")
    end

    δw = -1.0 * Rz(w, θ) \ Rθ(w, θ)

    return w, δw, flag, total_iter
end

w_solm, δw_solm, flag, total_iter = mehrotra(x0, p0)
@show total_iter

obj(w_soli[1:n])
obj(w_solm[1:n])

con(w_soli[1:n])
con(w_solm[1:n])

max_issues_repo_name: thowell/DirectMotionPlanning


text:
module Unwrap

export unwrap1d

# Unwrap two neighbouring values
function γ1d(a,b)
    d = a-b
    if abs(d) > π
        d-round(d/2π)*2π
    else
        d
    end
end

# 1d vector unwrapping, seems faster than unwrap() from DSP.jl
function unwrap1d(v::AbstractVector)
    u = Array{eltype(v)}(length(v))
    u[1] = v[1]
    for i = 2:length(v)
        u[i] = u[i-1] + γ1d(v[i],v[i-1])
    end
    u
end

function unwrap1d(M::AbstractMatrix, dim=1)
    m,n = size(M)
    U = Matrix{eltype(M)}(m, n)
    if dim == 1
        for j in 1:n
            U[:,j] = unwrap1d(M[:,j])
        end
    else
        for j in 1:n
            U[i,:] = unwrap1d(M[1,:])
        end
    end
    U
end

include("unwrap2d.jl")

end # module

ext: jl


text:
"""
    Checkpointer{I, T, P} <: AbstractOutputWriter

An output writer for checkpointing models to a JLD2 file from which models can be restored.
"""
mutable struct Checkpointer{I, T, P} <: AbstractOutputWriter
    iteration_interval :: I
         time_interval :: T
              previous :: Float64
                   dir :: String
                prefix :: String
            properties :: P
                 force :: Bool
               verbose :: Bool
end

"""
    Checkpointer(model; iteration_interval=nothing, time_interval=nothing, dir=".",
                 prefix="checkpoint", force=false, verbose=false,
                 properties = [:architecture, :boundary_conditions, :grid, :clock, :coriolis,
                               :buoyancy, :closure, :velocities, :tracers, :timestepper])

Construct a `Checkpointer` that checkpoints the model to a JLD2 file every so often as
specified by `iteration_interval` or `time_interval`. The `model.clock.iteration` is included
in the filename to distinguish between multiple checkpoint files.

Note that extra model `properties` can be safely specified, but removing crucial properties
such as `:velocities` will make restoring from the checkpoint impossible.

The checkpoint file is generated by serializing model properties to JLD2. However,
functions cannot be serialized to disk (at least not with JLD2). So if a model property
contains a reference somewhere in its hierarchy it will not be included in the checkpoint
file (and you will have to manually restore them).

Keyword arguments
=================
- `iteration_interval`: Save output every `n` model iterations.
- `time_interval`: Save output every `t` units of model clock time.
- `dir`: Directory to save output to. Default: "." (current working directory).
- `prefix`: Descriptive filename prefixed to all output files. Default: "checkpoint".
- `force`: Remove existing files if their filenames conflict. Default: `false`.
- `verbose`: Log what the output writer is doing with statistics on compute/write times
  and file sizes. Default: `false`.
- `properties`: List of model properties to checkpoint. Some are required.
"""
function Checkpointer(model; iteration_interval=nothing, time_interval=nothing,
                      dir=".", prefix="checkpoint", force=false, verbose=false,
                      properties = [:architecture, :grid, :clock, :coriolis,
                                    :buoyancy, :closure, :velocities, :tracers, :timestepper])

    validate_intervals(iteration_interval, time_interval)

    # Certain properties are required for `restore_from_checkpoint` to work.
    required_properties = (:grid, :architecture, :velocities, :tracers, :timestepper)
    for rp in required_properties
        if rp ∉ properties
            @warn "$rp is required for checkpointing. It will be added to checkpointed properties"
            push!(properties, rp)
        end
    end

    for p in properties
        p isa Symbol || error("Property $p to be checkpointed must be a Symbol.")
        p ∉ propertynames(model) && error("Cannot checkpoint $p, it is not a model property!")

        if has_reference(Function, getproperty(model, p)) && (p ∉ required_properties)
            @warn "model.$p contains a function somewhere in its hierarchy and will not be checkpointed."
            filter!(e -> e != p, properties)
        end
    end

    mkpath(dir)
    return Checkpointer(iteration_interval, time_interval, 0.0, dir, prefix, properties, force, verbose)
end

function write_output(model, c::Checkpointer)
    filepath = joinpath(c.dir, c.prefix * "_iteration" * string(model.clock.iteration) * ".jld2")
    c.verbose && @info "Checkpointing to file $filepath..."

    t1 = time_ns()
    jldopen(filepath, "w") do file
        file["checkpointed_properties"] = c.properties
        serializeproperties!(file, model, c.properties)
    end

    t2, sz = time_ns(), filesize(filepath)
    c.verbose && @info "Checkpointing done: time=$(prettytime((t2-t1)/1e9)), size=$(pretty_filesize(sz))"
end

# This is the default name used in the simulation.output_writers ordered dict.
defaultname(::Checkpointer, nelems) = :checkpointer

function restore_if_not_missing(file, address)
    if !ismissing(file[address])
        return file[address]
    else
        @warn "Checkpoint file does not contain $address. Returning missing. " *
              "You might need to restore $address manually."
        return missing
    end
end

function restore_field(file, address, arch, grid, loc, kwargs)
    field_address = file[address * "/location"]
    data = OffsetArray(convert_to_arch(arch, file[address * "/data"]), grid, loc)

    # Extract field name from address. We use 2:end so "tracers/T "gets extracted
    # as :T while "timestepper/Gⁿ/T" gets extracted as :Gⁿ/T (we don't want to
    # apply the same BCs on T and T tendencies).
    field_name = split(address, "/")[2:end] |> join |> Symbol

    # If the user specified a non-default boundary condition through the kwargs
    # in restore_from_checkpoint, then use them when restoring the field. Otherwise
    # restore the BCs from the checkpoint file as long as they're not missing.
    if :boundary_conditions in keys(kwargs) && field_name in keys(kwargs[:boundary_conditions])
        bcs = kwargs[:boundary_conditions][field_name]
    else
        bcs = restore_if_not_missing(file, address * "/boundary_conditions")
    end

    return Field(field_address, arch, grid, bcs, data)
end

const u_location = (Face, Cell, Cell)
const v_location = (Cell, Face, Cell)
const w_location = (Cell, Cell, Face)
const c_location = (Cell, Cell, Cell)

"""
    restore_from_checkpoint(filepath; kwargs=Dict())

Restore a model from the checkpoint file stored at `filepath`. `kwargs` can be passed
to the model constructor, which can be especially useful if you need to manually
restore forcing functions or boundary conditions that rely on functions.
"""
function restore_from_checkpoint(filepath; kwargs...)
    kwargs = length(kwargs) == 0 ? Dict{Symbol,Any}() : Dict{Symbol,Any}(kwargs)

    file = jldopen(filepath, "r")
    cps = file["checkpointed_properties"]

    if haskey(kwargs, :architecture)
        arch = kwargs[:architecture]
        filter!(p -> p ≠ :architecture, cps)
    else
        arch = file["architecture"]
    end

    grid = file["grid"]

    # Restore velocity fields
    kwargs[:velocities] =
        VelocityFields(arch, grid, u = restore_field(file, "velocities/u", arch, grid, u_location, kwargs),
                                   v = restore_field(file, "velocities/v", arch, grid, v_location, kwargs),
                                   w = restore_field(file, "velocities/w", arch, grid, w_location, kwargs))
    filter!(p -> p ≠ :velocities, cps) # pop :velocities from checkpointed properties

    # Restore tracer fields
    tracer_names = Tuple(Symbol.(keys(file["tracers"])))
    tracer_fields = Tuple(restore_field(file, "tracers/$c", arch, grid, c_location, kwargs) for c in tracer_names)
    tracer_fields_kwargs = NamedTuple{tracer_names}(tracer_fields)
    kwargs[:tracers] = TracerFields(arch, grid, tracer_names; tracer_fields_kwargs...)

    filter!(p -> p ≠ :tracers, cps) # pop :tracers from checkpointed properties

    # Restore time stepper tendency fields
    field_names = (:u, :v, :w, tracer_names...) # field names
    locs = (u_location, v_location, w_location, Tuple(c_location for c in tracer_names)...) # name locations

    G⁻_fields = Tuple(restore_field(file, "timestepper/G⁻/$(field_names[i])", arch, grid, locs[i], kwargs) for i = 1:length(field_names))
    Gⁿ_fields = Tuple(restore_field(file, "timestepper/Gⁿ/$(field_names[i])", arch, grid, locs[i], kwargs) for i = 1:length(field_names))

    G⁻_tendency_field_kwargs = NamedTuple{field_names}(G⁻_fields)
    Gⁿ_tendency_field_kwargs = NamedTuple{field_names}(Gⁿ_fields)

    # Restore time stepper
    kwargs[:timestepper] =
        QuasiAdamsBashforth2TimeStepper(eltype(grid), arch, grid, kwargs[:velocities], tracer_names;
                                        G⁻ = TendencyFields(arch, grid, tracer_names; G⁻_tendency_field_kwargs...),
                                        Gⁿ = TendencyFields(arch, grid, tracer_names; Gⁿ_tendency_field_kwargs...))

    filter!(p -> p ≠ :timestepper, cps) # pop :timestepper from checkpointed properties

    # Restore the remaining checkpointed properties
    for p in cps
        kwargs[p] = file["$p"]
    end

    close(file)

    model = IncompressibleModel(; kwargs...)

    return model
end


text:
# TODO: use parametric types with the type parameters corresponding
# to the lexer and parser

get_frontend(parser,lexer) = begin
    if parser=="lalr"
        if lexer === nothing
            throw(error("The LALR parser requires use of a lexer"))
        elseif lexer == "standard"
            return LALR_TraditionalLexer
        elseif lexer == "contextual"
            return LALR_ContextualLexer
        else
            throw(error("Unknown lexer: $lexer"))
        end
    else
        throw(error("Unknown parser: $parser"))
    end
end

#TODO: our transformers are structured differently
#==
_get_lexer_callbacks(transformer,terminals) = begin
    result = Dict()
    for terminal in terminals
        callback = getproperty(transformer, Symbol(terminal.name), nothing)
        if callback != nothing
            result[terminal.name] = callback
        end
    end
    return result
end
==#

struct PostLexConnector
    lexer
    postlexer
end

make_lexer_state(pc::PostLexConnector,text) = make_lexer_state(pc.lexer,text)

# TODO: use an iteration approach
lex(pc::PostLexConnector,lexer_state,parser_state) = begin
    i = lex(pc.lexer,lexer_state,parser_state)
    return process(pc.postlexer,i)
end

abstract type _ParserFrontend end
abstract type WithLexer <: _ParserFrontend end

_parse(pf::_ParserFrontend,start,input,args...) = begin
    if start === nothing
        start = pf.start
        if length(start) > 1
            throw(error("Lark initialized with more than 1 possible start rule. Must specify which start rule to parse: $start"))
        end
        start = start[]
    end
    p = parse(pf.parser,input, start,args...)
    if pf.after_transform !== nothing
        return transform(pf.after_transform,p)
    end
    return p
end

make_lexer(wl::WithLexer,text) = begin
    lexer = wl.lexer
    if wl.postlex != nothing
        lexer = PostLexConnector(wl.lexer,wl.postlex)
    end
    return LexerThread(lexer,text)
end

parse(wl::WithLexer,text;start=nothing) = begin
    _parse(wl,start,make_lexer(wl,text))
end

abstract type LALR_WithLexer <: WithLexer end

struct LALR_TraditionalLexer <: LALR_WithLexer
    lexer
    lexer_conf
    postlex
    parser
    start
    after_transform
end

LALR_TraditionalLexer(lexer_conf,parser_conf;options=nothing) = begin
    debug = options !== nothing ? options.debug : false
    LALR_TraditionalLexer(TraditionalLexer(lexer_conf), lexer_conf,
                          lexer_conf.postlex,
                          LALRParser(parser_conf,debug=debug),
                          parser_conf.start,
                          options !== nothing ? options.transformer : nothing
                          )
end

struct LALR_ContextualLexer <: LALR_WithLexer
    lexer
    lexer_conf
    postlex
    parser
    start
    after_transform
end

LALR_ContextualLexer(lexer_conf,parser_conf;options=nothing) = begin
    debug = options !== nothing ? options.debug : false
    parser = LALRParser(parser_conf,debug=debug)
    postlex = lexer_conf.postlex
    states = Dict([i=>collect(keys(t)) for (i,t) in enumerate(parser._parse_table.states)])
    always_accept = postlex !== nothing ? postlex.always_accept : ()
    lexer = ContextualLexer(lexer_conf,states,always_accept=always_accept)
    LALR_ContextualLexer(lexer,lexer_conf,postlex,parser,parser_conf.start,
                         options !== nothing ? options.transformer : nothing)
end

init_contextual_lexer(p::WithLexer,lexer_conf) = begin
    p.lexer_conf = lexer_conf
    states = Dict(idx=>collect(keys(t)) for (idx,t) in p.parser._parse_table.states)
    always_accept = if !isnothing(lexer_conf.postlex)
        lexer_conf.postlex.always_accept
    else
        ()
    end
    state_channel = Channel(0)
    p.lexer = ContextualLexer(lexer_conf.tokens,states,state_channel,
                              ignore=lexer_conf.ignore,
                              always_accept=always_accept,
                              user_callbacks=lexer_conf.callbacks)
end

# Note we don't have to pass the match method as multiple despatch
# will pick the right match to use
struct Earley <: WithLexer
    lexer_conf
    lexer
    parser
    Earley(lexer_conf,parser_conf,options) = begin
        x = new()
        init_traditional_lexer(x,lexer_conf)
        x.parser = earley.Parser(parser_conf,match,get_ambiguity_options(options)...)
        return x
    end
end

match(e::Earley,term,token) = term.name == token.type_


text:
#!/usr/bin/env bash
    #=
    exec julia --project="$(realpath $(dirname $(dirname $0)))" --color=yes --startup-file=no -e "include(popfirst!(ARGS))" \
    "${BASH_SOURCE[0]}" "$@"
    =#

include("distance.jl")
include("powers.jl")

"""
```julia
push_if_allowed!(C::AbstractArray{T}, w::T, d::Int)
```

Takes in an array and a word.  As long as the word does not mean that the distance is smaller than d, we add w to the array.  If we are successful in doing this, return true.  Otherwise, return false.  *This is a mutating function.*
"""
function push_if_allowed!(C::AbstractArray{T}, w::T, d::Int) where T <: AbstractWord
	isempty(C) && (push!(C, w); return true)

	for c in C
		if hamming_distance(c, w) < d
			return false
		end
	end

	push!(C, w)
	return true
end

"""
```julia
push_if_allowed!(C::AbstractArray{T}, C′::AbstractArray{T}, w::T, d::Int)
```

Takes in two arrays, A and B.  If w is allowed in C given distance d, push to C′.  If we are successful in doing this, return true.  Otherwise, return false.  *This is a mutating function.*
"""
function push_if_allowed!(C::AbstractVector{T}, C′::AbstractVector{T}, w::T, d::Int) where T <: AbstractWord
	isempty(C) && (push!(C′, w); return true)

	for c in C
		if hamming_distance(c, w) < d
			return false
		end
	end

	push!(C′, w)
	return true
end

"""
```julia
replace_if_allowed!(C::AbstractArray, d::Int, w, w′) -> Bool
```

Takes in an array and a word.  As long as the word does not mean that the distance is smaller than d, we replace a with b in the array.  Replaces and returns true if allowed; otherwise returns false.  *This is a mutating function.*
"""
function replace_if_allowed!(C::AbstractVector{T}, d::Int, ws::Pair) where T <: AbstractWord
	w, w′ = ws

	for c in C
		if hamming_distance(c, w′) < d
			return false
		end
	end

	replace!(C, ws)
	return true
end
function replace_if_allowed!(C::AbstractVector{T}, d::Int, w::T, w′::T) where T <: AbstractWord
	return replace_if_allowed!(C, d, Pair(w, w′))
end

_mutate_codeword(w::Word{N, T}, i::Int, a::T) where {T, N} =
	setindex!(w, a, i)

"""
```julia
mutate_codeword(w::NonStaticAbstractWord{N, T}, n::Int, i::Int, a::T) where {T, N} -> MVector{N, T}
mutate_codeword(w::Word{N, T}, n::Int, i::Int, a::T) where {T, N} -> MVector{N, T}
```

Mutates the word w, which is an `MVector` of length N, changing its iᵗʰ index to a.
"""
mutate_codeword(w::NonStaticAbstractWord{N, T}, i::Int, a::T) where {T, N} =
	_mutate_codeword(Word{N, T}(w), i, a)
mutate_codeword(w::Word{N, T}, n::Int, i::Int, a::T) where {T, N} =
	_mutate_codeword(w, i, a)

"""
```julia
get_all_words(Σ::Alphabet{N}, q::Int, n::Int) -> Codewords{M}
get_all_words(Σ::Alphabet{N}, n::Int) -> Codewords{M}
get_all_words(Σ::AbstractArray, q::Int, n::Int) -> Codewords{M}
get_all_words(Σ::AbstractArray, n::Int) -> Codewords{M}
get_all_words(q::Int, n::Int) -> Codewords{M}
```

Get the universe of *all* codewords of a given alphabet.  The alphabet will be uniquely generated if none is given.

Parameters:
  - `Σ::AbstractArray`: The alphabet allowed.
  - `q::Int`: The size of the alphabet.
  - `n::Int`: The (fixed) length of the words in the code.
  - `d::Int`: The minimum distance between words in the code.
  - `𝒰::AbstractArray`: The universe of all codewords of q many letters of block length n.

Returns:
  - `Codewords{M}`: An array of codewords, each of length `M`.  Each codewords is a tuple, and each character in said word is a symbol.

---

### Examples

```julia
julia> get_all_words(2, 2) # all words of block length 2 using 2 unique symbols
2×2 Array{Tuple{Symbol,Symbol},2}:
 (Symbol("##254"), Symbol("##254"))  (Symbol("##254"), Symbol("##253"))
 (Symbol("##253"), Symbol("##254"))  (Symbol("##253"), Symbol("##253"))
```
"""
get_all_words(Σ::Alphabet{N}, q::Int, n::Int) where {N} =
	collect(CodeUniverseIterator(UniverseParameters(Σ, q, n)))
get_all_words(Σ::Alphabet{N}, n::Int) where {N} =
	get_all_words(Σ, length(Set(Σ)), n) # if alphabet is given, then q is the length of that alphabet
get_all_words(Σ::AbstractArray, n::Int) =
	get_all_words(Alphabet(Σ), length(Set(Σ)), n) # if alphabet is given, then q is the length of that alphabet
get_all_words(q::Int, n::Int) =
	get_all_words(genalphabet(q), q, n) # generate symbols if no alphabet is given

"""
```julia
get_codewords_greedy(𝒰::UniverseParameters, d::Int) -> Codewords{M}
get_codewords_greedy(Σ::Alphabet{N}, q::Int, n::Int, d::Int) -> Codewords{M}
get_codewords_greedy(Σ::Alphabet{N}, n::Int, d::Int) -> Codewords{M}
get_codewords_greedy(q::Int, n::Int, d::Int) -> Codewords{M}
get_codewords_greedy(Σ::AbstractArray, q::Int, n::Int, d::Int) -> Codewords{M}
get_codewords_greedy(Σ::AbstractArray, n::Int, d::Int) -> Codewords{M}
get_codewords_greedy(Σ::AbstractArray, q::Int, n::Int, d::Int, 𝒰::AbstractArray) -> Codewords{M}
```

Search through the universe of all codewords and find a code of block length n and distance d, using the alphabet Σ.  The alphabet will be uniquely generated if none is given.  This uses a greedy algorithm, simply iterating through all words (see above) and choosing them if they fit in the code.  In some cases the greedy algorithm is the best, but in others it is very much not.

Parameters:
  - `𝒰::UniverseParameters`: The parameters of the universe of all codewords of q many letters of block length n.
  - `Σ::AbstractArray`: The alphabet allowed.
  - `q::Int`: The size of the alphabet.
  - `n::Int`: The (fixed) length of the words in the code.
  - `d::Int`: The minimum distance between words in the code.

Returns:
  - `Codewords{M}`: An array of codewords, each of length `M`.  Each codewords is a tuple, and each character in said word is a symbol.
"""
function get_codewords_greedy(𝒰::UniverseParameters, d::Int)
	C = eltype(𝒰)[]

	for wᵢ in CodeUniverseIterator(𝒰)
		push_if_allowed!(C, wᵢ, d)
	end

	return C
end
get_codewords_greedy(Σ::Alphabet{N}, q::Int, n::Int, d::Int) where {N} =
	get_codewords_greedy(UniverseParameters(Σ, q, n), d)
get_codewords_greedy(Σ::Alphabet{N}, n::Int, d::Int) where {N} =
	get_codewords_greedy(UniverseParameters(Σ, n), d)
get_codewords_greedy(q::Int, n::Int, d::Int) =
	get_codewords_greedy(UniverseParameters(q, n), d)
get_codewords_greedy(Σ::AbstractArray, q::Int, n::Int, d::Int) =
	get_codewords_greedy(Alphabet(Σ), q, n, d)
get_codewords_greedy(Σ::AbstractArray, n::Int, d::Int) =
	get_codewords_greedy(UniverseParameters(Σ, n), d)
get_codewords_greedy(Σ::AbstractArray, q::Int, n::Int, d::Int, 𝒰::AbstractArray) =
	get_codewords_greedy(Alphabet(Σ), q, n, d, 𝒰)

argmaxminima(A::AbstractArray; dims::Int) = getindex(argmin(A, dims = dims), argmax(argmin(A, dims = dims)))
maxminima(A::AbstractArray; dims::Int) = getindex(minimum(A, dims = dims), maximum(minimum(A, dims = dims)))
argminmaxima(A::AbstractArray; dims::Int) = getindex(argmax(A, dims = dims), argmin(argmax(A, dims = dims)))
minmaxima(A::AbstractArray; dims::Int) = getindex(maximum(A, dims = dims), minimum(maximum(A, dims = dims)))

"""
```julia
get_codewords_random(𝒰::UniverseParameters, d::Int; m::Int = 1000) -> Codewords{M}
get_codewords_random(Σ::Alphabet{N}, q::Int, n::Int, d::Int; m::Int=1000) -> Codewords{M}
get_codewords_random(Σ::Alphabet{N}, n::Int, d::Int; m::Int=1000)	-> Codewords{M}
get_codewords_random(q::Int, n::Int, d::Int; m::Int=1000) -> Codewords{M}
get_codewords_random(Σ::AbstractArray, q::Int, n::Int, d::Int; m::Int=1000) -> Codewords{M}
get_codewords_random(Σ::AbstractArray, n::Int, d::Int; m::Int=1000) -> Codewords{M}
get_codewords_random(Σ::AbstractArray, q::Int, n::Int, d::Int, 𝒰::AbstractArray; m::Int=1000) -> Codewords{M}
```

Search through the universe of all codewords at random and find a code of block length n and distance d, using the alphabet Σ.  The alphabet will be uniquely generated if none is given.  This is a cleverer algorithm than the greedy algorithm.  Increasing the `m` keyword argument arbitrarily _should_ produce a maximal code, as for each codeword it chooses, it collects a list of `m` many random words, and chooses the best one from that intermediate list.

Parameters:
  - `Σ::AbstractArray`: The alphabet allowed.
  - `q::Int`: The size of the alphabet.
  - `n::Int`: The (fixed) length of the words in the code.
  - `d::Int`: The minimum distance between words in the code.
  - `𝒰::AbstractArray`: The universe of all codewords of q many letters of block length n.

Returns:
  - `Codewords{M}`: An array of codewords, each of length `M`.  Each codewords is a tuple, and each character in said word is a symbol.
"""
function get_codewords_random(𝒰::UniverseParameters, d::Int; m::Int = 1000)
	C = eltype(𝒰)[]

	starting_word = rand(𝒰) # get a random word in the code start
	push!(C, starting_word)

	for _ in 1:length(𝒰)
		C′ = eltype(𝒰)[]
		# while length(C′) < m
		for _ in 1:m
			wᵣ = rand(𝒰)
			push_if_allowed!(C, C′, wᵣ, d) # if allowed in C, push to C′
			# wᵣ ∉ C′ && push_if_allowed!(C, C′, wᵣ, d)
		end
		isempty(C′) && break
		# [push_if_allowed!(C, C′, w, d) for _ in 1:m]
		distances = Int[hamming_distance(wᵢ, wⱼ) for wᵢ in C, wⱼ in C′]
		best_word = getindex(C′, getindex(argmaxminima(distances, dims = 1), 2))
		push!(C, best_word)
	end

	return C
end
get_codewords_random(Σ::Alphabet{N}, q::Int, n::Int, d::Int; m::Int=1000) where {N} =
	get_codewords_random(UniverseParameters(Σ, q, n), d, m=m)
get_codewords_random(Σ::Alphabet{N}, n::Int, d::Int; m::Int=1000) where {N} =
	get_codewords_random(UniverseParameters(Σ, n), d, m=m)
get_codewords_random(q::Int, n::Int, d::Int; m::Int=1000) =
	get_codewords_random(UniverseParameters(q, n), d, m=m)
get_codewords_random(Σ::AbstractArray, q::Int, n::Int, d::Int; m::Int=1000) =
	get_codewords_random(Alphabet(Σ), q, n, d, m=m)
get_codewords_random(Σ::AbstractArray, n::Int, d::Int; m::Int=1000) =
	get_codewords_random(UniverseParameters(Σ, n), d, m=m)
get_codewords_random(Σ::AbstractArray, q::Int, n::Int, d::Int, 𝒰::AbstractArray; m::Int=1000) =
	get_codewords_random(Alphabet(Σ), q, n, d, 𝒰, m=m)

# using Mmap
# function get_codewords_random_mmap(mmappath::AbstractString, 𝒰::UniverseParameters, d::Int)
# 	io = open(mmappath, "r+") # allow read and write
# 	# write(s, size(A,2))
# 	# read(s, Int)
# 	# close(s)
# 	# B = Mmap.mmap(io, BitArray, (25,30000))
# 	# Mmap.sync!(B);
# 	# close(io);
# 	# rm("mmap.bin")
# 	starting_word = rand(𝒰) # get a random word in the code start
# 	push!(C, starting_word)
#
# 	for _ in 1:length(𝒰)
# 		C′ = Tuple[]
# 		for _ in 1:m
# 			push_if_allowed!(C, C′, rand(𝒰), d) # if allowed in C, push to C′
# 		end
# 		isempty(C′) && break
# 		# [push_if_allowed!(C, C′, w, d) for _ in 1:m]
# 		distances = [hamming_distance(wᵢ, wⱼ) for wᵢ in C, wⱼ in C′]
# 		best_word = getindex(C′, getindex(argmaxminima(distances, dims = 1), 2))
# 		push!(C, best_word)
# 	end
#
# 	return C
# end

# get_codewords_random(𝒰::UniverseParameters, d::Int) = get_codewords_random(joinpath(tempdir(), "mmap.bin"), 𝒰, d)

"""
```julia
get_codewords(𝒰::UniverseParameters, d::Int; m::Int=10, m_random::Int = 1000) -> Codewords{M}
get_codewords(Σ::Alphabet{N}, q::Int, n::Int, d::Int; m::Int=10, m_random::Int=1000) -> Codewords{M}
get_codewords(q::Int, n::Int, d::Int; m::Int=10, m_random::Int=1000) -> Codewords{M}
get_codewords(q::Int, n::Int, d::Int; m::Int=10, m_random::Int=1000) -> Codewords{M}
get_codewords(Σ::AbstractArray, q::Int, n::Int, d::Int; m::Int=10, m_random::Int=1000) -> Codewords{M}
get_codewords(Σ::AbstractArray, n::Int, d::Int; m::Int=10, m_random::Int=1000) -> Codewords{M}
get_codewords(Σ::AbstractArray, q::Int, n::Int, d::Int, 𝒰::AbstractArray; m::Int=10, m_random::Int=1000) -> Codewords{M}
```

Use function `get_codewords_random` `m` many times (with `get_codewords_random(..., m = m_random)`), and `get_codewords_greedy`.  Return the code with the greatest number of words.  The alphabet will be uniquely generated if none is given.  You can omit Σ and 𝒰.  You can omit q if Σ is given.

Parameters:
  - `Σ::AbstractArray`: The alphabet allowed.
  - `q::Int`: The size of the alphabet.
  - `n::Int`: The (fixed) length of the words in the code.
  - `d::Int`: The minimum distance between words in the code.
  - `𝒰::AbstractArray`: The universe of all codewords of q many letters of block length n.
  - `m::Int` (kwarg): Try a random code m many times.
  - `m_random::Int` (kwarg): The number of possible words `get_codewords_random` chooses from for _each_ word it selects.

Returns:
  - `Codewords{M}`: An array of codewords, each of length `M`.  Each codewords is a tuple, and each character in said word is a symbol.

!!! note

    *If you are looking for a _maximal_ code, this is likely the function you need.*  Increasing `m` and `m_random` arbitrarily should ensure a maximal code—_however_, that computing power/time in not always possible, as it requires a lot of RAM to store certain codes in memeory.  Efforts are being made to make this process better by using memory-mapped filed instead of storing codewords in RAM, but this will make it much slower as well.  Help with this would be much appreciated.

---

### Examples

```julia
julia> get_codewords(["a", "b", "c"], 3, 2) # get codewords of block length 3 with distance 2.  Once again, are symbols for uniqueness
9-element Array{Tuple{Symbol,Symbol,Symbol},1}:
 (:a, :b, :a)
 (:c, :a, :b)
 (:b, :c, :c)
 (:b, :a, :a)
 (:c, :b, :c)
 (:a, :a, :c)
 (:a, :c, :b)
 (:c, :c, :a)
 (:b, :b, :b)
```
"""
function get_codewords(𝒰::UniverseParameters, d::Int; m::Int=10, m_random::Int = 1000)
	code_size = 0
	C = eltype(𝒰)[]

	for _ in 1:m
		random_code = get_codewords_random(𝒰, d; m = m_random)
		random_size = length(random_code)
		if random_size > code_size
			code_size = random_size
			C = random_code
		end
	end

	greedy_code = get_codewords_greedy(𝒰, d)
	greedy_size = length(greedy_code)
	if greedy_size > code_size
		code_size = greedy_size
		C = greedy_code
	end

	return C
end
get_codewords(Σ::Alphabet{N}, q::Int, n::Int, d::Int; m::Int=10, m_random::Int=1000) where {N} =
	get_codewords(UniverseParameters(Σ, q, n), d, m=m, m_random=m_random)
get_codewords(Σ::Alphabet{N}, n::Int, d::Int; m::Int=10, m_random::Int=1000) where {N} =
	get_codewords(UniverseParameters(Σ, n), d, m=m, m_random=m_random)
get_codewords(q::Int, n::Int, d::Int; m::Int=10, m_random::Int=1000) =
	get_codewords(UniverseParameters(q, n), d, m=m, m_random=m_random)
get_codewords(Σ::AbstractArray, q::Int, n::Int, d::Int; m::Int=10, m_random::Int=1000) =
	get_codewords(Alphabet(Σ), q, n, d, m=m, m_random=m_random)
get_codewords(Σ::AbstractArray, n::Int, d::Int; m::Int=10, m_random::Int=1000) =
	get_codewords(UniverseParameters(Σ, n), d, m=m, m_random=m_random)
get_codewords(Σ::AbstractArray, q::Int, n::Int, d::Int, 𝒰::AbstractArray; m::Int=10, m_random::Int=1000) =
	get_codewords(Alphabet(Σ), q, n, d, 𝒰, m=m, m_random=m_random)

"""
```julia
get_codewords(G::AbstractArray, m::Int) -> Codewords{M}
```

Get codewords of a code from the _generating matrix_ under a finite field of modulo `m`.  Precisely, computes all linear combinations of the rows of the generating matrix.

Parameters:
  - `G::AbstractArray`: A matrix of Ints which generates the code.
  - `m::Int`: The bounds of the finite field (i.e., the molulus you wish to work in).

Returns:
  - `Codewords{M}`: An array of codewords, each of length `M`.  Each codewords is a tuple, and each character in said word is a symbol.
"""
function get_codewords(G::AbstractArray, m::Int)
	codewords = Vector()
	rows = Vector(undef, size(G, 2))

	for i in 1:size(G, 1)
		rows[i] = [G[i, j] for j in 1:size(G, 2)]
	end

	for c in Base.Iterators.product([0:m-1 for _ in 1:size(G, 1)]...)
		word = Ref(c[1]) .* rows[1]

		for i in 2:size(G, 1)
			word = mod.(word .+ (Ref(c[i]) .* rows[i]), m)
		end

		push!(codewords, word)
	end

	return codewords
end

# function obtain_maximal_code(𝒰::UniverseParameters, d::Int)
# 	adj_matrix = Matrix{Int8}(undef, length(𝒰), length(𝒰))
#
# 	for u in CodeUniverseIterator(𝒰), u′ in CodeUniverseIterator(𝒰)
# 		distance = hamming_distance(u, u′)
# 	end
# end


text:
push!(LOAD_PATH, "src")

include("minperfhash.jl")


text:
###############################################################################
#
# Initialize WaveLengths
#
###############################################################################
"""
    create_wave_length(
                FT,
                sWLs = [collect(400.0:10.0: 650.1);
                        collect(655.0: 5.0: 770.1);
                        collect(780.0:25.0:2400.1)])

Create [`WaveLengths`](@ref) type struct, given
- `FT` Floating number type
"""
function create_wave_length(
            FT,
            sWLs = [collect(400.0:10.0: 650.1);
                    collect(655.0: 5.0: 770.1);
                    collect(780.0:25.0:2400.1)]
)
    sWL   = FT.(sWLs);
    optis = create_leaf_opticals(sWL, OPTI_2021);

    return WaveLengths{FT}(sWL=sWL, optis=optis)
end

max_stars_repo_name: bbuman/Land

max_issues_repo_name: bbuman/Land

max_forks_repo_name: bbuman/Land


text:
using ParaDwell
using Test

@testset "ParaDwell.jl" begin
    # These need data - need light dummy data for CI
    ParaDwell.DIVA(adminLevel=0);
    ParaDwell.DIVA(adminLevel=1);
    ParaDwell.DIVA(adminLevel=2,exclShapeBelowLen=100);
end


text:
export density

function convert_arguments(P::PlotFunc, d::KernelDensity.UnivariateKDE)
    ptype = plottype(P, Lines) # choose the more concrete one
    to_plotspec(ptype, convert_arguments(ptype, d.x, d.density))
end

function convert_arguments(P::PlotFunc, d::KernelDensity.BivariateKDE)
    ptype = plottype(P, Heatmap)
    to_plotspec(ptype, convert_arguments(ptype, d.x, d.y, d.density))
end

function trim_limits!(x, lims)
    len = column_length(x)
    for i in 1:ncols(x)
        xmin, xmax = extrema_nan(extract_view(x, 1:len, i))
        xmin′, xmax′ = lims[i]
        lims[i] = (max(xmin′, xmin), min(xmax′, xmax))
    end
end

function search_both(x, xlims)
    min, max = xlims
    x1 = searchsortedfirst(x, min)
    x2 = searchsortedlast(x, max)
    (x1, x2)
end

function trim_density(k::UnivariateKDE, lims::AbstractVector{<:Tuple})
    x1, x2 = search_both(k.x, lims[1])
    UnivariateKDE(k.x[x1:x2], k.density[x1:x2])
end

function trim_density(k::BivariateKDE, lims::AbstractVector{<:Tuple})
    x1, x2 = search_both(k.x, lims[1])
    y1, y2 = search_both(k.y, lims[2])
    BivariateKDE(k.x[x1:x2], k.y[y1:y2], k.density[x1:x2, y1:y2])
end

function trim_density(k::KernelDensity.AbstractKDE, x; xlims = (-Inf, Inf), ylims = (-Inf, Inf), trim = false)
    lims = Tuple[extrema(xlims), extrema(ylims)]
    trim && trim_limits!(x, lims)
    trim_density(k, lims)
end

function _density(x; xlims = (-Inf, Inf), ylims = (-Inf, Inf), trim = false, kwargs...)
    k = kde(x; kwargs...)
    trim_density(k, x; xlims = xlims, ylims = ylims, trim = trim)
end

function _density(x, y; xlims = (-Inf, Inf), ylims = (-Inf, Inf), trim = false, kwargs...)
    k = kde((x, y); kwargs...)
    trim_density(k, (x, y); xlims = xlims, ylims = ylims, trim = trim)
end

const density = Analysis(_density)

max_stars_repo_name: piever/StatMakie

max_issues_repo_name: piever/StatMakie

max_forks_repo_name: piever/StatMakie


text:
"""
Handles WebSockets communication logic.
"""
module Channels

using WebSockets, JSON, Logger

const ClientId =                          Int
const ChannelId =                         String
const ChannelClient =                     Dict{Symbol, Union{WebSockets.WebSocket,Vector{ChannelId}}}
const ChannelClientsCollection =          Dict{ClientId,ChannelClient} # { ws.id => { :client => ws, :channels => ["foo", "bar", "baz"] } }
const ChannelSubscriptionsCollection =    Dict{ChannelId,Vector{ClientId}}  # { "foo" => ["4", "12"] }
const MessagePayload =                    Union{Void,Dict}

mutable struct ChannelMessage
  channel::ChannelId
  client::ClientId
  message::String
  payload::MessagePayload
end

const CLIENTS       = ChannelClientsCollection()
const SUBSCRIPTIONS = ChannelSubscriptionsCollection()

"""
    subscribe(ws::WebSockets.WebSocket, channel::ChannelId) :: Void

Subscribes a web socket client `ws` to `channel`.
"""
function subscribe(ws::WebSockets.WebSocket, channel::ChannelId) :: Void
  if haskey(CLIENTS, ws.id)
    ! in(channel, CLIENTS[ws.id][:channels]) && push!(CLIENTS[ws.id][:channels], channel)
  else
    CLIENTS[ws.id] = Dict(
                          :client => ws,
                          :channels => ChannelId[channel]
                          )
  end

  push_subscription(ws.id, channel)

  nothing
end

"""
    unsubscribe(ws::WebSockets.WebSocket, channel::ChannelId) :: Void

Unsubscribes a web socket client `ws` from `channel`.
"""
function unsubscribe(ws::WebSockets.WebSocket, channel::ChannelId) :: Void
  if haskey(CLIENTS, ws.id)
    delete!(CLIENTS[ws.id][:channels], channel)
  end

  pop_subscription(ws.id, channel)

  nothing
end

"""
    unsubscribe_client(ws::WebSockets.WebSocket) :: Void

Unsubscribes a web socket client `ws` from all the channels.
"""
function unsubscribe_client(client::ClientId) :: Void
  unsubscribe_client(Channels.CLIENTS[client][:client])
end
function unsubscribe_client(ws::WebSockets.WebSocket) :: Void
  if haskey(CLIENTS, ws.id)
    for channel_id in CLIENTS[ws.id][:channels]
      pop_subscription(ws.id, channel_id)
    end

    delete!(CLIENTS, ws.id)
  end

  nothing
end

"""
    push_subscription(client::ClientId, channel::ChannelId) :: Void

Adds a new subscription for `client` to `channel`.
"""
function push_subscription(client::ClientId, channel::ChannelId) :: Void
  if haskey(SUBSCRIPTIONS, channel)
    ! in(client, SUBSCRIPTIONS[channel]) && push!(SUBSCRIPTIONS[channel], client)
  else
    SUBSCRIPTIONS[channel] = ClientId[client]
  end

  nothing
end

"""
    pop_subscription(client::ClientId, channel::ChannelId) :: Void

Removes the subscription of `client` to `channel`.
"""
function pop_subscription(client::ClientId, channel::ChannelId) :: Void
  if haskey(SUBSCRIPTIONS, channel)
    filter!(SUBSCRIPTIONS[channel]) do (client_id)
      client_id != client
    end
    isempty(SUBSCRIPTIONS[channel]) && delete!(SUBSCRIPTIONS, channel)
  end

  nothing
end

"""
    pop_subscription(client::ClientId) :: Void

Removes all subscriptions of `client`.
"""
function pop_subscription(channel::ChannelId) :: Void
  if haskey(SUBSCRIPTIONS, channel)
    delete!(SUBSCRIPTIONS, channel)
  end

  nothing
end

"""
    broadcast(channels::Union{ChannelId,Vector{ChannelId}}, msg::String) :: Void
    broadcast{U,T}(channels::Vector{ChannelId}, msg::String, payload::Dict{U,T}) :: Void

Pushes `msg` (and `payload`) to all the clients subscribed to the channels in `channels`.
"""
function broadcast(channels::Union{ChannelId,Vector{ChannelId}}, msg::String) :: Void
  isa(channels, Array) || (channels = ChannelId[channels])

  @parallel for channel in channels
    for client in SUBSCRIPTIONS[channel]
      ws_write_message(client, msg)
    end
  end

  nothing
end
function broadcast{U,T}(channels::Union{ChannelId,Vector{ChannelId}}, msg::String, payload::Dict{U,T}) :: Void
  isa(channels, Array) || (channels = ChannelId[channels])

  @parallel for channel in channels
    for client in SUBSCRIPTIONS[channel]
      try
        ws_write_message(client, ChannelMessage(channel, client, msg, payload) |> JSON.json)
      catch ex
        Logger.log(string(ex), :err)
        Logger.log("$(@__FILE__):$(@__LINE__)")

        rethrow(ex)
      end
    end
  end

  nothing
end

"""
    broadcast(msg::String) :: Void
    broadcast{U,T}(msg::String, payload::Dict{U,T}) :: Void

Pushes `msg` (and `payload`) to all the clients subscribed to all the channels.
"""
function broadcast(msg::String) :: Void
  broadcast(collect(keys(SUBSCRIPTIONS)), msg)
end
function broadcast{U,T}(msg::String, payload::Dict{U,T}) :: Void
  broadcast(collect(keys(SUBSCRIPTIONS)), msg, payload)
end

"""
  message(channel::ChannelId, msg::String) :: Void
  message{U,T}(channel::ChannelId, msg::String, payload::Dict{U,T}) :: Void

Pushes `msg` (and `payload`) to `channel`.
"""
function message(channel::ChannelId, msg::String) :: Void
  broadcast(ChannelId[channel], msg)
end
function message{U,T}(channel::ChannelId, msg::String, payload::Dict{U,T}) :: Void
  broadcast(ChannelId[channel], msg, payload)
end

"""
    ws_write_message(client::ClientId, msg::String) :: Void

Writes `msg` to web socket for `client`.
"""
function ws_write_message(client::ClientId, msg::String) :: Void
  write(CLIENTS[client][:client], msg)

  nothing
end

"""
    message(client::ChannelClient, msg::String) :: Void

Send message `msg` to `client`.
"""
function message(client::ChannelClient, msg::String) :: Void
  write(client[:client], msg)

  nothing
end

end


text:
#######################################################################
# Kernel Polynomial Method : momenta
#######################################################################
struct MomentaKPM{T,B<:Tuple}
    μlist::Vector{T}
    bandbracket::B
end
struct KPMBuilder{A,H<:AbstractMatrix,T,K,B}
    A::A
    h::H
    bandbracket::Tuple{B,B}
    order::Int
    missingket::Bool
    μlist::Vector{T}
    ket::K
    ket0::K
    ket1::K
end

function KPMBuilder(h::AbstractMatrix{Tv}, A = _defaultA(Tv); ket = missing, order = 10, bandrange = missing) where {Tv}
    eh = eltype(eltype(h))
    aA = eltype(eltype(A))
    μlist = zeros(promote_type(eh, aA), order + 1)
    bandbracket = bandbracketKPM(h, bandrange)
    missingket = ket === missing
    ket´ = missingket ? ketundef(h) : ket
    iscompatibleket(h, ket´) || throw(ArgumentError("ket is incompatible with Hamiltonian"))
    builder = KPMBuilder(A, h, bandbracket, order, missingket,
        μlist, ket´, similar(ket´), similar(ket´))
    return builder
end

ketundef(h::AbstractMatrix{T}) where {T<:Number} =
    Vector{T}(undef, size(h, 2))
ketundef(h::AbstractMatrix{S}) where {N,T,S<:SMatrix{N,N,T}} =
    Vector{SVector{N,T}}(undef, size(h, 2))

iscompatibleket(h::AbstractMatrix{T1}, ket::AbstractArray{T2}) where {T1,T2} =
    _iscompatibleket(T1, T2)
_iscompatibleket(::Type{T1}, ::Type{T2}) where {T1<:Real, T2<:Real} = true
_iscompatibleket(::Type{T1}, ::Type{T2}) where {T1<:Number, T2<:Complex} = true
_iscompatibleket(::Type{S1}, ::Type{S2}) where {N, S1<:SMatrix{N,N}, S2<:SVector{N}} =
    _iscompatibleket(eltype(S1), eltype(S2))
_iscompatibleket(::Type{S1}, ::Type{S2}) where {N, S1<:SMatrix{N,N}, S2<:SMatrix{N}} =
    _iscompatibleket(eltype(S1), eltype(S2))
_iscompatibleket(t1, t2) = false

function matrixKPM(h::Hamiltonian{<:Lattice,L}) where {L}
    iszero(L) ||
        throw(ArgumentError("Hamiltonian is defined on an infinite lattice. Convert it to a matrix first using `bloch(h, φs...)`"))
    return bloch(h)
end

matrixKPM(h::AbstractMatrix) = h
matrixKPM(A::UniformScaling) = A

"""
    momentaKPM(h::AbstractMatrix, A = I; ket = missing, order = 10, randomkets = 1, bandrange = missing)

Compute the Kernel Polynomial Method (KPM) momenta `μ_n = ⟨ket|T_n(h) A|ket⟩/⟨ket|ket⟩` where `T_n(x)`
is the Chebyshev polynomial of order `n`, for a given `ket::AbstractVector`, hamiltonian `h`, and
observable `A`. If `ket` is `missing`, momenta are computed by means of a stochastic trace
`μ_n = Tr[A T_n(h)] ≈ ∑ₐ⟨a|A T_n(h)|a⟩/N` over `N = randomkets` normalized random `|a⟩`.
Furthermore, the trace over a specific set of kets can also be computed; in this case
`ket::AbstractMatrix` must be a matrix where the columns are the kets involved in the calculation.

The order of the Chebyshev expansion is `order`. The `bandbrange = (ϵmin, ϵmax)` should completely encompass
the full bandwidth of `hamiltonian`. If `missing` it is computed automatically using `ArnoldiMethods` (must be loaded).

# Example
```
julia> h = LatticePresets.cubic() |> hamiltonian(hopping(1)) |> unitcell(region = RegionPresets.sphere(10));

julia> momentaKPM(bloch(h), bandrange = (-6,6))
Elsa.MomentaKPM{Float64}([0.9594929736144973, -0.005881595972403821, -0.4933354572913581, 0.00359537502632597, 0.09759451291347333, -0.0008081453185250322, -0.00896262538765363, 0.00048205637037715177, -0.0003705198310034668, 9.64901673962623e-20, 9.110915988898614e-18], (0.0, 6.030150753768845))
```
"""
momentaKPM(h::Hamiltonian, A = _defaultA(eltype(h)); kw...) = momentaKPM(matrixKPM(h), matrixKPM(A); kw...)

function momentaKPM(h::AbstractMatrix, A = _defaultA(eltype(h)); randomkets = 1, kw...)
   b = KPMBuilder(h, A; kw...)
    if b.missingket
        pmeter = Progress(b.order * randomkets, "Averaging moments: ")
        for n in 1:randomkets
            randomize!(b.ket)
            addmomentaKPM!(b, pmeter)
        end
        b.μlist ./= randomkets
    else
        pmeter = Progress(b.order, "Computing moments: ")
        addmomentaKPM!(b, pmeter)
    end
    return MomentaKPM(jackson!(b.μlist), b.bandbracket)
end

_defaultA(::Type{T}) where {T<:Number} = one(T) * I
_defaultA(::Type{S}) where {N,T,S<:SMatrix{N,N,T}} = one(T) * I

# This iterates bras <psi_n| = <psi_0|AT_n(h) instead of kets (faster CSC multiplication)
# In practice we iterate their conjugate |psi_n> = T_n(h') A'|psi_0>, and do the projection
# onto the start ket, |psi_0>
function addmomentaKPM!(b::KPMBuilder{<:AbstractMatrix,<:AbstractSparseMatrix}, pmeter)
    μlist, ket, ket0, ket1 = b.μlist, b.ket, b.ket0, b.ket1
    h´, A´, bandbracket = b.h', b.A', b.bandbracket
    order = length(μlist) - 1
    fill!(ket0, zero(eltype(ket0)))
    mul!(ket1, A´, ket)
    μlist[1] += proj(ket1, ket)
    ProgressMeter.next!(pmeter; showvalues = ())
    for n in 2:(order+1)
        ProgressMeter.next!(pmeter; showvalues = ())
        μ = iterateKPM!(ket0, ket1, ket, h´, bandbracket)
        μlist[n] += μ
        n + 1 > order + 1 && break
        ket0, ket1 = ket1, ket0
    end
    return μlist
end

function iterateKPM!(ket0::A, ket1::A, kini::A, adjh::Adjoint, (center, halfwidth)) where {S,A<:AbstractArray{S}}
    h = adjh.parent
    nzv = nonzeros(h)
    rv = rowvals(h)
    T = eltype(S)
    μ = zero(T)
    α = T(-2 * center / halfwidth)
    β = T(2 / halfwidth)
    for k in 1:size(ket0, 2)
        for col in 1:size(h, 2)
            @inbounds begin
                tmp = α * ket1[col, k] - ket0[col, k]
                for ptr in nzrange(h, col)
                    tmp += β * adjoint(nzv[ptr]) * ket1[rv[ptr],k]
                end
                ket0[col, k] = tmp
                μ += dot(tmp, kini[col, k])
            end
        end
    end
    return μ
end

function addmomentaKPM!(b::KPMBuilder{<:UniformScaling, <:AbstractSparseMatrix}, pmeter)
    μlist, ket, ket0, ket1 = b.μlist, b.ket, b.ket0, b.ket1
    h´, A, bandbracket = b.h', b.A, b.bandbracket
    order = length(μlist) - 1
    fill!(ket0, zero(eltype(ket0)))
    ket1 .= ket
    for n in 1:2:(order+1)
        ProgressMeter.next!(pmeter; showvalues = ())
        ProgressMeter.next!(pmeter; showvalues = ()) # twice because of 2-step
        μ, μ´ = iterateKPM!(ket0, ket1, h´, bandbracket)
        μlist[n] += μ
        n + 1 > order + 1 && break
        μlist[n + 1] += μ´
        ket0, ket1 = ket1, ket0
    end
    A.λ ≈ 1 || (μlist .*= A.λ)
    return μlist
end

function iterateKPM!(ket0::A, ket1::A, adjh::Adjoint, (center, halfwidth)) where {S,A<:AbstractArray{S}}
    h = adjh.parent
    nzv = nonzeros(h)
    rv = rowvals(h)
    T = eltype(S)
    μ = μ´ = zero(T)
    α = T(-2 * center / halfwidth)
    β = T(2 / halfwidth)
    for k in 1:size(ket0, 2)
        for col in 1:size(h, 2)
            @inbounds begin
                k1 = ket1[col, k]
                tmp = α * k1 - ket0[col, k]
                for ptr in nzrange(h, col)
                    tmp += β * adjoint(nzv[ptr]) * ket1[rv[ptr],k]
                end
                ket0[col, k] = tmp
                μ  += dot(k1, k1)
                μ´ += dot(tmp, k1)
            end
        end
    end
    return μ, μ´
end

# This is equivalent to tr(ket1'*ket2) for matrices, and ket1'*ket2 for vectors
proj(ket1, ket2) = dot(vec(ket1), vec(ket2))

function randomize!(v::AbstractVector{T}) where {T}
    v .= _randomize.(v)
    normalize!(v)
    return v
end
@inline _randomize(v::T) where {T<:Real} = randn(R)
@inline _randomize(v::T) where {R,T<:Complex{R}} = randn(R) + im * randn(R)
@inline _randomize(v::T) where {T<:SArray} = _randomize.(v)

function jackson!(μ::AbstractVector)
    order = length(μ) - 1
    @inbounds for n in eachindex(μ)
        μ[n] *= ((order - n + 1) * cos(π * n / (order + 1)) +
                sin(π * n / (order + 1)) * cot(π / (order + 1))) / (order + 1)
    end
    return μ
end

function bandbracketKPM(h, ::Missing)
    @warn "Computing spectrum bounds... Consider using the `bandrange` kwargs for faster performance."
    bandbracketKPM(h, bandrangeKPM(h))
end
bandbracketKPM(h, (ϵmin, ϵmax)::Tuple{T,T}, pad = float(T)(0.01)) where {T} = ((ϵmax + ϵmin) / 2, (ϵmax - ϵmin) / (2 - pad))

function bandrangeKPM(h::AbstractMatrix{T}) where {T}
    checkloaded(:ArnoldiMethod)
    R = real(T)
    decompl, _ = Main.ArnoldiMethod.partialschur(h, nev=1, tol=1e-4, which = Main.ArnoldiMethod.LR());
    decomps, _ = Main.ArnoldiMethod.partialschur(h, nev=1, tol=1e-4, which = Main.ArnoldiMethod.SR());
    ϵmax = R(real(decompl.eigenvalues[1]))
    ϵmin = R(real(decomps.eigenvalues[1]))
    @warn  "Computed bandrange = ($ϵmin, $ϵmax)"
    return (ϵmin, ϵmax)
end

#######################################################################
# Kernel Polynomial Method : observables
#######################################################################
"""
    dosKPM(h::AbstractMatrix; resolution = 2, kw...)

Compute, using the Kernel Polynomial Method (KPM), the local density of states `ρ(ϵ) =
⟨ket|δ(ϵ-h)|ket⟩/⟨ket|ket⟩` for a given `ket::AbstractVector` and hamiltonian `h`, or the
global density of states `ρ(ϵ) = Tr[δ(ϵ-h)]` if `ket` is `missing`.

If `ket` is an `AbstractMatrix` it evaluates the trace over the set of kets in `ket` (see
`momentaKPM` and its options `kw` for further details). The result is a tuple of energy
points `xk::Vector` and real `ρ::Vector` values (any imaginary part in ρ is dropped), where
the number of energy points `xk` is `order * resolution`, rounded to the closest integer.

    dosKPM(μ::MomentaKPM; resolution = 2)

Same as above with momenta `μ` as input.

    dosKPM(h::Hamiltonian; kw...)

Equivalent to `dosKPM(bloch(h); kw...)` for finite hamiltonians (zero dimensional).
"""
dosKPM(h::AbstractMatrix; resolution = 2, kw...) =
    dosKPM(momentaKPM(h; kw...), resolution = resolution)

dosKPM(μ::MomentaKPM; resolution = 2) = real.(densityKPM(μ; resolution = resolution))

dosKPM(h::Hamiltonian; kw...) = dosKPM(matrixKPM(h); kw...)

"""
    densityKPM(h::AbstractMatrix, A; resolution = 2, kw...)

Compute, using the Kernel Polynomial Method (KPM), the local spectral density of an operator
`A` `ρ_A(ϵ) = ⟨ket|A δ(ϵ-h)|ket⟩/⟨ket|ket⟩` for a given `ket::AbstractVector` and
hamiltonian `h`, or the global spectral density `ρ_A(ϵ) = Tr[A δ(ϵ-h)]` if `ket` is
`missing`. If `ket` is an `AbstractMatrix` it evaluates the trace over the set of kets in
`ket` (see `momentaKPM` and its options `kw` for further details). A tuple of energy points
`xk` and `ρ_A` values is returned where the number of energy points `xk` is `order *
resolution`, rounded to the closest integer.

    densityKPM(momenta::MomentaKPM; resolution = 2)

Same as above with the KPM momenta as input (see `momentaKPM`).

    densityKPM(h::Hamiltonian, A::Hamiltonian; kw...)

Equivalent to `densityKPM(bloch(h), bloch(A); kw...)` for finite Hamiltonians (zero dimensional).
"""
densityKPM(h::AbstractMatrix, A; resolution = 2, kw...) =
    densityKPM(momentaKPM(h, A; kw...); resolution = resolution)

densityKPM(h::Hamiltonian, A::Hamiltonian; kw...) = densityKPM(matrixKPM(h), matrixKPM(A); kw...)

function densityKPM(momenta::MomentaKPM{T}; resolution = 2) where {T}
    checkloaded(:FFTW)
    (center, halfwidth) = momenta.bandbracket
    numpoints = round(Int, length(momenta.μlist) * resolution)
    ρlist = zeros(T, numpoints)
    copyto!(ρlist, momenta.μlist)
    Main.FFTW.r2r!(ρlist, Main.FFTW.REDFT01, 1)  # DCT-III in FFTW
    xk = [cos(π * (k + 0.5) / numpoints) for k in 0:numpoints - 1]
    @. ρlist = center + halfwidth * ρlist / (π * sqrt(1.0 - xk^2))
    @. xk = center + halfwidth * xk
    return xk, ρlist
end

"""
    averageKPM(h::AbstractMatrix, A; kBT = 0, Ef = 0, kw...)

Compute, using the Kernel Polynomial Method (KPM), the thermal expectation value `<A> = Σ_k
f(E_k) <k|A|k> =  ∫dE f(E) Tr [A δ(E-H)] = Tr [A f(H)]` for a given hermitian operator `A`
and a hamiltonian `h` (see `momentaKPM` and its options `kw` for further details).
`f(E)` is the Fermi-Dirac distribution function, `kBT` is the temperature in energy
units and `Ef` the Fermi energy.

    averageKPM(μ::MomentaKPM, A; kBT = 0, Ef = 0)

Same as above with the KPM momenta as input (see `momentaKPM`).

    averageKPM(h::Hamiltonian, A::Hamiltonian; kw...)

Equivalent to `averageKPM(bloch(h), bloch(A); kw...)` for finite Hamiltonians (zero
dimensional).
"""
averageKPM(h::AbstractMatrix, A; kBT = 0, Ef = 0, kw...) = averageKPM(momentaKPM(h, A; kw...); kBT = kBT, Ef = Ef)

averageKPM(h::Hamiltonian, A::Hamiltonian; kw...) = averageKPM(matrixKPM(h), matrixKPM(A); kw...)

function averageKPM(momenta::MomentaKPM{T}; kBT = 0.0, Ef = 0.0) where {T}
    (center, halfwidth) = momenta.bandbracket
    order = length(momenta.μlist) - 1
    # if !iszero(kBT)
    #     @warn "Finite temperature requires numerical evaluation of the integrals"
    #     checkloaded(:QuadGK)
    # end
    average = sum(n -> momenta.μlist[n + 1] * fermicheby(n, Ef, kBT, center, halfwidth), 0:order)
    return average
end

# Pending issue: Unexpected behaviour with center != 0.
function fermicheby(n, Ef, kBT, center, halfwidth)
    kBT´ = kBT / halfwidth;
    Ef´ = (Ef - center) / halfwidth;
    T = typeof(Ef´)
    if kBT´ == 0
        int = n == 0 ? 0.5+asin(Ef´)/π : -2.0*sin(n*acos(Ef´))/(n*π)
    else
        throw(error("Finite temperature not yet implemented"))
        # η = 1e-10
        # int = Main.QuadGK.quadgk(E´ -> _intfermi(n, E´, Ef´, kBT´), -1.0+η, 1.0-η, atol= 1e-10, rtol=1e-10)[1]
    end
    return T(int)
end

_intfermi(n, E´, Ef´, kBT´) = fermifun(E´, Ef´, kBT´) * 2/(π*(1-E´^2)^(1/2)) * chebypol(n, E´) / (1+(n == 0 ? 1 : 0))

fermifun(E´, Ef´, kBT´) = kBT´ == 0 ? (E´<Ef´ ? 1.0 : 0.0) : (1/(1+exp((E´-Ef´)/(kBT´))))

function chebypol(m::Int, x::T) where {T<:Number}
    cheby0 = one(T)
    cheby1 = x
    if m == 0
        chebym = cheby0
    elseif m == 1
        chebym = cheby1
    else
        for i in 2:m
            chebym = 2x * cheby1 - cheby0
            cheby0, cheby1 = cheby1, chebym
        end
    end
    return chebym
end

max_issues_repo_name: pablosanjose/QBox

max_forks_repo_name: pablosanjose/QBox


text:
const davidon = let res_init=zeros(15), jac_init=zeros(15,15), x_init=fill(0.1, 15)

    function res(x, r)
        for i = 1:14
            r1 = 0.0
            ti = Float64(i)
            t = 1.0
            for j = 1:15
                r1 += t * x[j]
                t *= ti
            end
            r[i] = r1
        end
        r[15] = x[1] - 1.0
        return r
    end

    function jac(x, jac)
        for i = 1:14
            ti = Float64(i)
            t = 1.0
            for k = 1:15
                jac[i, k] = t
                t *= ti
            end
        end
        jac[15, 1] = 1.0
        for k = 2:15
            jac[15, k] = 0.0
        end
        return jac
    end

    f(;scale=1, verbose=false, print_steps=false) = testone("davidon", res, jac,
                                                   res_init, jac_init, x_init;
                                                   scale=scale, verbose=verbose,
                                                   print_steps=print_steps)
end


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

# Test getindex for skinny sparse matrix
function sparse_getindex_skinny_perf()
    seed = 1
    srand(seed)

    # matsize = (size(A,1), nnz(A))
    matsize = [
        (2^12,1  ),
        (2^12,2^6),
        (2^18,1  ),
        (2^18,2^9)]
    # indsize = (length(I), number of repetitions)
    indsize = [
        (1,   2^12),
        (2^8, 2^8 ),
        (2^16,2^4 )]

    c = 0 # counter
    for (m,nz) in matsize
        A = sprand(m,1,nz/m)
        for (n,p) in indsize
            c += 1
            I = rand(1:m,n)
            @timeit indexing(A,I,p) "sparse_getindex_skinny$c" ""
        end
    end
end

function indexing(A,I,p)
    J = [1]
    for k = 1:p
        A[I,J]
    end
    nothing
end

max_stars_repo_name: JuliaTagBot/Benchmarks

max_issues_repo_name: JuliaArchive/Benchmarks

max_forks_repo_name: JuliaArchive/Benchmarks


text:
using DrWatson
@quickactivate :Swalbe
using Plots, CUDA, CSV, DataFrames

println("Starting iteration through surface tensions and sliplengths")
df_area = DataFrame()
df_height = DataFrame()
for γ1 in [0.005]
    for δ1 in [1.0 0.5 0.25]
        println("Simulating with γ: $γ1 and δ: $δ1")
        itname = Symbol("ga_$(γ1)_sl_$(δ1)")
        sys = Swalbe.SysConst(Lx=300, Ly=300, γ=γ1, δ=δ1, n=3, m=2, hmin=0.07, Tmax=100000, tdump=1000000)
        h, diff = Swalbe.run_dropletrelax(sys, "GPU", radius=85)
        # rads = sqrt.(diff ./ π)
        # p = plot(rads/rads[1], xaxis=:log, xlabel="time", ylabel="R/R₀")
        # savefig(p, "../images/rads_$(γ1)_$(δ1).pdf")
        # pp = heatmap(Array(h), aspect_ratio=1, c=:viridis)
        # savefig(pp, "../images/shape_$(γ1)_$(δ1).png")
        df_area[!, itname] = diff
        df_height[!, itname] = vec(Array(h))
        # @tagsave(datadir("Relax_droplet", savename(itname, "bson")),
        #          @dict diff vec(h) itname)
        CUDA.reclaim()
    end
end

CSV.write("data/Relax_droplet/Area_evolution_smaller_g.csv", df_area)
CSV.write("data/Relax_droplet/Height_final_smaller_g.csv", df_height)

# ;gzip data/Relax_droplet/Area_evolution.csv
# ;gzip data/Relax_droplet/Height_final.csv


text:
using TimerOutputs

@params struct LinearElasticityResult
    comp
    u
end
Base.show(::IO, ::MIME{Symbol("text/plain")}, ::LinearElasticityResult) = println("TopOpt linear elasticity result")

function simulate(problem::StiffnessTopOptProblem, topology = ones(getncells(TopOptProblems.getdh(problem).grid)); round = true, hard = true, xmin = 0.001)
    if round
        if hard
            solver = FEASolver(Displacement, Direct, problem, xmin = 0.0)
        else
            solver = FEASolver(Displacement, Direct, problem, xmin = xmin)
        end
    else
        solver = FEASolver(Displacement, Direct, problem, xmin = xmin)
    end
    vars = solver.vars
    fill_vars!(vars, topology; round = round)
    solver(Val{true})
    comp = dot(solver.u, solver.globalinfo.f)
    return LinearElasticityResult(comp, copy(solver.u))
end

function fill_vars!(vars::Array, topology; round)
    if round
        vars .= Base.round.(topology)
    else
        copyto!(vars, topology)
    end
    return vars
end


text:
using JaCoP
using MathOptInterface
using ConstraintProgrammingExtensions

using Test

const MOI = MathOptInterface
const MOIT = MOI.Test
const MOIB = MOI.Bridges
const CP = ConstraintProgrammingExtensions
const COIT = CP.Test

@testset "JaCoP" begin
    include("MOI.jl")
end


text:
# AUTO GENERATED FILE - DO NOT EDIT

export dashiconpickeritem

"""
    dashiconpickeritem(;kwargs...)

A DashIconPickerItem component.

Keyword arguments:
- `id` (String; required): The ID used to identify this component in Dash callbacks.
- `icon` (String; required): The name of the icon to render. Example: "FaAdobe"
- `size` (Real; optional): The size of the icon
- `color` (String; optional): The color of the icon. Default: "#000"
- `className` (String; optional): CSS class for div
"""
function dashiconpickeritem(; kwargs...)
        available_props = Symbol[:id, :icon, :size, :color, :className]
        wild_props = Symbol[]
        return Component("dashiconpickeritem", "DashIconPickerItem", "dash_react_font_awesome_icon_picker", available_props, wild_props; kwargs...)
end


text:
using Pkg

PKGS = [
        "OhMyREPL",
        "CSV",
        "DataFrames",
        "Missings",
        "Parsers",
        "RData",
     ]

for pkg in PKGS
  if pkg ∉ keys(Pkg.installed())
    Pkg.add(pkg)
  end
end

using OhMyREPL

max_stars_repo_name: pedroagribeiro/dotfiles

max_issues_repo_name: pedroagribeiro/dotfiles

max_forks_repo_name: pedroagribeiro/dotfiles


text:
module Sym
using DiffRules, SpecialFunctions, NaNMath

include("engine/engine.jl")
include("types/types.jl")

end # module


text:
# Autogenerated wrapper script for MKL_jll for x86_64-apple-darwin
export libmkl_core, libmkl_rt

using IntelOpenMP_jll
JLLWrappers.@generate_wrapper_header("MKL")
JLLWrappers.@declare_library_product(libmkl_core, "@rpath/libmkl_core.1.dylib")
JLLWrappers.@declare_library_product(libmkl_rt, "@rpath/libmkl_rt.1.dylib")
function __init__()
    JLLWrappers.@generate_init_header(IntelOpenMP_jll)
    JLLWrappers.@init_library_product(
        libmkl_core,
        "lib/libmkl_core.1.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmkl_rt,
        "lib/libmkl_rt.1.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# ---------------------------------------------------------------------------- #
#
#   genSTLslice.jl
#
#   Generate STL slices directly from lattice frame representation
#
#   intrico
#   Summer 2018
#
#   Max Opgenoord
#
# ---------------------------------------------------------------------------- #

function genSTLslice( mesh::MeshF,  lat::Lattice, flname::String,
                      nprinter::SVector{3,Float64}, z::Float64,
                      actptsNods::Vector{Vector{SVector{3,Float64}}},
                      edgesNods::Vector{Vector{SVector{5,Int64}}};
                      n::Int64 = 8 )

    # z is defined from origin

    # open STL file
    fid = open( flname, "w" )

    # write header
    for jj in 1:80
        write(fid,' ')
    end
    write(fid,UInt32(0)) # number of faces is unknown now, so write 0

    # z location in part coordinates
    zpart = z * nprinter

    # find intersecting edges
    edgesInt = findEdgesInt( mesh, lat, nprinter, zpart )

    # generate the slices
    factot = [ 0 ]
    writeNodStrutSlice( mesh, lat, nprinter, zpart, z,
                        edgesInt, n, fid, factot,
                        actptsNods, edgesNods )

    # close STL
    close( fid )

    ## reopen to write the number of faces
    fid = open( flname, "a+" )
    seekstart(fid) # go back to top of file

    # rewrite header
    for jj in 1:80
        write(fid,' ')
    end

    # write correct number of faces
    write( fid, UInt32( factot[1] ) )

    # close STL
    close( fid )

end

function findEdgesInt( mesh::MeshF, lat::Lattice,
                       nprinter::SVector{3,Float64}, zpart::SVector{3,Float64} )

    ne = size(mesh.e,1)

    edgesInt = fill( false, ne )

    # loop over edges
    for ee in 1:ne

        if lat.ar[ee] < 1e-10
            continue
        end

        r = sqrt( lat.ar[ee] / π )

        inod1 = mesh.e[ee,1]
        inod2 = mesh.e[ee,2]
        nod1  = SVector{3,Float64}( mesh.p[inod1,:] )
        nod2  = SVector{3,Float64}( mesh.p[inod2,:] )
        ledge = norm( nod2 - nod1 )
        nvec  = SVector{3,Float64}( (nod2 - nod1) ./ ledge )

        # check if edge is intersecting or even remotely close to intersecting
        xint, d = findPlaneLineIntersec( zpart, nod1, nprinter, nvec )
        if d >= -2*r && d <= (ledge+2*r)
            # may be intersecting
            edgesInt[ee] = true
        end

    end

    return edgesInt

end

function writeNodStrutSlice( mesh::MeshF, lat::Lattice,
                             nprinter::SVector{3,Float64},
                             zpart::SVector{3,Float64}, z::Float64,
                             edgesInt::Vector{Bool},
                             n::Int64, fid::IOStream, factot::Vector{Int64},
                             actptsNods::Vector{Vector{SVector{3,Float64}}},
                             edgesNods::Vector{Vector{SVector{5,Int64}}} )

    zs = mesh.p * nprinter
    minz = minimum( zs )
    maxz = maximum( zs )

    φ = linspace( 0, 2*π, n+1 )
    xpts = Vector{ SVector{3,Float64} }( n+1 )

    edgWritten = fill( false, size(mesh.e,1) )

    # loop over nodes
    for nn in 1:mesh.n

        if !any(edgesInt[ mesh.n2e[nn] ])
            # no intersections
            continue
        end

        nedg = length(mesh.n2e[nn])

        # check if cutting through node
        intersecEdge = fill( false, nedg )
        intersecNodEdge = fill( false, length(edgesNods[nn]) )
        for kk in 1:length(edgesNods[nn])
            nod1 = mesh.p[nn,:] + actptsNods[nn][ edgesNods[nn][kk][1] ]
            nod2 = mesh.p[nn,:] + actptsNods[nn][ edgesNods[nn][kk][2] ]
            if (dot(nod1,nprinter) - z) * ( dot(nod2,nprinter) - z ) < 0.

                intersecNodEdge[kk] = true

                intersecEdge[ mesh.n2e[nn] .== edgesNods[nn][kk][3] ] = true
                intersecEdge[ mesh.n2e[nn] .== edgesNods[nn][kk][4] ] = true
            end
        end

        # write complicated node geometry
        writeNodSlice( mesh, lat, nn, edgesNods[nn], actptsNods[nn],
                       intersecEdge, intersecNodEdge, nprinter, zpart, z, n, fid, factot )

       if any(intersecEdge)
           # intersecting the node
           continue
       end

        # compute dist from node -- use pairings from edgesNods
        distNode = fill( 0., nedg )
        radmax = sqrt( maximum(lat.ar[ mesh.n2e[nn] ]) / π )
        for kk in 1:length(edgesNods[nn])

            inde1 = edgesNods[nn][kk][3]
            inde2 = edgesNods[nn][kk][4]
            inod1 = mesh.e[inde1,1] + mesh.e[inde1,2] - nn
            inod2 = mesh.e[inde2,1] + mesh.e[inde2,2] - nn
            vec1 = mesh.p[inod1,:] - mesh.p[nn,:]
            vec2 = mesh.p[inod2,:] - mesh.p[nn,:]

            arg = dot(vec1,vec2)/( norm(vec1)*norm(vec2) )
            arg = minimum( ( 1.0, arg) ) # guard for round-off errors
            arg = maximum( (-1.0, arg) ) # guard for round-off errors
            θ = acos( arg )

            # compute absolute distance between node center and intersection
            lv = radmax / tan( θ/2 ) # offset such that rods do not collide
            # l = sqrt( lv^2 + radmax^2 ) # absolute distance

            iedgloc1 = find(mesh.n2e[nn] .== edgesNods[nn][kk][3])[]
            iedgloc2 = find(mesh.n2e[nn] .== edgesNods[nn][kk][4])[]
            distNode[ iedgloc1 ] = max( distNode[ iedgloc1 ], lv )
            distNode[ iedgloc2 ] = max( distNode[ iedgloc2 ], lv )

        end

        # write slices for struts that do not collide with other struts
        for eloc in 1:nedg
            ee = mesh.n2e[nn][eloc]

            if intersecEdge[eloc] || edgWritten[ee]
                continue
            end

            rad = sqrt( lat.ar[ee] / π )

            inod1 = mesh.e[ee,1]
            inod2 = mesh.e[ee,2]
            nod1  = SVector{3,Float64}( mesh.p[nn,:] )
            nod2  = SVector{3,Float64}( mesh.p[inod1+inod2-nn,:] )
            ledge = norm( nod2 - nod1 )
            nvec  = SVector{3,Float64}( (nod2 - nod1) ./ ledge )

            nperp = cross( nvec, nprinter )
            nperp = nperp/norm(nperp)

            nzero = SVector{3}( nprinter[3], nprinter[1], nprinter[2] )
            if (abs( abs(dot(nprinter,nvec)) - 1.) ) < 1e-10
                nzero = cross(nzero,nprinter)
                nperp = cross(nprinter,nzero)
                nperp = nperp/norm(nperp)
            else
                nzero = cross(nvec,nperp)
            end
            nzero = nzero/norm(nzero)

            d = fill(0., 3)
            intsec0, d[1] = findPlaneLineIntersec( zpart, nod1, nprinter, nvec )

            # first line
            currpt = nod1 + nzero*rad
            intsec1, d[2] = findPlaneLineIntersec( zpart, currpt, nprinter, nvec )

            # second line
            currpt = nod1 - nzero*rad
            intsec2, d[3] = findPlaneLineIntersec( zpart, currpt, nprinter, nvec )

            if all( (d .< distNode[eloc]) .| (d .> 0.5*ledge) )
                # No intersections with first half of edge
                continue
            end

            edgWritten[ee] = true

            # generate points for STL
            vec1 = - rad * nperp
            vec2 = intsec1 - intsec0
            for jj in 1:n+1
                xpts[jj] = intsec0 + vec1 * cos(φ[jj]) + vec2 * sin(φ[jj])
            end

            writeSTLcontour( xpts, intsec0, nprinter, n, fid )
            factot[1] += n

        end

    end

end

function findPlaneLineIntersec( planept::SVector{3,Float64},
                                linept::SVector{3,Float64},
                                nPlane::SVector{3,Float64},
                                nLine::SVector{3,Float64} )

    d    = dot( planept-linept, nPlane ) / dot( nLine, nPlane )
    xloc = linept + d * nLine

    return xloc, d

end

function writeSTLcontour( xpts::Vector{SVector{3,Float64}},
                          intersec0::SVector{3,Float64},
                          nprinter::SVector{3,Float64},
                          n::Int64, fid::IOStream )

      for jj in 1:n

          xface = SVector{3,SVector{3,Float64}}( xpts[jj], xpts[jj+1], intersec0 )

          sterno.writeFacetSTLB( xface, nprinter, fid )

      end

end

function writeNodSlice( mesh::MeshF, lat::Lattice, nn::Int64,
                        edgesNods::Vector{SVector{5,Int64}},
                        actptsNods::Vector{SVector{3,Float64}},
                        intersecEdge::Vector{Bool}, intersecNodEdge::Vector{Bool},
                        nprinter::SVector{3,Float64},
                        zpart::SVector{3,Float64}, z::Float64,
                        n::Int64, fid::IOStream, factot::Vector{Int64} )

    nodvec = SVector{3}( mesh.p[nn,:] )

    Δn = ( z - dot(nodvec,nprinter) ) * nprinter

    Δnodvec = nodvec + Δn

    nint = count( intersecEdge )

    xpts   = Vector{ SVector{3,Float64} }( n+1 )
    xptscl = Vector{ SVector{3,Float64} }( 3 )
    xptsfl = Vector{ SVector{3,Float64} }( 4 )

    # loop over intersecting edges
    for kk in 1:length( mesh.n2e[nn] )

        e0 = mesh.n2e[nn][kk]

        if !intersecEdge[kk]
            continue
        end

        # find the other two neighboring edges
        e1 = 0
        e2 = 0
        enew = 0
        ihalf1 = 0
        ihalf2 = 0
        ct = 1
        for el in edgesNods
            if !intersecNodEdge[ct]
                ct += 1
                continue
            end
            if el[3] == e0
                enew = el[4]
            elseif el[4] == e0
                enew = el[3]
            else
                ct += 1
                continue
            end
            if e1 > 0
                e2 = enew
                ihalf2  = el[5]
            else
                e1 = enew
                ihalf1  = el[5]
            end
            ct += 1
        end

        # find which edge is active on which side by comparing ncut vectors
        nod0 = mesh.e[ e0, 1 ] + mesh.e[ e0, 2 ] - nn
        nod1 = mesh.e[ e1, 1 ] + mesh.e[ e1, 2 ] - nn
        nod2 = mesh.e[ e2, 1 ] + mesh.e[ e2, 2 ] - nn

        nvec0 = SVector{3}( mesh.p[nod0,:] - mesh.p[nn,:]); nvec0 = nvec0/norm(nvec0)
        nvec1 = SVector{3}( mesh.p[nod1,:] - mesh.p[nn,:]); nvec1 = nvec1/norm(nvec1)
        nvec2 = SVector{3}( mesh.p[nod2,:] - mesh.p[nn,:]); nvec2 = nvec2/norm(nvec2)

        ncut1 = nvec0 - nvec1; ncut1 = ncut1/norm(ncut1)
        ncut2 = nvec0 - nvec2; ncut2 = ncut2/norm(ncut2)
        nline1 = cross(ncut1,nprinter); nline1 = nline1/norm(nline1)
        nline2 = cross(ncut2,nprinter); nline2 = nline2/norm(nline2)

        # find correct orientation to make sure that nline is aligned with the halfvector of the edge (as projected from the edge)
        if dot( nline1 - dot(nline1,nvec1)*nvec1, actptsNods[ihalf1] ) < 0.
            nline1 = -nline1
        end
        if dot( nline2 - dot(nline2,nvec2)*nvec2, actptsNods[ihalf2] ) < 0.
            nline2 = -nline2
        end

        rad  = sqrt( lat.ar[ e0 ] / π )
        arg1 = dot(nvec0,nvec1)
        # arg = minimum( ( 1.0, arg) ) # guard for round-off errors
        # arg = maximum( (-1.0, arg) ) # guard for round-off errors
        θ1  = acos( arg1 )
        lv1 = rad / tan( θ1/2 ) # offset such that rods do not collide
        ℓ1  = sqrt( lv1^2 + rad^2 ) # absolute distance
        arg2 = dot(nvec0,nvec2)
        # arg = minimum( ( 1.0, arg) ) # guard for round-off errors
        # arg = maximum( (-1.0, arg) ) # guard for round-off errors
        θ2  = acos( arg2 )
        lv2 = rad / tan( θ2/2 ) # offset such that rods do not collide
        ℓ2  = sqrt( lv2^2 + rad^2 ) # absolute distance

        planeoff = dot( nprinter, Δn )
        cent1 = planeoff * cross( ncut1, cross(nprinter, ncut1 ) ) / norm( cross(ncut1, nprinter) )^2 # point of plane-plane intersection
        intersec1 = findIntersecEdge( nvec0, nline1, cent1, rad, ℓ1) + nodvec

        cent2 = planeoff * cross( ncut2, cross(nprinter, ncut2 ) ) / norm( cross(ncut2, nprinter) )^2 # point of plane-plane intersection
        intersec2 = findIntersecEdge( nvec0, nline2, cent2, rad, ℓ2) + nodvec

        # --- save this somewhere TODO for more efficient algorithm

        # --- write contour
        ledge = norm( mesh.p[nod0,:] - nodvec )
        if abs( dot(nvec0,nprinter) ) < 1e-10

            if dot( cross( intersec1 - nodvec, intersec2 - nodvec ), nprinter ) < 0.
                temp1 = intersec1
                intersec1 = intersec2
                intersec2 = temp1
            end
            # flat edge
            xptsfl[1] = intersec1
            xptsfl[2] = intersec1 + nvec0 * ( 0.5*ledge - dot(intersec1 - Δnodvec,nvec0) )
            xptsfl[3] = intersec2 + nvec0 * ( 0.5*ledge - dot(intersec2 - Δnodvec,nvec0) )
            xptsfl[4] = intersec2

            writeSTLcontour( xptsfl, Δnodvec, nprinter, 3, fid )
            factot[1] += 3

        else

            nperp = cross( nvec0, nprinter )
            nperp = nperp ./ norm(nperp)

            nzero = SVector{3}( nprinter[3], nprinter[1], nprinter[2] )
            if abs( abs(dot(nprinter,nvec0)) - 1.) < 1e-10
                nzero = cross(nzero,nprinter)
                nperp = cross(nprinter,nzero)
                nperp = nperp/norm(nperp)
            else
                nzero = cross(nvec0,nperp)
            end
            nzero = nzero/norm(nzero)

            intsec0, d0 = findPlaneLineIntersec( zpart, nodvec, nprinter, nvec0 )

            # first line
            currpt = nodvec + nzero*rad
            intsec1, = findPlaneLineIntersec( zpart, currpt, nprinter, nvec0 )

            # second line
            currpt = nodvec - nzero*rad
            intsec2, = findPlaneLineIntersec( zpart, currpt, nprinter, nvec0 )

            # find angle limits due to intersection
            vec1 = -rad * nperp
            vec2 =  intsec1 - intsec0
            if dot( cross(vec1,vec2), nprinter ) < 0.
                vec2 = -vec2
            end

            normvec1 = norm(vec1)
            normvec2 = norm(vec2)

            Δvec = intersec1 - intsec0
            φ1 = atan2( dot(vec2,Δvec)/normvec2^2, dot(vec1,Δvec)/normvec1^2 )

            Δvec = intersec2 - intsec0
            φ2 = atan2( dot(vec2,Δvec)/normvec2^2, dot(vec1,Δvec)/normvec1^2 )

            φmin = min(φ1,φ2)
            φmax = max(φ1,φ2)
            if (φmax - φmin) > π
                φ1 = φmax
                φmax = φmin
                φmin = φ1 - 2*π
            end

            φ = linspace( φmin, φmax, n+1 )
            for jj in 1:n+1
                xpts[jj] = intsec0 + vec1 * cos(φ[jj]) + vec2 * sin(φ[jj])
            end

            writeSTLcontour( xpts, intsec0, nprinter, n, fid )
            factot[1] += n

            # close the contour
            xptscl[1] = xpts[1]
            xptscl[3] = xpts[end]
            xptscl[2] = intsec0

            if d0 > 0.
                writeSTLcontour( xptscl, Δnodvec, nprinter, 2, fid )
                factot[1] += 2
            end

        end

    end

end

function findIntersecEdge( nvec1::SVector{3,Float64}, nline::SVector{3,Float64},
                           Δn::SVector{3,Float64}, r::Float64, ℓ::Float64 )

    # This uses the bisection method, so temporary.

    b1 = 0.
    b2 = 2*ℓ

    b = 0.
    rcomp = 0.
    while (b2 - b1)/b2 > 1e-5

        b = (b2 + b1)/2
        tempvec = b*nline + Δn
        rcomp = norm( tempvec - dot(tempvec,nvec1)*nvec1 )

        if rcomp > r
            b2 = b
        elseif rcomp < r
            b1 = b
        else
            break
        end

    end

    return (b*nline + Δn)

end


text:
abstract type PRModel <: ABCubicModel end

struct PRParam <: EoSParam
    a::PairParam{Float64}
    b::PairParam{Float64}
    Tc::SingleParam{Float64}
    Pc::SingleParam{Float64}
    Mw::SingleParam{Float64}
end

struct PR{T <: IdealModel,α,c,γ} <:PRModel
    components::Array{String,1}
    icomponents::UnitRange{Int}
    alpha::α
    mixing::γ
    translation::c
    params::PRParam
    idealmodel::T
    absolutetolerance::Float64
    references::Array{String,1}
end

@registermodel PR

export PR
function PR(components::Vector{String}; idealmodel=BasicIdeal,
    alpha = PRAlpha,
    mixing = vdW1fRule,
    activity=nothing,
    translation=NoTranslation,
    userlocations=String[],
    ideal_userlocations=String[],
    alpha_userlocations = String[],
    mixing_userlocations = String[],
    activity_userlocations = String[],
    translation_userlocations = String[],
     verbose=false)
    params = getparams(components, ["properties/critical.csv", "properties/molarmass.csv","SAFT/PCSAFT/PCSAFT_unlike.csv"]; userlocations=userlocations, verbose=verbose)
    k  = params["k"]
    pc = params["pc"]
    Mw = params["Mw"]
    Tc = params["Tc"]
    init_mixing = init_model(mixing,components,activity,mixing_userlocations,activity_userlocations,verbose)
    a,b = ab_premixing(PR,init_mixing,Tc,pc,k)
    init_idealmodel = init_model(idealmodel,components,ideal_userlocations,verbose)
    init_alpha = init_model(alpha,components,alpha_userlocations,verbose)
    init_translation = init_model(translation,components,translation_userlocations,verbose)
    icomponents = 1:length(components)
    packagedparams = PRParam(a,b,Tc,pc,Mw)
    references = String[]
    model = PR(components,icomponents,init_alpha,init_mixing,init_translation,packagedparams,init_idealmodel,1e-12,references)
    return model
end

function ab_consts(::Type{<:PRModel})
    return 0.457235,0.077796
end

function cubic_abp(model::PRModel, V, T, z)
    n = sum(z)
    āᾱ ,b̄, c̄ = cubic_ab(model,V,T,z,n)
    v = V/n+c̄
    _1 = one(b̄)
    denom = evalpoly(v,(-b̄*b̄,2*b̄,_1))
    p = R̄*T/(v-b̄) - āᾱ /denom
    return āᾱ, b̄, p
end

function cubic_poly(model::PRModel,p,T,z)
    a,b,c = cubic_ab(model,p,T,z)
    RT⁻¹ = 1/(R̄*T)
    A = a*p*RT⁻¹*RT⁻¹
    B = b*p*RT⁻¹
    k₀ = B*(B*(B+1.0)-A)
    k₁ = -B*(3*B+2.0) + A
    k₂ = B-1.0
    k₃ = 1.0
    return [k₀,k₁,k₂,k₃],c
end
#=
 (-B2-2(B2+B)+A)
 (-B2-2B2-2B+A)
 (-3B2-2B+A)
=#
function a_res(model::PRModel, V, T, z)
    n = sum(z)
    ā,b̄,c̄ = cubic_ab(model,V,T,z,n)
    Δ1 = 1+√2
    Δ2 = 1-√2
    ΔPRΔ = 2*√2
    RT⁻¹ = 1/(R̄*T)
    ρt = (V/n+c̄)^(-1) # translated density
    ρ  = n/V
    return -log(1+(c̄-b̄)*ρ) - ā*RT⁻¹*log((Δ1*b̄*ρt+1)/(Δ2*b̄*ρt+1))/(ΔPRΔ*b̄)

    #return -log(V-n*b̄) + āᾱ/(R̄*T*b̄*2^(3/2)) * log((2*V-2^(3/2)*b̄*n+2*b̄*n)/(2*V+2^(3/2)*b̄*n+2*b̄*n))
end

cubic_zc(::PRModel) = 0.3074


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

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

                     @variable(m, -23.82106505746807 <= q <= 20.87497636142748)

                     add_NL_constraint(m, :(1/(1 + exp(-(-0.7432570628691102 + 0.9028765283094948*1/(1 + exp(-(-0.11514299287103702 + 0.5608817277242992*$(x[1]) + 0.6700218006141032*$(x[2]) + -0.9693868545458724*$(x[3])))) + -0.9177729976456397*1/(1 + exp(-(-0.5568784315984976 + 0.281764667612443*$(x[1]) + 0.42150753955170206*$(x[2]) + -0.5829321104114529*$(x[3])))) + -0.6091404444372248*1/(1 + exp(-(-0.04388549286647825 + -0.9010617152865183*$(x[1]) + -0.7540101668519821*$(x[2]) + -0.4758905699629783*$(x[3])))) + 0.7075296204914552*1/(1 + exp(-(-0.511868403698823 + 0.5804725192481515*$(x[1]) + 0.7124288843879394*$(x[2]) + -0.06834010138653612*$(x[3])))) + -0.9115635329002854*1/(1 + exp(-(-0.9281357966735606 + -0.33421670762637223*$(x[1]) + -0.5412233562926474*$(x[2]) + -0.2936515077150612*$(x[3]))))))) + 1/(1 + exp(-(-0.2671327260424561 + -0.018499351907995898*1/(1 + exp(-(-0.11514299287103702 + 0.5608817277242992*$(x[1]) + 0.6700218006141032*$(x[2]) + -0.9693868545458724*$(x[3])))) + 0.12529725503678657*1/(1 + exp(-(-0.5568784315984976 + 0.281764667612443*$(x[1]) + 0.42150753955170206*$(x[2]) + -0.5829321104114529*$(x[3])))) + 0.9839396030143157*1/(1 + exp(-(-0.04388549286647825 + -0.9010617152865183*$(x[1]) + -0.7540101668519821*$(x[2]) + -0.4758905699629783*$(x[3])))) + 0.28500560463058955*1/(1 + exp(-(-0.511868403698823 + 0.5804725192481515*$(x[1]) + 0.7124288843879394*$(x[2]) + -0.06834010138653612*$(x[3])))) + 0.970334521625404*1/(1 + exp(-(-0.9281357966735606 + -0.33421670762637223*$(x[1]) + -0.5412233562926474*$(x[2]) + -0.2936515077150612*$(x[3]))))))) + 1/(1 + exp(-(0.40040910594674495 + -0.7454849743358074*1/(1 + exp(-(-0.11514299287103702 + 0.5608817277242992*$(x[1]) + 0.6700218006141032*$(x[2]) + -0.9693868545458724*$(x[3])))) + 0.7108986625228169*1/(1 + exp(-(-0.5568784315984976 + 0.281764667612443*$(x[1]) + 0.42150753955170206*$(x[2]) + -0.5829321104114529*$(x[3])))) + -0.30349624151173904*1/(1 + exp(-(-0.04388549286647825 + -0.9010617152865183*$(x[1]) + -0.7540101668519821*$(x[2]) + -0.4758905699629783*$(x[3])))) + -0.26057109207128715*1/(1 + exp(-(-0.511868403698823 + 0.5804725192481515*$(x[1]) + 0.7124288843879394*$(x[2]) + -0.06834010138653612*$(x[3])))) + 0.7858005505148937*1/(1 + exp(-(-0.9281357966735606 + -0.33421670762637223*$(x[1]) + -0.5412233562926474*$(x[2]) + -0.2936515077150612*$(x[3]))))))) + 1/(1 + exp(-(-0.7120753920031948 + -0.005433419997351141*1/(1 + exp(-(-0.11514299287103702 + 0.5608817277242992*$(x[1]) + 0.6700218006141032*$(x[2]) + -0.9693868545458724*$(x[3])))) + -0.8839514190404207*1/(1 + exp(-(-0.5568784315984976 + 0.281764667612443*$(x[1]) + 0.42150753955170206*$(x[2]) + -0.5829321104114529*$(x[3])))) + -0.2933551592941672*1/(1 + exp(-(-0.04388549286647825 + -0.9010617152865183*$(x[1]) + -0.7540101668519821*$(x[2]) + -0.4758905699629783*$(x[3])))) + -0.5418193939730438*1/(1 + exp(-(-0.511868403698823 + 0.5804725192481515*$(x[1]) + 0.7124288843879394*$(x[2]) + -0.06834010138653612*$(x[3])))) + 0.8978172795430157*1/(1 + exp(-(-0.9281357966735606 + -0.33421670762637223*$(x[1]) + -0.5412233562926474*$(x[2]) + -0.2936515077150612*$(x[3]))))))) + 1/(1 + exp(-(0.7950909223530527 + -0.6380481687930124*1/(1 + exp(-(-0.11514299287103702 + 0.5608817277242992*$(x[1]) + 0.6700218006141032*$(x[2]) + -0.9693868545458724*$(x[3])))) + -0.8467501205760102*1/(1 + exp(-(-0.5568784315984976 + 0.281764667612443*$(x[1]) + 0.42150753955170206*$(x[2]) + -0.5829321104114529*$(x[3])))) + 0.5001198683847332*1/(1 + exp(-(-0.04388549286647825 + -0.9010617152865183*$(x[1]) + -0.7540101668519821*$(x[2]) + -0.4758905699629783*$(x[3])))) + -0.03855272462979853*1/(1 + exp(-(-0.511868403698823 + 0.5804725192481515*$(x[1]) + 0.7124288843879394*$(x[2]) + -0.06834010138653612*$(x[3])))) + 0.3301524112226937*1/(1 + exp(-(-0.9281357966735606 + -0.33421670762637223*$(x[1]) + -0.5412233562926474*$(x[2]) + -0.2936515077150612*$(x[3]))))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions