text:
"""
    boykov_kolmogorov_impl(residual_graph, source, target, capacity_matrix)

Compute the max-flow/min-cut between `source` and `target` for `residual_graph`
using the Boykov-Kolmogorov algorithm.

Return the maximum flow in the network, the flow matrix and the partition
`{S,T}` in the form of a vector of 0's, 1's and 2's.

### References
- BOYKOV, Y.; KOLMOGOROV, V., 2004. An Experimental Comparison of
Min-Cut/Max-Flow Algorithms for Energy Minimization in Vision.

### Author
- Júlio Hoffimann Mendes (juliohm@stanford.edu)
"""
function boykov_kolmogorov_impl end
@traitfn function boykov_kolmogorov_impl(
    residual_graph::::IsDirected,          # the input graph
    source::Integer,                           # the source vertex
    target::Integer,                           # the target vertex
    capacity_matrix::AbstractMatrix    # edge flow capacities
    )

    T = eltype(capacity_matrix)
    U = eltype(residual_graph)

    n = nv(residual_graph)

    flow = 0
    flow_matrix = zeros(T, n, n)

    TREE = zeros(U, n)
    TREE[source] = U(1)
    TREE[target] = U(2)

    PARENT = zeros(U, n)

    A = [source,target]
    O = Vector{U}()

    while true
        # growth stage
        path = find_path!(residual_graph, source, target, flow_matrix, capacity_matrix, PARENT, TREE, A)

        isempty(path) && break

        # augmentation stage
        flow += augment!(path, flow_matrix, capacity_matrix, PARENT, TREE, O)

        # adoption stage
        adopt!(residual_graph, source, target, flow_matrix, capacity_matrix, PARENT, TREE, A, O)
    end

    return flow, flow_matrix, TREE
end

@traitfn function find_path!(
    residual_graph::::IsDirected,            # the input graph
    source::Integer,                       # the source vertex
    target::Integer,                       # the target vertex
    flow_matrix::AbstractMatrix,       # the current flow matrix
    capacity_matrix::AbstractMatrix,   # edge flow capacities
    PARENT::Vector,                     # parent table
    TREE::Vector,                       # tree table
    A::Vector                           # active set
    )
    T = eltype(residual_graph)
    tree_cap(p,q) = TREE[p] == one(T) ? capacity_matrix[p,q] - flow_matrix[p,q] :
    capacity_matrix[q,p] - flow_matrix[q,p]
    while !isempty(A)
        p = last(A)
        for q in neighbors(residual_graph, p)
            if tree_cap(p,q) > 0
                if TREE[q] == zero(T)
                    TREE[q] = TREE[p]
                    PARENT[q] = p
                    unshift!(A, q)
                end
                if TREE[q] ≠ zero(T) && TREE[q] ≠ TREE[p]
                    # p -> source
                    path_to_source = [p]
                    while PARENT[p] ≠ zero(T)
                        p = PARENT[p]
                        push!(path_to_source, p)
                    end

                    # q -> target
                    path_to_target = [q]
                    while PARENT[q] ≠ zero(T)
                        q = PARENT[q]
                        push!(path_to_target, q)
                    end

                    # source -> target
                    path = [reverse!(path_to_source); path_to_target]

                    if path[1] == source && path[end] == target
                        return path
                    elseif path[1] == target && path[end] == source
                        return reverse!(path)
                    end
                end
            end
        end
        pop!(A)
    end

    return Vector{T}()
end

function augment!(
    path::AbstractVector,                  # path from source to target
    flow_matrix::AbstractMatrix,       # the current flow matrix
    capacity_matrix::AbstractMatrix,   # edge flow capacities
    PARENT::Vector,                     # parent table
    TREE::Vector,                       # tree table
    O::Vector                           # orphan set
    )

    T = eltype(path)
    # bottleneck capacity
    Δ = Inf
    for i=1:length(path)-1
        p, q = path[i:i+1]
        cap = capacity_matrix[p,q] - flow_matrix[p,q]
        cap < Δ && (Δ = cap)
    end

    # update residual graph
    for i=1:length(path)-1
        p, q = path[i:i+1]
        flow_matrix[p,q] += Δ
        flow_matrix[q,p] -= Δ

        # create orphans
        if flow_matrix[p,q] == capacity_matrix[p,q]
            if TREE[p] == TREE[q] == one(T)
                PARENT[q] = zero(T)
                unshift!(O, q)
            end
            if TREE[p] == TREE[q] == 2
                PARENT[p] = zero(T)
                unshift!(O, p)
            end
        end
    end

    return Δ
end

@traitfn function adopt!(
    residual_graph::::IsDirected,               # the input graph
    source::Integer,                             # the source vertex
    target::Integer,                             # the target vertex
    flow_matrix::AbstractMatrix,       # the current flow matrix
    capacity_matrix::AbstractMatrix,   # edge flow capacities
    PARENT::Vector,                     # parent table
    TREE::Vector,                       # tree table
    A::Vector,                             # active set
    O::Vector                           # orphan set
    )

    T = eltype(residual_graph)
    tree_cap(p,q) = TREE[p] == 1 ? capacity_matrix[p,q] - flow_matrix[p,q] :
    capacity_matrix[q,p] - flow_matrix[q,p]
    while !isempty(O)
        p = pop!(O)
        # try to find parent that is not an orphan
        parent_found = false
        for q in neighbors(residual_graph, p)
            if TREE[q] == TREE[p] && tree_cap(q,p) > 0
                # check if "origin" is either source or target
                o = q
                while PARENT[o] ≠ 0
                    o = PARENT[o]
                end
                if o == source || o == target
                    parent_found = true
                    PARENT[p] = q
                    break
                end
            end
        end

        if !parent_found
            # scan all neighbors and make the orphan a free node
            for q in neighbors(residual_graph, p)
                if TREE[q] == TREE[p]
                    if tree_cap(q,p) > 0
                        unshift!(A, q)
                    end
                    if PARENT[q] == p
                        PARENT[q] = zero(T)
                        unshift!(O, q)
                    end
                end
            end

            TREE[p] = zero(T)
            B = setdiff(A, p)
            resize!(A, length(B))[:] = B
        end
    end
end

max_stars_repo_name: AsileBcd/JohnsonCycles

max_issues_repo_name: AsileBcd/JohnsonCycles

max_forks_repo_name: AsileBcd/JohnsonCycles


text:
function StatsBase.coefnames(m::MixedModel)
    Xtrm = m.feterm
    return invpermute!(copy(Xtrm.cnames), Xtrm.piv)
end

"""
cond(m::MixedModel)

Return a vector of condition numbers of the λ matrices for the random-effects terms
"""
LinearAlgebra.cond(m::MixedModel) = cond.(m.λ)

function StatsBase.dof(m::MixedModel)
    return m.feterm.rank + length(m.parmap) + dispersion_parameter(m)
end

"""
    dof_residual(m::MixedModel)

  Return the residual degrees of freedom of the model.

!!! note
  The residual degrees of freedom for mixed-effects models is not clearly defined due to partial pooling.
  The classical `nobs(m) - dof(m)` fails to capture the extra freedom granted by the random effects, but
  `nobs(m) - nranef(m)` would overestimate the freedom granted by the random effects. `nobs(m) - sum(leverage(m))`
  provides a nice balance based on the relative influence of each observation, but is computationally
  expensive for large models. This problem is also fundamentally related to [long-standing debates](https://bbolker.github.io/mixedmodels-misc/glmmFAQ.html#why-doesnt-lme4-display-denominator-degrees-of-freedomp-values-what-other-options-do-i-have)
  about the appropriate treatment of the denominator degrees of freedom for ``F``-tests.
  In the future, MixedModels.jl may provide additional methods allowing the user to choose the computation
  to use.

!!! warning
   Currently, the residual degrees of freedom is computed as `nobs(m) - dof(m)`, but this may change in
   the future without being considered a breaking change because there is no canonical definition of the
   residual degrees of freedom in a mixed-effects model.
"""
function StatsBase.dof_residual(m::MixedModel)
    # a better estimate might be nobs(m) - sum(leverage(m))
    # this version subtracts the number of variance parameters which isn't really a dimensional
    # and doesn't even agree with the definition for linear models
    return nobs(m) - dof(m)
end

"""
    issingular(m::MixedModel, θ=m.θ)

Test whether the model `m` is singular if the parameter vector is `θ`.

Equality comparisons are used b/c small non-negative θ values are replaced by 0 in `fit!`.
"""
issingular(m::MixedModel, θ=m.θ) = any(lowerbd(m) .== θ)

# FIXME: better to base this on m.optsum.returnvalue
StatsBase.isfitted(m::MixedModel) = m.optsum.feval > 0

StatsBase.meanresponse(m::MixedModel) = mean(m.y)

"""
    modelmatrix(m::MixedModel)

Returns the model matrix `X` for the fixed-effects parameters, as returned by [`coef`](@ref).

This is always the full model matrix in the original column order and from a field in the model
struct.  It should be copied if it is to be modified.
"""
StatsBase.modelmatrix(m::MixedModel) = m.feterm.x

StatsBase.nobs(m::MixedModel) = length(m.y)

StatsBase.predict(m::MixedModel) = fitted(m)

function retbl(mat, trm)
    return merge(
        NamedTuple{(fname(trm),)}((trm.levels,)),
        columntable(Tables.table(transpose(mat); header=Symbol.(trm.cnames))),
    )
end

"""
    raneftables(m::MixedModel; uscale = false)

Return the conditional means of the random effects as a `NamedTuple` of columntables
"""
function raneftables(m::MixedModel{T}; uscale=false) where {T}
    return NamedTuple{fnames(m)}((map(retbl, ranef(m; uscale=uscale), m.reterms)...,))
end

StatsBase.residuals(m::MixedModel) = response(m) .- fitted(m)

"""
    response(m::MixedModel)

Return the response vector for the model.

For a linear mixed model this is a `view` of the last column of the `XyMat` field.
For a generalized linear mixed model this is the `m.resp.y` field.
In either case it should be copied if it is to be modified.
"""
StatsBase.response(m::MixedModel) = m.y

function StatsBase.responsename(m::MixedModel)
    cnm = coefnames(m.formula.lhs)
    return isa(cnm, Vector{String}) ? first(cnm) : cnm
end

function σs(m::MixedModel)
    σ = dispersion(m)
    return NamedTuple{fnames(m)}(((σs(t, σ) for t in m.reterms)...,))
end

function σρs(m::MixedModel)
    σ = dispersion(m)
    return NamedTuple{fnames(m)}(((σρs(t, σ) for t in m.reterms)...,))
end

"""
    size(m::MixedModel)

Returns the size of a mixed model as a tuple of length four:
the number of observations, the number of (non-singular) fixed-effects parameters,
the number of conditional modes (random effects), the number of grouping variables
"""
function Base.size(m::MixedModel)
    dd = m.dims
    return dd.n, dd.p, sum(size.(m.reterms, 2)), dd.nretrms
end

"""
    vcov(m::MixedModel; corr=false)

Returns the variance-covariance matrix of the fixed effects.
If `corr` is `true`, the correlation of the fixed effects is returned instead.
"""
function StatsBase.vcov(m::MixedModel; corr=false)
    Xtrm = m isa GeneralizedLinearMixedModel ? m.LMM.feterm : m.feterm
    iperm = invperm(Xtrm.piv)
    p = length(iperm)
    r = Xtrm.rank
    Linv = inv(feL(m))
    T = eltype(Linv)
    permvcov = dispersion(m, true) * (Linv'Linv)
    if p == Xtrm.rank
        vv = permvcov[iperm, iperm]
    else
        covmat = fill(zero(T) / zero(T), (p, p))
        for j in 1:r, i in 1:r
            covmat[i, j] = permvcov[i, j]
        end
        vv = covmat[iperm, iperm]
    end

    return corr ? StatsBase.cov2cor!(vv, stderror(m)) : vv
end

StatsModels.formula(m::MixedModel) = m.formula

ext: jl


text:
## warning
function warn_is_registered(env::AbstractEnv)
    is_registered = NestedEnvironments.__REGISTERED_ENVS.__envs |> Map(__env -> typeof(__env) == typeof(env)) |> collect |> any
    return is_registered ? nothing : error("Unregistered env: $(typeof(env))")
end

## functions
# env names
function Base.names(env::AbstractEnv)
    return [name for name in fieldnames(typeof(env)) if typeof(getfield(env, name)) <: AbstractEnv]
end
# extend size for NamedTuple
function Base.size(x0::NamedTuple)
    env_names = keys(x0)
    if env_names == []
        return size(x0)
    else
        env_sizes = env_names |> Map(name -> size(x0[name]))
        return (; zip(env_names, env_sizes)...)  # NamedTuple
    end
end
function flatten_length(x0::NamedTuple)
    env_names = keys(x0)
    if env_names == []
        return prod(size(x0))
    else
        return env_names |> Map(name -> flatten_length(x0[name])) |> sum
    end
end
flatten_length(x0) = prod(size(x0))

# get the index
function _index(env::AbstractEnv, x0, _range)
    env_names = names(env)
    if env_names == []
        return _range
    else
        env_accumulated_flatten_lengths = env_names |> Map(name -> flatten_length(x0[name])) |> Scan(+) |> collect
        range_first = first(_range)
        env_ranges_tmp = (range_first-1) .+ [0, env_accumulated_flatten_lengths...] |> Consecutive(length(env_accumulated_flatten_lengths); step=1)
        env_ranges = zip(env_ranges_tmp...) |> MapSplat((x, y) -> x+1:y)
        env_indices = zip(env_names, env_ranges) |> MapSplat((name, range) -> _index(getfield(env, name), x0[name], range))
        return (; zip(env_names, env_indices)...)  # NamedTuple
    end
end
# get env_index and env_isze
function _preprocess(env::AbstractEnv, x0)
    env_size_nt = size(x0)
    env_flatten_length = flatten_length(x0)
    env_index_nt = _index(env, x0, 1:env_flatten_length)
    return env_index_nt, env_size_nt
end

## main functions
# transform readable to raw (flatten)
function _raw(env::AbstractEnv, x)
    env_names = names(env)
    if env_names == []
        return x
    else
        _x = env_names |> MapCat(name -> _raw(getfield(env, name), x[name])) |> collect
        return _x
    end
end
# transform raw (flatten) to readable
function _readable(_x, env_index_nt, env_size_nt)
    if typeof(env_index_nt) <: AbstractRange
        if env_size_nt == ()
            return _x[env_index_nt][1]  # scalar
        else
            return reshape(_x[env_index_nt], env_size_nt...)
        end
    else
        index_names = keys(env_index_nt)
        processed_values = index_names |> Map(name -> _readable(_x, env_index_nt[name], env_size_nt[name]))
        return (; zip(index_names, processed_values)...)
    end
end


text:
"""
One step of Francis transformation Transform hessenberg matrix
H is in upper Hessenberg form.
shifts is a vector of shift values.
With each shift value, its conjugate complex value is implicitly used as shift
All shift values should be distinct.
"""
function transform_Hess!(A::AbstractMatrix{T}, ilo::Integer, ihi::Integer, Q::AbstractM, s::Vector{S}, maxchase::Integer, iwindow::Integer) where {T<:Union{Real,Complex}, S<:Union{Real,Complex}}

  !isa(Q, AbstractMatrix) || size(A, 2) == size(Q, 2) || error("A and Q have incompatible sizes")

  n = size(A, 1)
  z = zero(T)
  m = counteigs(s)
  realcase = T <: Real
  minn = ifelse(realcase, 2, 1)
  #println("def_sub ilo = $ilo ihi = $ihi")

  # supress further calculations if A is already quasi-diagonal
  isize = ihi - ilo + 1
  if isize <= 1 || isize == minn && discriminant(A, ilo, ihi) < z
    m = 0
    maxchase = 0
  end

  # Any shifts provided
  if m > 0

    # First step: shift bulges blocking the undisturbed insertion
    if isize > 4
      i0, m0 = lastbulge(A, ilo, ihi, m + 1)
      if i0 > 0
        chase!(A, ilo, ihi, Q, m + 2 - i0, iwindow)
      end
    end

    # Second step: compute  prod(A-s(k))e1
    pe = prode1(A, s, ilo)

    # Third step: set in new upper left bulge
    set_in!(A, ilo, Q, pe)
  end

  # Forth step: chase all remaining bulks, oldest first
  if maxchase > 0
    chase!(A, ilo, ihi, Q, maxchase, iwindow)
  end
  ilo, ihi, Q
end

function transform_Hess!(A::AbstractMatrix{T}, Q::AbstractM, s::Vector{S}, maxchase::Integer, iwindow::Integer = 0) where {T<:Real, S<:Union{Real,Complex}}

  transform_Hess!(A, 1, size(A, 1), Q, s, maxchase, iwindow)
end

function transform_Hess!(A::AbstractMatrix{T}, s::Vector{S}, maxchase::Integer, iwindow::Integer = 0) where {T<:Real, S<:Union{Real,Complex}}

  #Q = LinearAlgebra.Rotation(LinearAlgebra.Givens{T}[])
  Q = Matrix{T}(I, size(A, 2), size(A, 2))
  transform_Hess!(A, Q, s, maxchase, iwindow)
end

# calculate product over k of (A - s[k]*I) * e1
function prode1(A::AbstractMatrix{T}, s::Vector, ilo::Integer) where {T<:Real}
  n = size(A, 1)
  pe = zeros(T, n)
  pe[ilo] = one(T)
  k = j = 1
  while k <= length(s)
    sr = real(s[k])
    si = imag(s[k])
    if si == 0
      pe = A * pe - sr * pe
      j += 1
    else
      sr = 2sr
      s2 = abs2(s[k])
      pe = A * ( A * pe - sr * pe) + s2 * pe

      if k < length(s) && imag(s[k+1]) == -si
        k += 1
      end
      j += 2
    end
    k += 1
  end
  pe[1:ilo-1] = 0
  pe
end

function prode1(A::AbstractMatrix{T}, s::Vector, ilo::Integer) where {T<:Complex}
  n = size(A, 1)
  pe = zeros(T, n)
  pe[ilo] = one(T)
  k = j = 1
  for sk in s
    pe = A * pe - sk * pe
  end
  pe
end

function set_in!(A::AbstractMatrix, ilo::Integer, Q::AbstractM, pe::Vector)
  m2 = findlast(x -> x != 0, pe)
  for k = m2:-1:ilo+1
    G, r = givens(pe, k-1, k)
    pe[k-1] = r
    pe[k] = 0
    A_mul_B!(G, A)
    A_mul_Bc!(A, G)
    A_mul_Bc!(Q, G)
  end
end

"""
Count number of eigenvalues.
Each non-real eigenvalue is counted twice, if not paired with its conjugate.
"""
function counteig(s::Vector{Complex{T}}) where T
  z = zero(T)
  m = length(s)
  count = 0
  k = 1
  while k <= m
    sk = s[k]
    si = imag(sk)
    if si == z
      count += 1
    else
      if k < m && imag(s[k+1]) == -si
        k += 1
      end
      count += 2
    end
    k += 1
  end
  count
end

"""
Chase all bulges towards right lower corner
"""
function chase!(A, ilo::Integer, ihi::Integer, Q, maxchase::Integer, iwindow::Integer)
  m1, m2 = size(A)
  n2 = min(ihi, m2)
  n1 = min(ihi, m1)

  m = i0 = n2 - 1 # column to start with

  while i0 >= ilo && m > 0
    i0, m = lastbulge(A, ilo, ihi, i0 - 1)

    if m > 0
      i1 = i0 + m > ihi -iwindow ? n2 - 2 : min(i0 + maxchase - 1, n2 - 2)
      for  i = i0:min(i0 + maxchase - 1, n2 - 2)
        r = A[i+1,i]
        for k = min(i+m+1,n1):-1:i+2
          aki = A[k,i]
          if aki != 0
            rp = hypot(r, aki)
            c = r / rp
            s = aki / rp
            A[i+1,i] = r = rp
            A[k,i] = 0
            for l = i+1:m2
              A[i+1,l], A[k,l] = A[i+1,l] * c' + A[k,l] * s', A[k,l] * c - A[i+1,l] * s
            end
            for l = 1:m1
              A[l,i+1], A[l,k] = A[l,i+1] * c + A[l,k] * s, A[l,k] * c' - A[l,i+1] * s'
            end
            for l = 1:m1
              Q[l,i+1], Q[l,k] = Q[l,i+1] * c + Q[l,k] * s, Q[l,k] * c' - Q[l,i+1] * s'
            end
          elseif r != 0
            rp = abs(r)
            c = r / rp
            A[i+1,i] = r = rp
            for l = i+1:m2
              A[i+1,l] *= c'; A[k,l] *= c
            end
            for l = 1:m1
              A[l,i+1] *= c; A[l,k] *= c'
            end
            for l = 1:m1
              Q[l,i+1] = Q[l,i+1] * c; Q[l,k] = Q[l,k] * c'
            end
          end
        end
      end
    end
  end
  A, Q
end

@inline function givens_improved(A::AbstractMatrix{T}, i1::Int, i2::Int, j::Int) where T<:Number
  rcs = Vector{T}(3)
  rcs[1] = hypot(A[i1,j], A[i2,j])
  rcs[2] = A[i1,j] / rcs[1]
  rcs[3] = A[i2,j] / rcs[1]
  rcs
end

"""
Find last bulge start column and size. Start search in column i.
"""
function lastbulge(A::AbstractMatrix{T}, ilo::Integer, ihi::Integer, i::Integer) where {T<:Number}
  n, m = size(A)
  n == m || error("need square matrix")
  #i <= n - 1 || error("column number must be less than n")

  maxpos = 2
  while i >= ilo
    m = findlast(A[i+2:ihi,i])
    if m == 0
      if maxpos > 2
        break
      end
    else
      maxpos = max(maxpos, m + i + 1)
    end
    i -= 1
  end
  m = maxpos - i - 2
  m == 0 ? 0 : i + 1, m
end

max_stars_repo_name: KlausC/LinAlgBigFloat

max_issues_repo_name: KlausC/LinAlgBigFloat

max_forks_repo_name: KlausC/LinAlgBigFloat


text:
# functions for reading/writing .csv files with optimization parameters

"""
   read_parameters(file; dict=true)

Reads the optimization parameters in a csv file and returns either an ordered
dictionary or a named tuple of the optimization parameters.
"""
function read_parameters(file; dict=true)

    # read file
    (data, header) = DelimitedFiles.readdlm(file, ',', header=true, comments=true, skipblanks=true)

    # keep only lines that have variable names
    idx = findall((x)->!isempty(strip(x)), data[:,1])
    data = data[idx, :]

    # get ordering of each column of data
    icol_name, icol_x0, icol_lb, icol_ub, icol_scale, icol_dv, icol_description = get_column_ordering(header)

    # parse each column
    nparam = size(data, 1)
    names = Symbol.(strip.(data[:, icol_name]))
    initial_values = parse_value.(data[:, icol_x0])
    lower_bounds = parse_value.(data[:, icol_lb])
    upper_bounds = parse_value.(data[:, icol_ub])
    scaling = icol_scale == nothing ? nothing : parse_value.(data[:, icol_scale])
    design_variable = parse_value.(data[:, icol_dv])
    description = strip.(data[:, icol_description])

    # replace missing lower bound values with -Inf
    lower_bounds[findall(isnothing.(lower_bounds))] .= -Inf

    # replace missing upper bound values with Inf
    upper_bounds[findall(isnothing.(upper_bounds))] .= Inf

    # replace missing scaling values with `1.0`
    scaling[findall(isnothing.(scaling))] .= 1.0

    # replace missing design variable flags with `false`
    design_variable[findall(isnothing.(design_variable))] .= false

    # create a named tuple of optimization parameters
    optparams = NamedTuple{Tuple(names)}(Tuple(OptimizationParameter(initial_values[i],
        lb=lower_bounds[i], ub=upper_bounds[i], scaling=scaling[i], dv=design_variable[i],
        description=description[i]) for i = 1:nparam))

    # return a dictionary if specified
    if dict
        return named_tuple_to_dict(optparams)
    else
        return optparams
    end

end

"""
   write_parameters([inputfile, ]outputfile, optparams)

Writes the optimization parameters in `optparams` to a csv file.  If given, the
file `inputfile` is used as a template for the organization of the resulting file.
"""
function write_parameters(inputfile, outputfile, optparams)

    # read input file, skip nothing
    data = DelimitedFiles.readdlm(inputfile, ',', Any, skipstart=0, comments=false, skipblanks=false)
    header = data[1,:]

    # get position of each column
    icol_name, icol_x0, icol_lb, icol_ub, icol_scale, icol_dv, icol_description = get_column_ordering(header)

    # get input file names indexed by line number
    inputnames = Symbol.(strip.(data[:,icol_name]))

    # loop through each optimization parameter overwriting data from input file
    for (name, optparam) in pairs(optparams)
        # find line for this optimization parameter
        irow = findfirst(inputnames .== name)

        if isnothing(irow)
            error("Template input file does not have parameter: `$name`")
        else
            # overwrite initial value
            data[irow,icol_x0] = write_value(optparam.x0)

            # overwrite lower bound
            if all(optparam.lb .== -Inf)
                # clear if not manually set to -Inf
                if any(parse_value(data[irow, icol_lb]) .!= -Inf)
                    data[irow,icol_lb] = ""
                end
            else
                data[irow,icol_lb] = write_value(optparam.lb)
            end

            # overwrite upper bound
            if all(optparam.ub .== Inf)
                # clear if not manually set to Inf
                if any(parse_value(data[irow, icol_ub]) .!= Inf)
                    data[irow,icol_ub] = ""
                end
            else
                data[irow,icol_ub] = write_value(optparam.ub)
            end

            # overwrite scaling
            if icol_scale != nothing
                if all(optparam.scaling .== 1.0)
                    # clear if not manually set to 1.0
                    if any(parse_value(data[irow, icol_scale]) .!= 1.0)
                        data[irow, icol_scale] = ""
                    end
                else
                    data[irow,icol_scale] = write_value(optparam.scaling)
                end
            end

            # overwrite design variable flag
            if all(optparam.dv .== false)
                # clear if not manually set to false
                if any(parse_value(data[irow, icol_dv]) .!= false)
                    data[irow, icol_dv] = ""
                end
            else
                data[irow,icol_dv] = write_value(optparam.dv)
            end

            # overwrite descriptions
            data[irow,icol_description] = optparam.description

        end
    end

    # finally write output file
    DelimitedFiles.writedlm(outputfile, data, ',', quotes=true)

    return nothing
end

function write_parameters(outputfile, optparams; descriptions=nothing)

    data = fill("", length(optparams)+1, 7)

    data[1,:] .= ["Parameter", "Initial Value(s)", "Lower Bound(s)", "Upper Bound(s)",
        "Scaling", "Design Variable(s)?", "Description"]

    irow = 1

    # loop through each optimization parameter overwriting data from input file
    for (name, optparam) in pairs(optparams)

        irow += 1

        # parameter name
        data[irow,1] = string(name)

        # initial value(s)
        data[irow,2] = write_value(optparam.x0)

        # lower bound(s)
        data[irow,3] = write_value(optparam.lb)

        # upper bound(s)
        data[irow,4] = write_value(optparam.ub)

        # scaling
        data[irow,5] = write_value(optparam.scaling)

        # design variable flag(s)
        data[irow,6] = write_value(optparam.dv)

        # descriptions
        data[irow,7] = optparam.description

    end

    # finally write output file
    DelimitedFiles.writedlm(outputfile, data, ',', quotes=true)

    return nothing
end

"""
    write_value(value)

Returns a string representing `value` using a set of custom rules.
"""
write_value

# by default writing a value writes a string
write_value(value) = string(value)

# for a string, we need to add quotes
write_value(value::AbstractString) = "\""*string(value)*"\""

# a value of nothing returns an empty string
write_value(value::Nothing) = ""

# an array prints a generic empty array if it is empty
function write_value(value::AbstractArray)
    if isempty(value)
        result = "[]"
    else
        result = string(value)
    end
    return result
end

"""
    get_column_ordering(header)

Returns the column ordering given the file header
"""
function get_column_ordering(header)

    # get position of each column
    header = lowercase.(strip.(header))[:]
    icol_name = findfirst(header .== "parameter")
    icol_x0 = findfirst(header .== "initial value(s)")
    icol_lb = findfirst(header .== "lower bound(s)")
    icol_ub = findfirst(header .== "upper bound(s)")
    icol_scale = findfirst(header .== "scaling")
    icol_dv = findfirst(header .== "design variable(s)?")
    icol_description = findfirst(header .== "description")

    # if any of the names can't be resolved, use the default ordering
    if any([icol_name, icol_x0, icol_lb, icol_ub, icol_scale, icol_dv, icol_description] .== nothing)
        # scaling column doesn't have to exist
        if length(header) == 6
            icol_name = 1
            icol_x0 = 2
            icol_lb = 3
            icol_ub = 4
            icol_scale = nothing
            icol_dv = 5
            icol_description = 6
        else
            icol_name = 1
            icol_x0 = 2
            icol_lb = 3
            icol_ub = 4
            icol_scale = 5
            icol_dv = 6
            icol_description = 7
        end
    end

    return icol_name, icol_x0, icol_lb, icol_ub, icol_scale, icol_dv, icol_description
end

"""
    parse_value(value)

Returns `value` if `value` is not a string.

If `value` is a string this function checks if `value` should have been parsed
as a boolean and returns the corresponding boolean.

Otherwise `value` is parsed as a Julia expression, with a few customizations:
- [] parses to Float64[]
- The unicode characters “ and ” are converted to "
"""
parse_value(value) = value

# for a string, special parsing is needed
function parse_value(str::AbstractString)
    str = strip(str)
    # check first for boolean input
    if lowercase(str) == "true"
        value = true
    elseif lowercase(str) == "false"
        value = false
    else # if not boolean parse as a Julia expression
        str = replace(str, "[]"=>"Float64[]")
        str = replace(str, "“"=>"\"")
        str = replace(str, "”"=>"\"")
        # str = ascii(str)
        # str = replace(str, r"(\"[^\"]*)([^\"]*(\"|$))" => s"FixedSizeString(\1\2)")
        value = eval(Meta.parse(str))
    end

    return value
end


text:
function test4()
  filename = "test4.asc"
  # exectuablepath = null string will not run LTspice.exe.  Test parsing only.
  sim = LTspiceSimulation(filename,executablepath="")
  @test LTspice.does_circuitfilearray_file_match(sim)
  show(IOBuffer(),sim)

  asteps = [1.0,2.0]
  bsteps = [10.0,15.0,20.0,25.0]

  @test stepvalues(sim) == (asteps,bsteps)
  @test measurementnames(sim) == ("sum","sump1000")
  @test stepnames(sim) == ("a","b")

  @static if Sys.iswindows()
      @test circuitpath(sim) == filename
  end

  @test parametervalues(sim) == []
  @test measurementvalues(sim)[1,1,1] == 11.0
  @test executablepath(sim) == ""
  @test haskey(sim,"sum") == true
  @test length(sim) == 16

  for i in eachindex(asteps)
      for j in eachindex(bsteps)
          @test measurementvalues(sim)[i,j,1] == asteps[i]+bsteps[j]
          @test sim["sum"][i,j] == asteps[i]+bsteps[j]
          @test measurementvalues(sim)[i,j,2] == asteps[i]+bsteps[j]+1000
          @test sim["sump1000"][i,j] == asteps[i]+bsteps[j]+1000
      end
  end

  for line in perlineiterator(sim)
      @test (line[1]+line[2] == line[3])
      @test (line[3] + 1000 == line[4])
  end

  for line in perlineiterator(sim,resultnames=["sum"],steporder=["b","a"])
      @test (line[1]+line[2] == line[3])
  end

  show(IOBuffer(),sim)
end
test4()


text:
module Collections

# Import the necessary packages.
using ..InfiniteOpt
import JuMP
const JuMPC = JuMP.Containers

include("VectorTuple.jl")
include("DualDict.jl")

end # end module


text:
using JuLIP, Test
using JuLIP.Testing
using JuLIP.FIO

h3("Testing single `Atoms` <-> `Dict`")
at = bulk(:Cu, cubic=true) * 3
set_pbc!(at, (true, false, true))
rattle!(at, 0.1)
D = write_dict(at)
at1 = Atoms(D)
println(@test(at == at1))
at2 = read_dict(D)
println(@test(at == at2))

h3("Test JSON fio")
fn = tempname()
save_dict(fn, D)
D1 = load_dict(fn)
# D1 == D  => this will be false so don't test it!
at3 = read_dict(D1)
println(@test at3 == at1)

h3("Test array of Atoms <-> Dict")
ats = [ (bulk(:Cu) * rand(2:4)) for n = 1:5 ]
Ds = Dict("ats" => write_dict.(ats))
ats1 = read_dict.(Ds["ats"])
println(@test ats1 == ats)

h3("Test JSON fio for array")
fn = tempname()
save_dict(fn, Ds)
Ds1 = load_dict(fn)
ats2 = read_dict.(Ds1["ats"])
println(@test ats == ats2)


text:
using ErlangTerm
import JSON, JSONWebTokens

"Serialize `data`, sha-256 encoded with `key` if it is not empty."
function serializek(data, key::AbstractString)
    if isempty(key)
        serialize(data)
    else
        JSONWebTokens.encode(
            JSONWebTokens.HS256(key),
            data
            |> serialize
            |> JSON.json
        ) |> Vector{UInt8}
    end
end

"Deserialize a `binary`, sha-256 decoded with `key` if it is not empty."
function deserializek(binary::Vector{UInt8}, key::AbstractString)
    if isempty(key)
        deserialize(binary)
    else
        j = JSONWebTokens.decode(
            JSONWebTokens.HS256(key),
            String(binary)
        )
        deserialize(convert(Vector{UInt8}, j))
    end
end

max_stars_repo_name: pbayer/erjulix

max_issues_repo_name: pbayer/erjulix

max_forks_repo_name: pbayer/erjulix


text:
#=
verifier.jl
this script implements functionality to independently verify a solution's correctness.
=#

using LightGraphs
using MetaGraphs
using Statistics

export VerifySolution, BasicStats

#=verifies if a given solution is valid per the LPBasic requirements
#this includes checking if:
    minimum and maximum frequency requirements are met

#return value:
    a tuple: (true, "Valid") - indicating that the solution is valid
             (false, error_string) - indicating that the solution is not valid and the reason why
=#
"""
    VerifySolution(instance::ProblemInstance, solution::Array{Bool,1})
Verify if a given solution is feasible per the LPBasic requirements.

Also check for construction errors in the instance and solution itself
Return a tuple (true,"Valid") or (false, <error_string>)
"""
function VerifySolution(instance::ProblemInstance, solution::Array{Bool,1})

    #we supose the instance was already validated?
    VerifyInstance(instance)

    #solution and instance match in number of lines?
    if(length(instance.lines) != length(solution))
        return false, "Number of lines in instance and solution do not match"
    end

    #minimum and maximum frequency requirements are met:
    #for each edge in the network,
    for e in edges(instance.network)
        #its coverage frequency
        f_min = get_prop(instance.network, e, :f_min)
        f_max = get_prop(instance.network, e, :f_max)

        coverage_freq = sum( if (solution[i]) IsEdgeOnLine(instance, e, i) * instance.lines[i].frequency else 0 end for i in 1:length(solution))

        if coverage_freq < f_min
            return false, "Coverage frequency is too low on edge " * repr(e) * ": " * string(coverage_freq) * " < " * string(f_min)
        elseif coverage_freq > f_max
            return false, "Coverage frequency too high on edge " * repr(e)* ": " * string(coverage_freq) * " > " + string(f_max)
        end
    end

    return true, "Valid"
end

"""
    CalculateCost(instance::ProblemInstance, solution::Array{Bool,1})::Float64
Calculate total line cost for given solution on given instance.

A solution for this method is an array of booleans indicating accepting or rejecting a line with its basic frequency
"""
function CalculateCost(instance::ProblemInstance, solution::Array{Bool,1})::Float64
    cost = Float64(0.0)
    for i in 1:length(instance.lines)
        if solution[i]
            cost += instance.lines[i].cost * instance.lines[i].frequency
        end
    end
    return cost
end

"""
    CalculateCost(instance::ProblemInstance, solution::Array{Bool,1})::Float64
Calculate total line cost for given solution on given instance.

A solution for this method is an array of integers representing each line's taken frequency
"""
function CalculateCost(instance::ProblemInstance, solution::Array{T,1} where T<:Integer)
    println("Calculate cost Integer")
    return sum(instance.lines[i].cost * solution[i] for i in 1:length(instance.lines))
end

#to do: implement some extra validation criteria, such as:

#returns a dictionary with relevant basic information about the instance-solution pair
"""
    BasicStats(instance::ProblemInstance, solution::Array{Bool,1})
Boolean/binary variant of the BasicStats function.

In this solution, each line has 0 frequency or its base-frequency
"""
function BasicStats(instance::ProblemInstance, solution::Array{Bool,1})

    return BasicStats(instance, collect(Int64, if solution[i] instance.lines[i].frequency else 0 end for i in 1:length(solution)))

end

"""
    BasicStats(instance::ProblemInstance, solution::Array{T,1} where {T<:Integer})
Populate and return a dictionary with basic information about the instance and solution.

The returned dictinary has the following key/value pairings:
- QntLines: the quantity of lines considered in the instance
- QntStops: the quantity of stops/nodes in the transportation network/graph
- QntEdges: the quantity of conections/edges in the transportation network/graph
- TakenLines: the quantity of lines used in the solution(taken frequency > 0)
- AverageFrequency: the average frequency of service in all conections/edges
- StandardDeviation: the standard deviation for the average above
- Variance: the variance for the average above
- Median: the median value for the average above
- LowestFrequency: the lowest frequency of service among all conections/edges
- HighestFrequency: the highest frequency of service among all conections/edges
- LineCost: the sum of costs for taking all lines in the solution with the given frequency
"""
function BasicStats(instance::ProblemInstance, solution::Array{T,1} where {T<:Integer})

    stats = Dict{Symbol, Any}()
    edge_frequencies = CalculateFrequencyArray(instance, solution)

    #it is maybe better to just do one loop and calculate these simultaneously
    stats[:QntLines] = length(instance.lines)
    stats[:QntStops] = nv(instance.network)
    stats[:QntEdges] = ne(instance.network)
    stats[:TakenLines] = count(i -> (i > 0), solution)
    stats[:AverageFrequency] = mean(edge_frequencies)
    stats[:StandardDeviation] = std(edge_frequencies, mean = stats[:AverageFrequency])
    stats[:Variance] = var(edge_frequencies, mean = stats[:AverageFrequency])
    stats[:Median] = median(edge_frequencies)
    stats[:LowestFrequency] = minimum(edge_frequencies)
    stats[:HighestFrequency] = maximum(edge_frequencies)
    stats[:LineCost] = CalculateCost(instance, solution)

    return stats

end

function CalculateFrequencyArray(instance::ProblemInstance, solution::Array{Bool,1})
    return CalculateFrequencyArray(instance, [if solution[i] instance.lines[i].frequency else 0 end] for i in 1:length(solution))
end

function CalculateFrequencyArray(instance::ProblemInstance, solution::Array{T,1} where {T<:Integer})
    edge_frequencies = Array{Int64, 1}()

    for e in edges(instance.network)
        push!(edge_frequencies, sum( solution[i] * IsEdgeOnLine(instance, e, i) for i in 1:length(solution)))
    end

    return edge_frequencies

end


text:
# Ref [Parisi-Rapuano '85](http://www.sciencedirect.com/science/article/pii/0370269385906707)
mutable struct ParisiRapuano <: AbstractRNG
    myrand::UInt32
    ira::Vector{UInt32}
    ip::UInt8
    ip1::UInt8
    ip2::UInt8
    ip3::UInt8
end

touint32(x::UInt64) = UInt32(x & 0x00000000ffffffff)

function ParisiRapuano(seed::Integer)
    myrand = convert(UInt32, seed)
    ip = 128
    ip1 = ip - 24
    ip2 = ip - 55
    ip3 = ip - 65
    ip3 =1
    ira = Vector{UInt32}(256)

    for i=ip3:ip-1
        y = myrand * UInt64(16807)
        myrand = touint32(xor(y, 0x7fffffff) + (y >> 31))
        if xor(myrand, 0x80000000) != 0
            myrand = (myrand & 0x7fffffff) + one(UInt32)
        end
        ira[i+1] = myrand
    end
    return ParisiRapuano(myrand, ira, ip, ip1, ip2, ip3)
end

ParisiRapuano() = ParisiRapuano(make_seed()) # should be greater then 0

function make_seed()
    seed = UInt64(0)
    while seed == 0
        seed = rand(RandomDevice(), UInt32)
    end
    return seed
end

srand(r::ParisiRapuano, seed) = deepcopy!(r, ParisiRapuano(seed))

function deepcopy!(dest::T, source::T) where T
    for name in fieldnames(T)
        setfield!(dest, name, deepcopy(getfield(source, name)))
    end
    return dest
end

function rand(r::ParisiRapuano, ::Type{Base.Random.CloseOpen})
    # assert( 1<= r.ip+1 <= 256)
	# assert( 1<= r.ip1+1 <= 256)
	# assert( 1<= r.ip2+1 <= 256)
	# assert( 1<= r.ip3+1 <= 256)
    uno = UInt8(1)
    x = xor(r.ira[r.ip1+1] + r.ira[r.ip2+1], r.ira[r.ip3+1])
    r.ira[r.ip+1] = x
    r.ip = (r.ip + uno) % UInt8;
    r.ip1 = (r.ip1 + uno) % UInt8;
    r.ip2 = (r.ip2 + uno) % UInt8;
    r.ip3 = (r.ip3 + uno) % UInt8
    return 2.3283064365e-10 * x
end

rand(r::ParisiRapuano, ::Type{Base.Random.Close1Open2}) = rand(r, Base.Random.CloseOpen) + 1.
rand(r::ParisiRapuano, ::Type{Float64}) = rand(r, Base.Random.CloseOpen)

max_stars_repo_name: JuliaTagBot/LittleScienceTools

max_issues_repo_name: JuliaTagBot/LittleScienceTools

max_forks_repo_name: JuliaTagBot/LittleScienceTools


text:
module PkgButlerEngine

import Mustache
import Pkg

function configure_pkg(path::AbstractString; channel=:auto)
    channel in (:auto, :stable, :dev) || error("Invalid value for channel.")

    path_for_butler_workflows_folder = joinpath(path, ".github", "workflows")

    path_for_main_butler_workflow = joinpath(path_for_butler_workflows_folder, "jlpkgbutler-butler-workflow.yml")
    path_for_main_butler_dev_workflow = joinpath(path_for_butler_workflows_folder, "jlpkgbutler-butler-dev-workflow.yml")

    mkpath(path_for_butler_workflows_folder)

    if channel==:auto
        channel = isfile(path_for_main_butler_dev_workflow) ? :dev : :stable
    end

    if channel==:stable
        if isfile(path_for_main_butler_dev_workflow)
            rm(path_for_main_butler_dev_workflow, force=true)
        end

        cp(joinpath(@__DIR__, "..", "templates", "jlpkgbutler-butler-workflow.yml"), path_for_main_butler_workflow, force=true)
    elseif channel==:dev
        if isfile(path_for_main_butler_workflow)
            rm(path_for_main_butler_workflow, force=true)
        end

        cp(joinpath(@__DIR__, "..", "templates", "jlpkgbutler-butler-dev-workflow.yml"), path_for_main_butler_dev_workflow, force=true)
    end
end

function cp_with_mustache(src, dest, vals)
    content = read(src, String)

    open(dest, "w") do file
        Mustache.render(file, content, vals, tags= ("\$[[", "]]"))
    end
end

function ensure_project_has_julia_compat(path)
    proj_file = isfile(joinpath(path, "JuliaProject.toml")) ? joinpath(path, "JuliaProject.toml") : joinpath(path, "Project.toml")

    pkg_toml_content = Pkg.TOML.parsefile(proj_file)

    if !haskey(pkg_toml_content, "compat")
        pkg_toml_content["compat"] = Dict{String,String}()
    end

    if !haskey(pkg_toml_content["compat"], "julia")
        pkg_toml_content["compat"]["julia"] = "1"

        open(proj_file, "w") do f
            Pkg.TOML.print(f, pkg_toml_content)
        end
    end
end

function ensure_project_uses_new_enough_documenter(path)
    doc_proj_file = isfile(joinpath(path, "docs", "JuliaProject.toml")) ? joinpath(path, "docs", "JuliaProject.toml") : joinpath(path, "docs", "Project.toml")

    if isfile(doc_proj_file)
        pkg_toml_content = Pkg.TOML.parsefile(doc_proj_file)

        if haskey(pkg_toml_content, "compat") && haskey(pkg_toml_content["compat"], "Documenter")
            documenter_compat_bound = pkg_toml_content["compat"]["Documenter"]
            version_bound = Pkg.Types.semver_spec(documenter_compat_bound)

            # This is the list of versions that don't work
            invalid_versions = Pkg.Types.semver_spec("0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,0.10,0.11,0.12,0.13,0.14,0.15,0.16,0.17,0.18,0.19,0.20,0.21,0.22,0.23")

            if !isempty(intersect(invalid_versions, version_bound))
                # TODO This is a bit crude. Ideally we would try to see whether a version >= 0.24 is listed in the compat section, and if so,
                # only remove the offending versions.
                pkg_toml_content["compat"]["Documenter"] = "~0.24"

                open(doc_proj_file, "w") do f
                    Pkg.TOML.print(f, pkg_toml_content)
                end

            end
        end

    end
end

function construct_version_matrix(path)
    proj_file = isfile(joinpath(path, "JuliaProject.toml")) ? joinpath(path, "JuliaProject.toml") : joinpath(path, "Project.toml")

    pkg_toml_content = Pkg.TOML.parsefile(proj_file)

    julia_compat_bound = pkg_toml_content["compat"]["julia"]

    version_spec = Pkg.Types.semver_spec(julia_compat_bound)

    versions = [v"1.0.5", v"1.1.1", v"1.2.0"]

    compat_versions = filter(i->i in version_spec, versions)

    return join(string.(compat_versions), ", ")
end

function update_pkg(path::AbstractString)
    configure_pkg(path)

    path_for_butler_workflows_folder = joinpath(path, ".github", "workflows")

    path_for_ci_master_workflow = joinpath(path_for_butler_workflows_folder, "jlpkgbutler-ci-master-workflow.yml")
    path_for_ci_pr_workflow = joinpath(path_for_butler_workflows_folder, "jlpkgbutler-ci-pr-workflow.yml")
    path_for_docdeploy_workflow = joinpath(path_for_butler_workflows_folder, "jlpkgbutler-docdeploy-workflow.yml")

    path_for_docs_make_file = joinpath(path, "docs", "make.jl")

    mkpath(path_for_butler_workflows_folder)

    ensure_project_has_julia_compat(path)

    view_vals = Dict{String, Any}()
    view_vals["JL_VERSION_MATRIX"] = construct_version_matrix(path)

    cp_with_mustache(joinpath(@__DIR__, "..", "templates", "jlpkgbutler-ci-master-workflow.yml"), path_for_ci_master_workflow, view_vals)
    cp_with_mustache(joinpath(@__DIR__, "..", "templates", "jlpkgbutler-ci-pr-workflow.yml"), path_for_ci_pr_workflow, view_vals)

    if isfile(path_for_docs_make_file)
        cp(joinpath(@__DIR__, "..", "templates", "jlpkgbutler-docdeploy-workflow.yml"), path_for_docdeploy_workflow, force=true)
    else isfile(path_for_docdeploy_workflow)
        rm(path_for_docdeploy_workflow, force=true)
    end

    ensure_project_uses_new_enough_documenter(path)
end

end # module


text:
using MultithreadingBenchmarks
using Printf
using InteractiveUtils

# ---------------------------------------------------------------
println("------ JULIA VERSIONINFO ------")
InteractiveUtils.versioninfo()

println("------ JULIA CPU INFO ------")
@show Sys.CPU_THREADS
Sys.cpu_summary()

# ---------------------------------------------------------------

const BENCHDIR = joinpath(dirname(pathof(MultithreadingBenchmarks)), "..", "bench")

const NUM_DATAPOINTS = 11  # Keep this odd, so we get a point at 50%

MultithreadingBenchmarks.perform_scaling_experiment(
    bench_file = "$BENCHDIR/simple_independent.jl",
    num_datapoints = NUM_DATAPOINTS,  # keep this odd, so we get a point at 50%
    nqueries = 1000,
    num_ops = 1_000_000_00,#0,
    plot_series_name = "1000 queries x 1e8 ops",
    )

MultithreadingBenchmarks.perform_scaling_experiment(
    bench_file = "$BENCHDIR/all_tasks_allocating.jl",
    num_datapoints = 5,  # Use fewer datapoints since it's expensive and clear
    nqueries = 1000,
    num_ops = 1_000_000,
    plot_series_name = "all tasks alloc garbage: 1000 queries x 1e6 ops",
    )

MultithreadingBenchmarks.perform_scaling_experiment(
     bench_file = "$BENCHDIR/all_tasks_compiling.jl",
     num_datapoints = 5,  # Use fewer datapoints since it's expensive and clear
     nqueries = 5,
     num_ops = 1_0,
     plot_series_name = "all tasks alloc garbage: 5 queries x 1e1 ops",
     )

# TODO: This benchmark is still very experimental
MultithreadingBenchmarks.perform_scaling_experiment(
    bench_file = "$BENCHDIR/onethread_alloc_garbage.jl",
    num_datapoints = NUM_DATAPOINTS,  # Unused -- we're setting it manually below:
    # Manually set the nthreads to test -- Skip nthreads=1 because one thread is occupied
    # allocating garbage and won't participate in the test. The test would deadlike for
    # nthreads=1, so we use nthreads=2 instead. Effectively, all these values are 1 higher
    # than are actually being tested, so we effectively test from [1, NUM_CPUS-1].
    # (NOTE: that we then also need to subtract 1 from all the nthreads in the results)
    nthreads_to_test = [2, Int.(round.(range(0, stop=MultithreadingBenchmarks.NUM_CORES, length=NUM_DATAPOINTS)))[2:end]...],
    # Subtract 1 from nthreads in result since 1 CPU is occupied allocating garbage.
    preprocess_results = results->(results.nthreads .-= 1; results),

    nqueries = 1000,
    num_ops = 1_000_000_00,#0,
    plot_series_name = "one task allocs garbage: 1000 queries x 1e8 ops",
    )


text:
using Oceananigans.TimeSteppers: ab2_step_field!

import Oceananigans.TimeSteppers: ab2_step!

combine_events(free_surface_event,
               velocities_events,
               tracer_events) = MultiEvent(tuple(free_surface_event, velocities_events..., tracer_events...))

combine_events(::Nothing, velocities_events, tracer_events) =
    MultiEvent(tuple(velocities_events..., tracer_events...))

combine_events(::Nothing, ::Nothing, tracer_events) = MultiEvent(tuple(tracer_events...))

function ab2_step!(model::HydrostaticFreeSurfaceModel, Δt, χ)

    workgroup, worksize = work_layout(model.grid, :xyz)

    barrier = Event(device(model.architecture))

    # Launch velocity update kernels

    velocities_events = []

    for name in (:u, :v)
        model.velocities isa PrescribedVelocityFields && break

        Gⁿ = model.timestepper.Gⁿ[name]
        G⁻ = model.timestepper.G⁻[name]
        velocity_field = model.velocities[name]

        event = launch!(model.architecture, model.grid, :xyz, ab2_step_field!,
                        velocity_field, Δt, χ, Gⁿ, G⁻,
                        dependencies=barrier)

        push!(velocities_events, event)
    end

    # Tracer update kernels

    tracer_events = []

    for name in tracernames(model.tracers)
        Gⁿ = model.timestepper.Gⁿ[name]
        G⁻ = model.timestepper.G⁻[name]
        tracer_field = model.tracers[name]

        event = launch!(model.architecture, model.grid, :xyz, ab2_step_field!,
                        tracer_field, Δt, χ, Gⁿ, G⁻,
                        dependencies=barrier)

        push!(tracer_events, event)
    end

    velocities_update = MultiEvent(Tuple(velocities_events))

    # Free surface update

    free_surface_event = ab2_step_free_surface!(model.free_surface, velocities_update, model, Δt, χ)

    prognostic_field_events = combine_events(free_surface_event, velocities_events, tracer_events)

    wait(device(model.architecture), prognostic_field_events)

    return nothing
end


text:
bonds(lattice::AbstractLattice, k::Int) = ((s, n) for s in sites(lattice) for n in neighbors(lattice, s, k))
bonds(lattice::AbstractLattice) = bonds(lattice, 1)


text:
a = BigInt(20)
f = factorial
c = Int(f(2a) / (f(a) * f(a)))
println(c)

max_stars_repo_name: sujimodern/project-euler-jl

max_issues_repo_name: sujimodern/project-euler-jl

max_forks_repo_name: sujimodern/project-euler-jl


text:
const GITLAB_GROUP = "UnofficialJuliaMirror"
const GITLAB_BOT_USERNAME = "UnofficialJuliaMirrorBot"
const GITLAB_BOT_PERSONAL_ACCESS_TOKEN = ENV["GITLAB_BOT_PERSONAL_ACCESS_TOKEN"]


text:
using GridapFSI
using Test

@testset "GridapFSI.jl" begin
  @testset "FSIDrivers.jl" begin include("FSIDriversTest.jl") end

  @testset "Analytical.jl" begin
    main(problemName="analytical",is_vtk=true)
    main(problemName="analytical",strategy="linearElasticity")
    #main(problemName="analytical",strategy="neoHookean")
    end
  @testset "ElasticFlag.jl" begin
    main()
    main(coupling="weak")
    main(strategy="biharmonic")
    main(coupling="weak",strategy="biharmonic",is_vtk=true,alpha_m_weight="volume")
  end
end


text:
struct NN
    g
end

function NN(nchars::Int, ntags::Int)
    T = Float32
    x = Node()
    embeds = embeddings(T, nchars, 20, init_w=Normal(0,0.05))
    h = lookup(embeds, x)
    for i = 1:2
        h = Conv1D(T,5,20,20,2,1)(h)
        h = relu(h)
    end
    h = Linear(T,20,ntags)(h)
    Graph(input=x, output=h)
end

function (nn::NN)(x, y=nothing)
    if y == nothing
        z = nn.g(x, false)
        argmax(z.data, 1)
    else
        z = nn.g(x, true)
        softmax_crossentropy(y, z)
    end
end


text:
"""
    HorizonLength(::Type{<:Union{MDP,POMDP})
    HorizonLength(::Union{MDP,POMDP})

Check whether MDP is Finite or Infinite Horizon and return corresponding struct (FiniteHorizon or InfiniteHorizon).
"""
abstract type HorizonLength end

"If HorizonLength(m::Union{MDP,POMDP}) == FiniteHorizon(), horizon(m) should be implemented and return an integer"
struct FiniteHorizon <: HorizonLength end
"HorizonLength(m::Union{MDP,POMDP}) == InfiniteHorizon() indicates that horizon(m) should not be called."
struct InfiniteHorizon <: HorizonLength end

HorizonLength(m::Union{MDP,POMDP}) = HorizonLength(typeof(m))
HorizonLength(::Type{<:Union{MDP,POMDP}}) = InfiniteHorizon()

"""
Return the number of *steps* that will be taken in the (PO)MDP, given it is Finite Horizon.

A simulation of a (PO)MDP with `horizon(m) == d` should contain *d+1* states and *d* actions and rewards.
"""
function horizon end

"""
    stage(m::Union{MDP,POMDP}, ss)::Int
    stage(m::Union{MDP,POMDP}, o)::Int
    stage(d)::Int

Considering a variable or distribution containing its stage assignment, return the number of its stage.
"""
function stage end

"""
    stage_states(m::Union{MDP,POMDP}, stage::Int)

Create (PO)MDP's states for given stage.
"""
function stage_states end

"""
    stage_stateindex(m::Union{MDP,POMDP}, ss}::Int

Compute the index of the given state in the corresponding stage.
"""
function stage_stateindex end

"""
    ordered_stage_states(w::FHWrapper, stage::Int)

Return an AbstractVector of states from given stage ordered according to stage_stateindex(mdp, s).
"""
function ordered_stage_states end

"""
    stage_observations(m::Union{MDP,POMDP}, stage::Int)

Create (PO)MDP's observations for given stage.
"""
function stage_observations end

"""
    stage_obsindex(m::Union{MDP,POMDP}, o::stage::Int)

Compute the index of the given observation in the corresponding stage.
"""
function stage_obsindex end

"""
    ordered_stage_observations(w::FHWrapper, stage::Int)

Return an AbstractVector of observations from given stage ordered according to stage_obsindex(w,o).
"""
function ordered_stage_observations end


text:
const marketOpenUS = TimeOfDay("9:30")
const marketCloseUS = TimeOfDay("16")

#############################################
# Exchange
#############################################

@enum Exchange NOEXCHANGE EDGX EDGA BZX BYX
exchange(s::AbstractString) = eval(symbol(s))
exchange(e::Exchange) = string(e)

const EXCH_LATENCY_NANOS = [Int(x * 1e3) for x in [199, 200, 100, 100]]
const EXCH_TAKE_FEES = [0.0035, 0.0, 0.0032, 0.0]
const EXCH_PROVIDE_FEES = [-0.0025, 0.0010, -0.0027, 0.0008]
getLatency(exch::Exchange) = (Int(exch) > 0 ? EXCH_LATENCY_NANOS[Int(exch)] : 0)
getTakeFee(exch::Exchange) = (Int(exch) > 0 ? EXCH_TAKE_FEES[Int(exch)] : 0.)
getProvideFee(exch::Exchange) = (Int(exch) > 0 ? EXCH_PROVIDE_FEES[Int(exch)] : 0.)
getFee(exch::Exchange, istake::Bool) = (istake ? getTakeFee : getProvideFee)(exch)
getPriceAdjustmentFromFee(exch::Exchange, buy::Bool, istake::Bool) = (buy ? 1 : -1) * getFee(exch, istake)

# Base.isless(e1::Exchange, e2::Exchange) = getTakeFee(e1) < getTakeFee(e2)
const EXCH_SORTED_BY_TAKE_FEE = sort([EDGX,EDGA,BZX,BYX], by=getTakeFee)

####################################################
# IDs
####################################################

typealias OID UInt64
typealias UID UInt64

iszero(i::UInt64) = i==0
# OID() = OID(0)
# UID() = UID(0)
# iszero(oid::OID) = oid == 0
# iszero(uid::UID) = uid == 0

type OID_Generator
  lastoid::OID
end
const OID_GENERATOR = OID_Generator(1e15)
function generateOID()
  OID_GENERATOR.lastoid += 1
  OID_GENERATOR.lastoid
end

#############################################
# Ticker
#############################################

immutable Ticker
  symbol::String

  Ticker(symbol::String) = new(strip(symbol))
  Ticker(byteArray::AbstractVector{UInt8}) = Ticker(String(byteArray))
  Ticker{T<:FixedLengthSymbol}(sym::T) = Ticker(sym.data)
end

Base.String(ticker::Ticker) = ticker.symbol
Base.hash(ticker::Ticker) = hash(ticker.symbol)
# @createIOMethods Ticker
Base.convert(::Type{Ticker}, s::String) = Ticker(s)
Base.show(io::IO, ticker::Ticker) = show(io, ticker.symbol)

for op in (:<, :>, :(==), :<=, :>=, :(Base.isless))
  @eval $op(t1::Ticker, t2::Ticker) = $op(t1.symbol, t2.symbol)
  @eval $op(t1::Ticker, t2::AbstractString) = $op(t1.symbol, t2)
  @eval $op(t1::AbstractString, t2::Ticker) = $op(t1, t2.symbol)
end


text:
P_referencia = 10^-3
P = 1
dB = 10*log10(abs(P/P_referencia))

max_stars_repo_name: PedroMoraesFranco/fisica-computacional

max_issues_repo_name: PedroMoraesFranco/fisica-computacional

max_forks_repo_name: PedroMoraesFranco/fisica-computacional


text:
################################################################################
# This set of functions provides various genetic utilities.
################################################################################
#
# Required external packages.
#
# using Statistics, Distributions

export map_function, inverse_map_function
export hardy_weinberg_test, xlinked_hardy_weinberg_test

"""
Calculate recombination fractions based on Haldane's or Kosambi's formula.
"""
function map_function(d::T where T <: Real, choice::String)
#
  one_d = one(typeof(d))
  if choice == "Haldane"
    theta = (one_d - exp(-2d)) / 2
  elseif choice == "Kosambi"
    theta = tanh(2d) / 2
  else
  throw(ArgumentError(
    "Only Haldane's or Kosambi's map function is allowed.\n \n"))
  end
  return theta
end # function map_function

"""
Calculates genetic map distances based on Haldane's or Kosambi's formula.
"""
function inverse_map_function(theta::T where T <: Real, choice::String)
#
  one_th = one(typeof(theta))
  if choice == "Haldane"
    if theta >= one_th / 2
      d = Inf
    else
      d = -(log(one_th - 2 * theta)) / 2
    end
  elseif choice == "Kosambi"
    if theta >= one_th / 2
      d = Inf
    else
      d = log((one_th + 2 * theta) / (one_th - 2 * theta)) / 4
    end
  else
    throw(ArgumentError(
      "Only Haldane's or Kosambi's map function is allowed.\n \n"))
  end
  return d
end # function inverse_map_function

"""
Test for Hardy-Weinberg equilibrium at a SNP.
The genotype vector conveys for each person his or her number of
reference alleles. Thus, all genotypes belong to {0, 1, 2, Missing}.
The pvalue of the chi-square statistic is returned.
"""
function hardy_weinberg_test(genotype::Vector)
#
  (expected, observed) = (zeros(3), zeros(3))
  #
  # Tally observations per cell.
  #
  n = 0
  for i = 1:length(genotype)
    if !ismissing(genotype[i])
      n = n + 1
      g = round(Int, genotype[i])
      j = g + 1
      observed[j] = observed[j] + 1.0
    end
  end
  if n == 0
    return 1.0
  end
  #
  # Compute allele frequencies.
  #
  p = (2.0 * observed[1] + observed[2]) / (2 * n)
  if p <= 0.0 || p >= 1.0
    return 1.0
  else
    q = 1.0 - p
    #
    # Sum over cells of the contingency table.
    #
    expected[1] = n * p^2
    expected[2] = n * 2.0 * p * q
    expected[3] = n * q^2
    test = 0.0
    for i = 1:3
    test = test + (observed[i] - expected[i])^2 / expected[i]
  end
  return ccdf(Chisq(1), test)
  end
end # function hardy_weinberg_test

"""
Test for Hardy-Weinberg equilibrium at an xlinked SNP.
The genotype vector conveys for each person his or her number of
reference alleles. Thus, all genotypes belong to {0, 1, 2, Missing}.
The p-value of the chi-square statistic is returned.
"""
function xlinked_hardy_weinberg_test(genotype::Vector,
  male::BitArray{1})

  (expected_female, observed_female) = (zeros(3), zeros(3))
  (expected_male, observed_male) = (zeros(2), zeros(2))
  #
  # Tally observations per cell and count members of each sex.
  #
  (females, males) = (0, 0)
  for i = 1:length(genotype)
    if !ismissing(genotype[i])
      g = round(Int, genotype[i])
      if male[i]
        males = males + 1
        j = min(g + 1, 2)
        observed_male[j] = observed_male[j] + 1.0
      else
        females = females + 1
        j = g + 1
        observed_female[j] = observed_female[j] + 1.0
      end
    end
  end
  if females + males == 0
    return 1.0
  end
  #
  # Compute allele frequencies.
  #
  p = 2.0 * observed_female[1] + observed_female[2] + observed_male[1]
  p = p / (2 * females + males)
  if p <= 0.0 || p >= 1.0
    return 1.0
  end
  q = 1.0 - p
  #
  # Sum over female cells of the contingency table.
  #
  expected_female[1] = females * p^2
  expected_female[2] = females * 2.0 * p * q
  expected_female[3] = females * q^2
  test = 0.0
  for i = 1:3
    term = (observed_female[i] - expected_female[i])^2 / expected_female[i]
    test = test + term
  end
  #
  # Sum over male cells of the contingency table.
  #
  expected_male[1] = males * p
  expected_male[2] = males * q
  for i = 1:2
    term = (observed_male[i] - expected_male[i])^2 / expected_male[i]
    test = test + term
  end
  return ccdf(Chisq(2), test)
end # function xlinked_hardy_weinberg_test

# Testing Routines
#
# d = rand()
# choice = "Haldane";
# theta = map_function(d, choice)
# di = inverse_map_function(theta, choice);
# println(di);
# d = rand()
# choice = "Kosambi";
# theta = map_function(d, choice)
# di= inverse_map_function(theta, choice);
# println(di);
# p = 1 / 3;
# n = 100;
# x = zeros(n);
# for i = 1:n
#   u = rand(2)
#   if u[1] <= p; x[i] = x[i] + 1.0; end
#   if u[2] <= p; x[i] = x[i] + 1.0; end
# end
# pvalue = hardy_weinberg_test(x);
# println("pvalue = ",pvalue);
# male = falses(n);
# male[1:div(n, 2)] .= true;
# for i = 1:div(n, 2)
#   x[i] = 0.0
#   u = rand(1)
#   if u[1] <= p; x[i] = x[i] + 1.0; end
# end
# pvalue = xlinked_hardy_weinberg_test(x, male);
# println("pvalue = ",pvalue);

max_stars_repo_name: OpenMendel/MendelBase

max_issues_repo_name: OpenMendel/MendelBase

max_forks_repo_name: OpenMendel/MendelBase


text:
"The average reflection coefficients"
function reflection_coefficients(ωs::Union{T,AbstractVector{T}}, psource::PlaneSource{T,Dim}, material::Material{Dim,Halfspace{T,Dim}}; kws...) where {T<:Number,Dim}

    Rs = map(ωs) do ω
        k_effs = wavenumbers(ω, psource.medium, material.species; kws...)
        w_effs = [
            WaveMode(ω, k_eff, psource, material; kws...)
        for k_eff in k_effs]

        reflection_coefficient(ω, w_effs, psource, material; kws...)
    end

    return Rs
end

"Calculates the reflection coefficient from each wave in waves and then sums the results."
reflection_coefficient(ω::T, waves::Vector{EffectivePlaneWaveMode{T,Dim}}, psource::PlaneSource{T,Dim}, material::Material{Dim,Halfspace{T,Dim}}; kws...) where {T<:Number,Dim} =
    sum(w -> reflection_coefficient(ω, w, psource, material; kws...), waves)

"Pairs each ω in ωs with each wave in waves to calculate the refleciton coefficients: reflection_coefficient(ω, wave)"
reflection_coefficients(ωs::AbstractVector{T}, waves::Vector{EffectivePlaneWaveMode{T,Dim}}, psource::PlaneSource{T,Dim}, material::Material{Dim,Halfspace{T,Dim}}; kws...) where {T<:Number,Dim} =
[ reflection_coefficient(ωs[i], waves[i], psource, material; kws...) for i in eachindex(ωs)]

max_stars_repo_name: arturgower/EffectiveWaves


text:
"""
    Estimator

Abstract base type for all estimators.
"""
abstract type Estimator end

# these two enable array broadcasting of all methods!
Base.size(::Estimator) = ()

Base.getindex(estimator::Estimator, i) = estimator

design(estimator::Estimator) = try return estimator.design catch error("not implemented") end

"""
    estimate{T<:Integer}(estimator::Estimator, x1::T, x2::T)

Estimate the response rate from observed x1 and x2.

# Parameters

| Parameter    | Description |
| -----------: | :---------- |
| estimator    | any BinaryTwoStageDeisgnEstimator object |
| x1           | stage-one responses |
| x2           | stage-two responses |

# Return Value

Real, estimated response rate.

# Examples
```julia-repl
julia> ss = SimpleSampleSpace(10:25, 100, n2min = 5)
julia> interimsamplesize(ss)
julia> design = getoptimaldesign(15, params, solver = Gurobi.GurobiSolver())
julia> est = MLE(design)
julia> estimate(est, 0, 0)
```
"""
estimate(estimator::Estimator, x1::T, x2::T) where {T<:Integer} = error("not implemented")

"""
    p{T1<:Integer, T2<:Real}(estimator::Estimator, x1::T1, x2::T1, p0::T2)

Compute the p value after observing (x1, x2) for null hypothesis H0: p <= p0 with
respect to ordering induced by `estimator`.

# Parameters

| Parameter    | Description |
| -----------: | :---------- |
| estimator    | any BinaryTwoStageDeisgnEstimator object |
| x1           | stage-one responses |
| x2           | stage-two responses |
| p0           | upper boundary of null hypothesis |

# Return Value

Real, p value.

# Examples
```julia-repl
julia> ss = SimpleSampleSpace(10:25, 100, n2min = 5)
julia> interimsamplesize(ss)
julia> design = getoptimaldesign(15, params, solver = Gurobi.GurobiSolver())
julia> est = MLE(design)
julia> p(est, 0, 0, .2)
```
"""
function pvalue(
  estimator::Estimator, x1::T1, x2::T1, p0::T2
) where {T1<:Integer, T2<:Real}

  checkx1x2(x1, x2, design(estimator))
  supp          = support(design(estimator))
  estimates     = estimate.(estimator, supp[:, 1], supp[:, 2])
  estimateobs   = estimate(estimator, x1, x2)
  ismoreextreme = estimates .>= estimateobs
  return pdf.(design(estimator), supp[ismoreextreme, 1], supp[ismoreextreme, 2], p0) |> sum

end

"""
    bias{T<:Real}(estimator::Estimator, p::T)

Bias of `estimator` given response rate `p`.

# Parameters

| Parameter    | Description |
| -----------: | :---------- |
| estimator    | any BinaryTwoStageDeisgnEstimator object |
| p0           | upper boundary of null hypothesis |

# Return Value

Real, bias given `p`.

# Examples
```julia-repl
julia> ss = SimpleSampleSpace(10:25, 100, n2min = 5)
julia> interimsamplesize(ss)
julia> design = getoptimaldesign(15, params, solver = Gurobi.GurobiSolver())
julia> est = MLE(design)
julia> bias(est, .3)
```
"""
function bias(estimator::Estimator, p::T) where {T<:Real}

  supp      = support(design(estimator))
  estimates = estimate.(estimator, supp[:, 1], supp[:, 2])
  probs     = pdf.(design(estimator), supp[:, 1], supp[:, 2], p)
  return dot(probs, estimates) - p

end

"""
    rmse{T<:Real}(estimator::Estimator, p::T)

Root mean squared error of `estimator` given response rate `p`.

# Parameters

| Parameter    | Description |
| -----------: | :---------- |
| estimator    | any BinaryTwoStageDeisgnEstimator object |
| p0           | upper boundary of null hypothesis |

# Return Value

Real, RMSE given `p`.

# Examples
```julia-repl
julia> ss = SimpleSampleSpace(10:25, 100, n2min = 5)
julia> interimsamplesize(ss)
julia> design = getoptimaldesign(15, params, solver = Gurobi.GurobiSolver())
julia> est = MLE(design)
julia> rmse(est, .3)
```
"""
function rmse(estimator::Estimator, p::T) where {T<:Real}

  supp     = support(design(estimator))
  return sqrt( sum(
    pdf.(design(estimator), supp[:, 1], supp[:, 2], p) .*
      (estimate.(estimator, supp[:, 1], supp[:, 2]) .- p).^2
  ) )

end

function mae(estimator::Estimator, p::T) where {T<:Real}

    supp     = support(design(estimator))
    return sum(
      pdf.(design(estimator), supp[:, 1], supp[:, 2], p) .*
        abs.(estimate.(estimator, supp[:, 1], supp[:, 2]) .- p)
    )

  end

"""
    incompatibleoutcomes(estimator::Estimator)

Find outcomes where the induced p value implies different decision than the
underlying design

# Parameters

| Parameter    | Description |
| -----------: | :---------- |
| estimator    | any BinaryTwoStageDeisgnEstimator object |

# Return Value

Array with respective outcomes

# Examples
```julia-repl
julia> ss = SimpleSampleSpace(10:25, 100, n2min = 5)
julia> interimsamplesize(ss)
julia> design = getoptimaldesign(15, params, solver = Gurobi.GurobiSolver())
julia> est = MLE(design)
julia> incompatibleoutcomes(est)
```
"""
function incompatibleoutcomes(estimator::Estimator)

    supp = support(design(estimator))
    res  = []
    p0   = estimator |> design |> parameters |> null
    a    = estimator |> design |> parameters |> mtoer
    for i in 1:size(supp, 1)
        x1, x2 = supp[i, :]
        if (pvalue(estimator, x1, x2, p0) <= a) != test(design(estimator), x1, x2)
            push!(res, supp[i, :])
        end
    end
    return hcat(res...)'

end

max_stars_repo_name: kkmann/binaryTwoStageDesigns

max_issues_repo_name: kkmann/binaryTwoStageDesigns

max_forks_repo_name: kkmann/binaryTwoStageDesigns


text:
# Test for Robust Optimization
# Version: 1.0
# Author: Edward J. Xu, edxu96@outlook.com
# Date: March 31th, 2019
# 0. Set working path, install packages and use module -----------------------------------------------------------------
push!(LOAD_PATH, "$(homedir())/Desktop/Production Planning, Robust Optimization, DTU02435")
cd("$(homedir())/Desktop/Production Planning, Robust Optimization, DTU02435")
using ExcelReaders
using JuMP
using GLPKMathProgInterface
# 1. Data Process ------------------------------------------------------------------------------------------------------
file = openxl("data.xlsx")
vec_timeProdDev = readxl(file, "timeProdDev!B2:B11")
vec_costInvest = readxl(file, "costInvest!B2:B5")
vec_costProd = readxl(file, "costProd!B2:B11")
vec_demand = readxl(file, "demand!B2:B11")
mat_prodOnMachine = readxl(file, "prodOnMachine!B2:E11")
vec_timeAvai = readxl(file, "timeAvai!B2:B5")
vec_timeProd = readxl(file, "timeProd!B2:B11")
vec_gammaCap = zeros(4)
for m = 1:4
    vec_gammaCap[m] = 0.0 * sum([mat_prodOnMachine[p, m]] for p = 1: 10)[1]
end
# 2. Robust Optim ------------------------------------------------------------------------------------------------------
J = 300
model = Model(solver = GLPKSolverMIP())
@variable(model, vec_y[1: 4] >= 0, Int)
@variable(model, mat_x[1: 4, 1: 10] >= 0)
@variable(model, beta)
@variable(model, vec_P[1: 10])
@variable(model, vec_Q[1: 10])
@variable(model, vec_theta1[1: 10])
@variable(model, vec_theta2[1: 10])
@objective(model, Min, sum([vec_y[m] * vec_costInvest[m]] for m = 1: 4)[1] +
    sum([mat_x[m, p] * vec_costProd[p]] for m = 1: 4, p = 1: 10)[1] + beta)
# Constraints
@constraint(model, beta >= J * sum(vec_P) + J * sum(vec_theta1))
@constraint(model, [p = 1: 10], - vec_Q[p] <= vec_theta1[p])
@constraint(model, [p = 1: 10], vec_theta1[p] <= vec_Q[p])

@constraint(model, [p = 1: 10], sum([mat_x[m, p]] for m = 1: 4)[1] + (vec_P[p] + vec_theta1[p]) >= vec_demand[p])
@constraint(model, [p = 1: 10], vec_P[p] + vec_theta1[p] >= 0)

@constraint(model, [m = 1: 4], sum([vec_timeProd[p] * mat_x[m, p] + vec_theta2[p]] for p = 1: 10)[1] <= vec_timeAvai[m] * vec_y[m])
@constraint(model, [m = 1: 4, p = 1: 10], - vec_timeProdDev[p] * mat_x[m, p] <= vec_theta2[p])
@constraint(model, [m = 1: 4, p = 1: 10], vec_theta2[p] <= vec_timeProdDev[p] * mat_x[m, p])
@constraint(model, [m = 1: 4, p = 1: 10], mat_x[m, p] <= mat_prodOnMachine[p, m] * vec_y[m] * 1000000)
solve(model)
obj = getobjectivevalue(model)
vec_result_y = getvalue(vec_y)
mat_result_x = getvalue(mat_x)

# 2. Optimization ------------------------------------------------------------------------------------------------------
model = Model(solver = GLPKSolverMIP())
@variable(model, vec_y[1: 4] >= 0, Int)
@variable(model, mat_x[1: 4, 1: 10] >= 0)
@objective(model, Min, sum([vec_y[m] * vec_costInvest[m]] for m = 1: 4)[1] +
    sum([mat_x[m, p] * vec_costProd[p]] for m = 1: 4, p = 1: 10)[1])
@constraint(model, [p = 1: 10], sum([mat_x[m, p]] for m = 1: 4)[1] >= vec_demand[p])
@constraint(model, [m = 1: 4, p = 1: 10], mat_x[m, p] <= mat_prodOnMachine[p, m] * vec_y[m] * 1000000)
@variable(model, vec_lambda[1: 4] >= 0)
@variable(model, mat_mu[1: 4, 1: 10] >= 0)
# @variable(model, vec_z[1: 4] >= 0)
@constraint(model, [m = 1: 4], sum([vec_timeProd[p] * mat_x[m, p]] for p = 1: 10)[1] + vec_gammaCap[m] * vec_lambda[m] +
    sum([mat_mu[m, p]] for p = 1: 10)[1] <= vec_timeAvai[m] * vec_y[m])
@constraint(model, [m = 1: 4, p = 1: 10], vec_lambda[m] + mat_mu[m, p] >= vec_timeProdDev[p] * mat_x[m, p])  # vec_z[m]
# @constraint(model, [m = 1: 4, p = 1: 10], - vec_z[m] <= mat_x[m, p])
# @constraint(model, [m = 1: 4, p = 1: 10], mat_x[m, p] <= vec_z[m])
@constraint(model, vec_y .== [15, 0, 15, 3])
solve(model)
obj = getobjectivevalue(model)
vec_result_y = getvalue(vec_y)
mat_result_x = getvalue(mat_x)
vec_result_prod = zeros(10)
for i = 1: 10
    vec_result_prod[i] = sum(mat_result_x[:,i])
end


text:
"""
    ReshapeLayer(dims)

Reshapes the passed array to have a size of `(dims..., :)`

## Arguments

* `dims`: The new dimensions of the array (excluding the last dimension).

## Inputs

* `x`: AbstractArray of any shape which can be reshaped in `(dims..., size(x, ndims(x)))`

## Returns

* AbstractArray of size `(dims..., size(x, ndims(x)))`
* Empty `NamedTuple()`
"""
struct ReshapeLayer{N} <: AbstractExplicitLayer
    dims::NTuple{N, Int}
end

@inline function (r::ReshapeLayer)(x::AbstractArray, ps, st::NamedTuple)
    return reshape(x, r.dims..., size(x, ndims(x))), st
end

function Base.show(io::IO, r::ReshapeLayer)
    print(io, "ReshapeLayer(output_dims = (", join(r.dims, ", "), ", :))")
end

"""
    FlattenLayer()

Flattens the passed array into a matrix.

## Inputs

* `x`: AbstractArray

## Returns

* AbstractMatrix of size `(:, size(x, ndims(x)))`
* Empty `NamedTuple()`
"""
struct FlattenLayer <: AbstractExplicitLayer end

@inline function (f::FlattenLayer)(x::AbstractArray{T, N}, ps, st::NamedTuple) where {T, N}
    return reshape(x, :, size(x, N)), st
end

"""
    SelectDim(dim, i)

Return a view of all the data of the input `x` where the index for dimension `dim` equals `i`. Equivalent to `view(x,:,:,...,i,:,:,...)` where `i` is in position `d`.

## Arguments

* `dim`: Dimension for indexing
* `i`: Index for dimension `dim`

## Inputs

* `x`: AbstractArray that can be indexed with `view(x,:,:,...,i,:,:,...)`

## Returns

* `view(x,:,:,...,i,:,:,...)` where `i` is in position `d`
* Empty `NamedTuple()`
"""
struct SelectDim{dim, index} <: AbstractExplicitLayer end

SelectDim(dim, index) = SelectDim{Val(dim), Val(index)}()

@inline function (s::SelectDim{dim, index})(x, ps, st::NamedTuple) where {dim, index}
    selectdim(x, get_known(dim), get_known(index)), st
end

function Base.show(io::IO, s::SelectDim{dim, index}) where {dim, index}
    print(io, "SelectDim(dim = ", get_known(dim), ", index = ", get_known(index), ")")
end

"""
    NoOpLayer()

As the name suggests does nothing but allows pretty printing of layers. Whatever input is passed is returned.
"""
struct NoOpLayer <: AbstractExplicitLayer end

@inline (noop::NoOpLayer)(x, ps, st::NamedTuple) = x, st

"""
    WrappedFunction(f)

Wraps a stateless and parameter less function. Might be used when a function is added to `Chain`. For example, `Chain(x -> relu.(x))` would not work and the right thing to do would be `Chain((x, ps, st) -> (relu.(x), st))`. An easier thing to do would be `Chain(WrappedFunction(Base.Fix1(broadcast, relu)))`

## Arguments

* `f::Function`: A stateless and parameterless function

## Inputs

* `x`: s.t `hasmethod(f, (typeof(x),))` is `true`

## Returns

* Output of `f(x)`
* Empty `NamedTuple()`
"""
struct WrappedFunction{F} <: AbstractExplicitLayer
    func::F
end

(wf::WrappedFunction)(x, ps, st::NamedTuple) = wf.func(x), st

function Base.show(io::IO, w::WrappedFunction)
    return print(io, "WrappedFunction(", w.func, ")")
end

"""
    ActivationFunction(f)

Broadcast `f` on the input but fallback to CUDNN for Backward Pass. Internally calls [`Lux.applyactivation`](@ref)

## Arguments

* `f`: Activation function

## Inputs

* `x`: Any array type s.t. `f` can be broadcasted over it

## Returns

* Broadcasted Activation `f.(x)`
* Empty `NamedTuple()`
"""
struct ActivationFunction{F} <: AbstractExplicitLayer
    func::F
end

(af::ActivationFunction)(x, ps, st::NamedTuple) = applyactivation(af.func, x), st

function Base.show(io::IO, af::ActivationFunction)
    return print(io, "ActivationFunction(", af.func, ")")
end

"""
    SkipConnection(layer, connection)

Create a skip connection which consists of a layer or [`Chain`](@ref) of consecutive layers and a shortcut connection linking the block's input to the output through a user-supplied 2-argument callable. The first argument to the callable will be propagated through the given `layer` while the second is the unchanged, "skipped" input.

The simplest "ResNet"-type connection is just `SkipConnection(layer, +)`.

## Arguments

* `layer`: Layer or `Chain` of layers to be applied to the input
* `connection`: A 2-argument function that takes `layer(input)` and the input

## Inputs

* `x`: Will be passed directly to `layer`

## Returns

* Output of `connection(layer(input), input)`
* Updated state of `layer`

## Parameters

* Parameters of `layer`

## States

* States of `layer`

See [`Parallel`](@ref) for a more general implementation.
"""
struct SkipConnection{T <: AbstractExplicitLayer, F} <:
       AbstractExplicitContainerLayer{(:layers,)}
    layers::T
    connection::F
end

@inline function (skip::SkipConnection)(x, ps::Union{ComponentArray, NamedTuple},
                                        st::NamedTuple)
    mx, st = skip.layers(x, ps, st)
    return skip.connection(mx, x), st
end

"""
    Parallel(connection, layers...)

Create a layer which passes an input to each path in `layers`, before reducing the output with `connection`.

## Arguments

* `layers`: A list of `N` Lux layers
* `connection`: An `N`-argument function that is called after passing the input through each layer. If `connection = nothing`, we return a tuple `Parallel(nothing, f, g)(x, y) = (f(x), g(y))`

## Inputs

* `x`: if `x` is not a tuple, then return is computed as `connection([l(x) for l in layers]...)`. Else one is passed to each layer, thus `Parallel(+, f, g)(x, y) = f(x) + g(y)`.

## Returns

* See the Inputs section for how the output is computed
* Updated state of the `layers`

## Parameters

* Parameters of each `layer` wrapped in a NamedTuple with `fields = layer_1, layer_2, ..., layer_N`

## States

* States of each `layer` wrapped in a NamedTuple with `fields = layer_1, layer_2, ..., layer_N`

See also [`SkipConnection`](@ref) which is `Parallel` with one identity.
"""
struct Parallel{F, T <: NamedTuple} <: AbstractExplicitContainerLayer{(:layers,)}
    connection::F
    layers::T
end

function Parallel(connection, layers...)
    names = ntuple(i -> Symbol("layer_$i"), length(layers))
    return Parallel(connection, NamedTuple{names}(layers))
end

function (m::Parallel)(x, ps::Union{ComponentArray, NamedTuple}, st::NamedTuple)
    return applyparallel(m.layers, m.connection, x, ps, st)
end

@generated function applyparallel(layers::NamedTuple{names}, connection::C, x::T,
                                  ps::Union{ComponentArray, NamedTuple},
                                  st::NamedTuple) where {names, C, T}
    N = length(names)
    y_symbols = [gensym() for _ in 1:(N + 1)]
    st_symbols = [gensym() for _ in 1:N]
    getinput(i) = T <: Tuple ? :(x[$i]) : :x
    calls = []
    append!(calls,
            [:(($(y_symbols[i]), $(st_symbols[i])) = layers[$i]($(getinput(i)),
                                                                ps.$(names[i]),
                                                                st.$(names[i])))
             for
             i in 1:N])
    push!(calls, :(st = NamedTuple{$names}((($(Tuple(st_symbols)...),)))))
    if C == Nothing
        push!(calls, :($(y_symbols[N + 1]) = tuple($(Tuple(y_symbols[1:N])...))))
    else
        push!(calls, :($(y_symbols[N + 1]) = connection($(Tuple(y_symbols[1:N])...))))
    end
    push!(calls, :(return $(y_symbols[N + 1]), st))
    return Expr(:block, calls...)
end

Base.keys(m::Parallel) = Base.keys(getfield(m, :layers))

"""
    BranchLayer(layers...)

Takes an input `x` and passes it through all the `layers` and returns a tuple of the outputs.

## Arguments

* `layers`: A list of `N` Lux layers

## Inputs

* `x`: Will be directly passed to each of the `layers`

## Returns

* Tuple: `(layer_1(x), layer_2(x), ..., layer_N(x))`
* Updated state of the `layers`

## Parameters

* Parameters of each `layer` wrapped in a NamedTuple with `fields = layer_1, layer_2, ..., layer_N`

## States

* States of each `layer` wrapped in a NamedTuple with `fields = layer_1, layer_2, ..., layer_N`

## Comparison with [`Parallel`](@ref)

This is slightly different from `Parallel(nothing, layers...)`

* If the input is a tuple, `Parallel` will pass each element individually to each layer

* `BranchLayer` essentially assumes 1 input comes in and is branched out into `N` outputs

## Example

An easy way to replicate an input to an NTuple is to do

```julia
l = BranchLayer(
    NoOpLayer(),
    NoOpLayer(),
    NoOpLayer(),
)
```
"""
struct BranchLayer{T <: NamedTuple} <: AbstractExplicitContainerLayer{(:layers,)}
    layers::T
end

function BranchLayer(layers...)
    names = ntuple(i -> Symbol("layer_$i"), length(layers))
    return BranchLayer(NamedTuple{names}(layers))
end

function (m::BranchLayer)(x, ps::Union{ComponentArray, NamedTuple}, st::NamedTuple)
    applybranching(m.layers, x, ps, st)
end

@generated function applybranching(layers::NamedTuple{names}, x,
                                   ps::Union{ComponentArray, NamedTuple},
                                   st::NamedTuple) where {names}
    N = length(names)
    y_symbols = [gensym() for _ in 1:N]
    st_symbols = [gensym() for _ in 1:N]
    calls = []
    append!(calls,
            [:(($(y_symbols[i]), $(st_symbols[i])) = layers[$i](x, ps.$(names[i]),
                                                                st.$(names[i])))
             for i in 1:N])
    push!(calls, :(st = NamedTuple{$names}((($(Tuple(st_symbols)...),)))))
    push!(calls, :(return tuple($(Tuple(y_symbols)...)), st))
    return Expr(:block, calls...)
end

Base.keys(m::BranchLayer) = Base.keys(getfield(m, :layers))

"""
    PairwiseFusion(connection, layers...)

```
x1 --> layer1 --> y1
                  |
                  |--> connection --> layer2 --> y2
                  |                              |
                  x2                             |--> connection --> layer3 --> y3
                                                 |                              |
                                                 x3                             |--> connection --> y4
                                                                                |
                                                                                x4
```

## Arguments

* `connection`: Takes 2 inputs and combines them
* `layers`: [`AbstractExplicitLayer`](@ref)s

## Inputs

Layer behaves differently based on input type:
1. If the input `x` is a tuple of length `N + 1`, then the `layers` must be a tuple of length `N`. The computation is as follows

```julia
y = x[1]
for i in 1:N
    y = connection(x[i + 1], layers[i](y))
end
```

2. Any other kind of input

```julia
y = x
for i in 1:N
    y = connection(x, layers[i](y))
end
```

## Returns

* See Inputs section for how the return value is computed
* Updated model state for all the contained layers

## Parameters

* Parameters of each `layer` wrapped in a NamedTuple with `fields = layer_1, layer_2, ..., layer_N`

## States

* States of each `layer` wrapped in a NamedTuple with `fields = layer_1, layer_2, ..., layer_N`
"""
struct PairwiseFusion{F, T <: NamedTuple} <: AbstractExplicitContainerLayer{(:layers,)}
    connection::F
    layers::T
end

function PairwiseFusion(connection, layers...)
    names = ntuple(i -> Symbol("layer_$i"), length(layers))
    return PairwiseFusion(connection, NamedTuple{names}(layers))
end

function (m::PairwiseFusion)(x, ps::Union{ComponentArray, NamedTuple}, st::NamedTuple)
    return applypairwisefusion(m.layers, m.connection, x, ps, st)
end

@generated function applypairwisefusion(layers::NamedTuple{names}, connection::C, x::T,
                                        ps::Union{ComponentArray, NamedTuple},
                                        st::NamedTuple) where {names, C, T}
    N = length(names)
    y_symbols = [gensym() for _ in 1:(N + 1)]
    st_symbols = [gensym() for _ in 1:N]
    getinput(i) = T <: Tuple ? :(x[$i]) : :x
    calls = [:($(y_symbols[N + 1]) = $(getinput(1)))]
    append!(calls,
            [:(($(y_symbols[i]), $(st_symbols[i])) = layers[$i]($(y_symbols[N + 1]),
                                                                ps.$(names[i]),
                                                                st.$(names[i]));
               $(y_symbols[N + 1]) = connection($(y_symbols[i]), $(getinput(i + 1))))
             for i in 1:N])
    push!(calls, :(st = NamedTuple{$names}((($(Tuple(st_symbols)...),)))))
    push!(calls, :(return $(y_symbols[N + 1]), st))
    return Expr(:block, calls...)
end

Base.keys(m::PairwiseFusion) = Base.keys(getfield(m, :layers))

"""
    Chain(layers...; disable_optimizations::Bool = false)

Collects multiple layers / functions to be called in sequence on a given input.

## Arguments

* `layers`: A list of `N` Lux layers

## Keyword Arguments

* `disable_optimizations`: Prevents any structural optimization

## Inputs

Input `x` is passed sequentially to each layer, and must conform to the input requirements of the internal layers.

## Returns

* Output after sequentially applying all the layers to `x`
* Updated model states

## Parameters

* Parameters of each `layer` wrapped in a NamedTuple with `fields = layer_1, layer_2, ..., layer_N`

## States

* States of each `layer` wrapped in a NamedTuple with `fields = layer_1, layer_2, ..., layer_N`

## Optimizations

Performs a few optimizations to generate reasonable architectures. Can be disabled using keyword argument `disable_optimizations`.
* All sublayers are recursively optimized.
* If a function `f` is passed as a layer and it doesn't take 3 inputs, it is converted to a WrappedFunction(`f`) which takes only one input.
* If the layer is a Chain, it is flattened.
* [`NoOpLayer`](@ref)s are removed.
* If there is only 1 layer (left after optimizations), then it is returned without the `Chain` wrapper.
* If there are no layers (left after optimizations), a [`NoOpLayer`](@ref) is returned.

## Example

```julia
c = Chain(
    Dense(2, 3, relu),
    BatchNorm(3),
    Dense(3, 2)
)
```
"""
struct Chain{T} <: AbstractExplicitContainerLayer{(:layers,)}
    layers::T
    function Chain(xs...; disable_optimizations::Bool=false)
        xs = disable_optimizations ? xs : flatten_model(xs)
        length(xs) == 0 && return NoOpLayer()
        length(xs) == 1 && return first(xs)
        names = ntuple(i -> Symbol("layer_$i"), length(xs))
        layers = NamedTuple{names}(xs)
        return new{typeof(layers)}(layers)
    end
    function Chain(xs::AbstractVector; disable_optimizations::Bool=false)
        Chain(xs...; disable_optimizations)
    end
end

function flatten_model(layers::Union{AbstractVector, Tuple})
    new_layers = []
    for l in layers
        f = flatten_model(l)
        if f isa Tuple || f isa AbstractVector
            append!(new_layers, f)
        elseif f isa Function
            if !hasmethod(f, (Any, Union{ComponentArray, NamedTuple}, NamedTuple))
                if f === identity
                    continue
                else
                    push!(new_layers, WrappedFunction(f))
                end
            else
                push!(new_layers, f)
            end
        elseif f isa Chain
            append!(new_layers, f.layers)
        elseif f isa NoOpLayer
            continue
        else
            push!(new_layers, f)
        end
    end
    return layers isa AbstractVector ? new_layers : Tuple(new_layers)
end

flatten_model(x) = x

function (c::Chain)(x, ps::Union{ComponentArray, NamedTuple}, st::NamedTuple)
    applychain(c.layers, x, ps, st)
end

@generated function applychain(layers::NamedTuple{fields}, x,
                               ps::Union{ComponentArray, NamedTuple},
                               st::NamedTuple{fields}) where {fields}
    N = length(fields)
    x_symbols = [gensym() for _ in 1:N]
    st_symbols = [gensym() for _ in 1:N]
    calls = [:(($(x_symbols[1]), $(st_symbols[1])) = layers[1](x, ps.layer_1, st.layer_1))]
    append!(calls,
            [:(($(x_symbols[i]), $(st_symbols[i])) = layers[$i]($(x_symbols[i - 1]),
                                                                ps.$(fields[i]),
                                                                st.$(fields[i])))
             for i in 2:N])
    push!(calls, :(st = NamedTuple{$fields}((($(Tuple(st_symbols)...),)))))
    push!(calls, :(return $(x_symbols[N]), st))
    return Expr(:block, calls...)
end

Base.keys(m::Chain) = Base.keys(getfield(m, :layers))

"""
    Dense(in_dims => out_dims, activation=identity; init_weight=glorot_uniform, init_bias=zeros32, bias::Bool=true)

Create a traditional fully connected layer, whose forward pass is given by: `y = activation.(weight * x .+ bias)`

## Arguments

* `in_dims`: number of input dimensions
* `out_dims`: number of output dimensions
* `activation`: activation function

## Keyword Arguments

* `init_weight`: initializer for the weight matrix (`weight = init_weight(rng, out_dims, in_dims)`)
* `init_bias`: initializer for the bias vector (ignored if `bias=false`)
* `bias`: whether to include a bias vector

## Input

* `x` must be a Matrix of size `in_dims × B` or a Vector of length `in_dims`

## Returns

* Matrix of size `out_dims × B` or a Vector of length `out_dims`
* Empty `NamedTuple()`

## Parameters

* `weight`: Weight Matrix of size `out_dims × in_dims`
* `bias`: Bias of size `out_dims × 1` (present if `bias=true`)
"""
struct Dense{bias, F1, F2, F3} <: AbstractExplicitLayer
    activation::F1
    in_dims::Int
    out_dims::Int
    init_weight::F2
    init_bias::F3
end

function Base.show(io::IO, d::Dense{bias}) where {bias}
    print(io, "Dense($(d.in_dims) => $(d.out_dims)")
    (d.activation == identity) || print(io, ", $(d.activation)")
    bias || print(io, ", bias=false")
    return print(io, ")")
end

function Dense(mapping::Pair{<:Int, <:Int}, activation=identity; init_weight=glorot_uniform,
               init_bias=zeros32, bias::Bool=true)
    return Dense(first(mapping), last(mapping), activation; init_weight=init_weight,
                 init_bias=init_bias, bias=bias)
end

function Dense(in_dims::Int, out_dims::Int, activation=identity; init_weight=glorot_uniform,
               init_bias=zeros32, bias::Bool=true)
    activation = NNlib.fast_act(activation)
    return Dense{bias, typeof(activation), typeof(init_weight), typeof(init_bias)}(activation,
                                                                                in_dims,
                                                                                out_dims,
                                                                                init_weight,
                                                                                init_bias)
end

function initialparameters(rng::AbstractRNG, d::Dense{bias}) where {bias}
    if bias
        return (weight=d.init_weight(rng, d.out_dims, d.in_dims),
                bias=d.init_bias(rng, d.out_dims, 1))
    else
        return (weight=d.init_weight(rng, d.out_dims, d.in_dims),)
    end
end

function parameterlength(d::Dense{bias}) where {bias}
    return bias ? d.out_dims * (d.in_dims + 1) : d.out_dims * d.in_dims
end
statelength(d::Dense) = 0

@inline function (d::Dense{false})(x::AbstractVecOrMat,
                                   ps::Union{ComponentArray, NamedTuple},
                                   st::NamedTuple)
    return applyactivation(d.activation, ps.weight * x), st
end

@inline function (d::Dense{false, typeof(identity)})(x::AbstractVecOrMat,
                                                     ps::Union{ComponentArray, NamedTuple},
                                                     st::NamedTuple)
    return ps.weight * x, st
end

@inline function (d::Dense{false})(x::AbstractArray,
                                   ps::Union{ComponentArray, NamedTuple},
                                   st::NamedTuple)
    sz = size(x)
    x_reshaped = reshape(x, sz[1], :)
    return reshape(applyactivation(d.activation, ps.weight * x_reshaped), d.out_dims,
                   sz[2:end]...), st
end

@inline function (d::Dense{false, typeof(identity)})(x::AbstractArray,
                                                     ps::Union{ComponentArray, NamedTuple},
                                                     st::NamedTuple)
    sz = size(x)
    x_reshaped = reshape(x, sz[1], :)
    return reshape(ps.weight * x_reshaped, d.out_dims, sz[2:end]...), st
end

@inline function (d::Dense{true})(x::AbstractVector, ps::Union{ComponentArray, NamedTuple},
                                  st::NamedTuple)
    return applyactivation(d.activation, elementwise_add(ps.weight * x, vec(ps.bias))), st
end

@inline function (d::Dense{true, typeof(identity)})(x::AbstractVector,
                                                    ps::Union{ComponentArray, NamedTuple},
                                                    st::NamedTuple)
    return elementwise_add(ps.weight * x, vec(ps.bias)), st
end

@inline function (d::Dense{true})(x::AbstractMatrix, ps::Union{ComponentArray, NamedTuple},
                                  st::NamedTuple)
    return applyactivation(d.activation, elementwise_add(ps.weight * x, ps.bias)), st
end

@inline function (d::Dense{true, typeof(identity)})(x::AbstractMatrix,
                                                    ps::Union{ComponentArray, NamedTuple},
                                                    st::NamedTuple)
    return elementwise_add(ps.weight * x, ps.bias), st
end

@inline function (d::Dense{true})(x::AbstractArray, ps::Union{ComponentArray, NamedTuple},
                                  st::NamedTuple)
    sz = size(x)
    x_reshaped = reshape(x, sz[1], :)
    return (reshape(applyactivation(d.activation,
                                    elementwise_add(ps.weight * x_reshaped, ps.bias)),
                    d.out_dims, sz[2:end]...), st)
end

@inline function (d::Dense{true, typeof(identity)})(x::AbstractArray,
                                                    ps::Union{ComponentArray, NamedTuple},
                                                    st::NamedTuple)
    sz = size(x)
    x_reshaped = reshape(x, sz[1], :)
    return (reshape(elementwise_add(ps.weight * x_reshaped, ps.bias), d.out_dims,
                    sz[2:end]...), st)
end

"""
    Scale(dims, activation=identity; init_weight=ones32, init_bias=zeros32, bias::Bool=true)

Create a Sparsely Connected Layer with a very specific structure (only Diagonal Elements are non-zero). The forward pass is given by: `y = activation.(weight .* x .+ bias)`

## Arguments

* `dims`: size of the learnable scale and bias parameters.
* `activation`: activation function

## Keyword Arguments

* `init_weight`: initializer for the weight matrix (`weight = init_weight(rng, out_dims, in_dims)`)
* `init_bias`: initializer for the bias vector (ignored if `bias=false`)
* `bias`: whether to include a bias vector

## Input

* `x` must be an Array of size `(dims..., B)` or `(dims...[0], ..., dims[k])` for `k ≤ size(dims)`

## Returns

* Array of size `(dims..., B)` or `(dims...[0], ..., dims[k])` for `k ≤ size(dims)`
* Empty `NamedTuple()`

## Parameters

* `weight`: Weight Array of size `(dims...)`
* `bias`: Bias of size `(dims...)`

!!! compat "Lux 0.4.3"
    `Scale` with multiple dimensions requires at least Lux 0.4.3.
"""
struct Scale{bias, F1, D, F2, F3} <: AbstractExplicitLayer
    activation::F1
    dims::D
    init_weight::F2
    init_bias::F3
end

function Base.show(io::IO, d::Scale)
    print(io, "Scale($(d.dims)")
    (d.activation == identity) || print(io, ", $(d.activation)")
    return print(io, ")")
end

function Scale(dims::Tuple{Vararg{Integer}}, activation=identity;
               init_weight=glorot_uniform,
               init_bias=zeros32, bias::Bool=true)
    activation = NNlib.fast_act(activation)
    return Scale{bias, typeof(activation), typeof(dims), typeof(init_weight),
                 typeof(init_bias)}(activation, dims, init_weight, init_bias)
end

function Scale(s1::Integer, s23::Integer...; _act=identity, kw...)
    Scale(tuple(s1, s23...), _act; kw...)
end
Scale(size_act...; kw...) = Scale(size_act[1:(end - 1)]...; _act=size_act[end], kw...)

function initialparameters(rng::AbstractRNG, d::Scale{true})
    return (weight=d.init_weight(rng, d.dims...), bias=d.init_bias(rng, d.dims...))
end
function initialparameters(rng::AbstractRNG, d::Scale{false})
    (weight=d.init_weight(rng, d.dims...),)
end

parameterlength(d::Scale{bias}) where {bias} = (1 + bias) * prod(d.dims)
statelength(d::Scale) = 0

function (d::Scale{true})(x::AbstractArray, ps::Union{ComponentArray, NamedTuple},
                          st::NamedTuple)
    return applyactivation(d.activation,
                           elementwise_add(elementwise_mul(ps.weight, x), ps.bias)), st
end

function (d::Scale{true, typeof(identity)})(x::AbstractArray,
                                            ps::Union{ComponentArray, NamedTuple},
                                            st::NamedTuple)
    return elementwise_add(elementwise_mul(ps.weight, x), ps.bias), st
end

function (d::Scale{false})(x::AbstractArray, ps::Union{ComponentArray, NamedTuple},
                           st::NamedTuple)
    return applyactivation(d.activation, elementwise_mul(ps.weight, x)), st
end

function (d::Scale{false, typeof(identity)})(x::AbstractArray,
                                             ps::Union{ComponentArray, NamedTuple},
                                             st::NamedTuple)
    return elementwise_mul(ps.weight, x), st
end