text:
#=
Copyright 2021 BlackRock, Inc.

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

    http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
=#

"""
    PiecewiseQuadratic()
    PiecewiseQuadratic(f::BoundedQuadratic)
    PiecewiseQuadratic(f_list::Vector{BoundedQuadratic}[; simplify_result=false])

Represents piecewise quadratic function, where each piece is a `BoundedQuadratic`.

The fields represent:
- `f_list`: a `Vector` of `BoundedQuadratic` functions

# Example
```jldoctest
julia> left = BoundedQuadratic(-Inf, 0., 0., -1., 0.)
BoundedQuadratic: f(x) = - 1.00000 x , ∀x ∈ [-Inf, 0.00000]

julia> right = BoundedQuadratic(0., Inf, 0., 1., 0.)
BoundedQuadratic: f(x) = + 1.00000 x , ∀x ∈ [0.00000, Inf]

julia> pwq = PiecewiseQuadratic([left, right])
Piecewise quadratic function:
BoundedQuadratic: f(x) = - 1.00000 x , ∀x ∈ [-Inf, 0.00000]
BoundedQuadratic: f(x) = + 1.00000 x , ∀x ∈ [0.00000, Inf]

```
"""
mutable struct PiecewiseQuadratic
    f_list::Vector{BoundedQuadratic}

    function PiecewiseQuadratic(f_list::Vector{BoundedQuadratic}; simplify_result=false)
        if simplify_result
            f_list = simplify(PiecewiseQuadratic(f_list; simplify_result=false)).f_list
        end
        return new(f_list)
    end
end

PiecewiseQuadratic() = PiecewiseQuadratic(BoundedQuadratic[])
PiecewiseQuadratic(f::BoundedQuadratic) = PiecewiseQuadratic([f])

"""
    indicator(dom::Interval)
    indicator(lb::Real, ub::Real)

Construct a PiecewiseQuadratic that is `0` on the interval and `Inf` everywhere else.

# Example
```jldoctest
julia> indicator(-5, 5)
Piecewise quadratic function:
BoundedQuadratic: f(x) = 0, ∀x ∈ [-5.00000, 5.00000]

```
"""
indicator(dom::Interval) = PiecewiseQuadratic(BoundedQuadratic(dom, 0, 0, 0))
indicator(lb::Real, ub::Real) = indicator(Interval(lb, ub))

function Base.show(io::IO, f::PiecewiseQuadratic)
    @printf(io, "Piecewise quadratic function:\n")
    for fi in f
        print(io, fi)
    end
end

"""
    zero(::Type{PiecewiseQuadratic})

Construct an empty PiecewiseQuadratic with no constraints.

# Example
```jldoctest
julia> zero(PiecewiseQuadratic)
Piecewise quadratic function:
BoundedQuadratic: f(x) = 0, ∀x ∈ ℝ

```
"""
zero(::Type{PiecewiseQuadratic}) = indicator(Interval())
copy(f::PiecewiseQuadratic) = PiecewiseQuadratic([copy(fi) for fi in f])
function copy!(out::PiecewiseQuadratic, f::PiecewiseQuadratic)
    out.f_list = f.f_list
    return out
end

#####
##### Indexing
#####

getindex(f::PiecewiseQuadratic, idx) = getindex(f.f_list, idx)
setindex!(f::PiecewiseQuadratic, val, idx) = setindex!(f.f_list, val, idx)
lastindex(f::PiecewiseQuadratic) = lastindex(f.f_list)
iterate(f::PiecewiseQuadratic) = iterate(f.f_list)
iterate(f::PiecewiseQuadratic, i::Int) = iterate(f.f_list, i::Int)
length(f::PiecewiseQuadratic) = length(f.f_list)
pop!(f::PiecewiseQuadratic) = pop!(f.f_list)

#####
##### Boolean
#####

"""
    isempty(f::PiecewiseQuadratic)

Return `true` if the PiecewiseQuadratic `f` is empty (`f_list` empty).
"""
isempty(f::PiecewiseQuadratic) = length(f) == 0

function ≈(f::PiecewiseQuadratic, g::PiecewiseQuadratic)
    if length(f) != length(g)
        return false
    end
    isequal = true
    for (fi, gi) in zip(f, g)
        isequal &= fi ≈ gi
    end
    return isequal
end

"""
    is_convex(f::PiecewiseQuadratic)

Return `true` if `f` is convex.

A `PiecewiseQuadratic` is convex if for all `i`:
* `f_i` is convex
* `f_i.ub == f_{i+1}.lb`
* `derivative(f_i)(f_i.ub) <= derivative(f_{i+1})(f_{i+1}.lb)`

# Example
```jldoctest
julia> left = BoundedQuadratic(-Inf, 0., 0., -1., 0.)
BoundedQuadratic: f(x) = - 1.00000 x , ∀x ∈ [-Inf, 0.00000]

julia> right = BoundedQuadratic(0., Inf, 0., 1., 0.)
BoundedQuadratic: f(x) = + 1.00000 x , ∀x ∈ [0.00000, Inf]

julia> pwq = PiecewiseQuadratic([left, right])
Piecewise quadratic function:
BoundedQuadratic: f(x) = - 1.00000 x , ∀x ∈ [-Inf, 0.00000]
BoundedQuadratic: f(x) = + 1.00000 x , ∀x ∈ [0.00000, Inf]

julia> is_convex(pwq)
true

```
"""
function is_convex(f::PiecewiseQuadratic)
    n = length(f.f_list)
    if n <= 1
        return is_convex(f.f_list[1])
    end
    for i in 1:(n - 1)
        left = f.f_list[i]
        right = f.f_list[i + 1]
        left_deriv = derivative(left, left.ub)
        right_deriv = derivative(right, right.lb)

        # not convex if not continuous or f_left's derivative is greater than f_right's
        if !continuous_and_overlapping(left, right) ||
           left_deriv - right_deriv > ϵ ||
           !is_convex(left)
            return false
        end
    end
    return true
end

#####
##### Operations
#####

"""
    (f::PiecewiseQuadratic)(x::Real)

Evaluate `f(x)` if `x` is in the domain of `f`, else return `Inf`.
"""
function (f::PiecewiseQuadratic)(x::Float64)
    for fi in f
        fi_x = fi(x)
        if fi_x < Inf
            return fi_x
        end
    end
    return Inf
end

+(f::PiecewiseQuadratic, g::PiecewiseQuadratic) = sum([f, g])
+(f::BoundedQuadratic, g::PiecewiseQuadratic) = g + PiecewiseQuadratic(f)
+(f::PiecewiseQuadratic, g::BoundedQuadratic) = g + f

-(f::PiecewiseQuadratic) = PiecewiseQuadratic(.-(f))

*(f::PiecewiseQuadratic, α::Real) = PiecewiseQuadratic(α .* f)
*(α::Real, f::PiecewiseQuadratic) = f * α

mult_scalar!(f::PiecewiseQuadratic, α::Real) = mult_scalar!(f, α, f)
function mult_scalar!(f::PiecewiseQuadratic, α::Real, out::PiecewiseQuadratic)
    @assert length(f) == length(out)
    return copy!(out, α * f)
end

"""
    scale(f::PiecewiseQuadratic, α::Real)

Return a new `PiecewiseQuadratic` that has been scaled by `α`. That is, given `f(x)` and `α`, returns `f(αx)`.

Note: this operation requires scaling the domain.

See also: [`scale!`](@ref)

# Example
```jldoctest
julia> f = PiecewiseQuadratic([BoundedQuadratic(-1., 0., 0., -1., 0.),
                               BoundedQuadratic(0., 1., 0., 1., 0.),
                               BoundedQuadratic(1., 5., 1., 1., 1.)])
Piecewise quadratic function:
BoundedQuadratic: f(x) = - 1.00000 x , ∀x ∈ [-1.00000, 0.00000]
BoundedQuadratic: f(x) = + 1.00000 x , ∀x ∈ [0.00000, 1.00000]
BoundedQuadratic: f(x) = 1.00000 x² + 1.00000 x + 1.00000, ∀x ∈ [1.00000, 5.00000]

julia> scale(f, 5)
Piecewise quadratic function:
BoundedQuadratic: f(x) = - 5.00000 x , ∀x ∈ [-0.20000, 0.00000]
BoundedQuadratic: f(x) = + 5.00000 x , ∀x ∈ [0.00000, 0.20000]
BoundedQuadratic: f(x) = 25.00000 x² + 5.00000 x + 1.00000, ∀x ∈ [0.20000, 1.00000]

```
"""
scale(f::PiecewiseQuadratic, α::Real) = PiecewiseQuadratic(scale.(f, α))

"""
    scale!(f::PiecewiseQuadratic, α::Real)
    scale!(f::PiecewiseQuadratic, α::Real, out::PiecewiseQuadratic)

Scale `f` inplace.

See also: [`scale`](@ref)
"""
scale!(f::PiecewiseQuadratic, α::Real) = scale!.(f, α)
scale!(f::PiecewiseQuadratic, α::Real, out::PiecewiseQuadratic) = copy!(out, scale(f, α))

"""
    perspective(f::PiecewiseQuadratic, α::Real)

Return the perspective function of `f`. That is, given `f(x)` and `α`, return `α * f(x / α)`.

Note: that this operation requires scaling of the domain.

See also: [`perspective!`](@ref)

# Example
```jldoctest
julia> f = PiecewiseQuadratic([BoundedQuadratic(-1., 0., 0., -5., 0.),
                               BoundedQuadratic(0., 2., 0., 2., 0.)])
Piecewise quadratic function:
BoundedQuadratic: f(x) = - 5.00000 x , ∀x ∈ [-1.00000, 0.00000]
BoundedQuadratic: f(x) = + 2.00000 x , ∀x ∈ [0.00000, 2.00000]

julia> persp = perspective(f, 5.)
Piecewise quadratic function:
BoundedQuadratic: f(x) = - 5.00000 x , ∀x ∈ [-5.00000, 0.00000]
BoundedQuadratic: f(x) = + 2.00000 x , ∀x ∈ [0.00000, 10.00000]

```
"""
perspective(f::PiecewiseQuadratic, α::Real) = PiecewiseQuadratic(perspective.(f, α))

"""
    perspective!(f::PiecewiseQuadratic, α::Real)
    perspective!(f::PiecewiseQuadratic, α::Real, out::PiecewiseQuadratic)

Shift perspective of `f` inplace.

See also: [`perspective`](@ref)
"""
perspective!(f::PiecewiseQuadratic, α::Real) = perspective!.(f, α)
function perspective!(f::PiecewiseQuadratic, α::Real, out::PiecewiseQuadratic)
    return copy!(out, perspective(f, α))
end

"""
    shift(f::PiecewiseQuadratic, δ::Real)

Return `f` shifted along the `x`-axis by `δ`.

Note: for `δ > 0`, this is a right shift. For `δ < 0`, this is a left shift.

See also: [`shift!`](@ref)

# Example
```jldoctest
julia> f = PiecewiseQuadratic([BoundedQuadratic(-1., 0., 0., -5., 0.),
                               BoundedQuadratic(0., 2., 0., 2., 0.)])
Piecewise quadratic function:
BoundedQuadratic: f(x) = - 5.00000 x , ∀x ∈ [-1.00000, 0.00000]
BoundedQuadratic: f(x) = + 2.00000 x , ∀x ∈ [0.00000, 2.00000]

julia> shift(f, 5.)
Piecewise quadratic function:
BoundedQuadratic: f(x) = - 5.00000 x + 25.00000, ∀x ∈ [4.00000, 5.00000]
BoundedQuadratic: f(x) = + 2.00000 x - 10.00000, ∀x ∈ [5.00000, 7.00000]

```
"""
shift(f::PiecewiseQuadratic, δ::Real) = PiecewiseQuadratic(shift.(f, δ))

"""
    shift!(f::PiecewiseQuadratic, δ::Real)
    shift!(f::PiecewiseQuadratic, δ::Real, out::PiecewiseQuadratic)

Shift `f` inplace along the `x`-axis by `δ`.

See also: [`shift`](@ref)
"""
shift!(f::PiecewiseQuadratic, δ::Real) = shift!.(f, δ)
shift!(f::PiecewiseQuadratic, δ::Real, out::PiecewiseQuadratic) = copy!(out, shift(f, δ))

"""
    tilt(f::PiecewiseQuadratic, α::Real)

Return `f` tilted by `α`. This shifts linear coefficient `q` by `α`.

See also: [`tilt!`](@ref)

# Example
```jldoctest
julia> f = PiecewiseQuadratic([BoundedQuadratic(-1., 0., 0., -5., 0.),
                               BoundedQuadratic(0., 2., 0., 2., 0.)])
Piecewise quadratic function:
BoundedQuadratic: f(x) = - 5.00000 x , ∀x ∈ [-1.00000, 0.00000]
BoundedQuadratic: f(x) = + 2.00000 x , ∀x ∈ [0.00000, 2.00000]

julia> tilt(f, 5.)
Piecewise quadratic function:
BoundedQuadratic: f(x) = 0, ∀x ∈ [-1.00000, 0.00000]
BoundedQuadratic: f(x) = + 7.00000 x , ∀x ∈ [0.00000, 2.00000]

```
"""
tilt(f::PiecewiseQuadratic, α::Real) = PiecewiseQuadratic(tilt.(f, α))

"""
    tilt!(f::PiecewiseQuadratic, α::Real)
    tilt!(f::PiecewiseQuadratic, α::Real, out::PiecewiseQuadratic)

Tilt `f` inplace.
"""
tilt!(f::PiecewiseQuadratic, α::Real) = tilt!.(f, α)
tilt!(f::PiecewiseQuadratic, α::Real, out::PiecewiseQuadratic) = copy!(out, tilt(f, α))

"""
     push!(f::PiecewiseQuadratic, g::BoundedQuadratic; simplify_result=false)

Add a BoundedQuadratic `g` to the right hand side of PiecewiseQuadratic `f`.
"""
function push!(f::PiecewiseQuadratic, g::BoundedQuadratic; simplify_result=false)
    push!(f.f_list, g)
    if simplify_result
        simplify(f; start_idx=length(f) - 1)
    end
end

"""
    append!(f::PiecewiseQuadratic, g::PiecewiseQuadratic; simplify_result=false)

Append a PiecewiseQuadratic `g` to the right hand side of PiecewiseQuadratic `f`.
"""
function append!(f::PiecewiseQuadratic, g::PiecewiseQuadratic; simplify_result=false)
    k = length(f)
    append!(f.f_list, g.f_list)
    if simplify_result && (k > 0)
        simplify(f; start_idx=k)
    end
end

"""
    reverse(f::PiecewiseQuadratic)

Return `f` reversed over the `y` axis. That is, given `f(x)`, return `f(-x)`.

See also: [`reverse`](@ref)
"""
reverse(f::PiecewiseQuadratic) = PiecewiseQuadratic(reverse(reverse.(f)))

"""
    reverse!(f::PiecewiseQuadratic)
    reverse!(f::PiecewiseQuadratic, out::PiecewiseQuadratic)

Reverse `f` inplace.

See also: [`reverse`](@ref)
"""
reverse!(f::PiecewiseQuadratic) = reverse!(f, f)
reverse!(f::PiecewiseQuadratic, out::PiecewiseQuadratic) = copy!(out, reverse(f))

"""
     simplify(f::PiecewiseQuadratic)

Return a simplified piecewise quadratic `f`.

A BoundedQuadratic `f_i` and and the current rightmost piece should be combined (to become the new rightmost piece) if
* Either of `f_i` or the rightmost piece is a point and they overlap at their lower and upper endpoints respectively
* `f_i` and the rightmost piece have the same coefficients and they correspond at their lower and upper endpoints respectively (they're just currently separate parts of the same function).

We also eliminate functions that are points at Inf or -Inf if they come up.

# Example
```jldoctest
julia> f = PiecewiseQuadratic([BoundedQuadratic(-1., 0., 0., -1., 0.),
                                      BoundedQuadratic(0., 0., 0., 0., 0.),
                                      BoundedQuadratic(0., 1., 0., 1., 0.),
                                      BoundedQuadratic(1., 2., 0., 1., 0.),
                                      BoundedQuadratic(3., 5., 1., 1., 1.)])
Piecewise quadratic function:
BoundedQuadratic: f(x) = - 1.00000 x , ∀x ∈ [-1.00000, 0.00000]
BoundedQuadratic: f(x) = 0, ∀x ∈ [0.00000, 0.00000]
BoundedQuadratic: f(x) = + 1.00000 x , ∀x ∈ [0.00000, 1.00000]
BoundedQuadratic: f(x) = + 1.00000 x , ∀x ∈ [1.00000, 2.00000]
BoundedQuadratic: f(x) = 1.00000 x² + 1.00000 x + 1.00000, ∀x ∈ [3.00000, 5.00000]

julia> simplify(f)
Piecewise quadratic function:
BoundedQuadratic: f(x) = - 1.00000 x , ∀x ∈ [-1.00000, 0.00000]
BoundedQuadratic: f(x) = + 1.00000 x , ∀x ∈ [0.00000, 2.00000]
BoundedQuadratic: f(x) = 1.00000 x² + 1.00000 x + 1.00000, ∀x ∈ [3.00000, 5.00000]

```
"""
function simplify(f::PiecewiseQuadratic)
    front, back = 1, length(f)
    if isempty(f[front])
        front += 1
    end
    # out = FixedMemoryPwq(back-front)
    out = PiecewiseQuadratic()
    push!(out, f[front])
    for (i, fi) in enumerate(f[(front + 1):end])
        if isempty(fi)
            continue
        end
        prev = out[end]

        val = fi(fi.lb)
        prev_val = prev(prev.ub)
        # if we have redundant points, keep the one with minimum value
        if is_point(prev) && is_point(fi) && prev.ub ≈ fi.lb
            # ignore right point if greater than the value of the left
            if prev_val ≳ val
                out[end] = fi
            end
        elseif continuous_and_overlapping(prev, fi) && (is_point(prev) || is_point(fi)) ||
               extend_dom(prev) ≈ extend_dom(fi)
            # take the coefficients of whichever function is not a point.
            # note: if we are combining and neither function is a point, both coefficient sets must be
            # the same, so it doesn't matter which we take.
            non_point = is_point(prev) ? fi : prev
            combined = restrict_dom(extend_dom(non_point), prev.lb, fi.ub)
            out[end] = combined
        elseif is_point(fi) && prev.ub ≈ fi.lb && val ≳ prev_val
            continue  # ignore
        elseif is_point(prev) && prev.ub ≈ fi.lb && prev_val ≳ val
            out[end] = fi  # replace
        else
            # no combination, add the current piece
            push!(out, fi)
        end
    end
    return out
end

"""
    sum(f_list::Vector{PiecewiseQuadratic})

Return the PiecewiseQuadratic sum of a list of PiecewiseQuadratics.

# Example
```jldoctest
julia> f = PiecewiseQuadratic([BoundedQuadratic(-1., 0., 0., -5., 0.),
                               BoundedQuadratic(0., 2., 0., 2., 0.)])
Piecewise quadratic function:
BoundedQuadratic: f(x) = - 5.00000 x , ∀x ∈ [-1.00000, 0.00000]
BoundedQuadratic: f(x) = + 2.00000 x , ∀x ∈ [0.00000, 2.00000]

julia> sum([f,f])
Piecewise quadratic function:
BoundedQuadratic: f(x) = - 10.00000 x , ∀x ∈ [-1.00000, 0.00000]
BoundedQuadratic: f(x) = + 4.00000 x , ∀x ∈ [0.00000, 2.00000]

julia> sum([f,-f])
Piecewise quadratic function:
BoundedQuadratic: f(x) = 0, ∀x ∈ [-1.00000, 0.00000]
BoundedQuadratic: f(x) = 0, ∀x ∈ [0.00000, 2.00000]

```
"""
function sum(f_list::Vector{PiecewiseQuadratic})
    work = PwqSumWorkspace(length(f_list))
    return _sum(f_list, work)
end

mutable struct PwqSumWorkspace
    f_active::Vector{BoundedQuadratic}
    active_idxs::Vector{Int64}
    active::BitArray
    ubs::Vector{Float64}
    next_fun_idxs::BitArray
    k::Int64

    function PwqSumWorkspace(k::Int64)
        f_active = Vector{BoundedQuadratic}(undef, k)
        active_idxs = Vector{Int64}(undef, k)
        active = BitArray(undef, k)
        ubs = Vector{Float64}(undef, k)
        next_fun_idxs = BitArray(undef, k)
        return new(f_active, active_idxs, active, ubs, next_fun_idxs, k)
    end
end

function _sum(f_list::Vector{PiecewiseQuadratic}, w::PwqSumWorkspace)
    k = length(f_list)
    @assert w.k == k
    @assert k > 0

    # Currently active functions
    for i in 1:k
        w.f_active[i] = f_list[i][1]
        w.active_idxs[i] = 1
        w.active[i] = length(f_list[i]) > 1
    end

    # Output function:
    g = PiecewiseQuadratic()
    sum_active, is_valid = _sum(w.f_active)
    if is_valid
        push!(g, sum(w.f_active))
    end

    # Loop to sum functions:
    while any(w.active)

        # Populate next_fun_idxs with all active functions with minimum upper bound:
        _get_min_ub!(w.f_active, w.active, w.ubs, w.next_fun_idxs)

        for i in 1:k
            if w.next_fun_idxs[i]
                w.active_idxs[i] += 1
                w.active[i] = w.active_idxs[i] < length(f_list[i])
                w.f_active[i] = f_list[i][w.active_idxs[i]]
            end
        end

        # Append the sum, if the domain is not empty:
        sum_active, is_valid = _sum(w.f_active)
        if is_valid
            push!(g, sum_active)
        end
    end
    return g
end

function _get_min_ub!(f_vec::Vector{BoundedQuadratic}, eligible::BitArray,
                      ubs::Vector{Float64}, idxs::BitArray)
    k = length(f_vec)
    ub_min = Inf
    for i in 1:k
        if eligible[i]
            ubs[i] = f_vec[i].ub
            ub_min = min(ubs[i], ub_min)
        end
    end
    for i in 1:k
        idxs[i] = eligible[i] && ubs[i] == ub_min
    end
end

function _get_min_ub(f_vec::Vector{BoundedQuadratic}, eligible::BitArray)
    k = length(f_vec)
    @assert k == length(eligible)
    ubs_and_idxs = [(fi.ub, i) for (i, fi) in enumerate(f_vec)][eligible]
    ub_min = minimum([ui[1] for ui in ubs_and_idxs])
    # Select all upper bounds equal to the minimum:
    idxs = [i for (ub, i) in ubs_and_idxs if (ub == ub_min)]
    return idxs
end

function _get_min_ub(f_vec::Vector{BoundedQuadratic}, eligible::Vector{Bool})
    return _get_min_ub(f_vec, convert(BitArray, eligible))
end

"""
    minimize(f::PiecewiseQuadratic)

Return the minimum `x` and `f(x)` of `f` over its domain.
"""
function minimize(f::PiecewiseQuadratic)
    v_star = Inf
    x_star = NaN
    for fi in f
        xi, vi = minimize(fi)
        if vi < v_star
            x_star = xi
            v_star = vi
        end
    end
    return x_star, v_star
end

"""
    prox(f::PiecewiseQuadratic, u::Float64[, ρ::Float64=1.0]; use_quadratic::Bool=true)

Return the proximal operator of `f`, `ρ` evaluated at `u`.

# Note: The proximal operator of `f`, `rho` is denoted:
```math
prox_{f, rho}(u) = argmin_{x ∈ dom(f)} f(x) + (rho / 2)||x - u||_2^2
```

See Section 6.2 of [arXiv:2103.05455](https://arxiv.org/abs/2103.05455) for more information.
"""
function prox(f::PiecewiseQuadratic, u::Float64, ρ::Float64=1.0; use_quadratic::Bool=true)
    ub_last = -Inf
    for fi in f
        p_adj = 2 * fi.p + (use_quadratic ? ρ : 0.0)
        lb = p_adj * fi.lb + fi.q
        ub = p_adj * fi.ub + fi.q
        if ub_last ≤ ρ * u ≤ lb
            return fi.lb
        elseif lb ≤ ρ * u ≤ ub
            return (ρ * u - fi.q) / p_adj
        end
        ub_last = ub
    end
    return f[end].ub
end

"""
    get_plot(f::PiecewiseQuadratic; N::Int64=1000)

Return `x`, `y` `N`-vectors for use with plotting libraries.

# Example
```julia
julia> using Plots

julia> f = PiecewiseQuadratic([
         BoundedQuadratic(-Inf, 3.0, 1.0, -3.0, 3.0),
         BoundedQuadratic(3.0, 4.0, 0.0, -1.0, 3.0),
         BoundedQuadratic(4.0, 6.0, 2.0, -20.0, 47.0),
         BoundedQuadratic(6.0, 7.5, 0.0, 1.0, -7.0),
         BoundedQuadratic(7.5, Inf, 0.0, 4.0, -29.0)
       ])
Piecewise quadratic function:
BoundedQuadratic: f(x) = 1.00000 x² - 3.00000 x + 3.00000, ∀x ∈ [-Inf, 3.00000]
BoundedQuadratic: f(x) = - 1.00000 x + 3.00000, ∀x ∈ [3.00000, 4.00000]
BoundedQuadratic: f(x) = 2.00000 x² - 20.00000 x + 47.00000, ∀x ∈ [4.00000, 6.00000]
BoundedQuadratic: f(x) = + 1.00000 x - 7.00000, ∀x ∈ [6.00000, 7.50000]
BoundedQuadratic: f(x) = + 4.00000 x - 29.00000, ∀x ∈ [7.50000, Inf]

julia> plot(get_plot(f); grid=false, linestyle=:dash, color=:black, label="piece-wise quadratic")
Plot{Plots.GRBackend() n=1}

julia> plot!(get_plot(simplify(envelope(f))); color=:blue, la=0.5, label="envelope")
Plot{Plots.GRBackend() n=2}

```
"""
function get_plot(f::PiecewiseQuadratic; N::Int64=1000)
    x = vcat([[fi.lb, fi.ub] for fi in f]...)
    x = x[isfinite.(x)]
    if isempty(x)
        x = [-5.0, 5]
    end
    x_min = minimum(x)
    x_max = maximum(x)
    span = x_max - x_min
    if span == 0.0
        span = 2.0
    end
    x_min -= 1 * span
    x_max += 1 * span
    x = vcat(range(x_min; stop=x_max, length=N), x)
    x = sort(x)
    x = [x0 for (x0, x1) in zip(x[1:(end - 1)], x[2:end]) if x0 != x1]
    y = f.(x)
    return x, y
end


text:
using ITensorGaussianMPS
using LinearAlgebra
using Test

@testset "ITensorGaussianMPS.jl" begin
  include("gmps.jl")
  include("electron.jl")
end

ext: jl


text:
# shows some data analysis using DataFrames
using CSV, DataFrames, GLM
dir = dirname(dirname(pathof(Econometrics)))
nerlove = CSV.read(dir*"/Examples/Data/nerlove.csv")
f = @formula(log(cost) ~ 1 + log(output) + log(labor) + log(fuel) + log(capital))
println()
println("using the GLM package")
@show lm(f, nerlove)
println()
println("using ols from this package")
ols(f, nerlove)
nothing

max_stars_repo_name: nilshg/Econometrics

max_issues_repo_name: nilshg/Econometrics

max_forks_repo_name: nilshg/Econometrics


text:
using ReadableRegex
using Test

@testset "ReadableRegex.jl" begin

    str1 = "Main Street 28"
    str2 = "Another Road 371"
    str3 = " Another Road 371"
    str4 = "Another Road 371a"
    str5 = "Another Road 371 "

    r = BEGIN * one_or_more(WORD) *
            zero_or_more(WHITESPACE * one_or_more(WORD)) *
            WHITESPACE * one_or_more(DIGIT) * END

    @test match(r, str1).match == str1
    @test match(r, str2).match == str2
    @test isnothing(match(r, str3))
    @test isnothing(match(r, str4))
    @test isnothing(match(r, str5))

    str6 = "a, bb, ccc"
    @test match(at_least(3, WORD), str6).match == "ccc"

    str7 = "a, bbb, ccccc"
    @test match(between(2, 4, WORD), str7).match == "bbb"
end

@testset "Lookarounds" begin
    str = "1a 2 b3 c4d"

    @test match(look_for(DIGIT, after = WORD), str).match == "3"
    @test match(look_for(DIGIT, not_before = WORD), str).match == "2"

    str2 = "123 for me, 456 for you"
    @test match(look_for(one_or_more(DIGIT) , before = " for me"), str2).match == "123"

    str3 = "a2 3"
    @test match(look_for(DIGIT, not_after = WORD), str3).match == "3"

    str4 = "1 2.0 .3 -.4 -5 60 700 800.9 +9000"

    reg = look_for(
            maybe(char_in("+-")) * one_or_more(DIGIT),
            not_after = ".",
            not_before = NON_SEPARATOR)

    @test collect(m.match for m in eachmatch(reg, str4)) == ["1", "-5", "60", "700", "+9000"]
end

@testset "Numbers" begin
    str = "1 2.0 .3 -.4 -5 60 700 800.9 +9000"

    matches = eachmatch(
        maybe(["-", "+"]) * maybe(zero_or_more(DIGIT) * ".") * one_or_more(DIGIT),
        str
    )

    @test all(split(str) .== (m.match for m in matches))
end

@testset "Range" begin
    str = "ab cd ef gh"

    matches = [m.match for m in eachmatch(exactly(2, 'a':'e'), str)]

    @test matches == ["ab", "cd"]
end

@testset "Compile" begin
    function matcher(str)
        reg = @compile BEGIN * maybe('a':'z') * exactly(1, DIGIT) * END
        match(reg, str)
    end

    @test !isnothing(matcher("a1"))
    @test !isnothing(matcher("3"))
    @test isnothing(matcher("2a"))
end

@testset "Capture and Reference" begin
    str = "pasta pesto"
    reg = "p" * capture(char_in("aeiou")) * "st" * reference(1)
    @test match(reg, str).match == "pasta"

    reg2 = "p" * capture(char_in("aeiou"), as = "vowel") * "st" * reference("vowel")
    @test match(reg, str).match == "pasta"
end

@testset "Unicode categories" begin
    str = "how is Jim doing?"
    @test match(UPPERCASE, str).match == "J"
    @test match(PUNCTUATION, str).match == "?"
end

@testset "multiarg char_in" begin
    str = "hakuna matata"

    reg = at_least(2, char_in("mat"))
    @test match(reg, str).match == "matata"

    reg2 = one_or_more(char_in('h':'n', 'a', 'u'))
    @test match(reg2, str).match == "hakuna"

    reg3 = at_least(2, char_not_in(WHITESPACE, 'a'))
    @test match(reg3, str).match == "kun"

    str = "-"
    reg = char_in("a-z") # this should give three chars, not a range
    @test match(reg, str).match == "-"
    reg2 = char_in('a':'z') # this should give a range
    @test match(reg2, str) === nothing
end

@testset "lazy and possessive quantifiers" begin

    str = "abcabcabc"

    @test match(maybe("abc"), str).match == "abc"
    @test match(lazy_maybe("abc"), str).match == ""

    @test match(at_least(2, "abc"), str).match == str
    @test match(lazy_at_least(2, "abc"), str).match == "abcabc"

    @test match(zero_or_more("abc"), str).match == str
    @test match(lazy_zero_or_more("abc"), str).match == ""

    @test match(one_or_more("abc"), str).match == str
    @test match(lazy_one_or_more("abc"), str).match == "abc"

    @test match(between(2, 3, "abc"), str).match == str
    @test match(lazy_between(2, 3, "abc"), str).match == "abcabc"

    str2 = "1234"

    @test match(one_or_more(DIGIT) * "4", str2).match == "1234"
    @test match(possessive_one_or_more(DIGIT) * "4", str2) === nothing

    @test match(zero_or_more(DIGIT) * "4", str2).match == "1234"
    @test match(possessive_zero_or_more(DIGIT) * "4", str2) === nothing

    @test match(between(1, 4, DIGIT) * "4", str2).match == "1234"
    @test match(possessive_between(1, 4, DIGIT) * "4", str2) === nothing

    @test match(at_least(2, DIGIT) * "4", str2).match == "1234"
    @test match(possessive_at_least(2, DIGIT) * "4", str2) === nothing

    # don't quite understand why this doesn't work, leave it out for now
    # str3 = "123"

    # @test match(maybe("123") * "3", str3).match == "3"
    # @test match(possessive_maybe("123") * "3", str3) === nothing
end

max_stars_repo_name: jkrumbiegel/ReadableRegex

max_issues_repo_name: jkrumbiegel/ReadableRegex

max_forks_repo_name: jkrumbiegel/ReadableRegex


text:
### A Pluto.jl notebook ###
# v0.12.15

using Markdown
using InteractiveUtils

# ╔═╡ 2bf8e4c2-339b-11eb-1d49-7f1b210a650a
using DelimitedFiles

# ╔═╡ 3e0c149a-3350-11eb-0b17-7bea5d1a28d4
md"# 2020 Day 1"

# ╔═╡ 18080fe2-339b-11eb-34f0-f393638506fc
md"## Part 1"

# ╔═╡ 4e084a4e-339b-11eb-3ffc-41118215f226
md"Find two elements that sum to 2020 and multiply them"

# ╔═╡ 24e1613c-339b-11eb-2976-59884053d1ee
input = readdlm("../data/day_1.txt", ' ', Int, '\n');

# ╔═╡ 67af36ec-339b-11eb-2bea-b1b1d06816d1
begin
	test_txt = """1721
	979
	366
	299
	675
	1456"""
	test_input = [parse(Int, x) for x in split(test_txt,"\n")]
end

# ╔═╡ c6762a40-339c-11eb-2313-a1266d053a73
17221 in test_input

# ╔═╡ bf6397c0-339b-11eb-346e-758d99ec2a5e
function find_sum_pairs(input::Array{Int},total::Int)
	for first_number in input
		second_number = total - first_number
		if second_number in input
			return [first_number, second_number]
		end
	end
end

# ╔═╡ d25fa824-339c-11eb-3cc2-75f419dc4740
begin
	test_pair = find_sum_pairs(test_input,2020)
	product_of_test_pair = test_pair[1] * test_pair[2]
	test_pair, product_of_test_pair
end

# ╔═╡ 0f8198aa-339d-11eb-0ae3-efa5f999b4f3
begin
	pair = find_sum_pairs(input,2020)
	product_of_pair = pair[1]*pair[2]
	pair, product_of_pair
end

# ╔═╡ 5f797e22-339d-11eb-1ff3-dda90f946838
md"## Part 2"

# ╔═╡ e3177eae-339e-11eb-3f9f-1f7230976832
md"Find the product of the triplet that sums to 2020"

# ╔═╡ 6b4a06e0-339d-11eb-0168-73ec316f0508
function find_sum_triplet(input::Array{Int},total::Int)
	for first_number in input
		remainder = total - first_number
		pair = find_sum_pairs(input, remainder)
		if pair != nothing
			return [first_number, pair[1], pair[2]]
		end
	end
end

# ╔═╡ 2fbb1bc0-339e-11eb-0917-376ee88ab59e
begin
	test_triplet = find_sum_triplet(test_input,2020)
	product_of_test_triplet = prod(test_triplet)
	test_triplet, product_of_test_triplet
end

# ╔═╡ 493f1fbc-339e-11eb-1d12-d1946e3f2f63
begin
	triplet = find_sum_triplet(input,2020)
	product_of_triplet = prod(triplet)
	triplet, product_of_triplet
end

# ╔═╡ Cell order:
# ╟─3e0c149a-3350-11eb-0b17-7bea5d1a28d4
# ╠═2bf8e4c2-339b-11eb-1d49-7f1b210a650a
# ╟─18080fe2-339b-11eb-34f0-f393638506fc
# ╟─4e084a4e-339b-11eb-3ffc-41118215f226
# ╠═24e1613c-339b-11eb-2976-59884053d1ee
# ╠═67af36ec-339b-11eb-2bea-b1b1d06816d1
# ╠═c6762a40-339c-11eb-2313-a1266d053a73
# ╠═bf6397c0-339b-11eb-346e-758d99ec2a5e
# ╠═d25fa824-339c-11eb-3cc2-75f419dc4740
# ╠═0f8198aa-339d-11eb-0ae3-efa5f999b4f3
# ╟─5f797e22-339d-11eb-1ff3-dda90f946838
# ╟─e3177eae-339e-11eb-3f9f-1f7230976832
# ╠═6b4a06e0-339d-11eb-0168-73ec316f0508
# ╠═2fbb1bc0-339e-11eb-0917-376ee88ab59e
# ╠═493f1fbc-339e-11eb-1d12-d1946e3f2f63

max_stars_repo_name: RensDimmendaal/advent-of-code-julia

max_issues_repo_name: RensDimmendaal/advent-of-code-julia

max_forks_repo_name: RensDimmendaal/advent-of-code-julia


text:
include("SLPDE.jl")
include("etdrk4.jl")
include("evaluate_lambda.jl")
include("sphcesaro.jl")

using PyPlot

srand(0)

n = 128
U0 = zeros(4n, 8n-1);
U0[1:n, 1:2n-1] = sphrandn(Float64, n, n)/n
n = size(U0, 1)

θ = [0;(0.5:n-0.5)/n;1]
φ = [(0:2n-2)*2/(2n-1);2]
x = [cospi(φ)*sinpi(θ) for θ in θ, φ in φ]
y = [sinpi(φ)*sinpi(θ) for θ in θ, φ in φ]
z = [cospi(θ) for θ in θ, φ in φ]

V = sph2fourier(U0)
F = zero(V)
A_mul_B!(F, FastTransforms.plan_synthesis(V), V)
F = [mean(F[1,:])*ones(1, size(F, 2)); F; mean(F[end,:])*ones(1, size(F, 2))]
F = [F F[:,1]]

function ETDRK4_free_energy(SLPDE::SemiLinearPDE{1, S}, T::Number, n::Int) where S
    U0 = SLPDE.U0[1]
    U = deepcopy(U0) # solution spherical harmonic coefficients
    V = zero(U) # solution Fourier coefficients
    F = zero(U) # solution function values
    NU = zero(U) # nonlinearity spherical harmonic coefficients
    NV = zero(U) # nonlinearity Fourier coefficients
    NF = zero(U) # nonlinearity function values

    # ETDRK4 stages
    A = zero(U)
    B = zero(U)
    C = zero(U)
    NA = zero(U)
    NB = zero(U)
    NC = zero(U)

    P = plan_sph2fourier(U; sketch = :none)
    Ps = FastTransforms.plan_synthesis(U)
    Pa = FastTransforms.plan_analysis(U)

    L = create_linear_operator(SLPDE.L, size(U, 1), size(U, 2)÷2 + 1)
    h = S(T)/n

    NE = NonlinearOperator(u -> 0.5*(u^2-1))

    ENERGY = zeros(S, n+1)
    ENERGY[1] = -0.5*sum(U.*L.*U) + sum(evaluate_nonlinear_operator!(NE, P, Ps, Pa, U, V, F, NF, NV, NU).^2)

    # Create all the matrices of diagonal scalings used in ETDRK4. These operate
    # component-wise on the spherical harmonic coefficients.

    Z = L*h
    ez = exp.(Z)
    ez2 = exp.(Z./2)
    h2ez2m1 = (h/2.0).*φ1.(Z./2)
    hezα = h.*expα.(Z)
    h2ezβ = 2h.*expβ.(Z)
    hezγ = h.*expγ.(Z)
    ez2u = zero(Z)

    for k=1:n
        ez2u .= ez2.*U
        evaluate_nonlinear_operator!(SLPDE.N, P, Ps, Pa, U, V, F, NF, NV, NU)
        A .= ez2u .+ h2ez2m1.*NU
        evaluate_nonlinear_operator!(SLPDE.N, P, Ps, Pa, A, V, F, NF, NV, NA)
        B .= ez2u .+ h2ez2m1.*NA
        evaluate_nonlinear_operator!(SLPDE.N, P, Ps, Pa, B, V, F, NF, NV, NB)
        C .= ez2.*A .+ h2ez2m1.*(2.0.*NB .- NU)
        evaluate_nonlinear_operator!(SLPDE.N, P, Ps, Pa, C, V, F, NF, NV, NC)
        U .= ez.*U .+ hezα.*NU .+ h2ezβ.*(NA .+ NB) .+ hezγ.*NC
        ENERGY[k+1] = -0.5*sum(U.*L.*U) + sum(evaluate_nonlinear_operator!(NE, P, Ps, Pa, U, V, F, NF, NV, NU).^2)
    end

    ENERGY
end

T = 64.0
n = round(Int, 10*T)
t = linspace(0, T, n+1)

SLPDE = SemiLinearPDE(LaplaceBeltrami(0.01), NonlinearOperator(u->u-u^3), U0)
local_energy = ETDRK4_free_energy(SLPDE, T, n)

SLPDE = SemiLinearPDE(NonlocalLaplaceBeltrami(0.01, -0.5, 1.0), NonlinearOperator(u->u-u^3), U0)
nonlocal_energy = ETDRK4_free_energy(SLPDE, T, n)

clf()
plot(t, local_energy, "-k", t, nonlocal_energy, "--k")
xlabel("\$t\$"); ylabel("\${\\cal E}(u)\$"); grid()
legend(["Local","Nonlocal"])
savefig("GL_free_energy.pdf")


text:
@testset "Examples" begin

    @testset "Quick reference example" begin
        # Input type: 3 channels 2D image
        invertex = inputvertex("in", 3, FluxConv{2}())

        # Mutable layers
        conv = mutable(Conv((3,3), 3 => 5, pad=(1,1)), invertex)
        batchnorm = mutable(BatchNorm(nout(conv), relu), conv)

        # Explore graph
        @test inputs(conv) == [invertex]
        @test outputs(conv) == [batchnorm]

        @test nin(conv) == [3]
        @test nout(conv) == 5

        @test layertype(conv) isa FluxConv{2}
        @test layertype(batchnorm) isa FluxBatchNorm

        # naming vertices is a good idea for debugging and logging purposes
        namedconv = mutable("namedconv", Conv((5,5), 3=>7, pad=(2,2)), invertex)

        @test name(namedconv) == "namedconv"

        # Concatenate activations
        conc = concat("conc", namedconv, batchnorm)

        @test nout(conc) == nout(namedconv) + nout(batchnorm)

        residualconv = mutable("residualconv", Conv((3,3), nout(conc) => nout(conc), pad=(1,1)), conc)

        # Elementwise addition. '>>' operation can be used to add metadata, such as a name in this case
        add = "add" >> conc + residualconv

        @test name(add) == "add"
        @test inputs(add) == [conc, residualconv]

        # Computation graph for evaluation
        graph = CompGraph(invertex, add)

        # Can be evaluated just like any function
        x = ones(Float32, 7, 7, nout(invertex), 2)
        @test size(graph(x)) == (7, 7, nout(add), 2) == (7 ,7, 12 ,2)

        # Graphs can be copied
        graphcopy = copy(graph)

        # Mutate number of neurons
        @test nout(add) == nout(residualconv) == nout(conv) + nout(namedconv) == 12
        Δnout(add, -3)
        @test nout(add) == nout(residualconv) == nout(conv) + nout(namedconv) == 9

        # Remove layer
        @test nv(graph) == 7
        remove!(batchnorm)
        @test nv(graph) == 6

        # Add layer
        insert!(residualconv, v -> mutable(BatchNorm(nout(v), relu), v))
        @test nv(graph) == 7

        # Change kernel size (and supply new padding)
        let Δsize = (-2, -2), pad = (1,1)
            mutate_weights(namedconv, KernelSizeAligned(Δsize, pad))
        end

        # Apply all mutations
        apply_mutation(graph)

        # Note: Parameters not changed yet...
        size(NaiveNASflux.weights(layer(namedconv))) == (5, 5, 3, 7)

        @test size(graph(x)) == (7, 7, nout(add), 2) == (7, 7, 9, 2)

        # ... because mutations are lazy by default so that no new layers are created until the graph is evaluated
        size(NaiveNASflux.weights(layer(namedconv))) == (3, 3, 3, 7)

        # Btw, the copy we made above is of course unaffected
        @test size(graphcopy(x)) == (7, 7, 12, 2)

    end

    @testset "Pruning xor example" begin
        using NaiveNASflux, Test
        import Flux: train!, mse
        import Random
        Random.seed!(666)
        niters = 100

        # First lets create a simple model
        # layerfun=ActivationContribution will wrap the layer and compute a pruning metric for it while the model trains
        mdense(in, outsize, act) = mutable(Dense(nout(in),outsize, act), in, layerfun=ActivationContribution)

        invertex = inputvertex("input", 2, FluxDense())
        layer1 = mdense(invertex, 32, relu)
        layer2 = mdense(layer1, 1, sigmoid)
        original = CompGraph(invertex, layer2)

        # Training params, nothing to see here
        opt = ADAM(0.1)
        loss(g) = (x, y) -> mse(g(x), y)

        # Training data: xor truth table: y = xor(x)
        x = Float32[0 0 1 1;
                    0 1 0 1]
        y = Float32[0 1 1 0]

        # Train the model
        for iter in 1:niters
            train!(loss(original), params(original), [(x,y)], opt)
        end
        @test loss(original)(x, y) < 0.001

        # Now, lets try three different ways to prune the network
        nprune = 16

        # Prune randomly selected neurons
        pruned_random = copy(original)
        Δnin(pruned_random.outputs[], -nprune)
        Δoutputs(pruned_random.outputs[], v -> rand(nout_org(v)))
        apply_mutation(pruned_random)

        # Prune the least valuable neurons according to the metric in ActivationContribution
        pruned_least = copy(original)
        Δnin(pruned_least.outputs[], -nprune)
        Δoutputs(pruned_least.outputs[], neuron_value)
        apply_mutation(pruned_least)

        # Prune the most valuable neurons according to the metric in ActivationContribution
        pruned_most = copy(original)
        Δnin(pruned_most.outputs[], -nprune)
        Δoutputs(pruned_most.outputs[], v -> -neuron_value(v))
        apply_mutation(pruned_most)

        # Can I have my free lunch now please?!
        @test loss(pruned_most)(x, y) > loss(pruned_random)(x, y) > loss(pruned_least)(x, y) > loss(original)(x, y)

        # The metric calculated by ActivationContribution is actually quite good (in this case).
        @test loss(pruned_least)(x, y) ≈ loss(original)(x, y) atol = 1e-5
    end

    @testset "Conv 2D xor example" begin
        using NaiveNASflux, Test
        import Flux: train!, logitbinarycrossentropy, glorot_uniform
        using Statistics
        import Random
        Random.seed!(666)
        niters = 20

        # Layers used in this example
        mconv(in, outsize, act; init=glorot_uniform) = mutable(Conv((3,3),nout(in)=>outsize, act, pad=(1,1), init=init), in)
        mavgpool(in, h, w) = mutable(MeanPool((h, w)), in)
        mdense(in, outsize, act) = mutable(Dense(nout(in),outsize, act), in)

        # Size of the input
        height = 4
        width = 4

        function model(nconv)
            invertex = inputvertex("in", 1, FluxConv{2}())
            l = invertex
            for i in 1:nconv
                l = mconv(l, 16, relu)
            end
            l = mavgpool(l, height, width)
            l = invariantvertex(x -> x[1,1,:,:], l)
            l = mdense(l, 2, identity)
            return CompGraph(invertex, l)
        end
        original = model(1)

        # Training params, nothing to see here
        opt_org = ADAM(0.01)
        loss(g) = (x, y) -> mean(logitbinarycrossentropy.(g(x), y))

        # Training data: 2D xor(-ish)
        # Class 1 are matrices A where A[n,m] ⊻ A[n+h÷2, m]) ⩓ A[n,m] ⊻ A[n, m+w÷2] ∀ n<h÷2, m<w÷2 is true, e.g. [1 0; 0, 1]
        function xy1(h,w)
            q1_3 = rand(Bool,h÷2,w÷2)
            q2_4 = .!q1_3
            return (x = Float32.(vcat(hcat(q2_4, q1_3), hcat(q1_3, q2_4))), y = Float32[0, 1])
        end
        xy0(h,w) = (x = rand(Float32[0,1], h,w), y = Float32[1, 0]) #Joke's on me when this generates false negatives :)
        # Generate 50% class 1 and 50% class 0 examples in one batch
        batch(h,w,batchsize) = mapfoldl(i -> i==0 ? xy0(h,w) : xy1(h,w), (ex1, ex2) -> (x = cat(ex1.x, ex2.x, dims=4), y = hcat(ex1.y,ex2.y)), (1:batchsize) .% 2)

        x_test, y_test = batch(height, width, 1024)
        startloss = loss(original)(x_test, y_test)

        # Train the model
        for iter in 1:niters
            train!(loss(original), params(original), [batch(height, width, 64)], opt_org)
        end
        # That didn't work so well...
        @test loss(original)(x_test, y_test) ≈ startloss atol=1e-1

        # Lets try three things:
        # 1. Just train the same model some more
        # 2. Add two more conv-layers to the already trained model
        # 3. Create a new model with three conv layers from scratch

        # Disclaimer: This experiment is intended to show usage of this library.
        # It is not meant to give evidence that method 2 is the better option.
        # Hyperparameters are tuned to strongly favor 2 in order to avoid sporadic failures

        # Add two layers after the conv layer
        add_layers = copy(original)
        function add2conv(in)
            l = mconv(in, nout(in), relu, init=idmapping)
            return mconv(l, nout(in), relu, init=idmapping)
        end
        insert!(vertices(add_layers)[2], add2conv)

        # New layers are initialized to identity mapping weights
        # We basically have the same model as before, just with more potential
        # Not guaranteed to be a good idea as it relies on existing layers to provide gradient diversity
        @test add_layers(x_test) == original(x_test)

        # Create a new model with three conv layers
        new_model = model(3)

        opt_add = deepcopy(opt_org)
        opt_new = ADAM(0.01)

        # Lets try again
        for iter in 1:niters
            b = batch(height, width, 64)
            train!(loss(original), params(original), [b], opt_org)
            train!(loss(add_layers), params(add_layers), [b], opt_add)
            train!(loss(new_model), params(new_model), [b], opt_new)
        end

        @test loss(add_layers)(x_test,y_test) < loss(new_model)(x_test,y_test) < loss(original)(x_test,y_test)
    end
end


text:
function main()
  #ch=Channel(1)
	f1() = (for i = 1:5 println("a"); put!(ch,1); sleep(1); end);
	f2() = (for i = 1:5 println("b: ", take!(ch)); sleep(0.1); end);
  f3() = (while true; println("c"); for i = 1:9999999 b=i^i; end; sleep(0); end);
  f4() = (while true; print("."); sleep(0.1); end);
	a = Task(f1);
	b = Task(f2);
  c = Task(f3);
  d = Task(f4);
	#schedule(a);
	#schedule(b);
  schedule(c);
  schedule(d);
	#yield()
  wait(c);
end

main()

max_stars_repo_name: Gilga/JuliaOpenGL

max_issues_repo_name: Gilga/JuliaOpenGL

max_forks_repo_name: Gilga/JuliaOpenGL


text:
module MatrixDepot
using GZip

import Base: show, search

export matrixdepot, @addgroup, @rmgroup

include("common.jl")        # main functions
include("higham.jl")          # test matrices
include("regu.jl")               # regularization test problem
include("graph.jl")             # adjacency matrices for graphs
include("data.jl")               # matrix data
include("download.jl")      # download data from the UF sparse matrix collection
include("datareader.jl")    # read matrix data
include("matrixmarket.jl")

const MY_DEPOT_DIR = joinpath(dirname(@__FILE__), "..", "myMatrixDepot")

if !isdir(MY_DEPOT_DIR)
    mkdir(MY_DEPOT_DIR)
    open(string(MY_DEPOT_DIR, "/group.jl"), "w") do f
        write(f, "usermatrixclass = \n Dict( \n \n \n );")
    end
    open(string(MY_DEPOT_DIR, "/generator.jl"), "w") do f
        write(f, "# include your matrix generators below ")
    end
end

files = Set(readdir(MY_DEPOT_DIR))
delete!(files, "generator.jl")
if isdir(MY_DEPOT_DIR)
    for file in files
        if split(file, '.')[2] == "jl"
            include("$(MY_DEPOT_DIR)/$(file)")
        end
    end
    include(string(MY_DEPOT_DIR, "/generator.jl"))
end

end # end module


text:
####################################################
# huawei 's  .conf wnejian
#
#####################################

Base.MainInclude.include("./ToPTformat/01-Apphw2PTFileTransfer.jl")

#########################################################################################
# 寻找制定文件夹下所有的原厂格式的原始数据文件，并转换成对应的RAW1~RAW4文件；
#　生成的文件名中的主机编号为　ＭＡｃｈｉｎｅＩD。
###########################
# rootdir = "/media/wang/705fc396-8c76-4812-9d0b-d17382d9dfc7/backup/t/app/fiberproc/sjz004/高邑元氏20191126/"
# machineID ="GYYS-001"

srcfile=ENV["splitedFileRootName"]
rootdir =  dirname(srcfile)
# machineID ="GYYS-001"
machineID =ENV["machineID"]
# dirname(srcfile)
startTime = ""
endTime = ""

#cfgfile = "../Jl-Aux/glbCfg/hw2PTCfg-BD-001.jl"
Main.julia_main(["",])

max_stars_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_issues_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_forks_repo_name: wangzhongtian/fiberAlarmAnalyseTool


text:
## TODO: add more types...
abstract type SymbolicObject <: Number end

struct Sym <: SymbolicObject
    py::Py
end

PythonCall.getpy(x::Sym) = x.py
PythonCall.ispy(x::Sym) = true

function Base.show(io::IO, s::SymbolicObject)
    out = PythonCall.pyrepr(String, getpy(s))
    out = replace(out, r"\*\*" => "^")
    print(io, out)
end

# allow dot access to underlying methods, properties
function Base.getproperty(x::Sym, a::Symbol)

    a == :py && return getfield(x,a)

    py = getpy(x)
    val = pygetattr(py, string(a), nothing)
    val === nothing && return nothing

    # hacky, but seems necessary for sympy.Poly, sympy.Curve, and others
    contains(string(pytype(val)), "ManagedProperties") && return (args...;kwargs...) -> ↑(val(unSym.(args)...; unSymkwargs(kwargs)...))

    meth = pygetattr(val, "__call__", nothing)
    meth == nothing && return Sym(val)

    return SymbolicCallable(meth)

end

# Call
(ex::Sym)(xs...) = subs(ex, xs...)

Sym(x::AbstractString) = symbols(x)
Sym(x::Symbol) = Sym(string(x))
Sym(x::Sym) = x
Sym(x::Complex{Bool}) = 1*real(x) + 1*imag(x)*IM
Sym(x::Complex) = Sym(real(x)) + Sym(imag(x))*IM
Sym(x::Number) = sympy.sympify(x) #Sym(sympy.pyconvert(Py, x))
Sym(x::PyArray) = Sym.(x)
Sym(x::PyDict) = Dict(Sym(k) => Sym(v) for (k,v) ∈ x)
Sym(x::PyIterable) = (Sym(u) for u ∈ x)
Sym(x::PyList) = Sym.(x)
Sym(x::PySet) = Set(Sym.(x))
#Sym(x::PyTable) =

Sym(x::Irrational{:π}) = sympy.pi
Sym(x::Irrational{:ℯ}) = sympy.E
Sym(x::Irrational{:e}) = sympy.E

struct SymbolicCallable
    val
end
function (v::SymbolicCallable)(args...; kwargs...)
    ↑(v.val(unSym.(args)...; unSymkwargs(kwargs)...))
end

struct SymFunction <: SymbolicObject
    py::Py
end
PythonCall.getpy(u::SymFunction) = u.py

function SymFunction(x::AbstractString; kwargs...)
    λ = __sympy__.Function(x; unSymkwargs(kwargs)...)
    SymFunction(λ)
end
(u::SymFunction)(xs...) = ↑(getpy(u)(↓(xs)...))

# Steal this idea from ModelingToolkit
# better than the **hacky** f'(0) stuff
"""
    Differential(x)

Use to find (partial) derivatives.

## Example
```
@syms x y u()
Dx = Differential(x)
Dx(u(x,y))  # resolves to diff(u(x,y),x)
Dx(u)       # will evaluate diff(u(x), x)
```
"""
struct Differential
    x::Sym
end
(∂::Differential)(u::Sym) = diff(u, ∂.x)
(∂::Differential)(u::SymFunction) = diff(u(∂.x), ∂.x)

max_stars_repo_name: jverzani/SymPyCall

max_issues_repo_name: jverzani/SymPyCall

max_forks_repo_name: jverzani/SymPyCall


text:
"""
    networks(dataset::MangalDataset, query::Pair...)

Returns networks that are part of a `MangalDataset`. Allows additional query parameters.
"""
function networks(dataset::MangalDataset, query::Pair...)
    return networks(Pair("dataset_id", dataset.id), query...)
end

"""
    network(name::AbstractString)

Returns a network of a given name.
"""
function network(name::AbstractString)
    return first(networks(Pair("name", name)))
end


text:
Ne = 800;
Ni = 200
adparam = SNN.ADEXParameter(;a = 6.0,
    b = 7.3,
    cm = 803.1,
    v_rest= -63.2,
    tau_m=19.7,
    tau_w=351.0,
    v_thresh=-39.2,
    delta_T=6.4,
    v_spike = -36.6,
    v_reset = -59.2,
    spike_delta = 16.3)

E = SNN.AD(;N = Ne, param=adparam)
I = SNN.AD(;N = Ni, param=adparam)

EE = SNN.SpikingSynapse(E, E, :v; σ = 0.5,  p = 0.8)
EI = SNN.SpikingSynapse(E, I, :v; σ = 0.5,  p = 0.8)
IE = SNN.SpikingSynapse(I, E, :v; σ = -1.0, p = 0.8)
II = SNN.SpikingSynapse(I, I, :v; σ = -1.0, p = 0.8)
P = [E, I]
C = [EE, EI, IE, II]

SNN.monitor([E, I], [:fire])
for t = 1:1000
    E.I = 795.6randn(Ne)
    I.I .= 200randn(Ni)
    SNN.sim!(P, C, 1ms)
end

max_stars_repo_name: russelljjarvis/backupSNN

max_issues_repo_name: russelljjarvis/backupSNN


text:
using LinearAlgebra: AbstractTriangular

# Matrix wrapper types that we know are square and are thus potentially invertible. For
# these we can use simpler definitions for `/` and `\`.
const SquareMatrix{T} = Union{Diagonal{T},AbstractTriangular{T}}

#####
##### `dot`
#####

function frule(::typeof(dot), x, y)
    return dot(x, y), Rule((Δx, Δy) -> sum(Δx * cast(y)) + sum(cast(x) * Δy))
end

function rrule(::typeof(dot), x, y)
    return dot(x, y), (Rule(ΔΩ -> ΔΩ * cast(y)), Rule(ΔΩ -> cast(x) * ΔΩ))
end

#####
##### `inv`
#####

function frule(::typeof(inv), x::AbstractArray)
    Ω = inv(x)
    m = @thunk(-Ω)
    return Ω, Rule(Δx -> m * Δx * Ω)
end

function rrule(::typeof(inv), x::AbstractArray)
    Ω = inv(x)
    m = @thunk(-Ω')
    return Ω, Rule(ΔΩ -> m * ΔΩ * Ω')
end

#####
##### `det`
#####

function frule(::typeof(det), x)
    Ω, m = det(x), @thunk(inv(x))
    return Ω, Rule(Δx -> Ω * tr(extern(m * Δx)))
end

function rrule(::typeof(det), x)
    Ω, m = det(x), @thunk(inv(x)')
    return Ω, Rule(ΔΩ -> Ω * ΔΩ * m)
end

#####
##### `logdet`
#####

function frule(::typeof(logdet), x)
    Ω, m = logdet(x), @thunk(inv(x))
    return Ω, Rule(Δx -> tr(extern(m * Δx)))
end

function rrule(::typeof(logdet), x)
    Ω, m = logdet(x), @thunk(inv(x)')
    return Ω, Rule(ΔΩ -> ΔΩ * m)
end

#####
##### `trace`
#####

frule(::typeof(tr), x) = (tr(x), Rule(Δx -> tr(extern(Δx))))

rrule(::typeof(tr), x) = (tr(x), Rule(ΔΩ -> Diagonal(fill(ΔΩ, size(x, 1)))))

#####
##### `*`
#####

function rrule(::typeof(*), A::AbstractMatrix{<:Real}, B::AbstractMatrix{<:Real})
    return A * B, (Rule(Ȳ -> Ȳ * B'), Rule(Ȳ -> A' * Ȳ))
end

#####
##### `/`
#####

function rrule(::typeof(/), A::AbstractMatrix{<:Real}, B::T) where T<:SquareMatrix{<:Real}
    Y = A / B
    S = T.name.wrapper
    ∂A = Rule(Ȳ -> Ȳ / B')
    ∂B = Rule(Ȳ -> S(-Y' * (Ȳ / B')))
    return Y, (∂A, ∂B)
end

function rrule(::typeof(/), A::AbstractVecOrMat{<:Real}, B::AbstractVecOrMat{<:Real})
    Aᵀ, dA = rrule(adjoint, A)
    Bᵀ, dB = rrule(adjoint, B)
    Cᵀ, (dBᵀ, dAᵀ) = rrule(\, Bᵀ, Aᵀ)
    C, dC = rrule(adjoint, Cᵀ)
    ∂A = Rule(dA∘dAᵀ∘dC)
    ∂B = Rule(dA∘dBᵀ∘dC)
    return C, (∂A, ∂B)
end

#####
##### `\`
#####

function rrule(::typeof(\), A::T, B::AbstractVecOrMat{<:Real}) where T<:SquareMatrix{<:Real}
    Y = A \ B
    S = T.name.wrapper
    ∂A = Rule(Ȳ -> S(-(A' \ Ȳ) * Y'))
    ∂B = Rule(Ȳ -> A' \ Ȳ)
    return Y, (∂A, ∂B)
end

function rrule(::typeof(\), A::AbstractVecOrMat{<:Real}, B::AbstractVecOrMat{<:Real})
    Y = A \ B
    ∂A = Rule() do Ȳ
        B̄ = A' \ Ȳ
        Ā = -B̄ * Y'
        _add!(Ā, (B - A * Y) * B̄' / A')
        _add!(Ā, A' \ Y * (Ȳ' - B̄'A))
        Ā
    end
    ∂B = Rule(Ȳ -> A' \ Ȳ)
    return Y, (∂A, ∂B)
end

#####
##### `norm`
#####

function rrule(::typeof(norm), A::AbstractArray{<:Real}, p::Real=2)
    y = norm(A, p)
    u = y^(1-p)
    ∂A = Rule(ȳ -> ȳ .* u .* abs.(A).^p ./ A)
    ∂p = Rule(ȳ -> ȳ * (u * sum(a->abs(a)^p * log(abs(a)), A) - y * log(y)) / p)
    return y, (∂A, ∂p)
end

function rrule(::typeof(norm), x::Real, p::Real=2)
    return norm(x, p), (Rule(ȳ -> ȳ * sign(x)), Rule(_ -> zero(x)))
end


text:
# K-Means

# Import Packages

using Distances
using Distributions
using PyPlot

pygui(true)

# Randomize Data

N = 100
K = 3
D = 2

randCenter = rand(K, D) * 100
dataDistributionMode = Distributions.Normal(0, 5)
randShifting = rand(dataDistributionMode, N, D)
data = zeros(N, D + K + 1)

# PyPlot.scatter(randCenter[:, 1], randCenter[:, 2])

for ii = 1: N

    global data

    k = Int64(floor(rand() * K) + 1)
    data[ii, D + 1] = k
    data[ii, 1: D] = randCenter[k, 1: D] .+ randShifting[ii, 1: D]

end

PyPlot.scatter(data[:, 1], data[:, 2])

# Simulate Artificial Selection after Observing Data Distribution

centerDistributionMode = Distributions.Normal(0, 10)
centerShifting = rand(centerDistributionMode, K, D)
center = randCenter[:, 1: D] .+ centerShifting[:, 1: D]

PyPlot.scatter(center[:, 1], center[:, 2])

# Run

iter = 100

for jj = 1: iter

    for kk = 1: N

        global data, center

        for mm = (D + 2): (D + K + 1)

            data[kk, mm] = Distances.euclidean(data[kk, 1: D], center[mm - D - 1, 1: D])

        end

        for nn = (D + 2): (D + K + 1)

            if data[kk, nn] == minimum(data[kk, (D + 2): (D + K + 1)])

                data[kk, D + 1] = nn - D - 1

            end

        end

    end

    centerNew = zeros(K, D + 1)

    for pp = 1: N

        kind = Int64(data[pp, D + 1])

        for qq = 1: K

            if kind == qq

                centerNew[qq, 1: D] += data[pp, 1: D]
                centerNew[qq, D + 1] += 1

            end

        end

    end

    centerNew[:, 1: D] = centerNew[:, 1: D] ./ centerNew[:, D + 1]
    center[:, 1: D] = centerNew[:, 1: D]

end

PyPlot.scatter(center[:, 1], center[:, 2])

max_stars_repo_name: mizu-bai/Try-Julia

max_issues_repo_name: mizu-bai/Data-Science-Demos

max_forks_repo_name: mizu-bai/Data-Science-Demos


text:
function softthresh(beta::AbstractVector{<:Real},
                    penalty::AbstractVector{<:Real})
    [max(bt - lam, 0.0) * sign(bt) for (bt, lam) in zip(beta, penalty)]
end


text:
using Base.Test
using BenchmarkTools
using AstronomicalTime

# using Revise

using PrecessionNutation

const TENTH_MICROARCSEC = deg2rad(0.1 / 3600e6)

@printf "Testing accuracy magnitude: %6.2f (0.1 μas)\n" log10(TENTH_MICROARCSEC)

@testset "PrecessionNutation" begin

    @testset "Validating PrecessionNutation against ERFA" begin

        @testset "IAU 2006 precession and IAU 2000A_R06 nutation CIO XYs" begin

            @testset "$year" for year in [2020] # [1980, 2000, 2020, 2050, 2100]
                ep = TTEpoch("$(year)-01-01T12:00:00.000")

                x, y, s = [precession_nutation_06(ep)...]
                x_erfa, y_erfa, s_erfa = [precession_nutation_erfa(ep)...]

                @printf "Differences in Year %4d: \n" year
                @printf "       PN.jl -      ERFA =      Diff | Rel. Err. | <= 0.1 μas \n"
                @printf "X: %9.2e - %9.2e = %9.2e | %9.3f | %s \n" x x_erfa (x-x_erfa) (x-x_erfa)/x_erfa norm(x-x_erfa) <= TENTH_MICROARCSEC
                @printf "Y: %9.2e - %9.2e = %9.2e | %9.3f | %s \n" y y_erfa (y-y_erfa) (y-y_erfa)/y_erfa norm(y-y_erfa) <= TENTH_MICROARCSEC
                @printf "s: %9.2e - %9.2e = %9.2e | %9.3f | %s \n" s s_erfa (s-s_erfa) (s-s_erfa)/s_erfa norm(s-s_erfa) <= TENTH_MICROARCSEC

                @test x ≈ x_erfa rtol=0 atol=TENTH_MICROARCSEC
                @test y ≈ y_erfa rtol=0 atol=TENTH_MICROARCSEC
                @test s ≈ s_erfa rtol=0 atol=TENTH_MICROARCSEC
            end

        end

    end

    @testset "Benchmarking against ERFA" begin
        ep = TTEpoch("2018-04-17T18:31:20.921")

        println("\n", repeat("=", 20))
        println("ERFA: \n")
        markERFA = @benchmark precession_nutation_erfa($ep)

        display(markERFA)
        println("\n")

        println("\n", repeat("=", 20))
        println("Variant 1: \n")

        markPN = @benchmark precession_nutation_06($ep)
        display(markPN)
        println("\n")

        judgement = judge(median(markPN), median(markERFA), time_tolerance = 0.5, memory_tolerance = 5.00)
        display(judgement)
        println("\n")

        @test !isregression(judgement)

    end
end


text:
"""
return the times when random assets will fail
according to total sim time (N), quantity (q) of
assets and proportion. It receives also random
"""
function get_dropping_times(seed,stabilisation_period,drop_proportion,q,N)
    Random.seed!(seed)
    #events
    k = Int(round(q * drop_proportion))
    # k = 30

    time_btwn_events = Int.(ceil.(first([ rand(Exponential(1/k),k) .* (N - (2 * stabilisation_period)* 1.7)]))) # Normalise time-btwn-events over the 1.7 of the observation period (assuming that we substract stabilisation at the beginning and at the end) to make sure it falls within the simulation time.

    # @show time_btwn_events
    event_times =  stabilisation_period .+ cumsum(time_btwn_events[1:k])

    #For testing
    #event_times = [30,50]#,70]

    log_info("Dropping times are: $event_times")
    return event_times
end

"""
Triggers random failures on active nodes. This function is coupled to the model as it detects nodes that are active and make sure disruption happen to some of those nodes. If decoupled, cleaner, but then requires to know in advance what paths are used at every time, to make sure that dropping nodes are always active.
"""
function trigger_random_node_drops!(model::ABM)
    #-1 pick node to remove
    #0 on_switch event
    #1remove from network
    #2in controller: update topology and paths
    #in switch detect path/port not available and ask controller

    if model.ticks in model.dropping_times
        #get ids of nodes that are part of active flows
        active_ids =  unique(vcat([ af[3]==f_S ? [af[2]] : af[3]==f_E ? [af[1]] : [af[1],af[2]] for af in get_state(model).active_flows if af[3]!=f_SE ]...))

        #ensure there is no replacement when selecting random node to drop
        active_always_up = setdiff(active_ids,unique([ last(pair) for pair in model.dropped_nodes]))

        #pick one random node
        dpn_id = !isempty(active_always_up) ? rand(active_always_up) : rand(get_live_snes(model))

        #for testing only
        # dpns_test = [5,9]#[3,10,13]
        # dpn_id = rand([ sid for sid in dpns_test if get_state(getindex(model,sid)).up  ])

        log_info(model.ticks,"Removing ntw node: $dpn_id...")
        dpn_ag = getindex(model,dpn_id)
        drop_node!(dpn_ag,model)

        dpn_ctl_ag = getindex(model,abs(get_control_agent(dpn_id,model)))
        schedule_event!(dpn_ctl_ag,CTL_Event(2),model.ticks+1,[dpn_id])
        schedule_event!(dpn_ctl_ag,CTL_Event(1),model.ticks+1,[-1*dpn_id])

    end

end

function hard_drop_node(model)
    #-1 pick node to remove
    #0 on_switch event
    #1remove from network
    #2in controller: update topology and paths
    #in switch detect path/port not available and ask controller

    dpn_ids = [3] # dropping node id
    dpt = 80 # dropping time

    if model.ticks == dpt

        for dpn_id in dpn_ids
            for nb in all_neighbors(model.ntw_graph,get_address(dpn_id,model.ntw_graph))
                link_down!(get_eid(nb,model),dpn_id,model)
            end
            #remove
            dpn_ag = getindex(model,dpn_id)
            #kill_agent!(dpn_ag,model)
            set_down!(dpn_ag)
            delete!(model.mapping_ctl_ntw,dpn_id)
            #remove_vertices!(model.ntw_graph,[get_address(i,model) for i in dpn_ids])
            model.ntw_graph = remove_vertex(model.ntw_graph,get_address(dpn_id,model.ntw_graph))
            update_addresses_removal!(dpn_id,model)
        end

    end

end


text:
import GEMPIC: get_charge, get_x, add_charge!, get_mass
import GEMPIC: inner_product, evaluate

"""
Simulation of 1d2v Vlasov-Maxwell with simple PIC method,
periodic boundary conditions, Weibel instability.
FEM with splines, degree 3 for B and 2 for E
"""

@testset " PIC VM 1D2V " begin
    delta_t = 0.05
    n_time_steps = 5
    beta = 0.0001
    initial_bfield = :cos

    k = [[1.25]]
    α = [0.0]
    σ = [[0.2, 0.005773502691896]]
    μ = [[0.0, 0.0]]

    nx = 32
    xmin = 0.0
    xmax = 5.02654824574

    n_particles = 100000
    sampling_case = :sobol
    symmetric = true
    splitting_case = :symplectic
    spline_degree = 3

    mesh = OneDGrid(xmin, xmax, nx)

    beta_cos_k(x) = beta * cos(2π * x / (xmax - xmin))
    beta_sin_k(x) = beta * sin(2π * x / (xmax - xmin))

    degree_smoother = spline_degree

    # Distribution function
    df = CosSumGaussian{1,2}(k, α, σ, μ)

    sampler = ParticleSampler{1,2}(sampling_case, symmetric, n_particles)

    # Initialize the particles : mass and charge set to 1.0 with one weight
    particle_group = ParticleGroup{1,2}(n_particles)

    # Initialize the field solver
    maxwell_solver = Maxwell1DFEM(mesh, spline_degree)

    kernel_smoother1 = ParticleMeshCoupling1D(mesh, n_particles, spline_degree, :galerkin)

    kernel_smoother0 = ParticleMeshCoupling1D(mesh, n_particles, spline_degree, :galerkin)

    # Initialize the arrays for the spline coefficients of the fields
    efield_dofs = [zeros(Float64, nx), zeros(Float64, nx)]
    bfield_dofs = zeros(Float64, nx)

    # Set the initial fields
    rho = zeros(Float64, nx)
    efield_poisson = zeros(Float64, nx)

    propagator = HamiltonianSplitting(
        maxwell_solver,
        kernel_smoother0,
        kernel_smoother1,
        particle_group,
        efield_dofs,
        bfield_dofs,
    )

    efield_dofs_n = propagator.e_dofs

    sample!(particle_group, sampler, df, mesh)

    # efield 1 by Poisson, rho is computed inside the function
    solve_poisson!(efield_dofs[1], particle_group, kernel_smoother0, maxwell_solver, rho)

    # bfield = beta*cos(kx): Use b = M{-1}(N_i,beta*cos(kx))

    if initial_bfield == :cos
        l2projection!(bfield_dofs, maxwell_solver, beta_cos_k, degree_smoother - 1)
    else
        l2projection!(bfield_dofs, maxwell_solver, beta_sin_k, degree_smoother - 1)
    end

    thdiag = TimeHistoryDiagnostics(
        particle_group, maxwell_solver, kernel_smoother0, kernel_smoother1
    )

    write_step!(
        thdiag,
        0.0,
        degree_smoother,
        efield_dofs,
        bfield_dofs,
        efield_dofs_n,
        efield_poisson,
    )

    for j in 1:n_time_steps # loop over time

        # Strang splitting
        strang_splitting!(propagator, delta_t, 1)

        # Diagnostics
        solve_poisson!(
            efield_poisson, particle_group, kernel_smoother0, maxwell_solver, rho
        )

        write_step!(
            thdiag,
            j * delta_t,
            degree_smoother,
            efield_dofs,
            bfield_dofs,
            efield_dofs,
            efield_poisson,
        )

        @test true
    end

    # Compute final rho
    fill!(rho, 0.0)
    for i_part in 1:(particle_group.n_particles)
        xi = get_x(particle_group, i_part)[1]
        wi = get_charge(particle_group, i_part)
        add_charge!(rho, kernel_smoother0, xi, wi)
    end

    @show thdiag.data
end


text:
# Plot text labels at arbitrary coordinate points
function PlotText(text::String, x::Real, y::Real, vertical::Bool = false,
                  pixel_offset::ImVec2 = ImVec2(0,0))
    LibCImPlot.PlotText(text, x, y, vertical, pixel_offset)
end

function PlotImage(
    user_texture_id,
    bounds_min::ImPlotPoint = (0,0),
    bounds_max::ImPlotPoint = (0,0),
    uv0::ImVec2 = (0,0),
    uv1::ImVec2=(1,1),
    tint_col::ImVec4=(1,1,1,1);
    label_id::String = ""
)
    LibCImPlot.PlotImage(label_id, user_texture_id, bounds_min, bounds_max, uv0, uv1, tint_col)
end

function PlotImage(
    label_id::String,
    user_texture_id,
    bounds_min::ImPlotPoint = (0,0),
    bounds_max::ImPlotPoint = (0,0),
    uv0::ImVec2 = (0,0),
    uv1::ImVec2=(1,1),
    tint_col::ImVec4=(1,1,1,1)
)
    LibCImPlot.PlotImage(label_id, user_texture_id, bounds_min, bounds_max, uv0, uv1, tint_col)
end


text:
function HyperLouvain(H::hypergraph,kmax::Int64,Ω::IntensityFunction,maxits::Int64=100,bigInt::Bool=true;α,verbose=true,scan_order="lexical", Z0 = collect(1:length(H.D)))
    HyperLouvain(H,Ω;α=α,kmax=kmax,maxits=maxits,bigInt=bigInt,verbose=verbose,scan_order=scan_order, Z0 = Z0)
end

function HyperLouvain(H::hypergraph,Ω::IntensityFunction;α,kmax=maximum(keys(H.E)),maxits::Int64=100,bigInt::Bool=true,verbose=true,scan_order="lexical", Z0 = collect(1:length(H.D)))
    """
    Basic step Louvain algorithm: iterate through nodes and greedily move
    nodes to adjacent clusters. Does not form supernodes and does not recurse.

    H: hypergraph
    Ω: group interaction function, as constructed by ΩFromDict(D)
    bigInt::Bool: whether to convert the degree sequence to an array of BigInt when evaluating the volume term in second_term_eval(). Recommended.
    return: Z::array{Int64, 1}: array of group labels.
    """
    Hyp, w = hyperedge_formatting(H)    # hyperedge to node list
    node2edges = EdgeMap(H)             # node to hyperedge list
    n = length(H.D)

    Neighbs = NeighborList(node2edges, Hyp)  # Store neighbors of each node

    Ẑ = Z0                              # All nodes start in singleton clusters
    Clusters = Vector{Vector{Int64}}()  # Also store as cluster list
    for v = 1:n
        push!(Clusters, Vector{Int64}())
        push!(Clusters[v],v)
    end

    iter = 0
    improving = true
    changemade = false

    # Data structures for fast local changes in volumes
    r = kmax
    V, μ, M = evalSums(Ẑ, H.D, r;constants=false, bigInt=bigInt)
    C = evalConstants(r)

    cuts = evalCuts(H, Ẑ, Ω)

    # initialize memoization of edge penalties
    Pn           = [log(Ω.ω(Ω.P(collect(1:i)), α)) for i in 1:kmax]
    edge2penalty = [w[e]*Pn[length(Hyp[e])] for e in 1:length(Hyp)]

    # Main Greedy Local Move Loop
    while improving && iter < maxits

        iter += 1
        if mod(iter,1) == 0;
            if verbose println("Louvain Iteration $iter, Q = $(round(Float64(modularity(H, Ẑ, Ω;α = α)), digits = 2))") end
        end
        improving = false

        scan = scan_order == "lexical" ? (1:n) : Random.shuffle(1:n)

        for i = scan                     # visit each node in turn

            Ci_ind = Ẑ[i]               # Cluster index for node i
            Ci = Clusters[Ci_ind]       # Indices of nodes in i's cluster
            Ni = Neighbs[i]             # Get the indices of i's neighbors
            NC = unique(Ẑ[Ni])          # Clusters adjacent to i (candidate moves)

            BestZ = Ẑ[i]                # The default is do nothing
            BestImprove = 0
            V_best = 0
            μ_best = 0
            M_best = 0
            pen_best = 0

            # Check if it's better to move to a nearby cluster, Cj
            for j = 1:length(NC)
                Cj_ind = NC[j]

                # Check improvement for move from cluster i to to cluster j
                if Cj_ind == Ci_ind
                    change = 0
                else
                    ΔV, Δμ, ΔM = increments(V, μ, M, [i], Cj_ind, H.D, Ẑ)
                    voldiff = sum(Ω.ω(Ω.aggregator(p),α)*ΔM[p]*C[p] for p in keys(ΔM))

                    cdiff, Δpen = CutDiff_OneNode(Hyp,w,node2edges,Ẑ,i,Cj_ind,edge2penalty,Ω;α=α) # NOTE: need to check on this
                    change =  cdiff - voldiff
#                     println(cdiff, " ", voldiff)
                end

                # Check if this is currently the best possible greedy move
                if change - BestImprove > 1e-8
                    V_best = ΔV
                    μ_best = Δμ
                    M_best = ΔM
                    pen_best = Δpen
                    BestImprove = change
                    BestZ = Cj_ind
                    improving = true
                end
            end

            # Move i to the best new cluster, if it strictly improves modularity
            if BestImprove > 1e-8

                changemade = true
                improving = true    # we have a reason to keep iterating!

                # update increments for volume computation
                V, μ, M = addIncrements(V, μ, M, V_best, μ_best, M_best)

                # update clustering
                ci_old = Ẑ[i]
                Ẑ[i] = BestZ

                # Remove i from its old cluster and add it to its new cluster
                Clusters[ci_old] = setdiff(Clusters[ci_old],i)
                push!(Clusters[BestZ],i)

                # update map from edge ID to penalty at edge
                for eid in keys(pen_best)
                    newpenalty = pen_best[eid]
                    edge2penalty[eid] = newpenalty
                end
            end
        end
    end

    if ~changemade
        if verbose println("No nodes moved clusters") end
    end
    Ẑ, Clusters = renumber(Ẑ,Clusters)
    return Ẑ

end

function compute_moddiff(edge2part,Cts,Hyp,w,node2edges,V::Array, μ::Array, M::Dict,I_::T, t::Int64, D::Vector{Int64}, Z::Vector{Int64},C::Dict,Ω::IntensityFunction;α) where T <: Union{Int64, Vector{Int64}}
    """
    NOTE: I_ can now be a Vector{Int64} of nodes, all of which are assumed to belong in the same cluster.
    """
    # increments due to proposal
    ΔV, Δμ, ΔM = increments(V, μ, M, I_, t, D, Z)

    voldiff = sum(Ω.ω(Ω.aggregator(p),α)*ΔM[p]*C[p] for p in keys(M))

    ΔC, Δe2p = CutDiff_Many(Cts, I_, t, Z, Hyp, w, node2edges,edge2part)
    cdiff = sum(ΔC[p]*log(Ω.ω(Ω.aggregator(p), α)) for p in keys(ΔC))
    mdiff = cdiff-voldiff

    return mdiff, ΔV, Δμ, ΔM, ΔC, Δe2p
end

function SuperNodeStep(H::hypergraph,Z::Vector{Int64},kmax::Int64,Ω,maxits::Int64=100,bigInt::Bool=true;α,verbose=true)
    """
    A Louvain step, but starting with all nodes in an arbitrary initial cluster
    assignment Z. Louvain only considers moving an entire cluster at once.
    Running this code with Z = collect(1:n) is equivalent to HypergraphLouvain

    H: hypergraph
    Z: initial cluster assignment as a length n vector
    Ω: group interaction function, as constructed by ΩFromDict(D)
    bigInt::Bool: whether to convert the degree sequence to an array of BigInt when evaluating the volume term in second_term_eval(). Recommended.
    return: Z::array{Int64, 1}: array of group labels.
    """
    Hyp, w = hyperedge_formatting(H)    # hyperedge to node list
    node2edges = EdgeMap(H)             # node to hyperedge list
    n = length(H.D)

    Ẑ = renumber(Z)     # Ensure clusters go from 1 to number unique clusters
    sn = maximum(Ẑ)     # number of supernodes

    Clusters = Vector{Vector{Int64}}()  # Also store as cluster list
    for j = 1:sn
        push!(Clusters, Vector{Int64}())
    end
    for v = 1:n
        push!(Clusters[Ẑ[v]],v) # put node v in cluster Z[v]
    end
    SuperNodes = deepcopy(Clusters)

    # Store indices of supernodes that neighbor each supernode
    tic = time()
    Neighbs = SuperNeighborList(Hyp, SuperNodes,n)
    toc = time()-tic
    if verbose
        @show toc
    end

    # Note: There's a very subtle difference between the initial clustering,
    # which defines a set of supernodes, and the current clustering that
    # is being built by Louvain.
    # Supernodes don't change--they move as a unit. Clusters change.

    iter = 0
    improving = true
    changemade = false

    # Data structures for fast local changes in volumes
    r = kmax
    V, μ, M = evalSums(Ẑ, H.D, r;constants=false, bigInt=bigInt);
    C = evalConstants(r)

    # Data structures for fast local changes in cuts
    Cts = evalCuts(Ẑ,H)                    # edge partition -> no. times it appears in Z
    edge2part = Vector{Vector{Int64}}()    # edge id -> partition type
    for e = 1:length(Hyp)
        edge = Hyp[e]
        push!(edge2part,partitionize(Ẑ[edge]))
    end
    t1 = 0
    # Main Greedy Local Move Loop
    while improving && iter < maxits

        iter += 1
        if mod(iter,1) == 0;
            if verbose println("Louvain Iteration $iter") end
        end
        improving = false

        for i = 1:sn          # visit each supernode in turn

            S = SuperNodes[i]       # Indices associated with this supernode
            rep = S[1]              # A "representative" node from this supernode
            Ci_ind = Ẑ[rep]         # Cluster associated with this supernode
            Ci = Clusters[Ci_ind]   # Indices of nodes in this cluster

            Ni = Neighbs[i]         # All supernodes that neighbor this supernode

            # Get all the clusters of neighboring supernodes
            NC = Vector{Int64}()
            for ni in Ni
                push!(NC,Ẑ[SuperNodes[ni][1]])
            end
            NC = unique(NC)
            # println("$i = i, length(NC) = $(length(NC))")

            BestZ = Ci_ind          # The default is to not move the cluster
            BestImprove = 0
            V_best = 0
            μ_best = 0
            M_best = 0
            C_best = 0
            e2p_best = 0

            # Check if it's better to move to a nearby cluster, Cj
            for j = 1:length(NC)
                Cj_ind = NC[j]

                # Check improvement for move from cluster i to to cluster j
                if Cj_ind == Ci_ind
                    change = 0
                else
                    tic = time()
                    change, ΔV, Δμ, ΔM, ΔC, Δe2p = compute_moddiff(edge2part,Cts,Hyp, w, node2edges,V, μ, M, S, Cj_ind, H.D, Ẑ, C, Ω; α=α)
                    t1 += time()-tic
                end

                # Check if this is currently the best possible greedy move to make
                if change - BestImprove > 1e-8
                    V_best = ΔV
                    μ_best = Δμ
                    M_best = ΔM
                    C_best = ΔC
                    e2p_best = Δe2p
                    BestImprove = change
                    BestZ = Cj_ind
                    improving = true
                end
            end

            # Move supernode i to the best new cluster, only if it strictly improves modularity
            if BestImprove > 1e-8

                changemade = true
                improving = true # we have a reason to keep iterating!

                # update increments for volume computation
                V, μ, M = addIncrements(V, μ, M, V_best, μ_best, M_best)

                # update clustering
                ci_old = Ẑ[rep] # index of the old cluster
                Ẑ[S] .= BestZ

                # Remove all nodes from supernode i from their old cluster.
                Clusters[ci_old] = setdiff(Clusters[ci_old],S)

                # Add them to their new cluster
                append!(Clusters[BestZ],S)

                # Update map from partition type to number of those partitions
                for part in keys(C_best)
                    Δcount = C_best[part]
                    Cts[part] = get(Cts,part,0) + Δcount
                end

                # update map from edge ID to partition type
                for eid in keys(e2p_best)
                    newpart = e2p_best[eid]
                    edge2part[eid] = newpart
                end
            end
        end
    end
    Ẑ, Clusters = renumber(Ẑ,Clusters)
    if verbose
        @show t1
    end
    return Ẑ, changemade
end

function SuperNodeLouvain(H::hypergraph,kmax::Int64,Ω,maxits::Int64=100,bigInt::Bool=true;α,verbose=true,scan_order="lexical", Z0 = nothing)
    SuperNodeLouvain(H,Ω;α=α,kmax=kmax,maxits=maxits,bigInt=bigInt,verbose=verbose,scan_order=scan_order, Z0 = copy(Z0))
end

function SuperNodeLouvain(H::hypergraph,Ω::IntensityFunction;α,kmax=maximum(keys(H.E)),maxits::Int64=100,bigInt::Bool=true,verbose=true,scan_order="lexical", Z0 = collect(1:length(H.D)))
    """
    Running Louvain and then the super-node louvain steps until no more
    progress is possible

    H: hypergraph
    Ω: group interaction function, as constructed by ΩFromDict(D)
    bigInt::Bool: whether to convert the degree sequence to an array of BigInt when evaluating the volume term in second_term_eval(). Recommended.
    return: Z::array{Int64, 1}: array of group labels.
    """

    phase = 1

    if verbose println("Faster SuperNode Louvain: Phase $phase") end
    Ẑ = HyperLouvain(H,kmax,Ω;α=α,verbose=verbose,scan_order=scan_order, Z0 = copy(Z0))

    n = length(Ẑ)

    changed = true

    while changed
        phase += 1
        if verbose println("SuperNode Louvain: Phase $phase") end
        Ẑ, changed = SuperNodeStep(H,Ẑ,kmax,Ω,maxits;α=α,verbose=verbose)
    end

    return Ẑ
end

function evalCuts(Z::Array{Int64,1}, H::hypergraph)
    """
    Gets map from partition type to the number of times that partition type
    shows up in edges for the clustering Z
    """
    C = Dict{Vector{Int64}, Int64}()
    for k in keys(H.E)
        Ek = H.E[k]
        for e in keys(Ek)
            p = partitionize(Z[e])
            C[p] = get(C, p, 0) + Ek[e]
        end
    end
    return C
end

function evalCuts(H::hypergraph, Z::Array{Int64,1}, Ω::IntensityFunction)

    C = Dict(p => 0 for p in Ω.range)
    for k in keys(H.E)
        Ek = H.E[k]
        for e in keys(Ek)
            p = Ω.P(Z[e])
            C[p] = get(C, p, 0) + Ek[e]
        end
    end
    return C
end

function CutDiff_Many(C, I, t, Z, Hyp, w, node2edges,edge2part)
    """
    C: Dict{Vector{Int64}, Int64}, the dict of current cut values as would be produced by evalCuts().
        C is keyed by partition vectors p.
        C[p] is the # of edges with label partition p.
        Note: Ω is NOT used in the calculation of C
    I: Vector{Int64}, the list of nodes to move, assumed to be in the same group
    t: Int64, the new group to which to move the nodes I
    Z: Vector{Int64}, the current group assignments
    Hyp: list of nodes in each hyperedge
    w:  weight of the hyperedge (number of times this set of nodes is assigned a hyperedge)
    node2edges: node ID to list of edges the node is in.
    edge2part: for hyperedge index e, returns the current partitionize vector for e for the clustering Z
    return: ΔC, a Dict{Vector{Int64}, Int64}, where ΔC[p] is the change in C[p] caused by moving the nodes in I to group t
    """

    E_id = Vector{Int64}()  # edge IDs that invole some node from I
    for i in I
        append!(E_id,node2edges[i])
    end
    E = unique(E_id)

    # This stores proposed changes to edge partitions
    Δe2p = Dict{Int64,Vector{Int64}}()

    # Store changes in the number of partition types in the clustering
    ΔC = Dict(p => 0 for p in keys(C))

    Z_prop = copy(Z)    # Store the proposed new clustering
    Z_prop[I] .= t

    for eid in E
        e = Hyp[eid]    # get nodes in edge
        we = w[eid]     # get the edge weight

        p_old = edge2part[eid]  # instead of: pold = partitionize(Z[e])
        p_new = partitionize(Z_prop[e])
        Δe2p[eid] = p_new

        # Update number of partition types for p_old and p_new
        ΔC[p_old] -= we
        ΔC[p_new] = get(ΔC,p_new,0) + we

    end
    return ΔC, Δe2p
end

function CutDiff_OneNode(H::Vector{Vector{Int64}},w::Array{Float64,1},node2edges::Vector{Vector{Int64}},Z::Array{Int64,1}, I::Int64,J::Int64,edge2penalty::Union{Vector{Float64}, Vector{BigFloat}},Ω::IntensityFunction; α)
    """
    CutDiff: Compute change in the first term of the modularity function
    resulting from moving a node I to cluster J.
    """
    orig = Z[I]
    He   = node2edges[I]
    Δpen = Dict{Int64,Float64}()    # Change in penalty for moving a node

    obj  = 0
    Z[I] = J
    for i in He
        e        = H[i]
        z        = Z[e]                      # set of clusters
        new_pen = w[i]*log(Ω.ω(Ω.P(z),α))   # new penalty if this move is made
        Δpen[i]  = new_pen
        obj     += Δpen[i]-edge2penalty[i]
    end
    Z[I] = orig

    return obj, Δpen
end

function SuperNeighborList(Hyp, SuperNodes,n)
    """
    For an initial clustering into supernodes, get a list of other supernodes
    adjacenct to each supernode
    """
    H = elist2incidence(Hyp,n)
    m = size(H,1)
    sn = length(SuperNodes)
    I = Vector{Int64}()
    J = Vector{Int64}()

    for s = 1:sn
        # For supernode j, get indices of all other supernodes i such that
        # some node in j shares a node with some node in i
        S = SuperNodes[s]

        # this is silly and slow. But also not a bottleneck, so I'm
        # not worrying about making it faster right now
        Edges = findall(x->x>0,vec(sum(H[:,S],dims = 2)))  # all hyperedges nodes from S touch

        for e in Edges
            push!(I,s)
            push!(J,e)
        end
    end

    # binary (unweighted) supernodes-by-edges incidence matrix
    Hnew = SparseArrays.sparse(J,I,ones(length(I)),m,sn)
    Neighbs = NeighborList(Hnew, SparseArrays.sparse(Hnew'))
    return Neighbs
end

max_stars_repo_name: PhilChodrow/HypergraphModularity

max_issues_repo_name: PhilChodrow/HypergraphModularity

max_forks_repo_name: PhilChodrow/HypergraphModularity


text:
using GroupsCore
using Test
import AbstractAlgebra
include("symmetric.jl")

include("conformance_test.jl")

include("cyclic.jl")
include("infinite_cyclic.jl")

@testset "GroupsCore.jl" begin

    include("test_notsatisfied.jl")

    @testset "Cyclic(12)" begin
        G = CyclicGroup(12)
        test_Group_interface(G)
        test_GroupElement_interface(rand(G, 2)...)
    end

    @testset "InfCyclic" begin
        G = InfCyclicGroup()
        test_Group_interface(G)
        test_GroupElement_interface(rand(G, 2)...)
    end

    @testset "Symmetric(5)" begin
        G = AbstractAlgebra.SymmetricGroup(5)
        test_Group_interface(G)
        test_GroupElement_interface(rand(G, 2)...)
    end

    include("extensions.jl")
end

@testset "GroupConstructions" begin

    @testset "DirectProduct" begin
        GH =
            let G = AbstractAlgebra.SymmetricGroup(3),
                H = AbstractAlgebra.SymmetricGroup(4)

                GroupsCore.Constructions.DirectProduct(G, H)
            end
        test_Group_interface(GH)
        test_GroupElement_interface(rand(GH, 2)...)
    end

    @testset "DirectPower" begin
        GGG = GroupsCore.Constructions.DirectPower{3}(
            AbstractAlgebra.SymmetricGroup(3),
        )
        test_Group_interface(GGG)
        test_GroupElement_interface(rand(GGG, 2)...)
    end
    @testset "WreathProduct" begin
        W = let G = AbstractAlgebra.SymmetricGroup(2),
                P = AbstractAlgebra.SymmetricGroup(4)

                GroupsCore.Constructions.WreathProduct(G, P)
            end
        test_Group_interface(W)
        test_GroupElement_interface(rand(W, 2)...)
    end
end


text:
using FiniteDifferences: add_tiny

@testset "Methods" begin

    @testset "add_tiny" begin
        @test add_tiny(convert(Float64, 5)) isa Float64
        @test add_tiny(convert(Float32, 5)) isa Float32
        @test add_tiny(convert(Float16, 5)) isa Float16

        @test add_tiny(convert(Int, 5)) isa Float64
        @test add_tiny(convert(UInt, 5)) isa Float64
        @test add_tiny(convert(Bool, 1)) isa Float64
    end

    # The different approaches to approximating the gradient to try.
    methods = [forward_fdm, backward_fdm, central_fdm]

    # The different floating point types to try and the associated required relative
    # tolerance.
    types = [Float32, Float64]

    # The different functions to evaluate (.f), their first derivative at 1 (.d1),
    # and second derivative at 1 (.d2).
    foos = [
        (f=sin, d1=cos(1), d2=-sin(1)),
        (f=exp, d1=exp(1), d2=exp(1)),
        (f=abs2, d1=2, d2=2),
        (f=x -> sqrt(x + 1), d1=0.5 / sqrt(2), d2=-0.25 / 2^(3/2)),
    ]

    # Test all combinations of the above settings, i.e. differentiate all functions using
    # all methods and data types.
    @testset "foo=$(foo.f), method=$m, type=$(T)" for foo in foos, m in methods, T in types
        @testset "method-order=$order" for order in [1, 2, 3, 4, 5]
            @test m(order, 0, adapt=2)(foo.f, T(1)) isa T
            @test m(order, 0, adapt=2)(foo.f, T(1)) == T(foo.f(1))
        end

        @test m(10, 1)(foo.f, T(1)) isa T
        @test m(10, 1)(foo.f, T(1)) ≈ T(foo.d1)

        @test m(10, 2)(foo.f, T(1)) isa T
        if T == Float64
            @test m(10, 2)(foo.f, T(1)) ≈ T(foo.d2)
        end
    end

    # Integration test to ensure that Integer-output functions can be tested.
    @testset "Integer output" begin
        @test isapprox(central_fdm(5, 1)(x -> 5, 0), 0; rtol=1e-12, atol=1e-12)
    end

    @testset "Adaptation improves estimate" begin
        @test forward_fdm(5, 1, adapt=0)(log, 0.001) ≈ 997.077814
        @test forward_fdm(5, 1, adapt=1)(log, 0.001) ≈ 1000
    end

    @testset "Limiting step size" begin
        @test !isfinite(central_fdm(5, 1)(abs, 0.001, max_step=0))
        @test central_fdm(5, 1)(abs, 0.001) ≈ 1.0
    end

    @testset "Accuracy at high orders, with high adapt" begin
        # Regression test against issues with precision during computation of _coeffs
        # see https://github.com/JuliaDiff/FiniteDifferences.jl/issues/64

        @test central_fdm(9, 5, adapt=4, condition=1)(exp, 1.0) ≈ exp(1) atol=1e-7

        poly(x) = 4x^3 + 3x^2 + 2x + 1
        @test central_fdm(9, 3, adapt=4)(poly, 1.0) ≈ 24 atol=1e-11
    end

    @testset "Printing FiniteDifferenceMethods" begin
        @test sprint(show, "text/plain", central_fdm(2, 1)) == """
            FiniteDifferenceMethod:
              order of method:       2
              order of derivative:   1
              grid:                  [-1, 1]
              coefficients:          [-0.5, 0.5]
            """
    end

    @testset "_is_symmetric" begin
        # Test odd grids:
        @test FiniteDifferences._is_symmetric([2, 1, 0, 1, 2])
        @test !FiniteDifferences._is_symmetric([2, 1, 0, 3, 2])
        @test !FiniteDifferences._is_symmetric([4, 1, 0, 1, 2])

        # Test even grids:
        @test FiniteDifferences._is_symmetric([2, 1, 1, 2])
        @test !FiniteDifferences._is_symmetric([2, 1, 3, 2])
        @test !FiniteDifferences._is_symmetric([4, 1, 1, 2])

        # Test zero at centre:
        @test FiniteDifferences._is_symmetric([2, 1, 4, 1, 2])
        @test !FiniteDifferences._is_symmetric([2, 1, 4, 1, 2], centre_zero=true)
        @test FiniteDifferences._is_symmetric([2, 1, 1, 2], centre_zero=true)

        # Test negation of a half:
        @test !FiniteDifferences._is_symmetric([2, 1, -1, -2])
        @test FiniteDifferences._is_symmetric([2, 1, -1, -2], negate_half=true)
        @test FiniteDifferences._is_symmetric([2, 1, 0, -1, -2], negate_half=true)
        @test FiniteDifferences._is_symmetric([2, 1, 4, -1, -2], negate_half=true)
        @test !FiniteDifferences._is_symmetric(
            [2, 1, 4, -1, -2],
            negate_half=true,
            centre_zero=true
        )

        # Test symmetry of `central_fdm`.
        for p in 2:10
            m = central_fdm(p, 1)
            @test FiniteDifferences._is_symmetric(m)
        end

        # Test asymmetry of `forward_fdm` and `backward_fdm`.
        for p in 2:10
            for f in [forward_fdm, backward_fdm]
                m = f(p, 1)
                @test !FiniteDifferences._is_symmetric(m)
            end
        end
    end

    @testset "extrapolate_fdm" begin
        # Also test an `Integer` argument as input.
        for x in [1, 1.0]
            for f in [forward_fdm, central_fdm, backward_fdm]
                estimate, _ = extrapolate_fdm(f(4, 3), exp, x, contract=0.8)
                @test estimate ≈ exp(1.0) atol=1e-7
            end
        end
    end
end


text:
#== # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   Rotations from an Earth-Fixed Inertial (ECI) reference frame to an
#   Earth-Fixed, Earth-Centered (ECEF) reference frame.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
#
#   [1] Vallado, D. A (2013). Fundamentals of Astrodynamics and Applications.
#       Microcosm Press, Hawthorn, CA, USA.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # ==#

export rECItoECEF

"""
    function rECItoECEF([T,] ECI, ECEF, JD_UTC::Number [, eop_data])

Compute the rotation from an Earth-Centered Inertial (`ECI`) reference frame to
an Earth-Centered, Earth-Fixed (`ECEF`) reference frame at the Julian Day [UTC]
`JD_UTC`. The rotation description that will be used is given by `T`, which can
be `DCM` or `Quaternion`. The ECI frame is selected by the input `ECI` and the
`ECEF` frame is selected by the input `ECEF`. The possible values are listed
below. The model used to compute the rotation is specified by the selection of
the origin and destination frames. Currently, there are two models supported:
IAU-76/FK5 and IAU-2006 with 2010 conventions (CIO approach only).

# Rotation description

The rotations that aligns the ECI with ECEF can be described by Direction Cosine
Matrices or Quaternions. This is selected by the parameter `T`. The possible
values are:

* `DCM`: The rotation will be described by a Direction Cosine Matrix.
* `Quaternion`: The rotation will be described by a Quaternion.

If no value is specified, then it falls back to `DCM`.

# Conversion model

The model that will be used to compute the rotation is automatically inferred
given the selection of the origin and destination frames. **Notice that mixing
IAU-76/FK5 and IAU-2006/2010 frames is not supported yet.**

# ECI Frame

The ECI frame is selected by the parameter `ECI`. The possible values are:

* `TEME()`: ECI will be selected as the True Equator Mean Equinox (TEME)
            reference frame.
* `TOD()`: ECI will be selected as the True of Date (TOD).
* `MOD()`: ECI will be selected as the Mean of Date (MOD).
* `J2000()`: ECI will be selected as the J2000 reference frame.
* `GCRF()`: ECI will be selected as the Geocentric Celestial Reference Frame
            (GCRF).
* `CIRS()`: ECEF will be selected as the Celestial Intermediate Reference System
            (CIRS).

# ECEF Frame

The ECEF frame is selected by the parameter `ECEF`. The possible values are:

* `ITRF()`: ECEF will be selected as the International Terrestrial Reference
            Frame (ITRF).
* `PEF()`: ECEF will be selected as the Pseudo-Earth Fixed (PEF) reference
           frame.
* `TIRS()`: ECEF will be selected as the Terrestrial Intermediate Reference
            System (TIRS).

# EOP Data

The conversion between the frames depends on EOP Data (see `get_iers_eop` and
`read_iers_eop`). If IAU-76/FK5 model is used, then the type of `eop_data` must
be `EOPData_IAU1980`. Otherwise, if IAU-2006/2010 model is used, then the type
of `eop_data` must be `EOPData_IAU2000A`. The following table shows the
requirements for EOP data given the selected frames.

|   Model       |   ECI   |  ECEF  |    EOP Data     |
|:--------------|:--------|:-------|:----------------|
| IAU-76/FK5    | `GCRF`  | `ITRF` | EOP IAU1980     |
| IAU-76/FK5    | `J2000` | `ITRF` | EOP IAU1980     |
| IAU-76/FK5    | `MOD`   | `ITRF` | EOP IAU1980     |
| IAU-76/FK5    | `TOD`   | `ITRF` | EOP IAU1980     |
| IAU-76/FK5    | `TEME`  | `ITRF` | EOP IAU1980     |
| IAU-76/FK5    | `GCRF`  | `PEF`  | EOP IAU1980     |
| IAU-76/FK5    | `J2000` | `PEF`  | Not required¹   |
| IAU-76/FK5    | `MOD`   | `PEF`  | Not required¹   |
| IAU-76/FK5    | `TOD`   | `PEF`  | Not required¹   |
| IAU-76/FK5    | `TEME`  | `PEF`  | Not required¹   |
| IAU-2006/2010 | `CIRS`  | `ITRF` | EOP IAU2000A    |
| IAU-2006/2010 | `GCRF`  | `ITRF` | EOP IAU2000A    |
| IAU-2006/2010 | `CIRS`  | `TIRS` | Not required¹   |
| IAU-2006/2010 | `GCRF`  | `TIRS` | Not required¹ ² |

`¹`: In this case, the Julian Time UTC will be assumed equal to Julian Time UT1
to compute the Greenwich Mean Sidereal Time. This is an approximation, but
should be sufficiently accurate for some applications. Notice that, if EOP Data
is provided, the Julian Day UT1 will be accurately computed.

`²`: In this case, the terms that account for the free-core nutation and time
dependent effects of the Celestial Intermediate Pole (CIP) position with respect
to the GCRF will not be available, reducing the precision.
The conversion between the frames depends on EOP Data (see `get_iers_eop` and
`read_iers_eop`). If IAU-76/FK5 model is used, then the type of `eop_data` must
be `EOPData_IAU1980`. The following table shows the requirements for EOP data
given the selected frames.

## MOD and TOD

In this function, if EOP corrections are not provided, then MOD and TOD frames
will be computed considering the original IAU-76/FK5 theory. Otherwise, the
corrected frame will be used.

# Returns

The rotation description represented by `T` that rotates the ECI reference frame
into alignment with the ECEF reference frame.

# Examples

```julia-repl
julia> eop_IAU1980 = get_iers_eop(:IAU1980);

julia> rECItoECEF(DCM, GCRF(), ITRF(), DatetoJD(1986, 06, 19, 21, 35, 0), eop_IAU1980)
3×3 StaticArrays.SArray{Tuple{3,3},Float64,2,9}:
 -0.619267    -0.78518     -0.000797313
  0.78518     -0.619267     0.00106478
 -0.00132979   3.33492e-5   0.999999

julia> rECItoECEF(GCRF(), ITRF(), DatetoJD(1986, 06, 19, 21, 35, 0), eop_IAU1980)
3×3 StaticArrays.SArray{Tuple{3,3},Float64,2,9}:
 -0.619267    -0.78518     -0.000797313
  0.78518     -0.619267     0.00106478
 -0.00132979   3.33492e-5   0.999999

julia> rECItoECEF(J2000(), PEF(), DatetoJD(1986, 06, 19, 21, 35, 0))
3×3 StaticArrays.SArray{Tuple{3,3},Float64,2,9}:
 -0.619271    -0.785177    -0.000796885
  0.785176    -0.619272     0.00106622
 -0.00133066   3.45854e-5   0.999999

julia> rECItoECEF(J2000(), PEF(), DatetoJD(1986, 06, 19, 21, 35, 0), eop_IAU1980)
3×3 StaticArrays.SArray{Tuple{3,3},Float64,2,9}:
 -0.619267    -0.78518     -0.000796879
  0.78518     -0.619267     0.00106623
 -0.00133066   3.45854e-5   0.999999

julia> rECItoECEF(Quaternion, GCRF(), ITRF(), DatetoJD(1986, 06, 19, 21, 35, 0), eop_IAU1980)
Quaternion{Float64}:
  + 0.4363098936462618 + 0.0005909969666939257.i - 0.00030510511316206974.j - 0.8997962182293519.k

julia> eop_IAU2000A = get_iers_eop(:IAU2000A);

julia> rECItoECEF(GCRF(), ITRF(), DatetoJD(1986, 06, 19, 21, 35, 0), eop_IAU2000A)
3×3 StaticArrays.SArray{Tuple{3,3},Float64,2,9}:
 -0.619267    -0.78518     -0.000797312
  0.78518     -0.619267     0.00106478
 -0.00132979   3.33502e-5   0.999999

julia> rECItoECEF(GCRF(), TIRS(), DatetoJD(1986, 06, 19, 21, 35, 0))
3×3 StaticArrays.SArray{Tuple{3,3},Float64,2,9}:
 -0.619271    -0.785177    -0.000796885
  0.785176    -0.619272     0.00106623
 -0.00133066   3.45884e-5   0.999999

julia> rECItoECEF(Quaternion, GCRF(), ITRF(), DatetoJD(1986, 06, 19, 21, 35, 0), eop_IAU2000A)
Quaternion{Float64}:
  + 0.4363098936309669 + 0.000590996988144556.i - 0.0003051056555230158.j - 0.8997962182365703.k
```
"""
@inline rECItoECEF(T_ECI::T_ECIs,
                   T_ECEF::T_ECEFs,
                   JD_UTC::Number,
                   eop_data::EOPData_IAU1980) =
    rECItoECEF(DCM, T_ECI, T_ECEF, JD_UTC, eop_data)

@inline rECItoECEF(T::T_ROT,
                   T_ECI::T_ECIs,
                   T_ECEF::T_ECEFs,
                   JD_UTC::Number,
                   eop_data::EOPData_IAU1980) =
    inv_rotation(rECEFtoECI(T, T_ECEF, T_ECI, JD_UTC, eop_data))

@inline rECItoECEF(T_ECI::T_ECIs_IAU_2006,
                   T_ECEF::T_ECEFs_IAU_2006,
                   JD_UTC::Number,
                   eop_data::EOPData_IAU2000A) =
    rECItoECEF(DCM, T_ECI, T_ECEF, JD_UTC, eop_data)

@inline rECItoECEF(T::T_ROT,
                   T_ECI::T_ECIs_IAU_2006,
                   T_ECEF::T_ECEFs_IAU_2006,
                   JD_UTC::Number,
                   eop_data::EOPData_IAU2000A) =
    inv_rotation(rECEFtoECI(T, T_ECEF, T_ECI, JD_UTC, eop_data))

# Specializations for those cases that EOP Data is not needed.
@inline rECItoECEF(T_ECI::Union{Type{Val{:J2000}},Type{Val{:TOD}},
                                Type{Val{:MOD}},Type{Val{:TEME}}},
                   T_ECEF::Type{Val{:PEF}}, JD_UTC::Number) =
    rECItoECEF(DCM, T_ECI, T_ECEF, JD_UTC)

@inline rECItoECEF(T::T_ROT,
                   T_ECI::Union{Type{Val{:J2000}},Type{Val{:TOD}},
                                Type{Val{:MOD}},Type{Val{:TEME}}},
                   T_ECEF::Type{Val{:PEF}}, JD_UTC::Number) =
    inv_rotation(rECEFtoECI(T, T_ECEF, T_ECI, JD_UTC))

@inline rECItoECEF(T_ECI::Union{Type{Val{:CIRS}},Type{Val{:GCRF}}},
                   T_ECEF::Type{Val{:TIRS}},
                   JD_UTC::Number) =
    rECItoECEF(DCM, T_ECI, T_ECEF, JD_UTC)

@inline rECItoECEF(T::T_ROT,
                   T_ECI::Union{Type{Val{:CIRS}},Type{Val{:GCRF}}},
                   T_ECEF::Type{Val{:TIRS}},
                   JD_UTC::Number) =
    inv_rotation(rECEFtoECI(T, T_ECEF, T_ECI, JD_UTC))


text:
# Welcome to your new Toolips server!
using Toolips
PUBLIC = "public"
IP = "127.0.0.1"
PORT = 8003
include("Home.jl")
include("Text_editor.jl")
function main()
        # Essentials
    global LOGGER = Logger()
    routes = []
    server_template = ServerTemplate(IP, PORT, routes , logger = LOGGER)
    server_template.add(Route("404", fn(four04)))
    server_template.add(Route("/", home))
    server_template.add(Route("/text-editor", te))
    server_template.add(Route("/text-editor/open", texteditsc))
    server_template.add(Route("/text-editor/save", fn(saver)))
    global TLSERVER = server_template.start()
    return(TLSERVER)
end
four04 = http -> write(http, "
404, not found
")
main()

max_issues_repo_name: emmettgb/ToolipsApp

max_forks_repo_name: emmettgb/ToolipsApp