text:
## Comrade Matrix
## The comrade matrix plays the role of the companion matrix for static basis polynomials

#
# generate the  matrix for an orhtogonal polynomial family
# should have eigvals(comrade_matrix(p)) ≈ roots(convert(Polynomial, p))
# This isn't exported
function comrade_matrix(p::P) where {P <: AbstractCCOP}
    U,V = comrade_pencil(p)
    U * inv(V)
end

# Generate the [comrade pencil](http://www.cecm.sfu.ca/personal/pborwein/MITACS/papers/LinMatPol.pdf)
# C₀, C₁ with λ C₁ - C₀ the linearization; and C₀ ⋅ C₁⁻¹ the companion matrix
# can be used with square matrix coefficients
function comrade_pencil(p::P, symmetrize=false) where {T, P <: AbstractCCOP{T}}
    comrade_pencil(p, αβγk(p)..., symmetrize)
end

function αβγk(p::P) where {T, P <: AbstractCCOP{T}}
    n = Polynomials.degree(p)
    𝐏 = _comrade_pencil_type(p)

    As = An.(𝐏, 0:n-1)
    Bs = Bn.(𝐏, 0:n-1)
    Cs = Cn.(𝐏, 1:n-1)
    kₙ, kₙ₋₁ = leading_term(𝐏{T}, n), leading_term(𝐏{T},n-1)

    # α = 1/A
    # β = -B/A
    # γ = C/A

    α = inv.(As[1:end-1])
    β = (-Bs ./ As)
    γ = (Cs ./ As[2:end])

    (α,β,γ,kₙ₋₁, kₙ)
end
_comrade_pencil_type(p::P) where {T, P <: AbstractCCOP{T}} =  ⟒(P){T}
_comrade_pencil_type(p::P) where {T, M<: AbstractMatrix{T}, P <: AbstractCCOP{M}} = ⟒(P){T}

# α = [α₀, α₁, ⋯, αₙ₋₂]
# β = [β₀, β₁, ⋯, βₙ₋₂, βₙ₋₁]
# γ = [γ₁, γ₂ ⋯, αₙ₋₁]
# kₙ₋₁/kₙ = aₙ₋₁
function comrade_pencil(p::P, α, β, γ, kₙ₋₁, kₙ, symmetrize=false) where {T, P <: Polynomials.AbstractPolynomial{T}}

    n = degree(p)
    R = eltype(one(T)/1)
    C₀ = diagm(2 => zeros(R, n-2),
               1 => γ,
               0 => β,
               -1 => α
               )

    C₀[end-1,end] *= kₙ * p[end]
    C₀[end,end] *= kₙ * p[end]

    for i ∈ 0:n-1
        C₀[1+i, end] -= kₙ₋₁ * p[i]
    end

    C₁ = diagm(0 => ones(R, n))
    C₁[end,end] = kₙ * p[end]

    C₀, C₁
end

function comrade_pencil(p::P, α, β, γ, kₙ₋₁, kₙ, symmetrize=false) where {T, M <: AbstractMatrix{T}, P <: Polynomials.AbstractPolynomial{M}}

    n = Polynomials.degree(p)
    m,m′ = size(p[0])
    @assert m == m′ # square
    𝐈 = I(m)

    R = eltype(one(T)/1)
    C₀ = zeros(R, n*m, n*m)
    C₁ = diagm(0 => ones(R, n*m))
    Δ = 1:m

    C₀[0*m .+ Δ, 0*m .+ Δ] .= β[1]
    C₀[1*m .+ Δ, 0*m .+ Δ] .= α[1]

    for i ∈ 1:n-2
        C₀[(i-1)*m .+ Δ, i*m .+ Δ] .= γ[i] * 𝐈
        C₀[i*m .+ Δ, i*m .+ Δ] .= β[1+i] * 𝐈
        C₀[(i+1)*m .+ Δ, i*m .+ Δ] .= α[1+i] * 𝐈
    end

    C₀[(n-2)*m .+ Δ, (n-1)*m .+ Δ] .= γ[end] * kₙ * p[end] * 𝐈
    C₀[(n-1)*m .+ Δ, (n-1)*m .+ Δ] .= β[end] * kₙ * p[end] * 𝐈

    for i ∈ 0:n-1
        C₀[i*m .+ Δ, (n-1)*m .+ Δ] .-= kₙ₋₁ * p[i]
    end

    C₁[(n-1)*m .+ Δ, (n-1)*m .+ Δ] .=  kₙ * p[end]

    if symmetrize
        H₀ = zeros(T, n*m, n*m)
        for k ∈ 1:n
            for j ∈ 1:k
                H₀[(k-j)*m .+ Δ ,(j-1)*m .+ Δ] .= p[k]
            end
        end
        return (C₀*inv(C₁)) * H₀, H₀
    else
        return C₀, C₁
    end
end

"""
    comrade_pencil_LU(p)

Return an LU decomposition of `λC₁ - C₀`, where `C₀` and `C₁` are
the pieces of the comrade pencil.

The returned value is a function. We have `L,Ũ = comrade_pencil_LU(p)(λ)`.
The upper-triangular part, `Ũ` is not quite `U`, as `U` is singluar at eigenvalues. Rather,
this holds up to floating point

```
C₀, C₁ = comrade_pencil(p)
L,Ũ = comrade_pencil_LU(p)(λ)
@assert det(Ũ) ≈ 1
Ũ[end,end] *= p(λ)
L*Ũ - (λ*C₁ - C₀) ≈ 0 # if atol=sqrt(eps())
```
"""
comrade_pencil_LU(p::P) where {P <: AbstractCCOP} =  comrade_pencil_LU(p, αβγk(p)...)

function comrade_pencil_LU(p::P, α, β, γ, kₙ₋₁, kₙ) where {T, P <: Polynomials.AbstractPolynomial{T}}
    n = degree(p)
    R = typeof(one(T)/1)
    # return a function of λ
    λ -> begin
        𝐏 = _comrade_pencil_type(p)
        ϕs = [basis(𝐏,i)(λ) for i ∈ 0:n-1]
        L = zeros(R, n, n)
        U = zeros(R, n, n)
        for i ∈ 1:n-1
            L[i,i] = one(R)
            L[i+1,i] = -ϕs[i]/ϕs[i+1]

            U[i,i+1] = -γ[i]
            U[i,i] = α[i] * ϕs[i+1]/ϕs[i]
        end
        L[end,end] = one(R)

        U[1,end] = kₙ₋₁ * p[0]
        for i ∈ 2:(n-2)
            j = 1+i
            U[i,end] = kₙ₋₁ * p[j-2] + ϕs[j-2]/ϕs[j-1] * U[i-1,end]
        end

        U[n-1,end] = kₙ₋₁ * p[end-2] + ϕs[end-2]/ϕs[end-1]*U[n-2,end] - kₙ*γ[end]*p[end]
        U[end,end] = ϕs[1]/ϕs[end] / prod(α) * 1 # p(λ)
        (L=L, U=U)
    end
end

function comrade_pencil_LU(p::P, α, β, γ, kₙ₋₁, kₙ) where {T, M <: AbstractMatrix{T}, P <: Polynomials.AbstractPolynomial{M}}

    m, m′ = size(p[0])
    @assert m == m′
    Iₘ = I(m)
    Δ = 1:m

    n = degree(p)
    R = typeof(one(T)/1)

    # return a function of λ
    λ -> begin
        𝐏 = _comrade_pencil_type(p)
        ϕs = [basis(𝐏,i)(λ) for i ∈ 0:n-1]
        L = zeros(R, n*m, n*m)
        Ũ = zeros(R, n*m, n*m)
        for i ∈ 1:n-1
            i′ = i - 1
            L[i′*m .+ Δ, i′*m .+ Δ] .= Iₘ
            L[(i′+1)*m .+ Δ, i′*m .+ Δ] .= -ϕs[i]/ϕs[i+1] * Iₘ

            Ũ[i′*m .+ Δ,(i′+1)*m .+ Δ] .= -γ[i] * Iₘ
            Ũ[i′*m .+ Δ, i′*m .+ Δ] .= α[i] * ϕs[i+1]/ϕs[i] * Iₘ
        end
        END = (n-1)*m .+ Δ
        L[END, END] .= Iₘ

        Ũ[Δ, END] .= kₙ₋₁ * p[0]
        for i ∈ 2:(n-2)
            j = 1+i
            i′ = i - 1
            Ũ[i′*m .+ Δ, END] .= kₙ₋₁ * p[j-2] + ϕs[j-2]/ϕs[j-1] * Ũ[(i′-1)*m .+ Δ,END]
        end

        Ũ[(n-1-1)*m .+ Δ,END] = kₙ₋₁ * p[end-2] +
            ϕs[end-2]/ϕs[end-1]*Ũ[(n-2-1)*m .+ Δ,END] - kₙ*γ[end]*p[end]

        Ũ[END,END] = ϕs[1]/ϕs[end] / prod(α) * Iₘ # use Ũ, not U, which uses p(λ)
        (L=L, U=Ũ)
    end

end

## -----

## Implement algorithm for eigenvalue of a comrade matrix from:
## Fast computation of eigenvalues of companion, comrade, and related matrices
## Jared L. Aurentz · Raf Vandebril · David S. Watkins

## Overlap with AMRVW.jl, but not really, as this paper uses different transforms:
## * a general 2x2 core transformation
## * a GaussTransform of the type `[1 0; m 1]`

abstract type AbstractTransform{T} end
# Ci CoreTransform [c 1; 1 0]
struct CoreTransform{T} <: AbstractTransform{T}
    g11::T
    g12::T
    g21::T
    g22::T
end
CoreTransform(m::Matrix) = CoreTransform(m[1,1], m[1,2], m[2,1], m[2,2])
LinearAlgebra.norm(c::CoreTransform) = maximum(abs, splat(c))
Base.isnan(c::CoreTransform) = isnan(c.g12) || isnan(c.g22) || isnan(c.g21) || isnan(c.g11)

# Gi GaussTransform [1 0, m 1]
struct GaussTransform{T} <: AbstractTransform{T}
    m::T
end

splat(c::CoreTransform) = (c.g11, c.g12, c.g21, c.g22)
splat(c::GaussTransform{T}) where {T}  = (one(T), zero(T), c.m, one(T))

function Base.inv(c::CoreTransform)
    b11,b12,b21,b22 = splat(c)
    Δ = b11 * b22 - b12 * b21
    CoreTransform(b22/Δ, -b12/Δ, -b21/Δ, b11/Δ)
end

Base.inv(g::GaussTransform) = GaussTransform(-g.m)

# overwrite R
# non-allocating would be better!
function LinearAlgebra.lmul!(c::AbstractTransform, R, i)
    ri = copy(R[i,:])
    c11,c12,c21,c22 = splat(c)
    R[i,:] = R[i,:] * c11 + R[i+1,:] * c12
    R[i+1,:] = ri * c21 + R[i+1,:] * c22
    nothing
end

## ----

# Helpful for developing
function Base.show(io::IO, c::AbstractTransform)
    c11,c12,c21,c22 = splat(c)
    show(IOContext(io, :compact => true), "text/plain", [c11 c12;
                                                         c21 c22])
    println(io)
end

## realize transform as a matrix
function realize(c::AbstractTransform{T}, i, n) where {T}
    A = diagm(0 => ones(T, n))
    a,b,c,d = splat(c)
    A[i,i] = a
    A[i,i+1] = b
    A[i+1,i] = c
    A[i+1,i+1] = d
    A
end

function Base.Matrix(Cs, B)
    n = length(Cs)
    m = size(B)[1]
    for i ∈ n:-1:1
        B = realize(Cs[i], i,m)*B
    end
    B
end

## ----- primary operations with CoreTransforms

# type 1
function turnover(Ci::CoreTransform,Ci1::CoreTransform,Gi::GaussTransform)
    a11,a12,a21,a22 = splat(Ci)
    b22,b23, b32, b33 = splat(Ci1)
    g21 = Gi.m

    â11 = a11 + a12 * (b22*g21)
    â12 = a12
    â21 = a21 + a22 * (b22 * g21)
    â22 = a22

    ĝ32 = b32 * g21 / â21

    b̂22 = b22
    b̂23 = b23
    b̂32 = b32 - (ĝ32 * a22) * b22
    b̂33 = b33 - (ĝ32 * a22) * b23

    GaussTransform(ĝ32), CoreTransform(â11,â12,â21,â22), CoreTransform(b̂22,b̂23,b̂32,b̂33)

end

# type 2
function turnover(Gi::GaussTransform, Gi1::GaussTransform, Ci::CoreTransform)
    a21 = Gi.m
    b32 = Gi1.m
    g11,g12,g21,g22 = Ci.g11,Ci.g12,Ci.g21,Ci.g22

    â11 = g11
    â12 = g12
    â21 = g21 + a21 * g11
    â22 = g22 + a21 * g12

    ĝ32 = b32*g21/â21
    b̂32 = b32*g22 - ĝ32*â22

    GaussTransform(ĝ32), CoreTransform(â11, â12, â21, â22), GaussTransform(b̂32)

end

# type 3
function turnover(Gi::GaussTransform, Gi1::GaussTransform, G′i::GaussTransform)
    a21 = Gi.m
    b32 = Gi1.m
    g21 = G′i.m

    â21 = a21 + g21
    ĝ32 = b32 * g21 / â21
    b̂32 = b32 - ĝ32

    GaussTransform(ĝ32), GaussTransform(â21), GaussTransform(b̂32)
end

function fuse(c::CoreTransform, g::GaussTransform)
    a11,a12,a21,a22 = splat(c)
    m = g.m
    CoreTransform(a11 + a12 * m, a12,
                  a21 + a22 * m, a22)
end

function fuse(g::GaussTransform, c::CoreTransform)
    m = g.m
    a11,a12,a21,a22 = splat(c)

    CoreTransform(a11,         a12,
                  a11*m + a21, a12*m + a22)
end

function fuse(g1::GaussTransform, g2::GaussTransform)
    GaussTransform(g1.m + g2.m)
end

# Rlᵢ -> l′ᵢR′
# modifies R; returns l′ᵢ
function passthrough!(R, l, i)
    m = l.m

    m′ = -m * R[i+1,i+1]/(R[i,i] + R[i,i+1]*m)

    for j ∈ 1:i+1
        R[j,i] += m * R[j,i+1] # R[:,i]   += m * R[:,i+1]  # R * Gᵢ
    end
    for j ∈ i:size(R)[2]
        R[i+1, j]  += m′ * R[i,j] # R[i+1,:] += m′ * R[i,:] # Ĝᵢ * R
    end

    R[i+1,i] = 0 # may accumulate round off

    GaussTransform(-m′)

end

## ---- initial shift

# get eigenvalues of lower 2x2 matrix of A = Cs * B
function lower_eigs(Cs, R::Matrix{T}, n=length(Cs)) where {T}

    cn11, cn12, cn21, cn22 = splat(Cs[n])

    if n >= 2
        cn1_11, cn1_12, cn1_21, cn1_22 = splat(Cs[n-1])

        u11 = cn1_21 * R[n-1,n] + cn1_22 * cn11 * R[n,n]
        u12 = cn1_21*R[n-1,n+1] + cn1_22 * cn11 * R[n,n+1] + cn1_22 * cn12 * R[n+1, n+1]
        u21 = cn21 * R[n,n]
        u22 = cn21 * R[n,n+1] + cn22 * R[n+1,n+1]

    else
        u11 = cn11 * R[1,1]
        u12 = cn11 * R[1,2] + cn12 * R[2,2]
        u21 = zero(T)
        u22 = cn22 * R[2,2]
    end

    ρ1, ρ2 = eigen_values(u11, u12, u21, u22)
    ρ1, ρ2
end

# We create bulge through Â L⁻¹ A L and chase it down through as sequence of
# unitary operations of the type G⁻¹ Â G, G a Gauss transform

L(Cs, R::Matrix{T}, random::Bool=false) where {T <: Complex} = L1(Cs, R, random)
L(Cs, R::Matrix{T}, random::Bool=false) where {T}  = L1L2(Cs, R, random)

function L1(Cs, R::Matrix{T}, random::Bool=false) where {T}
    random && return GaussTransform(rand(T))

    c11,c12,c21,c22 = splat(Cs[1])
    r11 = R[1,1]
    a11, a21 = c11*r11, c21*r11

    ρs = lower_eigs(Cs, R)
    ρ1 = norm(a21 - ρs[1]) < norm(a21 - ρs[2]) ? ρs[1] : ρs[2]

    m = a21 / (a11 - ρ1)

    (isinf(m) || isnan(m)) && return GaussTransform(rand(T))

    GaussTransform(m)
end

function L1L2(Cs,R::Matrix{T}, random::Bool=false) where {T}
    random && return (GaussTransform(rand(T)), GaussTransform(rand(T)))

    ρ1, ρ2 = lower_eigs(Cs, R)

    ρ1ρ2 = real(ρ1*ρ2)
    ρ1_ρ2 = real(ρ1 + ρ2)

    cn11,cn12,cn21,cn22 = splat(Cs[2])
    cn1_11,cn1_12, cn1_21, cn1_22 = splat(Cs[1])

    a11 = cn1_11 * R[1,1]
    a12 = cn1_11 * R[1,2] + cn1_12 * cn11 * R[2,2]
    a21 = cn1_21 * R[1,1]
    a22 = cn1_21 * R[1,2] + cn1_22 * cn11 * R[2,2]
    a32 = cn21 * R[2,2]

    x1 = a11^2 +  a12 * a21 - ρ1_ρ2*a11 + ρ1ρ2
    x2 = a21 * a11 + a22 * a21 - ρ1_ρ2*a21
    x3 = a32 * a21
    m1 = x2/x1
    m2 = x3/x2

    (isinf(m1) || isinf(m2) || isnan(m1) || isnan(m2)) && return (GaussTransform(rand(T)), GaussTransform(rand(T)))

    # return l1, l2 (L = L2 * L1 though)
    GaussTransform(m1), GaussTransform(m2)

end

# single-shift bulge chase
# updates Cs, R
function chase_bulge!(Cs, R, l1::GaussTransform)
    n = length(Cs)
    Cs[1] = fuse(inv(l1), Cs[1])
    for i ∈ 1:n-1
        l1 = passthrough!(R, l1, i)
        l1, Cs[i], Cs[i+1] = turnover(Cs[i],Cs[i+1], l1)
    end
    l1 = passthrough!(R, l1, n)
    Cs[n] = fuse(Cs[n], l1)
    nothing
end

# double-shift bulge chase
function chase_bulge!(Cs, R, l1l2::Tuple)
    n = length(Cs)

    l1, l2 = l1l2
    l1′, l2′ = inv(l1), inv(l2)
    l1′, Cs[1], l2′ = turnover(l1′, l2′, Cs[1])
    Cs[2] = fuse(l2′, Cs[2])

    for i in 1:n-2

        l2 = passthrough!(R, l2, i+1)
        l2, Cs[i+1], Cs[i+2] = turnover(Cs[i+1], Cs[i+2], l2)

        l1 = passthrough!(R, l1, i)
        l1, Cs[i], Cs[i+1] = turnover(Cs[i], Cs[i+1], l1)

        l2, l1, l1′ = turnover(l1′, l2, l1)

    end
    # knit in ends
    l2 = passthrough!(R, l2, n)
    Cs[n] = fuse(Cs[n], l2)

    l1 = passthrough!(R, l1, n-1)
    l1, Cs[n-1], Cs[n] = turnover(Cs[n-1], Cs[n], l1)

    l1 = fuse(l1′, l1)
    l1 = passthrough!(R, l1, n)
    Cs[n] = fuse(Cs[n], l1)

    nothing
end

# use AVW algorithm to find eigen value of matrix
# decomposed into Cs::Vector{GaussTransform}, B::UpperTriangular
function AVW_eigvals(Cs, B)

    max_ctr = 1
    ctr = 0
    while length(Cs) > 1

        # get L # random if no deflation
        if 0 < ctr < 15
            ctr += 1
            l = L(Cs, B)
        else
            l = L(Cs, B, true) # random=true
            ctr = 1
        end

        chase_bulge!(Cs, B, l)
        deflateQ!(Cs, B) && (ctr = 0)

        # check runaway
        max_ctr += 1
        (!isempty(Cs) && isnan(Cs[end])) && break
        max_ctr > 1000 && error("Failure to converge, too many steps")

    end

    !isempty(Cs) && lmul!(Cs[1], B, 1)

    _eigvals(B)

end

# modifies Cs and B on deflation
# return true if deflation
# we deflate when off-diagonal entry < ϵ
# paper would have |aₙ₊₁,ₙ| ≤ eps() * (|aₙ,ₙ| + |aₙ₊₁,ₙ₊₁|)
# but we just use eps()^(2/3) which seems to work better for Float64
function deflateQ!(Cs, B::Matrix{T}) where {T}
    # simple case
    # is aₙ₊₁,ₙ < eps(T)(|ann| + |an+1,n+1|)
    n = length(Cs)

    c_11, c_12, c_21, c_22 = splat(Cs[n])
    c1_11, c1_12, c1_21, c1_22 = n > 1 ? splat(Cs[n-1]) : (zero(T), zero(T), zero(T), zero(T))

    Rn_1n = n > 1 ? B[n-1,n] : zero(T)

    an1_n = c_21 * B[n,n]
    ann = c1_21 * Rn_1n + c1_22 * c_11 * B[n,n]
    an1n1 = c_21 * B[n,n+1] + c_22 * B[n+1,n+1]

    ϵ = cbrt(eps(real(T))^2)

    #if norm(an1_n) <=  ϵ * (norm(ann) + norm(an1n1))
    if norm(an1_n) <=  ϵ
        n = length(Cs)
        c = pop!(Cs)
        lmul!(c, B, n)
        return true
    end

    !(T <: Real) && return false
    n == 1 && return false

    c2_11, c2_12, c2_21, c2_22 = n > 2 ? splat(Cs[n-2]) : (zero(T), zero(T), zero(T), zero(T))
    Rn_1n_1 = n > 1 ? B[n-1, n-1] : zero(T)
    Rn_2n_1 = n > 2 ? B[n-2, n-1] : zero(T)
    ann_1 = c1_21 * Rn_1n_1
    an_1n_1 = c1_11 * c2_22 * Rn_1n_1 + c2_21 * Rn_2n_1

    # if norm(ann_1) <=  eps(real(T)) * (norm(an_1n_1) + norm(ann))
    if norm(ann_1) <=  ϵ
        A = Matrix(Cs, B)
        n = length(Cs)
        c = pop!(Cs)
        lmul!(c, B, n)
        c = pop!(Cs)
        lmul!(c, B, n-1)
        return true
    end

    return false
end

# once algorithm is complete, this grabs eigenvalues from B
function _eigvals(B::Matrix{T}) where {T}
    m = size(B)[1]
    es = Vector{complex(T)}(undef, m)
    i = 1

    #ϵ = 1000 * eps(real(T))
    ϵ = cbrt(eps(real(T))^2)
    while i < m
        #if abs(B[i+1,i]) <= ϵ * (norm(B[i,i]) + norm(B[i+1,i+1]))
        if abs(B[i+1,i]) <= ϵ
            es[i] = B[i,i]
            i += 1
        else
            r1, r2 = eigen_values(B[i,i], B[i,i+1], B[i+1,i], B[i+1,i+1])
            es[i] = r1
            es[i+1] = r2
            i += 2
        end
    end

    i == m && (es[end] = B[end,end])

    LinearAlgebra.sorteig!(es)
    es
end

# one step improves accuracy
function newton_refinement(λs, p::P)  where {P <: AbstractCCOP}
    p′ = derivative(p)
    λs .- p.(λs) ./ p′.(λs)
end

## -----
## from AMRVW
# altrnative to eigvals that works with other types
function  eigen_values(a11::T, a12::T, a21::T, a22::T) where {T <: Real}
    b = (a11 + a22) / 2
    c = a11 * a22 - a12 * a21

    e1r, e1i, e2r, e2i = qdrtc(b,c) #qdrtc(one(b), b, c)
    return  complex(e1r, e1i), complex(e2r, e2i)

end

# from `modified_quadratic.f90`
function  eigen_values(a11::S, a12::S, a21::S, a22::S) where {S <: Complex}

    tr = a11 + a22
    detm = a11 * a22 - a21 * a12
    disc::S = sqrt(tr * tr - 4.0 * detm)

    u::S = abs(tr + disc) > abs(tr - disc) ? tr + disc : tr - disc
    if iszero(u)
        return zero(S), zero(S)
    else
        e1::S = u / 2.0
        e2::S = detm / e1
        return e1, e2
    end

end

## Kahan quadratic equation with fma
##  https://people.eecs.berkeley.edu/~wkahan/Qdrtcs.pdf

## solve ax^2 - 2bx + c
function qdrtc(a::T, b::T, c::T) where {T <: Real}
    # z1, z2 roots of ax^2 - 2bx + c
    d = discr(a,b,c)  # (b^2 - a*c), as 2 removes 4

    if d <= 0
        r = b/a  # real
        s = sqrt(-d)/a #imag
        return (r,s,r,-s)
    else
        r = sqrt(d) * (sign(b) + iszero(b)) + b
        return (r/a, zero(T), c/r, zero(T))
    end
end

## more work could be done here.
function discr(a::T,b::T,c::T) where {T}
    pie = 3.0 # depends on 53 or 64 bit...
    d = b*b - a*c
    e = b*b + a*c

    pie*abs(d) > e && return d

    p = b*b
    dp = muladd(b,b,-p)
    q = a*c
    dq = muladd(a,c,-q)

    (p-q) + (dp - dq)
end

function qdrtc(b::T, c::T) where {T <: Real}
    d = b*b - c
    if d <= 0
        r = b  # real
        s = sqrt(-d) #imag
        return (r,s,r,-s)
    else
        r = sqrt(d) * (sign(b) + iszero(b)) + b
        return (r, zero(T), c/r, zero(T))
    end
end

## -----
# In SpecialPolynomials we have P_{n+1} = (An*x + Bn) * P_n + Cn P_{n-1}
# In paper, we have
# αₖ pₖ = (z  - βₖ) ⋅ pₖ₋₁ - γₖ ⋅ pₖ₋₂
## Decompose p::P into Cs, B where
## * Cs -- a vector of Gauss Transfroms
## * B -- an Upper triangular matrix
## Matrix(Cs, B) is comrade matrix
function comrade_decomposition(P::Type{<:AbstractCCOP}, ps)

    n = length(ps) - 1
    As = An.(P, 0:n-1)
    Bs = Bn.(P, 0:n-1)
    Cs = Cn.(P, 1:n-1)
    kₙ, kₙ₋₁ = leading_term(P, n), leading_term(P,n-1)

    p̂s = copy(ps) ./ 1

    # αs = [1/A0, 1/A1, ..., 1/An-2]
    # βs = [B0/A0, ..., Bn-1/An-1]
    # γs = [C1/A1, ..., Cn-1/An-1]
    # cs [ 0 ... 0 γ[end] β[end]] - ps[1:end-1]/ps[end] / An-1

    aₙ = An(P,n)
    βs = -Bs ./ As
    βn = pop!(βs)
    γs = [Cs[i]/As[i+1] for i ∈ 1:length(Cs)]
    γn = pop!(γs)
    pop!(As)
    αs = 1 ./ As

    p̂n = pop!(p̂s)
    p̂s ./= p̂n
    p̂s *= -kₙ₋₁ / kₙ
    p̂s[end-1] += γn #* kₙ
    p̂s[end] += βn #* kₙ

    comrade_decomposition(p̂s, αs, βs, γs)

end

# The comrade matrix is the following (though the paper flips it for Gaussian elimination)
# [β₁ γ₂ .  .  c₀;
#  α₁ β₂ γ₃ .  c₁;
#  .  α₂ β₃ γ₄ c₂;
#  .  .  α₃ β₄ c₃;
#  .  .  .  α₄ c₄]
# returns Cs, B with eigvals(C1 ⋅ C2 ⋅ ⋯ ⋅ Cn * B) = roots(P(p))
function comrade_decomposition(cs::Vector{T}, αs, βs, γs) where {T}

    n = length(cs)
    R = zeros(T, n, n)

    CTs = Vector{CoreTransform{T}}(undef, n-1)

    ĉ₀, ĉ₁ = cs[end+1-1], cs[end+1-2]

    for i ∈ 1:n-2
        λ = ĉ₀ / αs[end+1-i]
        ĉ₀ = ĉ₁ - λ * βs[end+1-i]
        ĉ₁ = cs[end+1 - (i+2)] - λ * γs[end+1-i]
        CTs[i] = CoreTransform(λ, one(T), one(T), zero(T))
    end
    λ = ĉ₀ / αs[1]
    CTs[end] = CoreTransform(λ, one(T), one(T), zero(T))
    ĉ₀ = ĉ₁ - λ * βs[1]

    for i ∈ 1:n-2
        R[i, i]   = αs[end+1-i]
        R[i, i+1] = βs[end+1-i]
        R[i, i+2] = γs[end+1-i]
    end

    R[n-1, n-1] = αs[1]
    R[n-1, n] = βs[1]
    R[n,n] = ĉ₀

    CTs, R
end

# find roots of P(p) without converting into standard basis
# using comrade_matrix
"""
    roots(p::P) where {P <: AbstractCCOP}

When `P` is a classical orthogonal family, finds the roots of `p` using the comrade matrix; otherwise falls back to finding the roots after conversion to the `Polynomial` type.

The eigenvalue of the comrade matrix are identified using an algorithm from

*Fast computation of eigenvalues of companion, comrade, and related matrices*
by Jared L. Aurentz, Raf Vandebril, and David S. Watkins
DOI 10.1007/s10543-013-0449-x
BIT Numer Math (2014) 54:7–30

The paper presents both a single- and double-shift version. The paper
has this caveat: "The double-shift code was less successful. It too is
very fast, but it is too unstable in its current form. We cannot
recommend it for accurate computation of roots of high-degree
polynomials. It could find use as an extremely fast method to get a
rough estimate of the spectrum"

"""
function Polynomials.roots(p::P) where {P <: AbstractCCOP}
    Cs, B = comrade_decomposition(P, coeffs(p))
    λs = AVW_eigvals(Cs, B)

    # hacky check on whether roots are good enough
    ns = p.(λs) ./ derivative(p).(λs)
    δ = maximum(abs,ns)
    if δ  >= sqrt(eps())
        rts = eigvals(Polynomials.companion(p))
        if maximum(abs, p.(λs)) > maximum(abs, p.(rts))
            return complex.(rts)
        end
    end
    return newton_refinement(λs, p)
end

## -----

## test root quality
## could presumably be much more efficient
## This test backwards stability
function a_posteriori_check(λs, p::P) where {P <: AbstractCCOP}
    ps = coeffs(p)
    n = length(ps) - 1
    v(λ) = [basis(P, i)(λ) for i ∈ n-1:-1:0]
    αₙ = 1/An(P,  n - 1)
    cₙ = ps[end-1]/ps[end] / αₙ
    p = P(ps)
    A = comrade_matrix(p)
    Aₒₒ = norm(A, Inf)
    [norm(αₙ * p(λ) / cₙ) /  (Aₒₒ * norm(v(λ), Inf)) for λ ∈ λs]
end


text:
f(x::Float64, y::Float64) = 2x + y;

try
    f(2.0, 3)
catch y
    println(y)
end

f(x::Number, y::Number) = 2x - y;

println(f(2.0, 3))

ext: jl


text:
my_f(x,y)=2*x+y
my_f2(x,y)=x+y

max_stars_repo_name: giandopal/MyExample

max_issues_repo_name: giandopal/MyExample

max_forks_repo_name: giandopal/MyExample


text:
#------------------------------------------------------------------------------
# Author: Edward Tate <edward.tate@erlang-solutions.com>
#------------------------------------------------------------------------------

export handle_response

#------------------------------------------------------------------------------
# Response Handlers
#------------------------------------------------------------------------------

function handle_response (s::TCPSocket, x::Any)
    info("[response] Response $x")
    reply_close(s, x)
end

#=
When a list of bytes is provided as a response, it is automatically assumed to
be streamed out. It is assumed that each byte array is one 'chunk' terminated
by a 4 byte 0 header. This should change soon.
=#

function handle_response (s::TCPSocket, bytes::Array{Uint8})
    info("[response] Forming streaming response")
    reply_info(s, :stream, [])
    reply_stream_close(s, bytes)
end

max_issues_repo_name: esl/FsServer

max_forks_repo_name: esl/FsServer


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 945c2bf1-d7dc-42c9-93d7-fd754f8fb1d7
html"""

Section 2.9

 Optimization

Lecture Video

"""

# ╔═╡ 400ebe26-0dea-4cf2-8744-6c73a45cd33e
begin
    import Pkg
    Pkg.activate(mktempdir())
    Pkg.add([
        Pkg.PackageSpec(name="PlutoUI", version="0.7"),
        Pkg.PackageSpec(name="Plots", version="1"),
        Pkg.PackageSpec(name="Optim", version="1"),
        Pkg.PackageSpec(name="JuMP", version="0.21"),
        Pkg.PackageSpec(name="Ipopt", version="0.6"),
        Pkg.PackageSpec(name="ForwardDiff", version="0.10"),
    ])
    using PlutoUI, Plots, Statistics, Optim, JuMP, Ipopt, ForwardDiff
end

# ╔═╡ b8d66df5-f593-40b4-8c46-3b638f9cc3e1
TableOfContents(title="📚 Table of Contents", aside=true)

# ╔═╡ 77253dd5-a2c8-4cf5-890a-5c8420c395b7
md"""
# Julia concepts

- Named tuples
"""

# ╔═╡ dccbd53d-33ed-4d37-9d2c-da76e090d5dd
md"""
# Line Fitting Many Ways
"""

# ╔═╡ 235919f2-a4b0-4bb5-870c-82809a170195
md"""
Last lecture we did some line fitting ("regression"). Let's see how we can actually solve that problem.
"""

# ╔═╡ 2ed86f33-bced-413c-9a8d-c6e49bfe5afb
md"""
# Exploratory Data Analysis
"""

# ╔═╡ a78aff8d-5ac0-4915-92de-ffefdb08f88e
md"""
Let's start off by making some noisy data.
"""

# ╔═╡ 0e43a6d3-7198-422b-b50c-b9caeaa53074
md"""
n = $(@bind n Slider(3:10:200, show_value=true))
"""

# ╔═╡ f8c98995-2152-4d45-996a-a0532a233719
x = sort((rand( -10:100, n)))

# ╔═╡ 8a5f1fdc-3cef-4c02-a73f-e5975b57b15a
y = 5/9 .* x  .- 17.7777777  .+  5 .* randn.() #  same as y =  5/9 .* (x .- 32)

# ╔═╡ 647093eb-a7e3-4175-8091-29c33407e5c9
begin
	plot(x,y, m=:c, mc=:red, legend=false, ls=:dash)
	xlabel!("°F")
	ylabel!("°C")
	# plot!( x, (x.-30)./2) Dave's cool approximation
end

# ╔═╡ cdc25782-65a8-43c5-8090-c1241b798b1a
md"""
# Least Squares fitting to a straight line
"""

# ╔═╡ 9ec4dd43-c95a-4f11-b844-fd6ccc93bb68
md"""
Suppose we are given data $x_i$ and measurements $y_i$. **Least-squares fitting** a straight line means finding the best ``m`` (slope) and ``b`` (intercept) that minimize
the "error" (distance from the data) in a least-squares sense:

$$\min_{m,b} \sum  ( (b + m x_i) - y_i)^2$$
"""

# ╔═╡ 9276b315-27b2-4b01-8fc8-4ebbba58d080
md"""
# Direct Formulas
"""

# ╔═╡ d22fd4bd-acfe-4e27-a484-3c2d6138f44e
md"""
## The Statistician's formula
"""

# ╔═╡ da0d208b-7d30-470a-b180-4cbfa98298e7
begin
	m = cov(x,y)/var(x) # same as (x.-mean(x))⋅(y.-mean(y))/sum(abs2,x.-mean(x))
	b = mean(y) - m * mean(x)
	(b=b, m=m)
end

# ╔═╡ 6cf233a7-9b8b-47aa-a3ad-2440d001af73
md"""
### Julia: Named Tuples
"""

# ╔═╡ 613c3e5f-bbdd-4cf9-b30f-69e2c42ae0ec
nt = (first=1, next=2, last=3.1) # kind of handy

# ╔═╡ 4cce580b-0032-419c-b386-e470b084ab96
typeof( nt )

# ╔═╡ 5503b4de-0b53-4223-8ce0-5e014be3f7ab
plot!(x -> m*x + b, lw=3, alpha=0.7)

# ╔═╡ 05e512ca-3123-48d9-9c11-5d6e9d90ef95
md"""
## The Linear Algebraist's Formula
"""

# ╔═╡ 939900b4-5327-43b4-883f-740c173c0db4
md"""
This is even shorter, but you need to know linear algebra. But it also generalizes.
"""

# ╔═╡ e0b4c2a9-a68b-47af-bf9c-f1a9f0256fd4
[one.(x) x]\y  # even shorter but you need to know linear algebra, but generalizes

# ╔═╡ 6d25e38e-c18a-48b3-8b12-b670f5a5180f
md"""
# Optimization Methods

Since the problem is an optimization problem, we can use optimization software to obtain an answer.  This is overkill for lines, but generalizes to so many nonlinear situations, including neural networks as in machine learning.
"""

# ╔═╡ f291c0cb-51ee-4b30-9e07-e7cf374f809e
md"""
## Optim.jl: A package written entirely in Julia for optimization
"""

# ╔═╡ aa06a447-d6c5-48ee-9864-c1f431fe5e4b
md"""
[Optim.jl Documentation](https://julianlsolvers.github.io/Optim.jl/stable/#)
"""

# ╔═╡ d3edfb26-7258-45a3-a88c-60831338df1f
md"""
We can  ask software to just solve the problem

$$\min_{b,m} \sum_{i=1}^n  [ (b + m x_i) - y_i]^2$$

or
``\min_{b,m}``  `loss(b,m)`

"""

# ╔═╡ 372b304a-3f57-4bec-88df-3d51ded57d5c
loss( (b, m) ) = sum( (b + m*x[i] - y[i])^2  for i=1:n )

# ╔═╡ 13b9ff38-225d-4ec1-be7f-bf0e0f5b4076
result =  optimize(loss, [0.0,0.0] )  # optimize f with starting guess

# ╔═╡ 7bd9bb8f-36c5-4ae1-ba20-25732d7fef2e
result.minimizer

# ╔═╡ 0af48ea2-698e-4919-aa96-97c5f46a928b
md"""
### Functions of Functions and Computing Power

Optimization such as
``\min_{b,m}``  `loss(b, m)`
is an example of a fairly heavy function of a function.  By this we mean that the input is a *function* such as  `loss(b, m)` and the output is the location or value of a minimum, say.  By "heavy" we mean that typically a large amount of computing power is needed.

Not that many years ago, computers were not strong enough for realistic problems.  Modern-day machine learning and so much more is enabled because computers can now surround entire codes with optimization, or if the software is compatible, automatic differentiation.
"""

# ╔═╡ 10386ce6-82fd-46ea-a44a-6ba14c5b0cd9
md"""
## JuMP.jl: A popular modelling language for Optimization Problems

JuMP = Julia for Mathematical Programming
"""

# ╔═╡ b7d8f11d-91ce-4b3a-87a1-1aa162e198ff
let

	n = length(x)
	model = Model(Ipopt.Optimizer)

	@variable(model, b)
	@variable(model, m)

    @objective(model, Min, sum( (b + m*x[i] - y[i])^2 for i in 1:n) )

	#set_silent(model)
	optimize!(model)

	(b=getvalue(b), m=getvalue(m))
end

# ╔═╡ 5ca85768-a19e-4ddf-89a4-88dca599d7a7
md"""
# Gradients
"""

# ╔═╡ dd39b088-f59f-43fa-bce0-5076398238f9
md"""
The above optimization methods made no explicit mention of derivative or gradient information.  For simple problems, gradients can be hand calculated, but for many real problems this is impractical.
"""

# ╔═╡ 5f41acf0-22bd-4224-a65a-81bd656e1c07
md"""
### Hand Computation
"""

# ╔═╡ 84f3a912-031c-40ed-ae29-02bbcc7b4612
md"""
``  \frac{\partial}{\partial b}\sum_{i=1}^n  ( (b + m x_i) - y_i)^2
=   2\sum_{i=1}^n  ( (b + m x_i) - y_i) ``

``   \frac{\partial}{\partial m }\sum_{i=1}^n  ( (b + m x_i) - y_i)^2
=   2\sum_{i=1}^n  x_i( (b + m x_i) - y_i) ``
"""

# ╔═╡ 6f64ede7-612e-47b3-b3a4-d22a1992a98d
begin
	∇loss(b,m,i) = 2*(b+m*x[i]-y[i]) .* [1,x[i]] # ith summand
	∇loss(b,m) = sum(∇loss(b,m,i) for i=1:n)

end

# ╔═╡ 36300b71-5a96-4964-b661-93de5631cf07
md"""
### Finite Difference Evaluation
"""

# ╔═╡ ad578b33-4387-49f5-b39d-92e05fca4ea5
∇loss(.1,.3)

# ╔═╡ 67fd90d7-bb34-411f-89f1-a410e6fb29ba
begin # finite difference
	ϵ = .000000001
	([loss([.1+ϵ ,.3]);loss([.1 ,.3+ϵ])] .- loss([.1 ,.3])) ./ ϵ
end

# ╔═╡ 3e229e4a-a697-460e-b995-a4773a6aca70
md"""
### Automatic Differentiation (AutoDiff)
 We're all so good at calculus I suppose, but nothing like letting the computer do it.  For real problems, what you learned in calculus is impractical. Note: Autodiff is not finite differences.  It is not as problematic as finite differences is in that
with finite differences it can be hard to know which ϵ to use, etc.
"""

# ╔═╡ 7566cb7e-f5da-4b81-af07-bf2c86963333
∇loss(.1,.3) # hand computation

# ╔═╡ e6d9aafd-fbd6-4ec4-a4a1-740a4e889dc5
ForwardDiff.gradient( loss, [.1,.3])

# ╔═╡ c9417d90-a9cb-4655-a258-8a8898e5576a
md"""
# Gradient Descent can be difficult for complicated functions
"""

# ╔═╡ 6535280a-e0ce-4e13-86dd-165d5f06cfe7
let
	b,m = 0,0  # starting guess

	for i=1:25
		db,dm = ∇loss(b,m)

		# Getting a good step size can be really hard
		# I worked out the line search η by hand
		η = sum( (b+m*x[i]-y[i])*(db+dm*x[i]) for i=1:n)/sum( (db+dm*x[i])^2 for i=1:n)

		b,m  = (b,m) .- η .* (db,dm)
	end
	(b=b,m=m)

end

# ╔═╡ 4c285bc2-b3c2-4d20-a904-ecaa07795342
md"""
Hoping for

`(b=-18.1716, m=0.56601)`
"""

# ╔═╡ 592397eb-ec52-423b-925b-d8becb9eac8e
md"""
# Stochastic Gradient Descent

Pick one coordinate (or a few coordinates) to update at a time.
This is what works in machine learning.

"""

# ╔═╡ 7086950b-c8db-49d4-b095-15be91c73b56
let
	b, m  = 0.0, 0.0
	for t=1:10_000_000
	    η = .00002  # there seems to be an art to picking these steplengths

	    b,m  =  (b,m) .- η *∇loss(b,m, rand(1:n))

	end
   	(b=b,m=m)
end

# ╔═╡ 327514f1-8081-4a6c-8be4-8ffd52ed3c46
md"""
# Bells and Whistles for optim.jl
"""

# ╔═╡ 98e00b2d-0802-4160-8e5c-302be5226916
optimize(loss, [0.0,0.0], BFGS(),  autodiff=:forward)

# ╔═╡ ef165ca5-bf4f-465e-8e9a-df1aec2d7caa
optimize(loss, [0.0,0.0], BFGS() )

# ╔═╡ 0305b418-51bb-47bb-98fb-319fc26b94cf
optimize(loss, [0.0,0.0], GradientDescent() )

# ╔═╡ 304a3a6e-c8c3-48d8-a101-313b3aa062f2
optimize(loss, [0.0,0.0], GradientDescent(), autodiff=:forward )

# ╔═╡ Cell order:
# ╟─945c2bf1-d7dc-42c9-93d7-fd754f8fb1d7
# ╠═400ebe26-0dea-4cf2-8744-6c73a45cd33e
# ╠═b8d66df5-f593-40b4-8c46-3b638f9cc3e1
# ╟─77253dd5-a2c8-4cf5-890a-5c8420c395b7
# ╟─dccbd53d-33ed-4d37-9d2c-da76e090d5dd
# ╟─235919f2-a4b0-4bb5-870c-82809a170195
# ╟─2ed86f33-bced-413c-9a8d-c6e49bfe5afb
# ╟─a78aff8d-5ac0-4915-92de-ffefdb08f88e
# ╟─0e43a6d3-7198-422b-b50c-b9caeaa53074
# ╠═f8c98995-2152-4d45-996a-a0532a233719
# ╠═8a5f1fdc-3cef-4c02-a73f-e5975b57b15a
# ╠═647093eb-a7e3-4175-8091-29c33407e5c9
# ╟─cdc25782-65a8-43c5-8090-c1241b798b1a
# ╟─9ec4dd43-c95a-4f11-b844-fd6ccc93bb68
# ╟─9276b315-27b2-4b01-8fc8-4ebbba58d080
# ╟─d22fd4bd-acfe-4e27-a484-3c2d6138f44e
# ╠═da0d208b-7d30-470a-b180-4cbfa98298e7
# ╟─6cf233a7-9b8b-47aa-a3ad-2440d001af73
# ╠═613c3e5f-bbdd-4cf9-b30f-69e2c42ae0ec
# ╠═4cce580b-0032-419c-b386-e470b084ab96
# ╠═5503b4de-0b53-4223-8ce0-5e014be3f7ab
# ╟─05e512ca-3123-48d9-9c11-5d6e9d90ef95
# ╟─939900b4-5327-43b4-883f-740c173c0db4
# ╠═e0b4c2a9-a68b-47af-bf9c-f1a9f0256fd4
# ╟─6d25e38e-c18a-48b3-8b12-b670f5a5180f
# ╟─f291c0cb-51ee-4b30-9e07-e7cf374f809e
# ╟─aa06a447-d6c5-48ee-9864-c1f431fe5e4b
# ╟─d3edfb26-7258-45a3-a88c-60831338df1f
# ╠═372b304a-3f57-4bec-88df-3d51ded57d5c
# ╠═13b9ff38-225d-4ec1-be7f-bf0e0f5b4076
# ╠═7bd9bb8f-36c5-4ae1-ba20-25732d7fef2e
# ╟─0af48ea2-698e-4919-aa96-97c5f46a928b
# ╟─10386ce6-82fd-46ea-a44a-6ba14c5b0cd9
# ╠═b7d8f11d-91ce-4b3a-87a1-1aa162e198ff
# ╟─5ca85768-a19e-4ddf-89a4-88dca599d7a7
# ╟─dd39b088-f59f-43fa-bce0-5076398238f9
# ╟─5f41acf0-22bd-4224-a65a-81bd656e1c07
# ╟─84f3a912-031c-40ed-ae29-02bbcc7b4612
# ╠═6f64ede7-612e-47b3-b3a4-d22a1992a98d
# ╟─36300b71-5a96-4964-b661-93de5631cf07
# ╠═ad578b33-4387-49f5-b39d-92e05fca4ea5
# ╠═67fd90d7-bb34-411f-89f1-a410e6fb29ba
# ╟─3e229e4a-a697-460e-b995-a4773a6aca70
# ╠═7566cb7e-f5da-4b81-af07-bf2c86963333
# ╠═e6d9aafd-fbd6-4ec4-a4a1-740a4e889dc5
# ╟─c9417d90-a9cb-4655-a258-8a8898e5576a
# ╠═6535280a-e0ce-4e13-86dd-165d5f06cfe7
# ╟─4c285bc2-b3c2-4d20-a904-ecaa07795342
# ╟─592397eb-ec52-423b-925b-d8becb9eac8e
# ╠═7086950b-c8db-49d4-b095-15be91c73b56
# ╟─327514f1-8081-4a6c-8be4-8ffd52ed3c46
# ╠═98e00b2d-0802-4160-8e5c-302be5226916
# ╠═ef165ca5-bf4f-465e-8e9a-df1aec2d7caa
# ╠═0305b418-51bb-47bb-98fb-319fc26b94cf
# ╠═304a3a6e-c8c3-48d8-a101-313b3aa062f2


text:
using SyntheticGrids
using Base.Test
using JSON

@testset "SyntheticGrids" begin
    const SEED = 666
    const GENPATH = joinpath(dirname(@__FILE__), "..", "data", "GenData.json")
    const CENSUSPATH = joinpath(dirname(@__FILE__), "..", "data", "Census_data.dat")
    const GENCOORDPATH = joinpath(dirname(@__FILE__), "..", "data", "Generator_coord.dat")
    const GENDATAPATH = joinpath(dirname(@__FILE__), "..", "data", "Generator_data.dat")
    const DUMMYPATH = joinpath(dirname(@__FILE__), "..", "data", "dummy.json")
    const DUMMYPPC = joinpath(dirname(@__FILE__), "..", "data", "dummy.p")
    const GRIDPATH = joinpath(dirname(@__FILE__), "..", "data", "testgrid.json")

    grid = Grid(SEED)
    grid1 = Grid(SEED)

    @testset "Basic grid creation" begin
        place_loads_from_zips!(grid; latlim = (38, 40), longlim = (-89, -88))
        @test length(buses(grid)) == 130
        place_gens_from_data!(grid; latlim = (38, 40), longlim = (-89, -88))
        @test length(buses(grid)) == 147
        connect!(grid)
        @test sum(adjacency(grid)) == 358
        count = SyntheticGrids.count_bus_type(grid)
        @test count == (130, 17)
        create_lines!(grid)
        @test grid.trans_lines[1].capacity == 4900
        @test length(trans_lines(grid)) == 179
        @test test_connectivity(grid.bus_conn, false)
        @test total_links(grid.bus_conn) == 179
        @test mean_node_deg(grid.bus_conn) ≈ 2.435374149659
        @test mean_shortest_path(adjacency(grid)) ≈ 6.764218612615115
        @test mean_shortest_path(adjacency(grid), distance(buses(grid))) ≈ 133.1803034218734
        @test robustness_line(grid.bus_conn, 10) > 1
        @test robustness_node(grid.bus_conn, 10) > 1
        place_loads_from_zips!(
            grid1;
            latlim = (38, 40),
            longlim = (-89, -88),
        )
        place_gens_from_data!(
            grid1;
            latlim = (38, 40),
            longlim = (-89, -88),
        )
        connect!(grid1)
        create_lines!(grid1)
        @test grid == grid1

        grid2 = Grid()
        function simple_func(t::Tuple)
            return t[1] > -t[2]/2
        end
        place_loads_from_zips!(grid2, simple_func)
        @test length(buses(grid2)) == 8171
        place_gens_from_data!(grid2, simple_func)
        @test length(buses(grid2)) == 9957

        grid3 = Grid(true)
        @test size(grid3.bus_conn) == (137, 137)
        @test size(grid3.sub_conn) == (43, 43)
        add_gen!(grid3, (11., 11.), 100, [11], ["ki"], reconnect = false)
        add_load!(grid3, (19.,12.), 100, 100, 100, reconnect = true)
        add_substation!(
            grid3,
            (11., 18.),
            [11],
            1000,
            6000,
            500,
            Set(),
            [grid3.buses[139], grid3.buses[138]]
        )
        @test size(grid3.bus_conn) == (139, 139)
        @test size(grid3.sub_conn) == (44, 44)
    end

    @testset "Clustering grids" begin
        count = SyntheticGrids.count_bus_type(grid)
        cluster!(
            grid,
            round(Int, 0.3 * count[1]),
            round(Int, 0.05 * count[2]),
            round(Int, 0.5 * count[2])
        )
        @test length(substations(grid)) == 47
        @test grid.substations[1].population == 6193
        @test grid.substations[47].generation == 16.1
        @test test_connectivity(grid.sub_conn, false)
        @test total_links(grid.sub_conn) == 69
        @test mean_shortest_path(sub_connectivity(grid) .> 0) ≈ 4.107741059302851
        @test mean_shortest_path(
            (sub_connectivity(grid) .> 0),
            distance(substations(grid))
            ) ≈ 113.21283328057127
        @test robustness_line(grid.sub_conn, 10) > 1
        @test robustness_node(grid.sub_conn, 10) > 1
        count = SyntheticGrids.count_bus_type(grid1)
        cluster!(
            grid1,
            round(Int, 0.3 * count[1]),
            round(Int, 0.05 * count[2]),
            round(Int, 0.5 * count[2])
        )
        @test (grid == grid1)
    end

    @testset "Input and Output" begin
        @testset "Pandapower interface" begin
            pgrid = to_pandapower(grid, DUMMYPPC)
            @test length(pgrid[:trafo]) == 36
            pgrid = SyntheticGrids.load_pp_grid(DUMMYPPC)
            @test length(pgrid[:trafo]) == 36
        end

        @testset "Generator data" begin
            SyntheticGrids.prepare_gen_data(GENDATAPATH, GENCOORDPATH, DUMMYPATH)
            genjson = JSON.parsefile(GENPATH)
            dummyjson = JSON.parsefile(DUMMYPATH)
            @test genjson == dummyjson
        end

        @testset "Save and Load Grid" begin
            save(grid, GRIDPATH)
            lgrid = load_grid(GRIDPATH)
            @test adjacency(grid) == adjacency(lgrid)
            @test sub_connectivity(grid) == sub_connectivity(lgrid)
            lgrid = 0 # trying to solve issues with Windows locking permissions to GRIDPATH
        end

        rm(DUMMYPATH)
        rm(DUMMYPPC)
        rm(GRIDPATH)
    end
end


text:
using JuliaDB, Interact, TableWidgets, Tables, Blink
using Observables: @map

# Pass a table with the last column containing only strings
function label!(t; categories=[""])

    category_vec = last(Tables.columntable(t))
    category_obs = Observable{Any}(category_vec)
    options = @map union(&category_obs, categories)

    function myrow(r, i; format = TableWidgets.format)
        fields = propertynames(r)
        category = getproperty(r, last(fields))
        wdg = autocomplete(options, value = category) |> onchange
        on(wdg) do val
            category_obs[][i] = val
            category_obs[] = category_obs[]
        end

        node("tr",
            node("th", format(i)),
            (node("td", format(getproperty(r, field))) for field in fields[1:end-1])...,
            node("td", wdg))
    end

    TableWidgets.rendertable(t, row = myrow)
end

# Annotate your expenses / income from a set of categories (the set of categories
# increases every time you type one manually)
df = table((Place = ["home", "home", "work"], Amount = [-12.3, -1.2, 1400], Category = ["", "", ""]))
t = label!(df, categories = ["Food", "Coffe", "Salary"])

w = Window()
body!(w, t)


text:
using Convex, ECOS

function go()
    # A = [-1 -2  6  3 8;
    #      2  3 -2 -1 3]
    # x = Variable(2, Positive())

    A = [2 0 -2 0 3;
         0 5 -3 0 3;
         0 0 5 -1 8;
         0 -1 0 5 8]
    x = Variable(4, Positive())

    prob = maximize(sum(log((A' * x)[1:4])),
                    [sum(x) == 100,
                     (A' * x)[5] == 500
                     ])

    solve!(prob, ECOS.Optimizer)

    A' * round.(evaluate(x))
end

Int(prod(go()[1:4]))


text:
using Sugar: expr_type, typed_expr, getfunction, similar_expr
"""
Emit a LazyMethod call Expr
"""
function emit_call(m, fun, typ, args...)
    typs = map(x-> expr_type(m, x), args)
    m2 = LazyMethod(m, fun, typs)
    typed_expr(typ, :call, m2, args...)
end

"""
Types that are accessible only via getfield,
e.g. Tuples that are translated to structs - so in Julia
They'll be addressable via getindex, but in another language it must be getfield.
"""
supports_indexing(m::LazyMethod, T) = false

"""
Types that are only indexable via static indices.
"""
needs_static_index(m::LazyMethod, T) = true

supports_indices(m::LazyMethod, T, index_types) = false

to_index(m, t, idx::T) where T <: Integer = (idx - T(1))
function to_index(m, T, idx)
    int_type = expr_type(m, idx)
    emit_call(m, -, int_type, idx, int_type(1))
end

function rewrite_indices(m::LazyMethod, T, indices)
    map(indices) do idx
        if supports_indexing(m, T) && !(is_fixedsize_array(T) && isa(idx, Integer))
            return to_index(m, T, idx)
        else
            return c_fieldname(m, T, idx) # indexing not supported, use getfield
        end
    end
end
@generated function cl_getindex(x::T, idx::Integer) where T
    expr = Expr(:block)
    N = nfields(T)
    for i = 1:N
        name = fieldname(T, i)
        if i == 1
            push!(expr.args, :(elem = $getfield(x, $(QuoteNode(name)))))
        else
            push!(expr.args, :(if $i == idx
                elem = $getfield(x, $(QuoteNode(name)))
            end))
        end
    end
    push!(expr.args, :(return elem))
    expr
end
function index_expression(m, expr, args, types)
    mparent = m.parent
    is_setindex = getfunction(m) == setindex!
    isempty(types) && error("Indexing into nothing with func: $m")
    T = types[1] # type to index into
    isimmutable(T) && is_setindex && error("Can't use setindex! on immutable. Found: $T with $types")
    indices, val, index_types = if is_setindex
        args[3:end], args[2], types[3:end]
    else
        args[2:end], nothing, types[2:end]
    end

    if T <: Tuple{X} where X <: cli.Numbers
        # drop indexing into tuples of length 1. They get mapped to the element type
        # which is not very nice, but allows to map fixedsize arrays to tuples,
        # and still have fast Tuple{N, <: Numbers} types.
        return Sugar.rewrite_ast(mparent, expr.args[2])
    end
    if supports_indices(mparent, T, index_types)
        fixed_size_static_idx = is_fixedsize_array(T) && length(indices) == 1 && isa(first(indices), Integer)
        ret = if supports_indexing(mparent, T) && !fixed_size_static_idx
            indices = rewrite_indices(mparent, T, indices)
            indexing = typed_expr(expr.typ, :ref, args[1], indices...)
            if is_setindex
                if Sugar.expr_type(mparent, val) == eltype(T)
                    indexing = :($indexing = $val)
                else
                    expr.args[1] = m # leave getindex
                    return expr
                end
            end
            indexing
        else
            if (T <: Tuple || is_fixedsize_array(T)) # Julia inbuilds allowing getindex, but need to use getfield in C
                if all(x-> isa(x, Integer), indices)
                    indices = rewrite_indices(mparent, T, indices)
                    emit_call(mparent, getfield, expr.typ, args[1], indices...)
                else
                    emit_call(mparent, cl_getindex, expr.typ, args[1], indices...)
                end
            else
                expr.args[1] = m # leave getindex
                expr
            end
        end
        return ret
    else
        expr.args[1] = m # leave getindex
        return expr
    end
end

function emit_constructor(m::LazyMethod, realtype, args)
    Sugar.typed_expr(realtype, :curly, m, args...)
end

gpu_select_value(x::Bool, a, b) = x ? a : b

function Sugar.rewrite_function(method::Union{LazyMethod{:CL}, LazyMethod{:GL}}, expr)
    f, _types = method.signature
    types = Sugar.to_tuple(_types)
    li = method.parent
    F = typeof(f)
    # first give the backends a chance to do backend specific rewriting
    rewritten, expr = rewrite_backend_specific(method, f, types, expr)
    rewritten && return expr
    rewritten, expr = rewrite_fastmath(f, types, li, expr)
    rewritten && return expr
    # if nothing got rewritten, do the backend independant rewrites:
    if f in (getindex, setindex!)
        return index_expression(method, expr, expr.args[2:end], types)
    elseif f == getfield && length(types) == 2 && isa(expr.args[3], Integer)
        return emit_call(
            li, getfield,
            Sugar.expr_type(method, expr),
            expr.args[2], c_fieldname(method, types[1], expr.args[3])
        )
    # elseif f == Base.indexed_next && length(types) == 3 && isa(expr.args[3], Integer)
    #     # if we have a static indexed next, we unfold it into a a getindex directly
    #     expr.args = expr.args[1:3]
    #     types = types[1:2]
    #     rexpr = index_expression(method, expr, expr.args[2:end], types)
    #     return rexpr
    elseif f == convert && length(types) == 2
        if types[1] == Type{types[2]}
            return Sugar.rewrite_ast(li, expr.args[3]) # no convert needed
        else # But for now we just change a convert into a constructor call
            realtype = Sugar.unspecialized_type(Sugar.expr_type(li, expr.args[2]))
            target_type = Sugar.unspecialized_type(Sugar.expr_type(li, expr.args[3]))
            if is_native_type(li, realtype) && is_native_type(li, target_type)
                constr_m = LazyMethod(li, realtype)
                return emit_constructor(constr_m, realtype, expr.args[3:end])
            else
                return emit_call(li, convert, realtype, realtype, expr.args[3])
            end
        end
    # Constructors
    elseif f == tuple
        realtype = Sugar.expr_type(li, expr)
        args = expr.args[2:end]
        if isempty(args) && sizeof(realtype) == 0
            push!(args, 0f0) # there are no empty types, so if empty, insert default
        end
        constr_m = LazyMethod(li, realtype)
        return emit_constructor(constr_m, realtype, args)
    elseif f == broadcast
        fb = Sugar.resolve_func(li, expr.args[2]) # rewrite if necessary
        # most shared intrinsics broadcast over fixedsize arrays
        if (fb in functions) && all(t-> is_native_type(li, t), types[2:end])
            shift!(expr.args) # remove broadcast from call expression
            expr.args[1] = LazyMethod(li, fb, types[2:end])
            return expr
        end
        expr.args[1] = method
        # expr.args[2] = LazyMethod(li, fb, types[2:end])
        return expr
    # div is / in c like languages
    elseif f == div && length(types) == 2 && all(x-> x <: cli.Ints, types)
        expr.args[1] = LazyMethod(li, (/), types)
        return expr
    # elseif f == Base.cttz_int && length(types) == 1 && types[1] <: cli.Ints
    #     expr.args[1] = LazyMethod(li, (/), types)
    #     return expr
    elseif f == length && length(types) == 1 && is_fixedsize_array(types[1])
        return fixed_array_length(types[1])
    return expr
    # Base.^ is pow in C
    elseif f == (^) && length(types) == 2 && all(t-> t <: cli.Numbers, types)
        expr.args[1] = LazyMethod(li, cl_pow, types)
        return expr
    elseif f == rem && length(types) == 2 && all(t-> t <: cli.Ints, types)
        expr.args[1] = LazyMethod(li, %, types)
        return expr
    elseif f == rem && length(types) == 2 && all(t-> t <: cli.Floats, types)
        expr.args[1] = LazyMethod(li, cli.remainder, types)
        return expr
    elseif (f == (⊻)) && length(types) == 2 && all(t-> t <: cli.Numbers, types)
        expr.args[1] = LazyMethod(li, ^, types)
        return expr
    elseif (f == Base.select_value) && length(types) == 3 && first(types) == Bool
        expr.args[1] = LazyMethod(li, gpu_select_value, types)
        return expr
    elseif f == (*) && length(types) == 1 && all(t-> t <: cli.Numbers, types)
        m = LazyMethod(li, identity, types)
        push!(li, m)
        expr.args[1] = m # drop *(x)
        return expr
    end
    expr.args[1] = method
    return expr
end

# get the concrete function objects
const fast_ops = getfield.(Base.FastMath, collect(values(Base.FastMath.fast_op)))
const slow_ops = getfield.(Base.FastMath, collect(keys(Base.FastMath.fast_op)))

const fast_op_2_op = Dict(zip(fast_ops, slow_ops))

function rewrite_fastmath(op, types, li, expr)
    if haskey(fast_op_2_op, op)
        slow_op = fast_op_2_op[op]
        if Sugar.isintrinsic(li, slow_op, Base.to_tuple_type(types))
            # actually, most intrinsics in OpenCL are already "fast ops"
            m = LazyMethod(li, slow_op, types)
            push!(li, m)
            expr.args[1] = m
            return true, expr
        end
    end
    return false, expr
end


text:
function evaluate(obj::WeightedEvaluation, arg0::Vector{jdouble}, arg1::Vector{jdouble})
    return jcall(obj, "evaluate", jdouble, (Vector{jdouble}, Vector{jdouble}), arg0, arg1)
end

function evaluate(obj::WeightedEvaluation, arg0::Vector{jdouble}, arg1::Vector{jdouble}, arg2::jint, arg3::jint)
    return jcall(obj, "evaluate", jdouble, (Vector{jdouble}, Vector{jdouble}, jint, jint), arg0, arg1, arg2, arg3)
end


text:
#This script executes the Trust Region Policy Optimisation Algorithm from Schulmann

include("src/init.jl")

θa = init_θ_single_norm(params.structure_net_pl)
θv = init_θ_single_norm(params.structure_net_vf)
best_net = init_θ_single_norm(params.structure_net_pl)
res = zeros(params.n_iter)
prec = zeros(params.n_iter)
all_traj  = init_all_traj(params)

"""
Trust Region Policy Optimisation
"""
function trpo(params::Params, θa::Any, θv::Any, res::Array{Float64,1}, prec::Array{Float64,1}, all_traj::Traj_Data, best_net::Any)
    #Init Part:
    best_net = copy(θa)
    min_cost = Inf
    #------------------------
    #Initialize the Policy
    ĝ_a = copy(θa)
    opt_para_a = optimizers(θa,   Adam; lr = params.lr_pl)

    #Initialize the Value Function
    ĝ_v = copy(θv)
    opt_para_v = optimizers(θv,   Adam; lr = params.lr_vf)

    #Init Progress bar
    progress = Progress(params.n_iter,1)

    # Init rms
    rms_all = init_rms(params.s_init)
    rms_loc = init_rms(params.s_init)

    #saving old data to see differences
    a_t_old = all_traj.traj[1].a_t
    weights_old = zeros(length(net_to_vec(θa)),2)

    for i_iter = 1:params.n_iter
        ĝ_a *= 0
        ĝ_v *= 0

        #execute the parallel_episodes
        res[i_iter], prec[i_iter] = parallel_episodes(params,θa,θv,all_traj,rms_all)

        #Reorganice Data in minibatches
        mini_batch = build_minibatches(all_traj,params)

        #First refit the value function for new data batch:
        θv = refit_VF(params, θv, mini_batch, opt_para_v)

        #Now apply the gae for the new data based on updated value function
        gae_all!(all_traj,θv,params)

        #Reorganice Data in minibatches as advantages have changed...
        mini_batch = build_minibatches(all_traj,params)

        #Apply the gradient
        ĝ_a = parallel_gradient(params, θa, mini_batch)
        θa  = do_TRPO_update(ĝ_a, mini_batch, θa, params)

        if min_cost > res[i_iter]
            min_cost = res[i_iter]
            best_net = copy(θa)
            print(" new_best: ", round(min_cost, digits = 3))
        end

        #Visualization_stuff
        #-----------------------------------------------------------------------
        ProgressMeter.update!(progress,i_iter) #update progress bar
        if floor(i_iter/50) == ceil(i_iter/50)
            a_t_old ,weights_old = plot_stuff(params, all_traj, a_t_old,
            weights_old, i_iter,θa, θv, res, prec/params.n_traj, ĝ_a, ĝ_v)
        end
        #-----------------------------------------------------------------------

        @everywhere GC.gc(true)
    end
    return res
end

trpo(params, θa, θv, res, prec, all_traj, best_net)

max_stars_repo_name: NikonPic/REINFORCEMENT-JULIA

max_issues_repo_name: NikonPic/REINFORCEMENT-JULIA

max_forks_repo_name: NikonPic/REINFORCEMENT-JULIA


text:
"""
Sls fileformat
Sector list file

The link between this and the region file ARE is made by volume name

See EUROCONTROL NEST Manual for Sls fileformat description

Sls DataFrame Column Names
#   Field           Type        Size    Comment
1   :SECTORNAME     String      Max 19
2   :VOLUMENAME     String              Refer to name in volumne IN the ARE file (max 24 characters)
3   :VOLUMEBOTTOMLEVEL  Int64           In FL, overwrite level present in ARE file for that volume
4   :VOLUMETOPLEVEL Int64               In FL, overwrite level present in ARE file for that volume

Example SLS fileformat
CASTOR_UP + AX 000 295
CASTOR_UP + GU 245 265
TMA + SLP 000 095
TMA + SLX 000 125
ABRON + TOP 45 660
"""
module Sls

export read

using Format
using CSV
using DataFrames

const fileformat = Dict(1=>String, 2=>String, 3=>String, 4=>Int64,
5=>Int64)

const header_format = ["SECTORNAME", "VOLUMESIGN", "VOLUMENAME",
"VOLUMEBOTTOMLEVEL", "VOLUMETOPLEVEL"]

function read(file)
    df = CSV.read(file, delim=" ", types=fileformat, header=header_format,
    copycols=true)
    select!(df, Not(:VOLUMESIGN))
    return df
end

end # module


text:
@testset "Function Call" begin
s = """
function f(x; y = 1, z::Int = 4)
    println(z)
    x + y
end

f(1; y = 3)
"""
msgs = lintstr(s)
@test isempty(msgs)

s = """
function f{T}(x::T, y::T)
    a = Array{T, 1}()
    append!(a, [x, y])
    a
end
"""
msgs = lintstr(s)
@test isempty(msgs)

s = """
function f(::Type{Int}, x, y)
    x + y
end
"""
msgs = lintstr(s)
@test isempty(msgs)

s = """
function f(x, y, x)
    x + y
end
"""
msgs = lintstr(s)
@test msgs[1].code == :E331
@test contains(msgs[1].message, "duplicate argument")

s = """
function f{Int64}(x::Int64, y::Int64)
    x + y
end
"""
msgs = lintstr(s)
@test msgs[1].code == :E534
@test msgs[1].variable == "Int64"
@test contains(msgs[1].message, "introducing a new name for an implicit argument to the " *
    "function, use {T<:Int64}")

s = """
function f{T<:Int64}(x::T, y::T)
    x + y
end
"""
msgs = lintstr(s)
@test msgs[1].code == :E513
@test msgs[1].variable == "T <: Int64"
@test contains(msgs[1].message, "leaf type as a type constraint makes no sense")

s = """
function f{Int<:Real}(x::Int, y::Int)
    x + y
end
"""
msgs = lintstr(s)
@test msgs[1].code == :E536
@test contains(msgs[1].message, "use {T<:...} instead of a known type")

s = """
function f(x, args...)
    x + length(args)
end
"""
msgs = lintstr(s)
@test isempty(msgs)

s = """
function f(x, args..., bogus...)
    x + length(args) + length(bogus)
end
"""
msgs = lintstr(s)
@test msgs[1].code == :E413
@test contains(msgs[1].message, "positional ellipsis ... can only be the last argument")

s = """
function f(x=1, y, args...)
    x + length(args) + y
end
"""
msgs = lintstr(s)
@test msgs[1].code == :E411
@test contains(msgs[1].message, "non-default argument following default arguments")

s = """
function f(x, y; z, q=1)
    x + q + y
end
"""
msgs = lintstr(s)
@test msgs[1].code == :E423
@test contains(msgs[1].message, "named keyword argument must have a default")

s = """
function f(x, y; args..., z=1)
    x + length(args) + y
end
"""
msgs = lintstr(s)
@test msgs[1].code == :E412
@test contains(msgs[1].message, "named ellipsis ... can only be the last argument")

s = """
function f(x, args...; kwargs...)
    x + length(args) + length(kwargs)
end
"""
msgs = lintstr(s)
@test isempty(msgs)

s = """
function f(x::Array{Number,1})
    length(x)
end
"""
msgs = lintstr(s)
@test msgs[1].code == :E533
@test contains(msgs[1].message, "type parameters are invariant, try f{T<:Number}(x::T)...")

s = """
function f(x::Dict{Symbol,Number})
    length(x)
end
"""
msgs = lintstr(s)
@test msgs[1].code == :E533
@test contains(msgs[1].message, "type parameters are invariant, try f{T<:Number}(x::T)...")

s = """
function f(x; y = 1, z::Int = 0.1)
    x + y
end
"""
msgs = lintstr(s)
@test msgs[1].code == :E516
@test contains(msgs[1].message, "type assertion and default")

s = """
function f(x; y = 1, z::Int = error("You must provide z"))
    x + y + z
end
"""
msgs = lintstr(s)
@test isempty(msgs)

s = """
g(x) = -x
function f(arr::Array)
    map(g, arr)
end
"""
msgs = lintstr(s)
@test isempty(msgs)

s = """
g(; x=0, y=1) = x+y
function f()
    d = Dict{Symbol, Int}(:x=>6, :y=>4)
    g(; d...)
end
"""
msgs = lintstr(s)
@test isempty(msgs)

s = """
function f{T}(a::Vector{T})
    n = size(a)
    tmp  = Array(T, n)
    tmp2 = zeros(T, 1)
    tmp3 = zeros(1, 2, 3)
    @lintpragma("Info type a")
    @lintpragma("Info type n")
    @lintpragma("Info type tmp")
    @lintpragma("Info type T")
    @lintpragma("Info type tmp2")
    @lintpragma("Info type tmp3")
    tmp
end
"""
msgs = lintstr(s)
@test msgs[1].code == :I271
@test contains(msgs[1].message, "typeof(a) == Array")
@test msgs[2].code == :I271
@test contains(msgs[2].message, "typeof(n) == Tuple{Int64}")
if VERSION ≥ v"0.6-"
    @test msgs[3].code == :I271
    @test contains(msgs[3].message, "typeof(tmp) == Array")
end
@test msgs[4].code == :I271
@test_broken contains(msgs[4].message, "typeof(T) == Type")
if VERSION ≥ v"0.6-"
    @test msgs[5].code == :I271
    @test contains(msgs[5].message, "typeof(tmp2) == Array")
end
@test msgs[6].code == :I271
@test contains(msgs[6].message, "typeof(tmp3) == Array{Float64,3}")

s = """
function f(x)
    f = x
    f
end
"""
msgs = lintstr(s)
@test msgs[1].code == :W355
@test msgs[1].variable == "f"
@test contains(msgs[1].message, "conflicts with function name")

s = """
function f(x)
    local f = x
    f
end
"""
msgs = lintstr(s)
@test_broken isempty(msgs)

s = """
function f(x)
    x.somefunc(1)
end
"""
msgs = lintstr(s)
@test isempty(msgs)

s = """
function f(x=1)
    @lintpragma("Info type x")
    return x
end
"""
msgs = lintstr(s)
@test msgs[1].code == :I271
@test contains(msgs[1].message, "typeof(x) == Int")

s = """
function f(x::Int8=Int8(1))
    @lintpragma("Info type x")
    return x
end
"""
msgs = lintstr(s)
@test msgs[1].code == :I271
@test contains(msgs[1].message, "typeof(x) == Int8")

s = """
function f(c::Char)
    x = convert(Int, c)
    @lintpragma("Info type x")
    return x
end
"""
msgs = lintstr(s)
@test msgs[1].code == :I271
@test contains(msgs[1].message, "typeof(x) == Int")

s = """
function f(args...; dict...)
    @lintpragma("Info type args")
    @lintpragma("Info type dict")
    length(args) + length(dict)
end
"""
msgs = lintstr(s)
@test msgs[1].code == :I271
@test contains(msgs[1].message, "typeof(args) == Tuple")
@test msgs[2].code == :I271
@test contains(msgs[2].message, "typeof(dict) == Tuple")

s = """
function f(args::Float64...)
    @lintpragma("Info type args")
    length(args)
end
"""
msgs = lintstr(s)
@test msgs[1].code == :I271
@test contains(msgs[1].message, "typeof(args) == Tuple{Vararg{Float64")

s = """
function f(args::Float64...)
    if args[1] == 1
        0
    else
        length(args)
    end
end
"""
msgs = lintstr(s)
@test isempty(msgs)

s = """
function f(x::Symbol)
    if x == "blah"
        0
    else
        1
    end
end
"""
msgs = lintstr(s)
@test msgs[1].code == :W542
@test contains(msgs[1].message, "comparing apparently incompatible type")

s = """
function f(X::Int)
    Y = X + 1
    return Y
end
"""
msgs = lintstr(s)
@test isempty(msgs)

@test_broken isempty(lintstr("""
function f1(a::Float64, b::Float64)
  function f2(x::Float64, a=a)
    x + a^2
  end
  return f2(b)
end
"""))

s="""
function f1(a::Float64; b=string(a))
    string(a) * b
end
"""
msgs = lintstr(s)
@test isempty(msgs)

s = """
function f end
"""
msgs = lintstr(s)
@test isempty(msgs)

s = """
f(; x=1) = x
f(; 2)
"""
msgs = lintstr(s)
@test msgs[1].code == :E133
@test contains(msgs[1].message, "unknown keyword pattern")

s = """
function () end
"""
msgs = lintstr(s)
@test isempty(msgs)

s = """
function foo end
"""
msgs = lintstr(s)
@test isempty(msgs)

s = """
(==)(1, 1)
"""
msgs = lintstr(s)
@test isempty(msgs)

s="""
(==)(1)
"""
msgs = lintstr(s)
@test msgs[1].code == :E132
@test contains(msgs[1].message, "Lint does not understand argument")

s = """
(==)((1, 1)...)
"""
msgs = lintstr(s)
@test isempty(msgs)
end


text:
module LinearSolvers

using Gridap
using Gridap.Algebra
using IterativeSolvers
using IncompleteLU

import Gridap.Algebra: LinearSolver
import Gridap.Algebra: symbolic_setup
import Gridap.Algebra: numerical_setup
import Gridap.Algebra: numerical_setup!
import Gridap.Algebra: solve!

export GmresSolver
export symbolic_setup
export numerical_setup
export numerical_setup!
export solve!

struct GmresSolver <: LinearSolver
  preconditioner
  precond_kwargs::Dict
  verbose::Bool
  function GmresSolver(;
    verbose=false, preconditioner=IterativeSolvers.Identity, precond_kwargs...)
    new(preconditioner,precond_kwargs, verbose)
  end
end

struct GmresSymbolicSetup <: SymbolicSetup
  preconditioner
  precond_kwargs::Dict
  verbose::Bool
end

mutable struct GmresNumericalSetup{T<:AbstractMatrix} <: NumericalSetup
  A::T
  preconditioner
  precond_kwargs::Dict
  verbose::Bool
end

function symbolic_setup(s::GmresSolver, mat::AbstractMatrix)
  GmresSymbolicSetup(s.preconditioner, s.precond_kwargs, s.verbose)
end

function numerical_setup(ss::GmresSymbolicSetup,mat::AbstractMatrix)
  GmresNumericalSetup(mat, ss.preconditioner, ss.precond_kwargs, ss.verbose)
end

function numerical_setup!(ns::GmresNumericalSetup, mat::AbstractMatrix)
  ns.A = mat
end

function get_preconditioner(ns::GmresNumericalSetup)

end

function solve!(
  x::AbstractVector,ns::GmresNumericalSetup,b::AbstractVector)
  p=ns.preconditioner(ns.A; ns.precond_kwargs...)
  # initialize the solution to 0
  x .= 0.0
  gmres!(x, ns.A, b,verbose=ns.verbose, Pl=p, initially_zero=true)
end

end # module


text: module SparkMethodsForDegenerateLagrangianSystems end


text:
##### Beginning of file

@info("Importing the RepoSnapshots module...")

import RepoSnapshots

import TimeZones

@info("Reading config files...")

include(joinpath("config","preferences","branches.jl",))
include(joinpath("config","preferences","git-hosts.jl",))
include(joinpath("config","preferences","git-user.jl",))
include(joinpath("config","preferences","time-zone.jl",))

include(
    joinpath(
        "config", "repositories",
        "do-not-push-to-these-destinations.jl",
        )
    )
include(
    joinpath(
        "config", "repositories",
        "do-not-try-url-list.jl",
        )
    )
include(
    joinpath(
        "config", "repositories",
        "try-but-allow-failures-url-list.jl",
        )
    )

RepoSnapshots.CommandLine.run_snapshots_command_line!!(
    ;
    arguments = ARGS,
    git_user_name = GIT_USER_NAME,
    git_user_email = GIT_USER_EMAIL,
    src_provider = src_provider,
    dst_provider = dst_provider,
    include_branches = INCLUDE_BRANCHES,
    exclude_branches = EXCLUDE_BRANCHES,
    do_not_try_url_list = DO_NOT_TRY_URL_LIST,
    do_not_push_to_these_destinations = DO_NOT_PUSH_TO_THESE_DESTINATIONS,
    try_but_allow_failures_url_list = TRY_BUT_ALLOW_FAILURES_URL_LIST,
    time_zone = TIME_ZONE,
    )

##### End of file


text:
# I could enhance the assignment performance by constructing oid3d (not oind3d), which
# stores unique object IDs rather than the reference itself.  This array is different from
# oind3d in that it distinguishes objects joined by periodic boundary condition.  Then, I
# can retrieve objects from a map from this IDs to objects.  I could use oid3d in smoothing
# as well, by constructing an 8-vector oid3d_vxl, like pind3d_vxl and oind3d_vxl.
#
# The rationale for this trick is that assigning an Int to VecInt is faster than assigning a
# concrete Object (like Cuboid) to Vector{Object{3}}.  In my test,
#   setindex!(::VecInt, ::Int, <index>)
# was about twice as fast as
#   setindex!(::Vector{Object{3}}, ::Cuboid, <index>), though both took only a few nanoseconds.
#
# Currently the performance of assignment and smoothing seems satisfactory, so I will not
# pursue this extra optimization.  If I need to assign many objects, the situation may
# change and I may need to implement this optimization.

export create_param_array, create_oind_array, assign_param!, assign_val_shape!

# About the order of indices of paramKd
#
# (Bottom line: paramKd[i,j,k,v,w] uses a different indexing scheme than conventional
# indexing scheme F[w,i,j,k] for vector field F; note that the polarization indices v and w
# come last in paramKd, whereas the polarization index w comes first in F.  This is to speed
# up the material parameter assignment.  See below.)
#
# paramKd has 5 indices for K = 3: the first three are positional indices (i,j,k) and the
# last two are material parameter tensor component indices (v,w).  This mean paramKd[i,j,k,:,:]
# is the 3×3 material tensor.
#
# This is different from the conventional indexing scheme used for material parametetrs.
# For example, we usually write the xy-component of ε at a location (i,j,k) as ε_xy[i,j,k],
# where the tesnor component subscripts v = x and w = y comes earlier than the positional
# indices (i,j,k).  Similarly, the x-component of the E-field is usually written E_x[i,j,k],
# not E[i,j,k]_x.  (Also, the usual indexing scheme for vector fields is F[w,i,j,k].  This
# puts the x-, y-, z-componnets of F at the same location (i,j,k) as adjacent entries of
# the column vector F[:], which helps reducing the bandwidth of the matrix versions of local
# differential operators acting on F.)
#
# However, when assigning the material parameters in the code in this file, we often assign
# the same value to a block of paramKd.  This block is chosen differently for different
# material tensor components, because different tensor components are evaluated at different
# locations in Yee's grid.  Therefore, in this assignment, we first fix the material tensor
# components v and w, determine the range of (i,j,k) to which the same material parameter
# value to assign, and perform the assignment.
#
# Now, if we index paramKd as paramKd[v,w,i,j,k], v and w are the fastest-varing indices.
# Therefore, fixing them and assigning for a contiguous range of (i,j,k) does not actually
# assign in a contiguous memory block.  On the other hand, if we index paramKd as
# paramKd[i,j,k,v,w], for fixed v and w a contiguous range of (i,j,k) actually assigns in
# a more contiguous memory block.  (The entries with adjacent i for the same j,k,v,w are
# actually adjacent in the memory space.)  Therefore, the latter indexing scheme results in
# faster performance.  In my experiment of assigning 8 objects in 3D, the latter turns out
# to be about 30% faster.
#
# For this reason, I use the less conventional indexing scheme of paramKd[i,j,k,v,w].

# About the material tensors stored as paramKd
#
# Keyword argument ncmp is type parameter Kf elsewhere.
#
# paramKd[i,j,k,:,:] is a 3×3 tensor for ncmp = 3.  Suppose we are dealing with electric
# materials. Then this is the ε tensor.  It is important to note that not all the entries of
# this tensor are measured at the same physical locations.
#
# For example, if Ex, Ey, Ez are all tangential to primal grid planes (like in the standard
# Yee's grid construction), then εxx, εyy, εzz are defined at (i+½,j,k), (i,j+½,k), (i,j,k+½),
# respectively.  Also, six εvw's with v ≠ w are defined at (i,j,k).
#
# The size of paramKd in location dimensions (i.e., the length in the (i,j,k)-directions) is
# N.+1.  In order to simplify the assignment algorithm, paramKd is created with size N.+1 to
# match the sizes of the arrays created by create_oind_arrays.
#
# For length(N) = 2, the output paramKd is indexd as paramKd[i,j,v,w], where (i,j) is the
# grid cell location, and v and w are the row and column indices of the ncmp×ncmp tensorial
# material parameters (like the ε tensor and μ tensor).
create_param_array(N::SInt; ncmp::Int=3) = (s = (N.+1).data; zeros(ComplexF, s..., ncmp, ncmp))  # (i,j,...) element is ncmp×ncmp tensor

#  Create an array to store object indices.
create_oind_array(N::SInt) = zeros(ObjInd, (N.+1).data)

# Notes on ghost point transformation by boundary conditions
#
# - First of all, do not transform objects according to boundary conditions.  Periodization,
# for example, must be done by explicitly putting the translated object.  In other words,
# the solver periodizes the whatever composition inside the domain, but composing the space
# inside the domain is the user's responsibility.  (The user also must make sure that the
# surface normal is correctly calculated at the points on periodic boundaries.  This is
# achieved by making the object to extrude the periodic boundary sufficiently deep with the
# shape appearing on the opposite periodic boundary.  If the object extrudes the periodic
# boundary too shallow, with the extruded part is cut with a plane parallel to the domain
# boundary, then some boundary points within the object could be closer to the cut surface
# than the actual object surface, such that the boundary normal could be wrongly chosen as
# the surface normal.  It is safe, and also easy, to make the object to extrude a periodic
# boundary by the entire shape appearing on the opposite boundary, because then the boundary
# points cannot feel the effect of the termination of the object by periodic boundaries at
# all when calculating the surface normal.  This means that when mimicking a 2D problem by a
# 3D problem with one cell in the z-direction, putting prism-shaped objects with prism
# height of ∆z is not sufficient; it is safer to make the prism height 3∆z, as on the
# +z-boundary the object appearing on the opposite boundary has a thickness of ∆z, and on
# the -z-boundary the object appearing on the opposite boundary has also a thickness of ∆z.
# See example/usage_uni.jl.)
#
# - Transform ghost points back to the corresponding points inside the domain, and see which
# object is there.  Assign that object to the ghost points.
#
# - Ghost points are usually copied from different non-ghost points by translation (for
# the periodic BC) or reflection (for the symmetry BC), but NOT if their transformed points
# are themselves (e.g., primal ghost points on the symmetry BC).  For the latter case, the
# fields of interest on the ghost points are still inferrable by other means than spatial
# transformation (i.e., the fields are zero at primal ghost points on the symmetry BC).
# Therefore, at the equation solving step we don't need to keep the degrees of freedom for
# those ghost points.  Still, to smooth material parameterts at non-ghost points around
# ghost boundaries, we need objects assigned to ghost points (because ghost points are the
# corners of the voxels centered at those non-ghost points).  These ghost points are not
# copied from other non-ghost points by spatial transformation, so we need to assign objects
# to these ghost points.
#
# - When the ghost points are the copies of non-ghost points, both the fields and objects at
# the ghost points must be copied for consistency, not only one of them.  For example, if a
# periodic boundary condition is used, the primal grid point at the positive end must be
# copies of the primal grid point at the negative end, from which both the objects and
# fields must be copied to the ghost points.

# Notes on assignments and ghost points
#
# The goal is to "paint" grid points with a given object.  Initially I thought I would not
# need to paint ghost points independently, because I thought ghost points are by definition
# points whose properties are inferrable from other points.  This is not always the case.
# See the above "Notes on ghost point transformation by boundary conditions".
#
# If the ghost points are not exactly on the symmetry BC, we can recover the ghost point
# arrays from somewhere else.  The algorithm below prepares grid point indices such that we can
# assign objects the same way regardless of the kinds of boundary conditions.

# Notes on setting object indices instead of Objects themselves
#
# Why do we set up an array oindKd′ of object indices rather than an array of Objects
# themselves?  By setting up an array of object indices, we have to pass the maps from
# object indices to shapes (oind2shp), to parameter indices (oind2pind), and also a map from
# parameter indices to material parameters (pind2matprm) in order to get the shape and
# material parameter at each grid point.  If we set up an array of Objects, we would not
# need to pass all these maps, because an Object already contains the shape and material
# parameter.  Still, there are a few reasons to pass an array of object indices instead of
# an array of objects:
#
# - Dispatching the shape and material paramter from an Object turns out to be slow, because
# Objects have different concrete types depending on the type of the Shape it contains.
# This makes the multiple dispatch to kick in, which is slow.
#
# - Setting the argumenty types Base types (e.g., AbstractArray) rather custom types (e.g.,
# Object in the present context) increases the reusability of the code.  For example, there
# are cases where we want to pass a submatrix of the 3×3 material parameter tensor (rather
# than the material parameter tensor itself), because we want to treat the transverse and
# longitudinal components of the material parameter tensor separately (e.g., in the
# waveguide mode solver).  By taking the elementary types, this can be achieved without
# changing the assign_param! function, because we can create a vector of 2×2 submatrices of
# the material parameter tensors outside assign_param! and pass it.  If assign_param! took
# an array of Objects, handling such cases would be difficult because we created an Object
# with a 3×3 material parameter tensor in this case.

# Unlike smoothing.jl/smooth_param! and param.jl/create_paramop, the paramKd and oindKd′
# passed to assign_param! are usually (but not always; see benchmark/smoothing2d.jl) for
# different field types (e.g., paramKd for E-field and oindKd′ for the H-field).
#
# Explanation.  We set up oindKd′ such that it contains the voxel corner information to use
# in smoothing for determining whether paramKd at the voxel centers needs smoothing or not.
# For example, to determine whether paramKd for K = 3 defined at an Ew point needs smoothing
# or not, we take the voxel centered at the Ew point and examine the material parameter
# indices (pind) at the eight voxel corners.  If all the pind's at the eight voxel corners
# are the same, we assume that the voxel is filled with a uniform material parameter
# identified by the pind.  This means even though these voxel corners are the H-field points,
# we still need to examine the electric material properties at these voxel corners.
#
# Conversely, the Ew points, at which we are going to set paramKd, are the voxel corners of
# the voxel centered at the Hw points.  So, at the Ew points we want to set oindKd′ to use
# for smoothing the magnetic material parameterts.  This is why we need to pass paramKd and
# oindKd′ for different field types.  (Note that we assign paramKd and oindKd′ at the same
# locations, because we set them up point-by-point.)

# About Kf⏐₁
#
# The number of entries in the tuple oindKd′ is Kf⏐₁, which is Kf or 1.  (Kf is the dimesnion
# of the field and different from K, the dimension of space.)  Kf⏐₁ is the parameter used to
# control the assignment behavior:
# - For Kf⏐₁ == 1, the material parameter tensor is evaluated at the location specified by gt₀
# in each Yee cell, and the off-diagonal entries are assigned.
# - For Kf⏐₁ == Kf, the material parameter tensor is evaluated at the half-index shifted
# locations from gt₀ in each Yee cell, and at each half-index shifted location the
# corresponding diagonal entry is assigned.
#
# If Kf == 1, it is ambiguous which of the two assignment behaviors is taken, because Kf⏐₁ is
# 1 and Kf simultaneously.  In fact, a mixed behavior is chosen: the choice of the location
# to evaluate the material parameter tensor follows the behavior for Kf⏐₁ == 1 (i.e., does
# not shift from the location specified by gt₀), but the choice of tensor entries to assign
# follows the behavior for Kf⏐₁ == Kf (i.e., assign the diagonal entries, which actually
# correspond to the entire material parameter tensor (= scalar) in this case).

# Below, primed variables indicate that the variables may not be of the same material type
# as the non-primed variables.  Compare this with the usage of the prime in smoothing.jl
# that indicates voxel corner quantities instead of voxel center quantities.
function assign_param!(paramKd::AbsArrComplexF{K₊₂},  # electric (magnetic) parameter array to set; K₊₂ = K+2, where 2 is rank of material tensor
                       oindKd′::NTuple{Kf⏐₁,AbsArr{ObjInd,K}},  # object index arrays to set at electric (magnetic) parameter locations; Kf⏐₁ = Kf or 1
                       oind2shp::AbsVec{Shape{K,K²}},  # map from oind to shape; K² = K^2
                       oind2pind::AbsVec{ParamInd},  #  map from oind to electric (magnetic) material parameter index
                       pind2matprm::AbsVec{S²ComplexF{Kf,Kf²}},  # map from pind to electric (magnetic) material parameters; Kf² = Kf^2
                       gt₀::SVec{K,GridType},  # grid type of voxel corners; generated by ft2gt.(ft, boundft)
                       τl::Tuple2{NTuple{K,AbsVecReal}},  # field component locations transformed by boundary conditions
                       isbloch::SBool{K}  # boundary conditions
                       ) where {K,Kf,K²,Kf²,K₊₂,Kf⏐₁}
    @assert K²==K^2 && Kf²==Kf^2 && K₊₂==K+2 && (Kf⏐₁==Kf || Kf⏐₁==1)
    @assert size(paramKd,K+1)==size(paramKd,K+2)==Kf
    # `sub` is the subscripts for obj3d, pind3d, oind3d.
    #
    # Store circularly shifted subscripts in sub for Bloch boundary condition, such that
    # τl[sub[w]] is always sorted for all boundary conditions.  Sorted τl is necessary to
    # use findfirst and findlast in assign_param_obj!.
    #
    # Primal grid: ghost points are always at the positive end.
    # - Bloch: ghost points copy the values of non-ghost points at the negative end, so τl
    # must be circshifted to right in order to be sorted.
    # - symmetry: ghost points have own degrees of freedom at the positive end, so τl is
    # already sorted.
    # Therefore, the primal grid needs to be circshifted to right only for the Bloch BC.
    #
    # Dual grid: ghost points are always at the negative end.
    # - Bloch: ghost points copy the values of non-ghost points at the positive end, so τl
    # must be circshifted to left in order to be sorted.
    # - symmetry: ghost points copy the values of non-ghost points at the negative end, so τl is
    # already sorted.
    @inbounds M = length.(τl[nPR])  # NTuple{K,Int}: N.+1
    sub = (map((m,b)->circshift(1:m,b), M, isbloch.data),  # NTuple{K,VecInt}: primal grid
           map((m,b)->circshift(1:m,-b), M, isbloch.data))  # NTuple{K,VecInt}: dual grid

    # `psub` is the subscripts for paramKd.
    #
    # The subscripts for paramKd are a bit different from the subscripts for other arrays.
    # paramKd always have ghost cells at the positive end, regardless of the grid type.  (To
    # be more precise, only the first N portion of paramKd is used to construct the Maxwell
    # operator, and this used portion exclude the ghost points.)  Therefore, if the above
    # sub has ghost points' subscripts at the positive end, the psub below does not have to
    # be circshifted.  On the other hand, if the above sub has ghost points's subscripts at
    # the negative end, the psub below must be circshifted by +1 to push the ghost points to
    # the negative end.
    #
    # Bottom line: prepare subscripts such that circshifted paramKd and circshifted τl have
    # ghost points at the same indices.  Note that before circshift,
    # - paramKd has ghost cells at the positive end,
    # - the primal grid has ghost points at the positive end, and
    # - the dual grid has ghost points at the negative end.
    # This means that before circshift, primal paramKd and the primal τl have ghost points
    # at the same indices, whereas dual paramKd and the dual τl don't.  Therefore, primal
    # paramKd must be circshifted to right if the primal τl is circshifted to right.  (This
    # is to send paramKd's ghost points existing at the positive end to the negative end).
    # Similarly, dual paramKd must be circshifted to right if the dual τl is NOT circshifted
    # to left.  (This is also to send paramKd's ghost points existing at the positive end to
    # the negative end, where the ghost points of the non-circshifted dual τl exist).
    psub = (map((m,b)->circshift(1:m,b), M, isbloch.data),  # NTuple{K,VecInt}: primal grid
            map((m,b)->circshift(1:m,!b), M, isbloch.data))  # NTuple{K,VecInt}: dual grid

    # Perform assignment.
    for nw = 1:Kf⏐₁  # if Kf⏐₁ == Kf for Kf = 3, w = X̂, Ŷ, Ẑ
        # Set the grid types of the x-, y-, z-locations of Fw.
        #
        # I thought I could replace the condition Kf⏐₁==1 with isfield_ortho_shape(Kf,K),
        # because gt_cmp becomes gt₀ when isfield_ortho_shape(Kf,K) = true.  The replacement
        # is incorrect, however, because we need gt_cmp to become gt₀ even when the shape
        # and field subspaces coincide (e.g., K = Kf = 3) to assign the off-diagonal entries
        # of the material parameter tensors at voxel corners.
        gt_cmp = Kf⏐₁==1 ? gt₀ : gt_w(nw, gt₀)

        # Choose the circularly shifted subscripts to use.
        sub_cmp = t_ind(sub, gt_cmp)  # NTuple{K,VecInt}
        psub_cmp = t_ind(psub, gt_cmp)  # NTuple{K,VecInt}

        # Prepare the circularly shifted locations of the field components.
        τlcmp = view.(t_ind(τl,gt_cmp), sub_cmp)  # NTuple{K,VecFloat}: locations of Fw = Uw or Vw

        # Note that view() below is used to get a circularly shifted version of the
        # array, not a portion of.
        oind′_cmp = view(oindKd′[nw], sub_cmp...)  # circularly shifted nw-component of oindKd′

        # Below, we have four combinations depending on the values of Kf⏐₁ and Kf:
        # Kf⏐₁ is either ==Kf or ==1, and Kf is either ≥2 or ==1.
        #
        # The case of Kf⏐₁ == Kf covers Kf≥2 and Kf==1.  When Kf==1, the case
        # corresponds to Kf⏐₁==1 and Kf==1.  Therefore, the "if" case Kf⏐₁ == Kf
        # covers the two combinations:
        # - Kf⏐₁==Kf and Kf≥2,
        # - Kf⏐₁==1 and Kf==1.
        #
        # Therefore, the "else" case should cover the remaining two combinations:
        # - Kf⏐₁==Kf && Kf==1,
        # - Kf⏐₁==1 && Kf≥2.
        # However, the first combination is equivalent to Kf⏐₁==1 && Kf==1, which
        # was already covered in the "if" case.  Therefore, the "else" case
        # deals with only the case of Kf⏐₁==1 && Kf≥2.
        if Kf⏐₁ == Kf  # Kf⏐₁==Kf && Kf≥2, or Kf⏐₁==1 && Kf==1
            # Overwrite the (nw,nw)-diagonal entry of paramKd.
            param_cmp = view(paramKd, psub_cmp..., nw, nw)  # circularly shifted (nw,nw)-component of paramKd; ndims(param_cmp) = 3
            pind2matprm_w = [mp[nw,nw] for mp = pind2matprm]
            assign_param_cmp!(param_cmp, oind′_cmp, oind2shp, oind2pind, pind2matprm_w, τlcmp)
        else  # Kf⏐₁==1 && Kf≥2
            # Overwrite the off-diagonal entires of paramKd.
            param_cmp = view(paramKd, psub_cmp..., 1:Kf, 1:Kf)  # circularly shifted paramKd; ndims(param_cmp) = K+2
            assign_param_cmp!(param_cmp, oind′_cmp, oind2shp, oind2pind, pind2matprm, τlcmp)  # this does nothing for Kf = 1
        end
    end

    return nothing
end

# This is the innermost function that is customized for MaxwellBase (i.e., it takes arrays
# of specific element types).
#
# The goal of this function is to set all the arrays' entries, where the arrays are
# param_cmp and oind′_cmp at the same locations (τlcmp).  For each shape we iterate over
# grid points and test if each point is included in the shape.  Because there are so many
# grid points, the total time taken for these tests is somewhat substantial, so we don't
# want to test the inclusion of the same point again.  Therefore, for a point that is
# confirmed included in the shape, we want to set up all the arrays that can be set up.
# If the point is an Ew point, it is where the electric material properties ε_ww is set up
# in paramKd, but it is also where the object index for smoothing magnetic material
# parameters is set up in oindKd′, because this Ew point is the corners of the voxel whose
# center is at an Hw point, and we use the material parameter index at these voxel corners
# to determine whether the voxel centered at the Hw point is filled with a single magnetic
# material.
#
# For Kf⏐₁ == Kf.
assign_param_cmp!(param_cmp::AbsArrComplexF{K},  # output material parameter array (array of scalars)
                  oind′_cmp::AbsArr{ObjInd,K},  # output object index array
                  oind2shp::AbsVec{Shape{K,K²}},  # input map from oind to shape
                  oind2pind::AbsVec{ParamInd},  # input map from oind to pind
                  pind2matprm::AbsVecComplexF,  #  input map from pind to material parameter
                  τlcmp::NTuple{K,AbsVecReal}) where {K,K²} = # location of field components
    assign_param_cmp_impl!(param_cmp, oind′_cmp, oind2shp, oind2pind, pind2matprm, τlcmp)

# For Kf⏐₁ == 1.
# Below, note that the element type of pint2matprm, S²ComplexF{Kf,Kf²}, is always a matrix
# (i.e., rank-2 tensor), regardless of the value of Kf.  For example, if we are assigning μz
# in a 2D TE problem, μz is a type of S²ComplexF{1,1}.  This is why K₊₂ is K + 2 regardless
# of the value of Kf, where 2 is the rank of the material parameter tensor.
assign_param_cmp!(param_cmp::AbsArrComplexF{K₊₂},  # output material parameter array (array of tensors); M = K+2
                  oind′_cmp::AbsArr{ObjInd,K},  # output object index array
                  oind2shp::AbsVec{Shape{K,K²}},  # input map from oind to shape
                  oind2pind::AbsVec{ParamInd},  # input map from oind to pind
                  pind2matprm::AbsVec{S²ComplexF{Kf,Kf²}},  #  input map from pind to material parameter
                  τlcmp::NTuple{K,AbsVecReal}) where {K,Kf,K²,Kf²,K₊₂} = # location of field components
    assign_param_cmp_impl!(param_cmp, oind′_cmp, oind2shp, oind2pind, pind2matprm, τlcmp)

function assign_param_cmp_impl!(param_cmp, oind′_cmp, oind2shp, oind2pind, pind2matprm, τlcmp)
    for oind = 1:length(oind2shp)  # last index of oind2shp is index of last object; see object.jl/add!
        # Prepare values to set.
        @inbounds shp = oind2shp[oind]
        @inbounds pind = oind2pind[oind]
        @inbounds param = pind2matprm[pind]

        # Set the values to the arrays.
        #
        # Below, each of the two arrays oind′_cmp and param_cmp is individually set by
        # assign_val_shape!.  This means that for each point p and shape sh, the same test
        # p ϵ sh is repeated two times.  This is inefficient. An alternative implementation
        # is to create a function similar to assign_val_shape! that takes all these two
        # arrays, test p ϵ sh once, and set all the arrays simultaneously.  This may seem
        # more efficient, but this turns out to perform a lot of dynamic dispatches of
        # assign_val!, because two different realizations of assign_val! depending on the
        # output array type are called for every p ∈ sh in the for loop.
        #
        # A huge number of dynamic dispatches can be avoided by setting up each array
        # individually as below.  Then, the output array type of assign_val! is determined
        # at the interface of assign_val_shape!, so the dynamic dispatch of assign_val!
        # occurs only once per assign_val_shape!.
        assign_val_shape!(oind′_cmp, ObjInd(oind), shp, τlcmp)

        # Assign parameters at Yee's field points where the v-component E and H affects
        # the w=v components of D and B.  This will use assign_val(..., scalar, ...) for
        # setting diagonal entries of the material parameter tensor.
        assign_val_shape!(param_cmp, param, shp, τlcmp)
    end

    return nothing
end

# The only role of these wrappers assign_val_shape! is to limit the types of array and val
# to either (AbsArr{T,3}, T) or (AbsArr{T,5}, AbsMat{T}).  There is no performance benefit.
assign_val_shape!(array::AbsArr{T,K}, val::T, shape::Shape{K}, τlcmp::NTuple{K,AbsVecReal}) where {T,K} =
    assign_val_shape_impl!(array, val, shape, τlcmp)

assign_val_shape!(array::AbsArr{T,K₊₂}, val::SMatrix{Kf,Kf,T}, shape::Shape{K}, τlcmp::NTuple{K,AbsVecReal}) where {T,K,Kf,K₊₂} =
    assign_val_shape_impl!(array, val, shape, τlcmp)

# Given a shape, assign value at the points within the shape.
function assign_val_shape_impl!(array::AbsArr{T},
                                val::Union{T,SMatrix{Kf,Kf,T}},
                                shape::Shape{K},
                                τlcmp::NTuple{K,AbsVecReal}) where {T,K,Kf}
    # Set the location indices of object boundaries.
    @assert all(issorted.(τlcmp))
    bn, bp = bounds(shape)  # (SVec{K}, SVec{K})

    # Below, if shape is beyond the positive end of the domain, then findlast(l.≤b) = length(l)
    # but findfirst(l.≥b) = nothing.  To prevent any assignment, we have to replace nothing
    # with length(l)+1.  Then, the index range becomes length(l)+1:length(l), over which no
    # iteration occurs.  Similar consideration for shape beyond the negative end of the
    # domain.
    subn = map((l,b) -> (n = findfirst(l.≥b); n==nothing ? length(l)+1 : n), τlcmp, bn.data)  # NTuple{K,Int}
    subp = map((l,b) -> (n = findlast(l.≤b); n==nothing ? 0 : n), τlcmp, bp.data)  # NTuple{K,Int}
    CI = CartesianIndices(map((nᵢ,nₑ) -> nᵢ:nₑ, subn, subp))  # CartesianIndices{K}

    # Below, I think using the type assertion in (shape::Cuboid{K,K*K}).p achieves type
    # stability of p, because Cuboid{K,L}.p is a type of SMatrix{N,N,Float64,L}.
    if shape isa Cuboid{K,K*K} && (shape::Cuboid{K,K*K}).p == LinearAlgebra.I  # shape is Cartesian box
        assign_val!(array, val, CI)
    else  # shape is not Cartesian box
        for ci = CI
            pt = t_ind(τlcmp, ci)
            if pt ∈ shape
                assign_val!(array, val, ci)
            end  # if pt ∈ shape
        end  # for ci = ...
    end  # if shape isa ...

    return nothing
end

# Could be named Base.setindex!, but didn't want this to be exported, so named different.
function assign_val!(array::AbsArr{T,K}, scalar::T, ci::CartesianIndex{K}) where {T,K}
    @inbounds array[ci] = scalar
    return nothing
end

function assign_val!(array::AbsArr{T,K}, scalar::T, CI::CartesianIndices{K}) where {T,K}
    @inbounds array[CI] .= Ref(scalar)
    return nothing
end

# Similar to assign_val!(..., scalar, ...), but set the off-diagonal entries of a given
# tensor.
function assign_val!(array::AbsArr{T,K₊₂}, tensor::SMatrix{Kf,Kf,T}, ci::CartesianIndex{K}) where {T,K,Kf,K₊₂}  # K₊₂ = K+2
    # Below the main diagonal.
    for c = 1:Kf, r = c+1:Kf  # column- and row-indices
        @inbounds array[ci,r,c] = tensor[r,c]
    end

    # Above the main diagonal.
    for c = 2:Kf, r = 1:c-1  # column- and row-indices
        @inbounds array[ci,r,c] = tensor[r,c]
    end

    return nothing
end

function assign_val!(array::AbsArr{T,K₊₂}, tensor::SMatrix{Kf,Kf,T}, CI::CartesianIndices{K}) where {T,K,Kf,K₊₂}  # K₊₂ = K+2
    # Below the main diagonal.
    for c = 1:Kf, r = c+1:Kf, ci = CI  # column- and row-indices and location index
        @inbounds array[ci,r,c] = tensor[r,c]
    end

    # Above the main diagonal.
    for c = 2:Kf, r = 1:c-1, ci = CI  # column- and row-indices and location index
        @inbounds array[ci,r,c] = tensor[r,c]
    end

    return nothing
end