text:
println("# === Below Generated by ",PROGRAM_FILE," ===")
println("")

function arg_name(m)
  m.captures[3]
end

function arg_value(m)
  #if m.captures[2] == "*" && m.captures[1] == "char" return "String($(m.captures[3]))" end
  if m.captures[3] == "env"
      "penv"
  else
     m.captures[3]
  end
end

decl_arg_type(m) = decl_arg_type(m.captures[1], m.captures[2])
ccall_arg_type(m; r=false) = ccall_arg_type(m.captures[1], m.captures[2], r=r)

function decl_arg_type(t, s)
  if s == "*"
    if t == "char"
      return "AbstractString"
    elseif t == "void"
      return "Ptr{Nothing}"
    elseif t == "JNIEnv"
      return "Ptr{JNIEnv}"
    else
      return "Array{$t,1}"
#      return "Ptr{$t}"
    end
  elseif t == "jsize" #|| t == "jint" || t == "jlong" || t == "jshort" || t == "jbyte"
    return Integer
  elseif t == "jobject"
    return "jobject"
  elseif t == "jobjectArray"
    return "jobjectArray"
  end

  if t == "void"
    return "Nothing"
  end

  return s == "" ? t : "Array{$t,1}"
#  return s == "" ? t : "Ptr{$t}"
end

function ccall_arg_type(t, s; r=false)
  if s == "*"
    if t == "char"
      return "Cstring"
    elseif t == "void"
      return "Ptr{Nothing}"
    elseif t == "JNIEnv"
      return "Ptr{JNIEnv}"
    else
      return "Ptr{$t}"
    end
  end

  if t == "void"
    return "Nothing"
  end

  # No asterisk: type
  # If return type, Ptr
  # Else Array
  return s == "" ? t : r ? "Ptr{$t}" : "Array{$t,1}"
end

# julia_arg(m) = string(arg_name(m), "::", decl_arg_type(m))
function julia_arg(m)
    if arg_name(m) == "isCopy"
        "isCopy::Ptr{jboolean}"
    elseif arg_name(m) == "elems"
        string(arg_name(m), "::", "Ptr{$(m.captures[1])}")
    else
        string(arg_name(m), "::", decl_arg_type(m))
    end
end

function to_julia_fnname(fn)::String
  #lowercase(fn)
  stringbuilder = IOBuffer()
  inuppercase = false
  for c in fn
    if !inuppercase && isuppercase(c)
      print(stringbuilder, '_')
      inuppercase = true
    elseif inuppercase && islowercase(c)
      inuppercase = false
    end
    print(stringbuilder, lowercase(c))
  end
  String(take!(stringbuilder))[2:end]
end

for line in open(readlines, "Interfaces.jl", "r")
  # m: match comments
  # Example:
  # # jclass ( *GetSuperclass) (JNIEnv *env, jclass sub);
  # Group 1: Return type ((?:void|char|j\w+)) "jclass"
  # Group 2: Asterisk (\**) ""
  # Group 3: Function name (\w+) "GetSuperclass"
  # Group 4: Arguments \((.*)\) "JNIEnv *env, jclass sub"
  m = match(r"\# \s* (?:const\s*)? ((?:void|char|j\w+)) \s* (\**) \s* \( \s* \* (\w+) \s* \) \s* \((.*)\) \s* ;"x, line)
  if m === nothing continue end

  # Ignore functions with variable argument syntax
  # Only process vararg functions that end in A
  if occursin("...", m.captures[4]) continue end
  if occursin("va_list", m.captures[4]) continue end

  # Get return type
  rtype = ccall_arg_type(m.captures[1], m.captures[2], r=true)

  # Function name
  fname = m.captures[3]
  # Split arguments
  args = split(m.captures[4], ",")

  # mm: Analyze arguments
  # Example: [1] "JNIEnv *env", [2] "jclass sub"
  # Group 1: Argument type ((?:void|j\w+|char|JNI\w+|JavaVM))
  # Group 2: Asterisk (\**)
  # Group 3: Argument name (\w+)
  mm = map(x->match(r"^\s* (?:const\s+)? \s* ((?:void|j\w+|char|JNI\w+|JavaVM)) \s*? (\**) \s* (\w+) \s*$"x, x), args)

  julia_fnname = to_julia_fnname(fname)
  # skip the JNIEnv arg for julia since it is passed as a global to ccall
  julia_args = join(push!(map(julia_arg, mm)[2:end],"penv=jnienvptrs[]"), ", ")
  arg_types = join(map(ccall_arg_type, mm), ", ")
  arg_names = join(map(arg_value, mm), ", ")

  # Commented out export command
  # print("#export $fname\n")
  print("$(julia_fnname)($julia_args) =\n  ccall(jninativeinterfaceref[].$(fname), $rtype, ($arg_types,), $arg_names)\n\n")
end

# println("end")
println("")
println("# === Above Generated by ",PROGRAM_FILE," ===")


text:
type AtomicUnitSet
    α::Int
    atoms::Vector{AtomicUnit}

    AtomicUnitSet(α::Int) = new(α, AtomicUnit[])
end

typealias AUSet AtomicUnitSet

A⁽⁸⁾ = AtomicUnitSet(8)

function Base.show(io::IO, A::AtomicUnitSet)
    print(io, "AtomicUnitSet(α = $(A.α), atoms = $(A.atoms))")
end

function Base.push!(A::AtomicUnitSet, a::AtomicUnit)
    a.e ≤ A.α/2 ? push!(A.atoms, a) : error("α = $(A.α) ∴ a ∉ A^(α)")    # how to render?: A⁽ᵅ⁾ (A^(α))
end

type CovalentBonds
end

abstract MolecularUnit

type Molecule <: MolecularUnit
    A::AtomicUnitSet
    B::CovalentBonds
    state::Symbol    # :stable or :unstable

end

# molecular unit behaviour
phi{T<:Real}(x::T, d::Int, vH) = sum([a*x^i for i=1:d])
ϕ = phi

ext: jl


text:
include("bayesquad.jl")
include("logbayesquad.jl")

function check_kernel_parameters(l, σ)
    σ > 0 || throw(ArgumentError("σ should be positive"))
    l isa AbstractMatrix && (
        l isa LowerTriangular || throw(
            ArgumentError(
                "For l an AbstractMatrix, only LowerTriangular matrices are accepted"
            ),
        )
    )
    l isa Union{AbstractVector,Real} && (
        all(>=(0), l) || throw(ArgumentError("At least one lengthscale value is negative"))
    )
    return nothing
end

function get_kernel_params(k::Kernel; kwargs...)
    return (k, (; kwargs...))
end

function get_kernel_params(k::ScaledKernel; kwargs...)
    # return get_kernel_params(k.kernel; merge((;kwargs...), (σ=first(k.σ²),)))
    return get_kernel_params(k.kernel; kwargs..., σ=first(k.σ²))
end

function get_kernel_params(k::TransformedKernel; kwargs...)
    check_transform(k.transform)
    return get_kernel_params(k.kernel; kwargs..., l=transform_param(k.transform))
end

function check_transform(transform)
    transform isa Union{ScaleTransform,ARDTransform,LinearTransform} ||
        throw(ArgumentError("No lengthscale could be extracted from kernel,\n
        only ScaleTransform, ARDTransform and LinearTransform are allowed"))
    return nothing
end

function kernel(b::AbstractBQ)
    if b.l isa AbstractMatrix
        b.σ * (b.kernel ∘ LinearTransform(inv(b.l)))
    else
        return b.σ * (b.kernel ∘ ScaleTransform(inv.(b.l)))
    end
end


text:
function demoHMMburglar()

    # setup a grid representation of the rooom
    Gx = 5 # two dimensional grid size
    Gy = 5
    H = Gx * Gy # number of states on grid

    # matrix representing the possible states of the system
    st = reshape(1:H, Gx, Gy) # assign each grid point a state

    # make a deterministic state transition matrix HxH on a 2D grid:
    phgh = zeros(H, H) # transition from state j to state i
    for x = 1:Gx
        for y = 1:Gy
            # from the current state-cell (j coord in state transition matrix)
            # to next state-cell on (i coord in state transition matrix):
            #   the next row, same column
            if validgridposition(x + 1, y, Gx, Gy) # sample for x = 1, y = 2
                phgh[st[x + 1, y], st[x, y]] = 1 # 2,2=7  1,2=6
            end
            #   the previous row, same column
            if validgridposition(x - 1, y, Gx, Gy)
                phgh[st[x - 1, y], st[x, y]] = 1 # 0,2  1,2
            end
            #   the same row, next column
            if validgridposition(x, y + 1, Gx, Gy)
                phgh[st[x, y + 1], st[x, y]] = 1 # 1,3=11 1,2=6
            end
            #   the same row, previous column
            if validgridposition(x, y - 1, Gx, Gy)
                phgh[st[x, y - 1], st[x, y]] = 1 # 1,1=1 1,2=6
            end
        end
    end
    # conditional distribution from state transition matrix
    phghm = condp(phgh) # matrix with sum(phghm, 1) = 1 with phghm[i, j] = p(hg=i | hm=j)
    ph1=condp(ones(H,1)) # initialise probabilities for the states of the hidden variable at timestep 1
    pvgh=zeros(4,H) # initialise emision matrix

    pv1gh = 0.01 * ones(1,H); r = randperm(H); pv1gh[r[1:10]] = 0.9; # Creaks in 10 randomly chosen cells
    pv2gh = 0.01 * ones(1,H); r = randperm(H); pv2gh[r[1:10]] = 0.9; # Bumps in 10 randomly chosen cells

    PyPlot.figure()
    axc = PyPlot.subplot(2, 1, 1)
    axc[:set_title]("creaks layout")
    PyPlot.imshow(reshape(pv1gh, Gx, Gy), cmap="bone");
    axb = PyPlot.subplot(2, 1, 2)
    axb[:set_title]("bumps layout")
    PyPlot.imshow(reshape(pv2gh, Gx, Gy), cmap="bone");
    map([axb, axc]) do axesfig
        PyPlot.setp(axesfig[:get_xticklines](), visible=false)
        PyPlot.setp(axesfig[:get_xticklabels](), visible=false)
        PyPlot.setp(axesfig[:get_yticklines](), visible=false)
        PyPlot.setp(axesfig[:get_yticklabels](), visible=false)
    end

    # Form the joint distribution p(v|h)=p(v1|h)p(v2|h)
    # v = (v1, v2) and v1 and v2 are independent given h
    vv = zeros(4, 2)
    for i = 1:4
        pvgh[1, :] =  pv1gh .* pv2gh; vv[1, :] = [1 1]; # p(v1=1|h)*p(v2=1|h)
        pvgh[2, :] =  pv1gh.*(1-pv2gh); vv[2, :] = [1 2]; # p(v1=1|h)*p(v2=1|h)
        pvgh[3, :] =  (1-pv1gh).*pv2gh; vv[3, :] = [2 1]; # p(v1=1|h)*p(v2=1|h)
        pvgh[4, :] =  (1-pv1gh).*(1-pv2gh); vv[4, :] = [2 2]; # p(v1=1|h)*p(v2=1|h)
    end

    # draw some random samples:
    T=10
    h = zeros(Integer, 1, T) # holds the state value for the hidden variable at a specific timestep
    v = zeros(Integer, 1, T) # holds the state value for the visible variable at a specific timestep

    h[1]=randgen(ph1) # initialize the hidden variable @t=1 with a random state based on ph1 distribution
    v[1]=randgen(pvgh[:, h[1]]) # initialize the visible variable @t=1 with a random state based on pvgh( vg | h@t=1)

    for t=2:T
        h[t]=randgen(phghm[:, h[t-1]]) # set the hidden variable state @t based on h@t-1 using the transition matrix
        v[t]=randgen(pvgh[:, h[t]]) # set the visible variable state @t based on h@t using the emission matrix
    end

    # Perform inference based on the observed v:
    (alpha, loglik) = HMMforward(v, phghm, ph1, pvgh) # filtering
    phtgV1t = alpha # filtered posterior - infer the current hidden state p(ht | v_1:t)

    phtgV1T = HMMgamma(alpha, phghm) # Smoothed Burglar distribution
    maxstate, logprob = HMMviterbi(v, phghm, ph1, pvgh) # Most likely Burglar path

    PyPlot.figure()
    for t = 1:T
        axg = PyPlot.subplot(5, T, t); PyPlot.imshow(repmat(vv[v[t], :], 2, 1), cmap="bone");
        if t == 2 # used t == 1 or t == 2 for title alignment only
            axg[:set_title]("Creaks and Bumps")
        end
        # add Filtering data row of T images from the previous row offset
        axf = PyPlot.subplot(5, T, T+t);  PyPlot.imshow(reshape(phtgV1t[:, t], Gx, Gy), cmap="bone");
        if t == 1
            axf[:set_title]("Filtering")
        end
        # add Smoothing data row of T images from the previous row offset
        axs = PyPlot.subplot(5, T, 2*T+t);  PyPlot.imshow(reshape(phtgV1T[:, t], Gx, Gy), cmap="bone");
        if t == 1
            axs[:set_title]("Smoothing")
        end
        z=zeros(H,1); z[maxstate[t]]=1;
        # add Viterbi data row of T images from the previous row offset
        axv = PyPlot.subplot(5,T,3*T+t); PyPlot.imshow(reshape(z,Gx,Gy), cmap="bone")
        if t == 1
            axv[:set_title]("Viterbi")
        end
        z = zeros(H,1); z[h[t]] = 1;
        # add true data row of T images from the previous row offset
        axt = PyPlot.subplot(5,T,4*T+t); PyPlot.imshow(reshape(z,Gx,Gy), cmap="bone")
        if t == 2
            axt[:set_title]("True Burglar position")
        end
        map([axg, axf, axs, axv, axt]) do axesfig
            PyPlot.setp(axesfig[:get_xticklines](), visible=false)
            PyPlot.setp(axesfig[:get_xticklabels](), visible=false)
            PyPlot.setp(axesfig[:get_yticklines](), visible=false)
            PyPlot.setp(axesfig[:get_yticklabels](), visible=false)
        end
    end

end


text:
module One

str = "hello this is a string!"
str = 'hello this is a string!'
str = """hello this is a string!"""
str = '''hello this is a string!'''

type Point
  x::Float32
  y::Float32
  Point(x,y) = new(x,y)
end

x->begin
   c = A + B;
end

end


text:
module DiBitVectors

import Base: getindex, setindex!, size, length, checkbounds, push!, pop!, zero

"""
    DiBitVector <: AbstractVector{UInt8}

A bitvector whose elements are two bits wide, allowing
storage of integer values between 0 and 3. Optimized
for performance and memory savings.
"""
mutable struct DiBitVector <: AbstractVector{UInt8}
    data::Vector{UInt64}
    len::UInt

    function DiBitVector(n::Integer, v::Integer)
        if !(Int(v) in 0:3)
            throw(ArgumentError("v must be in 0:3"))
        end
        fv = (0x0000000000000000, 0x5555555555555555,
        0xaaaaaaaaaaaaaaaa, 0xffffffffffffffff)[v + 1]
        vec = Vector{UInt64}(undef, cld(n, 32))
        fill!(vec, fv)
        return new(vec, n % UInt)
    end
end

@inline checkbounds(D::DiBitVector, n::Integer) =  0 < n ≤ length(D.data) << 5 || throw(BoundsError(D, n))

DiBitVector(n::Integer) = DiBitVector(n, 0)
DiBitVector() = DiBitVector(0, 0)

@inline Base.length(x::DiBitVector) = x.len % Int
@inline Base.size(x::DiBitVector) = (length(x),)

@inline index(n::Integer) = ((n-1) >>> 5) + 1
@inline offset(n::Integer) = ((n-1) << 1) & 63

@inline function Base.getindex(x::DiBitVector, i::Int)
    @boundscheck checkbounds(x, i)
    return UInt8((@inbounds x.data[index(i)] >>> offset(i)) & 3)
end

@inline function unsafe_setindex!(x::DiBitVector, v::UInt64, i::Int)
    bits = @inbounds x.data[index(i)]
    bits &= ~(UInt(3) << offset(i))
    bits |= convert(UInt64, v) << offset(i)
    @inbounds x.data[index(i)] = bits
end

@inline function Base.setindex!(x::DiBitVector, v::Integer, i::Int)
    v & 3 == v || throw(DomainError("Can only contain 0:3 (tried $v)"))
    @boundscheck checkbounds(x, i)
    unsafe_setindex!(x, convert(UInt64, v), i)
end

@inline function Base.push!(x::DiBitVector, v::Integer)
    len = length(x)
    len == length(x.data) << 5 && push!(x.data, zero(UInt))
    x.len = (len + 1) % UInt
    x[len+1] = convert(UInt64, v)
    return x
end

@inline function Base.pop!(x::DiBitVector)
    x.len == 0 && throw(ArgumentError("array must be non-empty"))
    v = x[end]
    x.len = (x.len - 1) % UInt
    x.len == (length(x.data) -1) << 5 && pop!(x.data)
    return v
end

@inline zero(x::DiBitVector) = DiBitVector(x.len, 0)

export DiBitVector

end # module


text:
using LoopVectorization, Random, Test
T = Float32

@testset "ifelse (masks)" begin

    function promote_bool_store!(z, x, y)
        for i ∈ eachindex(x)
            z[i] = (x[i]*x[i] + y[i]*y[i]) < 1
        end
        z
    end
    function promote_bool_storeavx!(z, x, y)
        @avx for i ∈ eachindex(x)
            z[i] = (x[i]*x[i] + y[i]*y[i]) < 1
        end
        z
    end
    function promote_bool_store_avx!(z, x, y)
        @_avx for i ∈ eachindex(x)
            z[i] = (x[i]*x[i] + y[i]*y[i]) < 1
        end
        z
    end
    # @macroexpand @_avx for i ∈ eachindex(x)
            # z[i] = (x[i]*x[i] + y[i]*y[i]) < 1
        # end
    function promote_bool_storeavx2!(z, x, y)
        @avx for i ∈ eachindex(x)
            z[i] = (x[i]*x[i] + y[i]*y[i]) < 1 ? 1 : 0
        end
        z
    end
    function promote_bool_store_avx2!(z, x, y)
        @_avx for i ∈ eachindex(x)
            z[i] = (x[i]*x[i] + y[i]*y[i]) < 1 ? 1 : 0
        end
        z
    end

    function Bernoulli_logit(y::BitVector, α::AbstractVector{T}) where {T}
        t = zero(promote_type(Float64,T))
        @inbounds for i ∈ eachindex(α)
            invOmP = 1 + exp(α[i])
            nlogOmP = log(invOmP)
            nlogP = nlogOmP - α[i]
            t -= y[i] ? nlogP : nlogOmP
        end
        t
    end
    function Bernoulli_logitavx(y::BitVector, α::AbstractVector{T}) where {T}
        t = zero(T === Int32 ? Float32 : Float64)
        @avx for i ∈ eachindex(α)
            invOmP = 1 + exp(α[i])
            nlogOmP = log(invOmP)
            nlogP = nlogOmP - α[i]
            t -= y[i] ? nlogP : nlogOmP
        end
        t
    end
    function Bernoulli_logit_avx(y::BitVector, α::AbstractVector{T}) where {T}
        t = zero(T === Int32 ? Float32 : Float64)
        @_avx for i ∈ eachindex(α)
            invOmP = 1 + exp(α[i])
            nlogOmP = log(invOmP)
            nlogP = nlogOmP - α[i]
            t -= y[i] ? nlogP : nlogOmP
        end
        t
    end

    function addormul!(c, a, b)
        for i ∈ eachindex(c,a,b)
            c[i] = a[i] > b[i] ? a[i] + b[i] : a[i] * b[i]
        end
    end
    function addormul_avx!(c, a, b)
        @_avx for i ∈ eachindex(c,a,b)
            c[i] = a[i] > b[i] ? a[i] + b[i] : a[i] * b[i]
        end
    end
    function addormulavx!(c, a, b)
        @avx for i ∈ eachindex(c,a,b)
            c[i] = if a[i] > b[i]
                a[i] + b[i]
            else
                a[i] * b[i]
            end
        end
    end
    function addormulp1!(c, a, b)
        for i ∈ eachindex(c,a,b)
            c[i] = 1 + (a[i] > b[i] ? a[i] + b[i] : a[i] * b[i])
        end
    end
    function addormulp1_avx!(c, a, b)
        @_avx for i ∈ eachindex(c,a,b)
            c[i] = 1 + (a[i] > b[i] ? a[i] + b[i] : a[i] * b[i])
        end
    end
    function addormulp1avx!(c, a, b)
        @avx for i ∈ eachindex(c,a,b)
            a_greater_than_b = a[i] > b[i]
            c[i] = 1 + (a_greater_than_b ? a[i] + b[i] : a[i] * b[i])
        end
    end
    function addifelsemul_avx!(c, a, b)
        @_avx for i ∈ eachindex(c,a,b)
            c[i] = ifelse(a[i] > b[i], a[i] + b[i], a[i] * b[i])
        end
    end
    function addifelsemulavx!(c, a, b)
        @avx for i ∈ eachindex(c,a,b)
            c[i] = ifelse(a[i] > b[i], a[i] + b[i], a[i] * b[i])
        end
    end
    function addifelsemulp1_avx!(c, a, b)
        @_avx for i ∈ eachindex(c,a,b)
            c[i] = 1 + ifelse(a[i] > b[i], a[i] + b[i], a[i] * b[i])
        end
    end
    function addifelsemulp1avx!(c, a, b)
        @avx for i ∈ eachindex(c,a,b)
            c[i] = 1 + ifelse(a[i] > b[i], a[i] + b[i], a[i] * b[i])
        end
    end
    function ifelseoverwrite!(p)
        for i ∈ eachindex(p)
            p[i] = p[i] < 0.5 ? p[i]^2 : p[i]^3
        end
    end
    function ifelseoverwriteavx!(p)
        @avx for i ∈ eachindex(p)
            p[i] = p[i] < 0.5 ? p[i]^2 : p[i]^3
        end
    end

    function maybewriteand!(c, a, b)
        @inbounds for i ∈ eachindex(c,a,b)
            a[i] > b[i] && (c[i] = a[i] + b[i])
        end
    end
    function maybewriteand_avx!(c, a, b)
        @_avx for i ∈ eachindex(c,a,b)
            a[i] > b[i] && (c[i] = a[i] + b[i])
        end
    end
    function maybewriteandavx!(c, a, b)
        @avx for i ∈ eachindex(c,a,b)
            a[i] > b[i] && (c[i] = a[i] + b[i])
        end
    end
    function maybewriteor!(c, a, b)
        @inbounds for i ∈ eachindex(c,a,b)
            a[i] > b[i] || (c[i] = a[i] ^ b[i])
        end
    end
    function maybewriteor_avx!(c, a, b)
        @_avx for i ∈ eachindex(c,a,b)
            a[i] > b[i] || (c[i] = a[i] ^ b[i])
        end
    end
    function maybewriteoravx!(c, a, b)
        @avx for i ∈ eachindex(c,a,b)
            a[i] > b[i] || (c[i] = a[i] ^ b[i])
        end
    end
    function maybewriteor!(c::AbstractVector{<:Integer}, a, b)
        @inbounds for i ∈ eachindex(c,a,b)
            a[i] > b[i] || (c[i] = a[i] & b[i])
        end
    end
    function maybewriteor_avx!(c::AbstractVector{<:Integer}, a, b)
        @_avx for i ∈ eachindex(c,a,b)
            a[i] > b[i] || (c[i] = a[i] & b[i])
        end
    end
    function maybewriteoravx!(c::AbstractVector{<:Integer}, a, b)
        @avx for i ∈ eachindex(c,a,b)
            a[i] > b[i] || (c[i] = a[i] & b[i])
        end
    end

    function notacondload!(C, A, b)
        @inbounds for n ∈ 1:size(C,2), m ∈ 1:size(C,1)
            C[m,n] = A[m,n] * (b[n] > 0 ? b[n] : -b[n])
        end
    end
    function notacondloadavx!(C, A, b)
        @avx for n ∈ 1:size(C,2), m ∈ 1:size(C,1)
            C[m,n] = A[m,n] * (b[n] > 0 ? b[n] : -b[n])
        end
    end
    function condloadscalar!(C, A, c, b)
        @inbounds for n ∈ 1:size(C,2), m ∈ 1:size(C,1)
            C[m,n] = A[m,n] * (c[n] > 0 ? b[n] : 1) + c[n]
        end
    end
    function condloadscalaravx!(C, A, c, b)
        @avx for n ∈ 1:size(C,2), m ∈ 1:size(C,1)
            C[m,n] = A[m,n] * (c[n] > 0 ? b[n] : 1) + c[n]
        end
    end
    function maskedloadscalar!(C, A, b)
        @inbounds for n ∈ 1:size(C,2), m ∈ 1:size(C,1)
            C[m,n] = A[m,n] * (A[m,n] > 0 ? b[n] : 1)
        end
    end
    function maskedloadscalaravx!(C, A, b)
        @avx for n ∈ 1:size(C,2), m ∈ 1:size(C,1)
            C[m,n] = A[m,n] * (A[m,n] > 0 ? b[n] : 1)
        end
    end
    function AtmulBpos!(C, A, B)
        @inbounds for n ∈ 1:size(C,2), m ∈ 1:size(C,1)
            Cₘₙ = zero(eltype(C))
            @simd ivdep for k ∈ 1:size(A,1)
                Cₘₙ += A[k,m] * B[k,n]
            end
            C[m,n] > 0 && (C[m,n] = Cₘₙ)
        end
    end
    function AtmulBposavx!(C, A, B)
        @avx for n ∈ 1:size(C,2), m ∈ 1:size(C,1)
            Cₘₙ = zero(eltype(C))
            for k ∈ 1:size(A,1)
                Cₘₙ += A[k,m] * B[k,n]
            end
            C[m,n] > 0 && (C[m,n] = Cₘₙ)
        end
    end
    function AtmulBpos_avx!(C, A, B)
        @_avx for n ∈ 1:size(C,2), m ∈ 1:size(C,1)
            Cₘₙ = zero(eltype(C))
            for k ∈ 1:size(A,1)
                Cₘₙ += A[k,m] * B[k,n]
            end
            C[m,n] > 0 && (C[m,n] = Cₘₙ)
        end
    end
    function condstore!(x)
        @inbounds for i ∈ eachindex(x)
            x1 = 2*x[i]-100
            x2 = x1*x1
            x3 = x2 + x1
            x4 = x3
            x[i] = x1
            (x1 < -50) && (x[i] = x2)
            (x1 < 60) || (x[i] = x4)
        end
    end
    function condstore1avx!(x)
        @avx for i ∈ eachindex(x)
            x1 = 2*x[i]-100
            x2 = x1*x1
            x3 = x2 + x1
            x[i] = x1
            (x1 < -50) && (x[i] = x2)
            (x1 < 60) || (x[i] = x3)
        end
    end
    function condstore1_avx!(x)
        @_avx for i ∈ eachindex(x)
            x1 = 2*x[i]-100
            x2 = x1*x1
            x3 = x2 + x1
            x[i] = x1
            (x1 < -50) && (x[i] = x2)
            (x1 < 60) || (x[i] = x3)
        end
    end
    function condstore2avx!(x)
        @avx for i ∈ eachindex(x)
            x1 = 2*getindex(x, i)-100
            x2 = x1*x1
            x3 = x2 + x1
            setindex!(x, x1, i)
            (x1 < -50) && setindex!(x, x2, i)
            (x1 < 60) || setindex!(x, x3, i)
        end
    end
    function condstore2_avx!(x)
        @_avx for i ∈ eachindex(x)
            x1 = 2*getindex(x, i)-100
            x2 = x1*x1
            x3 = x2 + x1
            setindex!(x, x1, i)
            (x1 < -50) && setindex!(x, x2, i)
            (x1 < 60) || setindex!(x, x3, i)
        end
    end
    function andorassignment!(x, y, z)
        @inbounds for i ∈ eachindex(x, y, z)
            yᵢ = y[i]
            zᵢ = z[i]
            (yᵢ > 0.5) || (yᵢ *= 2)
            (zᵢ < 0.5) && (zᵢ *= 2)
            x[i] = yᵢ * zᵢ
        end
    end
    function andorassignmentavx!(x, y, z)
        @avx for i ∈ eachindex(x, y, z)
            yᵢ = y[i]
            zᵢ = z[i]
            (yᵢ > 0.5) || (yᵢ *= 2)
            (zᵢ < 0.5) && (zᵢ *= 2)
            x[i] = yᵢ * zᵢ
        end
    end
    function andorassignment_avx!(x, y, z)
        @avx for i ∈ eachindex(x, y, z)
            yᵢ = y[i]
            zᵢ = z[i]
            (yᵢ > 0.5) || (yᵢ *= 2)
            (zᵢ < 0.5) && (zᵢ *= 2)
            x[i] = yᵢ * zᵢ
        end
    end

    function twoifelses!(res, half, m, keep=nothing, final=true)
        𝒶𝓍j=axes(half,1)
        for j in 𝒶𝓍j
            𝓇𝒽𝓈 = if isnothing(keep)
                log(half[j]) + m[j]
            else
                res[j] + (log(half[j]) + m[j])
            end
            res[j] = isnothing(final) ? 𝓇𝒽𝓈 : exp(𝓇𝒽𝓈)
        end
        res
    end
    function twoifelses_avx!(res, half, m, keep=nothing, final=true)
        𝒶𝓍j=axes(half,1)
        @avx for j in 𝒶𝓍j
            𝓇𝒽𝓈 = if isnothing(keep)
                log(half[j]) + m[j]
            else
                res[j] + (log(half[j]) + m[j])
            end
            res[j] = if isnothing(final)
                𝓇𝒽𝓈
            else
                exp(𝓇𝒽𝓈)
            end
        end
        res
    end

    N = 117
    for T ∈ (Float32, Float64, Int32, Int64)
        @show T, @__LINE__
        if T <: Integer
            a = rand(-T(100):T(100), N); b = rand(-T(100):T(100), N);
        else
            a = rand(T, N); b = rand(T, N);
        end;
        c1 = similar(a); c2 = similar(a);

        promote_bool_store!(c1, a, b);
        promote_bool_storeavx!(c2, a, b);
        @test c1 == c2
        fill!(c2, -999999999); promote_bool_store_avx!(c2, a, b);
        @test c1 == c2
        fill!(c2, -999999999); promote_bool_storeavx2!(c2, a, b);
        @test c1 == c2
        fill!(c2, -999999999); promote_bool_store_avx2!(c2, a, b);
        @test c1 == c2

        fill!(c1,  999999999); addormul!(c1, a, b)
        fill!(c2, -999999999); addormul_avx!(c2, a, b)
        @test c1 ≈ c2
        fill!(c2, -999999999); addormulavx!(c2, a, b)
        @test c1 ≈ c2
        fill!(c2, -999999999); addifelsemul_avx!(c2, a, b)
        @test c1 ≈ c2
        fill!(c2, -999999999); addifelsemulavx!(c2, a, b)
        @test c1 ≈ c2
        addormulp1!(c1, a, b)
        addormulp1_avx!(c2, a, b)
        @test c1 ≈ c2
        fill!(c2, -999999999); addormulp1avx!(c2, a, b)
        @test c1 ≈ c2
        fill!(c2, -999999999); addifelsemulp1_avx!(c2, a, b)
        @test c1 ≈ c2
        fill!(c2, -999999999); addifelsemulp1avx!(c2, a, b)
        @test c1 ≈ c2

        fill!(c1, -999999999); maybewriteand!(c1, a, b)
        fill!(c2, -999999999); maybewriteand_avx!(c2, a, b)
        @test c1 ≈ c2
        fill!(c2, -999999999); maybewriteandavx!(c2, a, b)
        @test c1 ≈ c2

        fill!(c1, -999999999); maybewriteor!(c1, a, b)
        fill!(c2, -999999999); maybewriteor_avx!(c2, a, b)
        @test c1 ≈ c2
        fill!(c2, -999999999); maybewriteoravx!(c2, a, b)
        @test c1 ≈ c2

        andorassignment!(c1, a, b);
        andorassignmentavx!(c2, a, b);
        @test c1 ≈ c2
        fill!(c2, -999999999); andorassignment_avx!(c2, a, b);
        @test c1 ≈ c2

        a1 = copy(a); a2 = copy(a);
        ifelseoverwrite!(a1)
        ifelseoverwriteavx!(a2)
        @test a1 ≈ a2

        if T <: Union{Float32,Float64}
            a .*= 100;
        end;
        b1 = copy(a);
        b2 = copy(a);
        condstore!(b1)
        condstore1avx!(b2)
        @test b1 == b2
        copyto!(b2, a); condstore1_avx!(b2)
        @test b1 == b2
        copyto!(b2, a); condstore2avx!(b2)
        @test b1 == b2
        copyto!(b2, a); condstore2_avx!(b2)
        @test b1 == b2

        M, K, N = 83, 85, 79;
        if T <: Integer
            A = rand(T(-100):T(100), K, M);
            B = rand(T(-100):T(100), K, N);
            C1 = rand(T(-100):T(100), M, N);
            b = rand(T(-100):T(100), N);
            d = rand(T(-100):T(100), N);
        else
            A = randn(T, K, M);
            B = randn(T, K, N);
            C1 = randn(T, M, N);
            b = randn(T, N);
            d = randn(T, N);
        end;
        C2 = copy(C1); C3 = copy(C1);
        AtmulBpos!(C1, A, B)
        AtmulBposavx!(C2, A, B)
        AtmulBpos_avx!(C3, A, B)
        @test C1 ≈ C2
        @test C1 ≈ C3
        C1 = similar(B);
        C2 = similar(B);
        notacondload!(C1, B, b)
        notacondloadavx!(C2, B, b)
        @test C1 ≈ C2
        maskedloadscalar!(C1, B, b)
        maskedloadscalaravx!(C2, B, b)
        @test C1 ≈ C2
        condloadscalar!(C1, B, b, d)
        condloadscalaravx!(C2, B, b, d)
        @test C1 ≈ C2

        if T <: Integer
            half = rand(T(1):T(100), 7);
            m = rand(T(-10):T(10), 7);
        else
            half = rand(T, 7); m = rand(T, 7);
        end;
        if sizeof(T) == 4
            res1 = Vector{Float32}(undef, 7);
            res2 = Vector{Float32}(undef, 7);
        else
            res1 = Vector{Float64}(undef, 7);
            res2 = Vector{Float64}(undef, 7);
        end

        for keep ∈ (nothing,true), final ∈ (nothing,true)
            @test twoifelses!(res1, half, m) ≈ twoifelses_avx!(res2, half, m)
        end

    end

    a = rand(-10:10, 43);
    bit = a .> 0.5;
    t = Bernoulli_logit(bit, a);
    @test isapprox(t, Bernoulli_logitavx(bit, a), atol = Int === Int32 ? 0.1 : 0)
    @test isapprox(t, Bernoulli_logit_avx(bit, a), atol = Int === Int32 ? 0.1 : 0)
    a = rand(43);
    bit = a .> 0.5;
    t = Bernoulli_logit(bit, a);
    @test t ≈ Bernoulli_logitavx(bit, a)
    @test t ≈ Bernoulli_logit_avx(bit, a)

    ai = [rand(Bool) for _ in 1:71];
    bi = [rand(Bool) for _ in 1:71];
    # if LoopVectorization.VectorizationBase.AVX2 || Base.libllvm_version ≥ v"8" #FIXME Why doesn't this work on Travis Ivy Bridge Julia 1.1?
        @test (ai .& bi) == (@avx ai .& bi)
        @test (ai .| bi) == (@avx ai .| bi)
        @test (ai .⊻ bi) == (@avx ai .⊻ bi)
    # else
    #     @test_broken (ai .& bi) == (@avx ai .& bi)
    #     @test_broken (ai .| bi) == (@avx ai .| bi)
    #     @test_broken (ai .⊻ bi) == (@avx ai .⊻ bi)
    # end
    a = bitrand(127); b = bitrand(127);
    @test (a .& b) == (@avx a .& b)
    @test (a .| b) == (@avx a .| b)
    @test (a .⊻ b) == (@avx a .⊻ b)
end


text:
# Autogenerated wrapper script for H3_jll for x86_64-apple-darwin
export libh3

JLLWrappers.@generate_wrapper_header("H3")
JLLWrappers.@declare_library_product(libh3, "@rpath/libh3.1.dylib")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libh3,
        "lib/libh3.1.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# trying to replicate Bai and Ng 2002

using DynamicFactorModels

tuples1 = [(100, 100, 1), (100, 40, 1), (100, 60, 1), (200, 60, 1), (500, 60, 1), (1000, 60, 1), (2000, 60, 1), (4000, 60, 1), (4000, 100, 1), (8000, 60, 1), (8000, 100, 1), (50, 10, 1), (100, 10, 1), (100, 20, 1)]  # T, N, r combinations I can replicate
tuples2 = [(100, 100, 3), (100, 40, 3), (100, 60, 3), (200, 60, 3), (500, 60, 3), (1000, 60, 3), (2000, 60, 3), (4000, 60, 3), (4000, 100, 3), (8000, 60, 3), (8000, 100, 3), (50, 10, 3), (100, 10, 3), (100, 20, 3)]  # T, N, r combinations I can replicate
tuples3 = [(100, 100, 5), (100, 40, 5), (100, 60, 5), (200, 60, 5), (500, 60, 5), (1000, 60, 5), (2000, 60, 5), (4000, 60, 5), (4000, 100, 5), (8000, 60, 5), (8000, 100, 5), (50, 10, 5), (100, 10, 5), (100, 20, 5)]  # T, N, r combinations I can replicate
reduced_tuples1 = [(100, 40, 1), (100, 60, 1), (200, 60, 1), (500, 60, 1), (50, 10, 1), (100, 10, 1), (100, 20, 1)]  # T, N, r combinations I can replicate fast
reduced_tuples2 = [(100, 40, 3), (100, 60, 3), (200, 60, 3), (500, 60, 3), (50, 10, 3), (100, 10, 3), (100, 20, 3)]  # T, N, r combinations I can replicate fast
reduced_tuples3 = [(100, 40, 5), (100, 60, 5), (200, 60, 5), (500, 60, 5), (50, 10, 5), (100, 10, 5), (100, 20, 5)]  # T, N, r combinations I can replicate fast

function replicate_table(sample_tuples, num_reps=1)  # replicate main information criteria of tables 1,2,3
    criteria = ["PCp1", "PCp2", "PCp3", "ICp1", "ICp2", "ICp3"]
    results = zeros(length(sample_tuples), length(criteria))
    for criterion_idx in 1:length(criteria)
        criterion = criteria[criterion_idx]
        println("crunching criterion $criterion")
        for i in 1:length(sample_tuples)
            results_per_rep = zeros(num_reps)
            for rep in 1:num_reps
                println("tuple: ", sample_tuples[i], "\t repetition: ", rep)
                y, x, f, lambda, epsilon_x, epsilon_y = apply(factor_model_DGP, sample_tuples[i])
                x = normalize(x)

                w = reshape(ones(length(y)), (length(y), 1))  # reshape to Array{Float64, 2}
                #dfm = try
                dfm = DynamicFactorModel(y, w, x, criterion)
                #catch exception
                #    warn("caught error. Probably because of negative variances.")
                #    println(exception)
                #end
                results_per_rep[rep] = dfm.number_of_factors
            end
            results[i, criterion_idx] = mean(results_per_rep)
            println("result for tuple ", sample_tuples[i], " is: ", results[i, criterion_idx])
        end
    end
    results
end

@time results1 = replicate_table(reduced_tuples1); gc()
@time results2 = replicate_table(reduced_tuples2); gc()
@time results3 = replicate_table(reduced_tuples3); gc()
results = vcat(results1, results2, results3)


text:
using Tweaks.BaseTweaks

@testset "argmin" begin
    minelem = argmin(-5:5) do i
        i^2
    end
    @test minelem == 0
end

@testset "argmax" begin
    maxelem = argmax(-5:5) do i
        -i^2
    end
    @test maxelem == 0
end

@testset "totuple" begin
    @test totuple((1, 2, 3, 4, 5)) == (1, 2, 3, 4, 5)
    @test totuple(1:5) == (1, 2, 3, 4, 5)
    @test totuple(i for i in 1:5) == (1, 2, 3, 4, 5)
    @test totuple(collect(1:5)) == (1, 2, 3, 4, 5)
end

@testset "@tuple" begin
    @test @tuple((1, 2, 3, 4, 5)) == (1, 2, 3, 4, 5)
    @test @tuple(1:5) == (1, 2, 3, 4, 5)
    @test @tuple(i for i in 1:5) == (1, 2, 3, 4, 5)
    @test @tuple(collect(1:5)) == (1, 2, 3, 4, 5)
    n = 5
    @test @tuple(i for i in 1:n) == (1, 2, 3, 4, 5)
end


text:
module DriveShare
export get_dataset

using Requests: get, save
using DataFrames
import JSON

function get_dataset(version)

    dirprefix = Pkg.dir("DriveShare", "data")
    if !isdir(dirprefix)
        mkdir(dirprefix)
    end

    if version == :trace
        url = "http://driveshare.me/anonUTM"
        filepath = joinpath(dirprefix, "anonUTM.json")
        success(`wget --continue --timestamping --output-document=$filepath $url`)
        sleep(0.1)
        return JSON.parsefile(filepath, dicttype=Dict, use_mmap=true)
    elseif version == :highlevel
        url = "http://driveshare.me/highlevel"
        filepath = joinpath(dirprefix, "highlevel.csv")
        success(`wget --continue --timestamping --output-document=$filepath $url`)
        sleep(0.1)
        return readtable(filepath)
    end

	# if version == :trace
	# 	save(get("http://driveshare.me/anonUTM"), "anonUTM.json")
	# 	return JSON.parsefile("anonUTM.json", dicttype=Dict, use_mmap=true)
	# elseif version == :highlevel
	# 	save(get("http://driveshare.me/highlevel"), "highlevel.csv")
	# 	return readtable("highlevel.csv")
	# end
end

end


text:
module Polarisation

using StaticArrays
using LinearAlgebra

"Make horizontally polarised Jones vector"
function H()
   SVector(1.0, 0.0)
end

"Make horizontally orientated linear polariser Jones matrix"
function LP()
    @SMatrix [ 1.0  0.0 ;
               0.0  0.0 ]
end

"Arbitrary waveplate with phase ϕ; fast axis horizontal"
function WP(ϕ)
    @SMatrix [ exp(im*ϕ/2)  0.0 ;
               0.0          exp(-im*ϕ/2) ]
end

"rotation operator"
function rot(θ)
    @SMatrix [ cos(θ)  sin(θ) ;
              -sin(θ)  cos(θ) ]
end

"rotate Jonesmatrix J by θ"
function rotate(J, θ)
    rot(-θ)*J*rot(θ)
end

"Get Stokes parameters for input field E = (Ex, Ey)"
function Stokes(E; normalise=false)
    Ex = E[1]
    Ey = E[2]
    I = abs(Ex)^2 + abs(Ey)^2
    Q = abs(Ex)^2 - abs(Ey)^2
    U = 2*real(Ex*conj(Ey))
    V = -2*imag(Ex*conj(Ey))
    S = SVector(I, Q, U, V)
    if normalise
        S = S./S[1]
    end
    S
end

"Get normalised cartesian coordinates from Stokes parameters
 (for Poincare sphere) "
function cartesian(S)
    S[2:end]./S[1]
end

"Get polarization ellipse parameters from Stokes parameters"
function ellipse(S)
    I = S[1]
    Q = S[2]
    U = S[3]
    V = S[4]
    aL = sqrt(Q^2 + U^2)
    θ = angle(aL)/2
    A = sqrt((I + aL)/2)
    B = sqrt((abs(I - aL))/2)
    h = sign(V)
    A, B, θ, h
end

"Calculate ellipticity from Stokes parameters"
function ellipticity(S)
    A, B, θ, h = ellipse(S)
    r = A/B
    r > 1 ? 1/r : r
end

end


text:
module SoftRank

import ChainRulesCore: rrule, NoTangent

export rank, softrank

"""
    partition(x; ϵ=1e-9)

Compute the sizes of each partition of constant value obtained by isotonic regression of data `x`.
`ϵ` denotes the tolerance for what is considered equal in floating point terms.
"""
function partition(x; ϵ=1e-9)
	if length(x) == 0
		return eltype(x)[]
	end

	sizes = Int[1]
	for i in 2:length(x)
		if abs(x[i] - x[i-1]) > ϵ
			append!(sizes, 0)
		end
		sizes[end] += 1
	end

	return sizes
end

"""
    isotonic(x)

Isotonically regress on data vector `x`.
Returns the monotonically increasing fit.
"""
function isotonic(x)
	n = length(x)
	τ = collect(vec(1:n))
	c = ones(Float32,n)
	Σ = copy(x)
	y = copy(x)

	i = 1
	while i <= n
		k = τ[i] + 1
		if k == (n + 1)
			break
		end

		if y[i] > y[k]
			i = k
			continue
		end

		Σʸ = Σ[i]
		Σᶜ = c[i]
		# We are within an increasing subsequence
		while true
			y₀  = y[k]
			Σʸ += Σ[k]
			Σᶜ += c[k]

			k = τ[k] + 1
			if k == (n + 1) || y₀ > y[k]
				y[i] = Σʸ / Σᶜ
				Σ[i] = Σʸ
				c[i] = Σᶜ

				τ[i]   = k-1
				τ[k-1] = i
				if i > 1
					i = τ[i-1]
				end

				break
			end
		end
	end

	i = 1
	while i <= n
		k = τ[i]
		y[(i+1):k] .= y[i]
		i = k + 1
	end

	return y

end

"""
    ∇isotonic(x)

Compute the gradient of isotonic regression with respect to the inputs about the solution point.
"""
function ∇isotonic(soln, x⃗)
	y⃗ = zeros(Float32, size(x⃗))

	i₁ = 1
	for δ in partition(soln)
		i₂ = i₁ + δ
		y⃗[i₁:(i₂-1)] .= sum(x⃗[i₁:(i₂-1)])/δ
		i₁ = i₂
	end

	return y⃗
end

"""
    projection(x)

Regularize the ranking algorithm by euclidean projection onto the permutehedron.
"""
function projection(x)
	n = length(x)
    w = Float32.(reverse(range(1/n, 1; length=n))) #vec(collect(n:-1:1))
	ι = reverse(sortperm(x))
	s = x[ι]

	dual   = isotonic(s.-w)
	primal = s .- dual

	return primal[invperm(ι)]
end

"""
    rank(x)

Rank the items of vector `x` in ascending order.
"""
rank(x)             = invperm(sortperm(x))
"""
    softrank(x)

Regularization of the ranking algorithm.
Scaled euclidean projection onto the permutehedron.
"""
softrank(x; ϵ=Float32(1e-2)) = projection(x ./ ϵ)

function rrule(::typeof(softrank), x; ϵ=Float32(1e-2))
    # scale input
    x = x ./ ϵ

    # compute value
    n = length(x)
    w = Float32.(reverse(range(1/n, 1; length=n)))
    ι = reverse(sortperm(x))
    s = x[ι]

    dual   = isotonic(s.-w)
    primal = s .- dual

    ι¯¹ = invperm(ι)
    return primal[ι¯¹], (∇) -> (NoTangent(), (∇ .- ∇isotonic(dual, ∇[ι])[ι¯¹]) ./ ϵ)
end

# FIXME: topological flips are wrong!!

function hsplit!(r, i, δx, δy, ι, n)
    x₀ = (δx.hi + δx.lo) / 2

    lo = findall(r[1,:] .≤ x₀)
    hi = findall(r[1,:] .> x₀)

    @show lo, hi

    # left partition
    if length(lo) > 1
        n = vsplit!(r[:,lo], i[lo], (lo=δx.lo,hi=x₀), δy, ι, n)
    elseif length(lo) == 1
        ι[i[lo[1]]] = n
        n += 1
    end

    # right partition
    if length(hi) > 1
        n = vsplit!(r[:,hi], i[hi], (lo=x₀,hi=δx.hi), δy, ι, n)
    elseif length(hi) == 1
        ι[i[hi[1]]] = n
        n += 1
    end

    return n
end

function vsplit!(r, i, δx, δy, ι, n)
    y₀ = (δy.hi + δy.lo) / 2

    lo = findall(r[2,:] .≤ y₀)
    hi = findall(r[2,:] .> y₀)

    @show lo, hi

    # left partition
    if length(lo) > 1
        n = hsplit!(r[:,lo], i[lo], δx, (lo=δy.lo,hi=y₀), ι, n)
    elseif length(lo) == 1
        ι[i[lo[1]]] = n
        n += 1
    end

    # right partition
    if length(hi) > 1
        n = hsplit!(r[:,hi], i[hi], δx, (lo=y₀,hi=δy.hi), ι, n)
    elseif length(hi) == 1
        ι[i[hi[1]]] = n
        n += 1
    end

    return n
end

# NOTE: assumes 2D
#       assumes coordinates fall between -1 and +1
function hilbert(r)
    ι = zeros(Int, size(r,2))
    i = collect(1:size(r,2))
    n = hsplit!(r, i, (lo=-1.0,hi=+1.0), (lo=-1.0,hi=+1.0), ι, 1)

    return ι
end

end

max_stars_repo_name: nnoll/seqspace

max_issues_repo_name: nnoll/seqspace

max_forks_repo_name: nnoll/seqspace


text:
using DataFrames, ROOTDataFrames

N = 10000

#pre-allocate input data
bigdf = similar(

    #here we define TTree structure
    DataFrame(
        evt=Int64[],
        njet=Int64[],
        jet_pt=Vector{Float64}[],
        jet_eta=Vector{Float64}[],
        jet_id=Vector{Int64}[],

        nlep=Int64[],
        lep_pt=Vector{Float64}[],
        lep_eta=Vector{Float64}[],
        lep_id=Vector{Int64}[]

    ),
    N
)

#fill event by event
for iev=1:N
    nj = round(Int, rand()*10)
    bigdf[iev, :evt] = iev
    bigdf[iev, :njet] = nj

    #create jets with random vectors
    bigdf[iev, :jet_pt] = [rand()*200 for i=1:nj]
    bigdf[iev, :jet_eta] = [randn()*2 for i=1:nj]
    bigdf[iev, :jet_id] = [i for i=1:nj]

    nl = round(Int, rand()*2)
    bigdf[iev, :nlep] = nl
    bigdf[iev, :lep_pt] = [rand()*200 for i=1:nl]
    bigdf[iev, :lep_eta] = [randn()*2 for i=1:nl]
    bigdf[iev, :lep_id] = [i for i=1:nl]

end

#write output tree
dtf = TreeDataFrame(bigdf, "dyntree.root";
    #need to specify which branches correspond to dynamic branch lengths
    size_branches=Dict(
        :jet_pt=>:njet, :jet_eta=>:njet, :jet_id=>:njet,
        :lep_pt=>:nlep, :lep_eta=>:nlep, :lep_id=>:nlep,
    )
)


text:
using Dates
using SHA
using DataFrames
using CSV
using JSON

hashlog(ID) = (read(`docker logs $(ID)`) |> sha256 |> bytes2hex)
containerDestroyByID(ID) = read(`docker rm -f $(ID)`)

function servicePort(str)
	a = read(`docker service inspect $(str)`, String) |> JSON.parse
	return a[1]["Spec"]["EndpointSpec"]["Ports"][1]["PublishedPort"]
end

function stats()
	a = read(`docker stats --no-stream`, String)

	b = split(a, '\n')
	c = map(x -> split(x), b)
	containerID = map(x -> String(x[1]), c[2:end-1])
	containerStr = map(x -> String(x[2] |> x -> split(x, '.')[1]), c[2:end-1])
	NetIO = map(x -> x[10], c[2:end-1])

	function mparse(x)
		N = try
			parse(Float64, x)
		catch
			nothing
		end
	end

	function parseIO(x)
		testunit(x, unit) = split(x, unit)[1] |> mparse

		testunit(x, "MB")
		if ~isnothing(testunit(x,"MB"))
			return (testunit(x, "MB"), "MB")
		elseif ~isnothing(testunit(x,"kB"))
			return (testunit(x, "kB"), "kB")
		elseif ~isnothing(testunit(x,"B"))
			return (testunit(x, "B"), "B")
		else
			return (0, "unknown")
		end
	end

	function netIOtoBytes(val)
		if (val[2] == "B")
			return val[1]
		elseif (val[2] == "kB")
			return val[1] * 1024
		elseif (val[2] == "MB")
			return val[1] * 1024 * 1024
		elseif (val[2] == "GiB")
			return val[1] * 1024 * 1024 * 1024
		else
			return 0
		end
	end

	parsedNetIO = map(parseIO, NetIO)
	cleanNetIO = map(netIOtoBytes, parsedNetIO)
	return containerID, containerStr, cleanNetIO
end

function containerAgeByID(ID)
	c = read(`docker ps -f "ID=$(ID)" --format "{{.CreatedAt}}"`, String)
	if c == ""
		return nothing
	end
	d = split(c) |> x -> String.(x)
	t = Date(d[1], dateformat"YYYY-mm-dd") + Time(d[2], dateformat"HH:MM:SS")
	age = round(now() - t, Minute)
end

function updateActivity(containerActivity, ID, netIO)
	try
		netIO - containerActivity.value[ID]
	catch
		netIO
	end
end

function regenerate_service(containerStr::String, port::Int, delay::Int)
	sleep(delay)
	z = String.(split(containerStr, '-'))
    thisPort = servicePort(containerStr)
	if thisPort ∉ resolve_ports[z[1]]
		println("Port was checked out, regenerate")
		push!(resolve_ports[z[1]], port)
	else
		println("Container regenerated, but container was not checked out")
	end
end

function containerManager()
	containerID, containerStr, netIO = stats()
	age = map(containerAgeByID, containerID)
	activity = map((x,y) -> updateActivity(containerActivity,x,y), containerID, netIO)
	map((x,y) -> containerActivity.value[x] = y, containerID, netIO)

	taintedFlag = map(containerID, containerStr) do x,y
		z = String.(split(y, '-'))
		# hash = open(f->read(f, String), z[1])
        thisPort = servicePort(y)
		println(y, "  ", thisPort, "  ", resolve_ports[z[1]])
		# (hashlog(x) != hash) || (thisPort ∉ resolve_ports[z[1]])
		thisPort ∉ resolve_ports[z[1]]
	end

	log = DataFrame(
		t = now(),
		containerID = containerID,
		containerStr = containerStr,
		age = age,
		netIO = netIO,
		activity = activity,
		taintedFlag = taintedFlag
	)

	function destroyContainer(containerID, containerStr, age, taintedFlag, activity)
		if taintedFlag && (activity < 1.0)
			z = String.(split(containerStr, '-'))
			port = port_base[z[1]] + parse(Int, z[2])
			println("Destroy $(containerID)")
			println("Restore Service port: $(port)")
			containerDestroyByID(containerID)
			@async regenerate_service(containerStr, port, 300) # restore to pool w/300s delay
		end

		return nothing
	end

	map(destroyContainer, containerID, containerStr, age, taintedFlag, activity)
	log |> CSV.write("logfile.csv", append = true)
	show(log)
	println()

	return log
end


text:
# Copyright (c) 2020: Matthew Wilhelm & Matthew Stuber.
# This work is licensed under the Creative Commons Attribution-NonCommercial-
# ShareAlike 4.0 International License. To view a copy of this license, visit
# http://creativecommons.org/licenses/by-nc-sa/4.0/ or send a letter to Creative
# Commons, PO Box 1866, Mountain View, CA 94042, USA.
#############################################################################
# DynamicBounds.jl
# A package for compute bounds and relaxations of the solutions of
# parametric differential equations.
# See https://github.com/PSORLab/DynamicBoundsBase.jl
#############################################################################
# src/reexport_using.jl
# A copy of the reexport module that contains the `@reexport using foo: bar`
# syntax functionality contained in Reexport.jl master branch. Since no such
# tagged release is currently available.

#=
Copyright (c) 2014: Simon Kornblith.
Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
THE SOFTWARE.
=#

#############################################################################

macro reexport(ex)
    isa(ex, Expr) && (ex.head == :module ||
                      ex.head == :using ||
                      (ex.head == :toplevel &&
                       all(e->isa(e, Expr) && e.head == :using, ex.args))) ||
        error("@reexport: syntax error")

    if ex.head == :module
        modules = Any[ex.args[2]]
        ex = Expr(:toplevel, ex, :(using .$(ex.args[2])))
    elseif ex.head == :using && all(e->isa(e, Symbol), ex.args)
        modules = Any[ex.args[end]]
    elseif ex.head == :using && ex.args[1].head == :(:)
        symbols = [e.args[end] for e in ex.args[1].args[2:end]]
        return esc(Expr(:toplevel, ex, :(eval(Expr(:export, $symbols...)))))
    else
        modules = Any[e.args[end] for e in ex.args]
    end

    esc(Expr(:toplevel, ex,
             [:(eval(Expr(:export, names($mod)...))) for mod in modules]...))
end


text:
@testset "Mat.jl" begin
    begin
        mat = PetscMat()

        setSize!(mat, m_local=(Int32)(3), n_local=(Int32)(4))
        @test mat.sized

        setPreallocation!(mat, (Int32)[3, 3, 2], (Int32)[0,0,0])
        @test mat.preallocated

        # Test individual assignment
        mat[1,2] = 1.2

        assemble!(mat)

        @test size(mat) == (3,4)

        @test mat[1,2] == 1.2
    end

    begin
        mat = PetscMat()

        setSize!(mat, m_local=(Int32)(3), n_local=(Int32)(4))
        @test mat.sized

        setPreallocation!(mat, (Int32)[3, 3, 2], (Int32)[0,0,0])
        @test mat.preallocated

        # Test range assignment
        mat[1:2,2:3] = (Float64)[1 2; 3 4]

        assemble!(mat)

        @test size(mat) == (3,4)
        @test mat[1:2,2:3] == (Float64)[1 2; 3 4]
    end

    begin
        mat = PetscMat()

        setSize!(mat, m_local=(Int32)(3), n_local=(Int32)(4))
        @test mat.sized

        setPreallocation!(mat, (Int32)[3, 3, 2], (Int32)[0,0,0])
        @test mat.preallocated

        mat[[1,3],[2,3]] = (Float64)[1 2; 3 4]
        assemble!(mat)

        plusEquals!(mat, (Float64)[2 3; 4 5], (Int32)[1,3], (Int32)[2,3])
        assemble!(mat)

        @test mat[1:3,2:3] == (Float64)[3 5; 0 0; 7 9]
    end
end


text:
using NetLogo

NetLogo.init(joinpath(homedir(), "bin/netlogo/6.0.4/app"))
NetLogo.withmodel("models/Sample Models/Biology/Moths.nlogo") do w
  NetLogo.command(w, "setup")
  light = NetLogo.report(w, "one-of lights")
  @show light
  @show light["intensity"]
  Base.length(light)
  println(collect(light))
end


text:
using DataFrames
using SparseArrays
using Statistics

struct Dataset{T <: Number}
    ratings::SparseMatrixCSC{AbstractRating{T}, Int}
    preference::Preference{T}
    users::Int
    items::Int
end

function Dataset(df::DataFrame, users::Int, items::Int)
    @assert in(:user, names(df))
    @assert in(:item, names(df))
    @assert in(:rating, names(df))

    if users < maximum(df[:user]) || items < maximum(df[:item])
        throw(ArgumentError("users or items must satisfy maximum[df[:k]] >= k"))
    end

    preference = Preference(df[:rating])

    ratings = convert(df[:rating], preference)

    matriz = sparse(df[:user], df[:item], ratings, users, items)

    return Dataset(matriz, preference, users, items)
end

Dataset(df::DataFrame) = Dataset(df, maximum(df[:user]), maximum(df[:item]))

users(dataset::Dataset) = dataset.users
items(dataset::Dataset) = dataset.items

Base.size(dataset::Dataset) = (users(dataset), items(dataset))
Base.length(dataset::Dataset) = length(nonzeros(dataset.ratings))

function Base.Array(dataset::Dataset{T})::Matrix{Union{Missing, T}} where T <: Number
    matrix = Array{Union{Missing, T}}(missing, Persa.users(dataset), Persa.items(dataset))

    users = rowvals(dataset.ratings)
    ratings = nonzeros(dataset.ratings)
    for item = 1:items(dataset)
       for j in nzrange(dataset.ratings, item)
          user = users[j]
          rating = ratings[j]
          matrix[user, item] = value(rating)
       end
    end

    return matrix
end

checkuser(dataset::Dataset, user::Int) = (user >= 0 && user <= users(dataset)) ? user : error("invalid user id ($user)")
checkitem(dataset::Dataset, item::Int) = (item >= 0 && item <= items(dataset)) ? item : error("invalid item id ($item)")

function Base.getindex(dataset::Dataset, user::Int, item::Int)::AbstractRating
    checkuser(dataset, user)
    checkitem(dataset, item)

    dataset.ratings[user, item]
end

function Base.getindex(dataset::Dataset, i::Int)::Tuple{Int, Int, Rating}
    if i > length(dataset)
        throw(ArgumentError("index must satisfy 1 <= i <= length(dataset)"))
    end

    users = rowvals(dataset.ratings)
    ratings = nonzeros(dataset.ratings)

    for item = 1:items(dataset)
        if i in nzrange(dataset.ratings, item)
            user = users[i]
            rating = ratings[i]

            return (user, item, rating)
        end
    end
end

function Base.getindex(dataset::Dataset, user::Int, c::Colon)
    elements = collect(dataset.ratings[user, :])
    return [(idx, elements[idx]) for idx in findall(!isnan, elements)]
end

function Base.getindex(dataset::Dataset, c::Colon, item::Int)
    elements = collect(dataset.ratings[:, item])
    return [(idx, elements[idx]) for idx in findall(!isnan, elements)]
end

Base.iterate(dataset::Dataset, state = 1) = state > length(dataset) ? nothing : (dataset[state], state+1)

function Statistics.mean(dataset::Dataset)
    μ = 0
    total = 0
    for (u, v, r) in dataset
        if !isnan(r)
            μ += r
            total += 1
        end
    end
    return μ ./ total
end


text: # Functions and methods for cascade computation


text:
import Base: +, -, ==

immutable BoundingBox{T}
  i::T
  j::T
  h::T    # height
  w::T    # width
end

eps = 1e-6

BoundingBox() = BoundingBox(0,0,0,0)
BoundingBox(a, b, c, d) = BoundingBox(promote(a,b,c,d)...)

"""
Add bounding boxes to find BB of their union
"""
function +(bbA::BoundingBox, bbB::BoundingBox)
  i = min(bbA.i, bbB.i)
  j = min(bbA.j, bbB.j)
  h = max(bbA.h+bbA.i, bbB.h+bbB.i) - i
  w = max(bbA.w+bbA.j, bbB.w+bbB.j) - j
  return BoundingBox(i,j,h,w)
end

"""
Subtract bounding boxes to find BB of their intersection
"""
function -(bbA::BoundingBox, bbB::BoundingBox)
  i = max(bbA.i, bbB.i)
  j = max(bbA.j, bbB.j)
  h = min(bbA.h+bbA.i, bbB.h+bbB.i)-i
  w = min(bbA.w+bbA.j, bbB.w+bbB.j)-j
  if h <= 0 || w <= 0
    bb = BoundingBox(NaN, NaN, NaN, NaN)
  else
    bb = BoundingBox(i,j,h,w)
  end
  return bb
end

"""
Boolean if bounding boxes intersect
"""
function intersects(bbA::BoundingBox, bbB::BoundingBox)
  bb = bbA - bbB
  return !isnan(bb.i)
end

"""
Test if BoundingBoxes have the same origin and dimensions
"""
function ==(bbA::BoundingBox, bbB::BoundingBox)
  return bbA.i == bbB.i && bbA.j == bbB.j && bbA.w == bbB.w && bbA.h == bbB.h
end

"""
Convert bounding box object to a polygon point list (counter-clockwise)
"""
function bb_to_pts(r)
  return [r.i r.j;
          r.i+r.h r.j;
          r.i+r.h r.j+r.w;
          r.i r.j+r.w;
          r.i r.j];
end

"""
Return xmin, ymin, xmax, ymax of a bounding box (opposing corner definition)
"""
function minsandmax(bb)
  return bb.i, bb.j, bb.i+bb.h, bb.j+bb.w
end

"""
`FIND_MESH_BB` - Find bounding box around mesh

    BoundingBox(ilow, jlow, height, width) = find_mesh_bb(nodes)

* `nodes`: 2xN matrix of mesh nodes
* `BoundingBox`: smallest integer-valued rectangle containing all mesh nodes

"""
function find_mesh_bb(nodes)
    ilow = floor(Int64,minimum(nodes[:,1]))
    jlow = floor(Int64,minimum(nodes[:,2]))
    ihigh = ceil(Int64,maximum(nodes[:,1]))
    jhigh = ceil(Int64,maximum(nodes[:,2]))
    return BoundingBox(ilow, jlow, ihigh-ilow, jhigh-jlow)
end

"""
Snap bounding box to integer values (smallest rectangle containing original)
Returns BoundingBox of integers
"""
function snap_bb(bb)
  r = bb_to_pts(bb)
  i = floor(Int,bb.i)
  j = floor(Int,bb.j)
  h = ceil(Int,r[2,1]) - i
  w = ceil(Int,r[3,2]) - j
  return BoundingBox{Int}(i, j, h, w)
end

"""
Apply affine transform to points in a bounding box & find new bounding box
"""
function tform_bb(bb, tform)
  tform_pts = [bb_to_pts(bb) ones(size(bb_to_pts(bb),1),1)] * tform
  i = minimum(tform_pts[:,1])
  j = minimum(tform_pts[:,2])
  h = maximum(tform_pts[:,1])-i
  w = maximum(tform_pts[:,2])-j
  return BoundingBox(i, j, h, w)
end

"""
Convert Tuple for image size into a BoundingBox at (1,1)
"""
function sz_to_bb(sz)
  return BoundingBox(0, 0, sz[1], sz[2])
end

"""
Convert a BoundingBox to its snapped sizes
"""
function bb_to_sz(bb)
  bb = snap_bb(bb);
  return bb.h, bb.w
end

"""
Get bounding box offset
"""
function get_offset(bb::BoundingBox)
  return [bb.i, bb.j]
end

"""
Get height & width of bounding box
"""
function get_size(bb::BoundingBox)
  return [bb.h, bb.w]
end

"""
Get four-tuple of bounding box upper-left & lower-right in i,j coordinates
"""
function get_bounds(bb::BoundingBox)
  return (bb.i, bb.j, bb.i+bb.h, bb.j+bb.w)
end

"""
Get four-tuple of bounding box upper-left & dimensions in x,y coordinates

Convention used by Cairo package
"""
function get_rect(bb::BoundingBox)
  return (bb.j, bb.i, bb.w, bb.h)
end

"""
Return the product of the width & height of bounding box
"""
function get_area(bb::BoundingBox)
  return bb.w*bb.h
end

"""
Convert bounding box to tuple of ranges for easy array slicing
"""
function bb_to_slice(bb::BoundingBox{Int64})
  return (bb.i):(bb.i+bb.h-1), (bb.j):(bb.j+bb.w-1)
end

function bb_to_slice(bb::BoundingBox{Float64})
  return round(Int64, bb.i):round(Int64, bb.i+bb.h-1),
              round(Int64, bb.j):round(Int64, bb.j+bb.w-1)
end

"""
Convert tuple of ranges to bounding box
"""
function slice_to_bb(slice)
  return BoundingBox(slice[1][1], slice[2][1],
                      slice[1][end]-slice[1][1]+1, slice[2][end]-slice[2][1]+1)
end

"""
Shift bounding box by 2-element array
"""
function translate_bb(bb::BoundingBox, offset)
  return BoundingBox(bb.i + offset[1], bb.j + offset[2], bb.h, bb.w)
end

"""
Transform bounding box by scaling matrix and snap for nearest integer
"""
function scale_bb(bb::BoundingBox{Int64}, scale)
  tform = make_scale_matrix(scale)
  tbb = tform_bb(bb, tform)
  return snap_bb(tbb)
end

"""
Transform bounding box by scaling matrix
"""
function scale_bb(bb::BoundingBox{Float64}, scale)
  tform = make_scale_matrix(scale)
  return tform_bb(bb, tform)
end


text:
# test forest of graphs can solve with CSM, specifically #518

using IncrementalInference
using Statistics
using Test
using TensorCast

@testset "Test forest of orphaned graphs" begin

fg = initfg()
addVariable!(fg, :x0, ContinuousScalar)
addFactor!(fg, [:x0;], Prior(Normal(0,0.1)))
addVariable!(fg, :x1, ContinuousScalar)
addFactor!(fg, [:x0;:x1], LinearRelative(Normal(10,0.1)))
addVariable!(fg, :x2, ContinuousScalar)
addFactor!(fg, [:x1;:x2], LinearRelative(Normal(10,0.1)))

addVariable!(fg, :x10, ContinuousScalar)
addFactor!(fg, [:x10;], Prior(Normal()))
addVariable!(fg, :x11, ContinuousScalar)
addFactor!(fg, [:x10;:x11], LinearRelative(Normal(-10,1.0)))
addVariable!(fg, :x12, ContinuousScalar)
addFactor!(fg, [:x11;:x12], LinearRelative(Normal(-10,1.0)))

# plotDFG(fg)
# getSolverParams(fg).drawtree = true
# getSolverParams(fg).showtree = true
# solve factor graph with two orphaned components
vo = Symbol[:x12, :x2, :x0, :x11, :x1, :x10]
tree = solveTree!(fg, eliminationOrder=vo)

# test tree will have two different root nodes
@test getEliminationOrder(tree) == vo

@test getParent(tree, getClique(tree, :x1)) |> length == 0
@test getParent(tree, getClique(tree, :x10)) |> length == 0

@test getChildren(tree, getClique(tree, :x1)) |> length == 1
@test getChildren(tree, getClique(tree, :x10)) |> length == 1

@test getChildren(tree, getClique(tree, :x2)) |> length == 0
@test getChildren(tree, getClique(tree, :x12)) |> length == 0

## Test the numerical values are correct

pts_ = getBelief(fg, :x0) |> getPoints
@cast pts[i,j] := pts_[j][i]
@test  pts  |> mean |> abs < 1.0
pts_ = getBelief(fg, :x1) |> getPoints
@cast pts[i,j] := pts_[j][i]
@test (pts  |> mean) - 10 |> abs < 2.0
pts_ = getBelief(fg, :x2) |> getPoints
@cast pts[i,j] := pts_[j][i]
@test (pts  |> mean) - 20 |> abs < 3.0

pts_ = getBelief(fg, :x10) |> getPoints
@cast pts[i,j] := pts_[j][i]
@test  pts  |> mean |> abs < 2.0
pts_ = getBelief(fg, :x11) |> getPoints
@cast pts[i,j] := pts_[j][i]
@test (pts  |> mean) + 10 |> abs < 4.0
pts_ = getBelief(fg, :x12) |> getPoints
@cast pts[i,j] := pts_[j][i]
@test (pts  |> mean) + 20 |> abs < 5.0

# using RoMEPlotting
# Gadfly.set_default_plot_size(35cm, 25cm)
#
# plotKDE(fg, ls(fg))

end


text:
# ------------------------------------------------------------------------------
# bddc
# ------------------------------------------------------------------------------

"""
```julia
BDDC(operator, coarseoperator, smoother, prolongation)
```

BDDC preconditioner for finite element operators

# Arguments:
- `operator`:       finite element operator to precondition
- `injectiontype`:  type of injection into subassembled space to use

# Returns:
- BDDC preconditioner object
"""
mutable struct BDDC <: AbstractPreconditioner
    # data never changed
    operator::Operator
    injectiontype::BDDCInjectionType.BDDCInjectType

    # data empty until assembled
    primalnodes::AbstractArray{Int,1}
    subassemblednodes::AbstractArray{Int,1}
    interfacenodes::AbstractArray{Int,1}
    interiornodes::AbstractArray{Int,1}
    primalmodes::AbstractArray{Int,1}
    subassembledmodes::AbstractArray{Int,1}
    interfacemodes::AbstractArray{Int,1}
    interiormodes::AbstractArray{Int,1}
    subassembledinverse::AbstractArray{Float64,2}
    interiorinverse::AbstractArray{Float64,2}
    schur::AbstractArray{Float64,2}
    mixedmultiplicity::AbstractArray{Float64}
    JDT::AbstractArray{Float64,2}
    primalrowmodemap::AbstractArray{Float64,2}
    primalcolumnmodemap::AbstractArray{Float64,2}
    subassembledrowmodemap::AbstractArray{Float64,2}
    subassembledcolumnmodemap::AbstractArray{Float64,2}
    mixedrowmodemap::AbstractArray{Float64,2}
    mixedcolumnmodemap::AbstractArray{Float64,2}
    modepermutation::AbstractArray{Bool,2}
    jacobi::Jacobi

    # inner constructor
    BDDC(operator::Operator, injectiontype::BDDCInjectionType.BDDCInjectType) = (
        # validity check
        if operator.dimension == 1
            error("BDDC preconditioner only valid for 1 and 2 dimensional operators") # COV_EXCL_LINE
        end;

        # constructor
        new(operator, injectiontype)
    )
end

# printing
# COV_EXCL_START
function Base.show(io::IO, preconditioner::BDDC)
    if preconditioner.injectiontype == BDDCInjectionType.scaled
        print(io, "lumped ")
    elseif preconditioner.injectiontype == BDDCInjectionType.harmonic
        print(io, "Dirichlet ")
    end
    print(io, "BDDC preconditioner")
end
# COV_EXCL_STOP

# ------------------------------------------------------------------------------
# data for computing symbols
# ------------------------------------------------------------------------------

"""
```julia
getprimalnodes(bddc)
```

Compute or retrieve the primal nodes for a BDDC preconditioner

# Returns:
- Vector of primal nodes for BDDC preconditioner

# Example:
```jldoctest
# setup
mesh = Mesh2D(1.0, 1.0);
p = 4
diffusion = GalleryOperator("diffusion", p, p, mesh);
bddc = LumpedBDDC(diffusion)

# note: either syntax works
primalnodes = LFAToolkit.getprimalnodes(bddc);
primalnodes = bddc.primalnodes;

# verify
@assert primalnodes == [1, p, p^2 - p + 1, p^2]

# output

```
"""
function getprimalnodes(bddc::BDDC)
    # assemble if needed
    if !isdefined(bddc, :primalnodes)
        operator = bddc.operator
        primalnodes = []
        numbernodes = 0
        for input in operator.inputs
            if input.evaluationmodes[1] != EvaluationMode.quadratureweights
                primalnodes =
                    vcat(primalnodes..., [input.basis.primalnodes .+ numbernodes]...)
                numbernodes += input.basis.numbernodes
            end
        end

        # store
        bddc.primalnodes = primalnodes
    end

    # return
    return getfield(bddc, :primalnodes)
end

"""
```julia
getsubassemblednodes(bddc)
```

Compute or retrieve the subassembled nodes for a BDDC preconditioner

# Arguments:
- `preconditioner`:  preconditioner to compute sumassembled nodes

# Returns:
- Vector of subassembled nodes for BDDC preconditioner

# Example:
```jldoctest
# setup
mesh = Mesh2D(1.0, 1.0);
p = 4
diffusion = GalleryOperator("diffusion", p, p, mesh);
bddc = LumpedBDDC(diffusion)

# note: either syntax works
subassemblednodes = LFAToolkit.getsubassemblednodes(bddc);
subassemblednodes = bddc.subassemblednodes;

# verify
@assert subassemblednodes == setdiff(1:p^2, [1, p, p^2-p+1, p^2])

# output

```
"""
function getsubassemblednodes(bddc::BDDC)
    # assemble if needed
    if !isdefined(bddc, :subassemblednodes)
        numbernodes, _ = size(bddc.operator.elementmatrix)

        # store
        bddc.subassemblednodes = setdiff(1:numbernodes, bddc.primalnodes)
    end

    # return
    return getfield(bddc, :subassemblednodes)
end

"""
```julia
getinterfacenodes(bddc)
```

Compute or retrieve the interface nodes for a BDDC preconditioner

# Arguments:
- `preconditioner`:  preconditioner to compute interface nodes

# Returns:
- Vector of interface nodes for BDDC preconditioner

# Example:
```jldoctest
# setup
mesh = Mesh2D(1.0, 1.0);
p = 4
diffusion = GalleryOperator("diffusion", p, p, mesh);
bddc = LumpedBDDC(diffusion)

# note: either syntax works
interfacenodes = LFAToolkit.getinterfacenodes(bddc);
interfacenodes = bddc.interfacenodes;

# verify
trueinterfacenodes = [1:p..., p^2-p+1:p^2...]
for i = 1:p-2
    push!(trueinterfacenodes, i*p+1)
    push!(trueinterfacenodes, (i+1)*p)
end
trueinterfacenodes = sort(trueinterfacenodes)
trueinterfacenodes = setdiff(trueinterfacenodes, bddc.primalnodes)
@assert interfacenodes == trueinterfacenodes

# output

```
"""
function getinterfacenodes(bddc::BDDC)
    # assemble if needed
    if !isdefined(bddc, :interfacenodes)
        operator = bddc.operator
        interfacenodes = []
        numbernodes = 0
        for input in operator.inputs
            if input.evaluationmodes[1] != EvaluationMode.quadratureweights
                interfacenodes =
                    vcat(interfacenodes..., [input.basis.interfacenodes .+ numbernodes]...)
                numbernodes += input.basis.numbernodes
            end
        end
        interfacenodes = setdiff(interfacenodes, bddc.primalnodes)

        # store
        bddc.interfacenodes = interfacenodes
    end

    # return
    return getfield(bddc, :interfacenodes)
end

"""
```julia
getinteriornodes(bddc)
```

Compute or retrieve the interior nodes for a BDDC preconditioner

# Arguments:
- `preconditioner`:  preconditioner to compute interior nodes

# Returns:
- Vector of interior nodes for BDDC preconditioner

# Example:
```jldoctest
# setup
mesh = Mesh2D(1.0, 1.0);
p = 4
diffusion = GalleryOperator("diffusion", p, p, mesh);
bddc = LumpedBDDC(diffusion)

# note: either syntax works
interiornodes = LFAToolkit.getinteriornodes(bddc);
interiornodes = bddc.interiornodes;

# verify
trueinterfacenodes = [1:p..., p^2-p+1:p^2...]
for i = 1:p-2
    push!(trueinterfacenodes, i*p+1)
    push!(trueinterfacenodes, (i+1)*p)
end
trueinterfacenodes = sort(trueinterfacenodes)
trueinteriornodes = setdiff(bddc.subassemblednodes, trueinterfacenodes)
@assert interiornodes == trueinteriornodes

# output

```
"""
function getinteriornodes(bddc::BDDC)
    # assemble if needed
    if !isdefined(bddc, :interiornodes)
        # store
        bddc.interiornodes = setdiff(bddc.subassemblednodes, bddc.interfacenodes)
    end

    # return
    return getfield(bddc, :interiornodes)
end

"""
```julia
getprimalmodes(bddc)
```

Compute or retrieve the primal modes for a BDDC preconditioner

# Arguments:
- `preconditioner`:  preconditioner to compute primal nodes

# Returns:
- Vector of primal modes for BDDC preconditioner

# Example:
```jldoctest
# setup
mesh = Mesh2D(1.0, 1.0);
p = 4
diffusion = GalleryOperator("diffusion", p, p, mesh);
bddc = LumpedBDDC(diffusion)

# note: either syntax works
primalmodes = LFAToolkit.getprimalmodes(bddc);
primalmodes = bddc.primalmodes;

# verify
@assert primalmodes == [1]

# output

```
"""
function getprimalmodes(bddc::BDDC)
    # assemble if needed
    if !isdefined(bddc, :primalmodes)
        primalmodes::Array{Int,1} = []
        numberprimalnodes = maximum(size(bddc.primalnodes))
        numbermodes, _ = size(bddc.operator.rowmodemap)
        modemap = bddc.operator.rowmodemap[:, bddc.primalnodes] * ones(numberprimalnodes)
        for i = 1:numbermodes
            if modemap[i] > 0.0
                push!(primalmodes, i)
            end
        end

        # store
        bddc.primalmodes = primalmodes
    end

    # return
    return getfield(bddc, :primalmodes)
end

"""
```julia
getsubassembledmodes(bddc)
```

Compute or retrieve the subassembled modes for a BDDC preconditioner

# Arguments:
- `preconditioner`:  preconditioner to compute subassembled modes

# Returns:
- Vector of subassembled modes for BDDC preconditioner

# Example:
```jldoctest
# setup
mesh = Mesh2D(1.0, 1.0);
p = 4
diffusion = GalleryOperator("diffusion", p, p, mesh);
bddc = LumpedBDDC(diffusion)

# note: either syntax works
subassembledmodes = LFAToolkit.getsubassembledmodes(bddc);
subassembledmodes = bddc.subassembledmodes;

# verify
@assert subassembledmodes == setdiff(1:(p-1)^2, [1])

# output

```
"""
function getsubassembledmodes(bddc::BDDC)
    # assemble if needed
    if !isdefined(bddc, :subassembledmodes)
        numbermodes, _ = size(bddc.operator.rowmodemap)

        # store
        bddc.subassembledmodes = setdiff(1:numbermodes, bddc.primalmodes)
    end

    # return
    return getfield(bddc, :subassembledmodes)
end

"""
```julia
getinterfacemodes(bddc)
```

Compute or retrieve the interface modes for a BDDC preconditioner

# Arguments:
- `preconditioner`:  preconditioner to compute interface modes

# Returns:
- Vector of interface modes for BDDC preconditioner

# Example:
```jldoctest
# setup
mesh = Mesh2D(1.0, 1.0);
p = 4
diffusion = GalleryOperator("diffusion", p, p, mesh);
bddc = LumpedBDDC(diffusion)

# note: either syntax works
interfacemodes = LFAToolkit.getinterfacemodes(bddc);
interfacemodes = bddc.interfacemodes;

# verify
trueinterfacemodes = [1:p-1...]
for i = 1:p-2
    push!(trueinterfacemodes, i*(p-1)+1)
end
trueinterfacemodes = sort(trueinterfacemodes)
trueinterfacemodes = setdiff(trueinterfacemodes, bddc.primalmodes)
@assert interfacemodes == trueinterfacemodes

# output

```
"""
function getinterfacemodes(bddc::BDDC)
    # assemble if needed
    if !isdefined(bddc, :interfacemodes)
        interfacemodes::Array{Int,1} = []
        numberinterfacenodes = maximum(size(bddc.interfacenodes))
        numbermodes, _ = size(bddc.operator.rowmodemap)
        modemap =
            bddc.operator.rowmodemap[:, bddc.interfacenodes] * ones(numberinterfacenodes)
        for i = 1:numbermodes
            if modemap[i] > 0.0
                push!(interfacemodes, i)
            end
        end
        interfacemodes = setdiff(interfacemodes, bddc.primalmodes)

        # store
        bddc.interfacemodes = interfacemodes
    end

    # return
    return getfield(bddc, :interfacemodes)
end

"""
```julia
getinteriormodes(bddc)
```

Compute or retrieve the interior modes for a BDDC preconditioner

# Arguments:
- `preconditioner`:  preconditioner to compute interior modes

# Returns:
- Vector of interior modes for BDDC preconditioner

# Example:
```jldoctest
# setup
mesh = Mesh2D(1.0, 1.0);
p = 4
diffusion = GalleryOperator("diffusion", p, p, mesh);
bddc = LumpedBDDC(diffusion)

# note: either syntax works
interiormodes = LFAToolkit.getinteriormodes(bddc);
interiormodes = bddc.interiormodes;

# verify
@assert interiormodes == setdiff(bddc.subassembledmodes, bddc.interfacemodes)

# output

```
"""
function getinteriormodes(bddc::BDDC)
    # assemble if needed
    if !isdefined(bddc, :interiormodes)
        # store
        bddc.interiormodes = setdiff(bddc.subassembledmodes, bddc.interfacemodes)
    end

    # return
    return getfield(bddc, :interiormodes)
end

"""
```julia
getsubassembledinverse(bddc)
```

Compute or retrieve the solver for subdomain for a BDDC preconditioner

# Arguments:
- `preconditioner`:  preconditioner to compute subdomain solver

# Returns:
- Solver for subdomain

# Example:
```jldoctest
# setup
mesh = Mesh2D(1.0, 1.0);
p = 4
diffusion = GalleryOperator("diffusion", p, p, mesh);
bddc = LumpedBDDC(diffusion)

# note: either syntax works
subassembledinverse = LFAToolkit.getsubassembledinverse(bddc);
subassembledinverse = bddc.subassembledinverse;

# verify
@assert subassembledinverse ≈ Matrix(diffusion.elementmatrix[bddc.subassemblednodes, bddc.subassemblednodes])^-1

# output

```
"""
function getsubassembledinverse(bddc::BDDC)
    # assemble if needed
    if !isdefined(bddc, :subassembledinverse)
        elementmatrix = bddc.operator.elementmatrix
        subassembledmatrix =
            Matrix(elementmatrix[bddc.subassemblednodes, bddc.subassemblednodes])

        # store
        bddc.subassembledinverse = subassembledmatrix^-1
    end

    # return
    return getfield(bddc, :subassembledinverse)
end

"""
```julia
getinteriorinverse(bddc)
```

Compute or retrieve the solver for subdomain interior for a BDDC preconditioner

# Arguments:
- `preconditioner`:  preconditioner to compute subdomain interior solver

# Returns:
- Solver for subdomain interior

# Example:
```jldoctest
# setup
mesh = Mesh2D(1.0, 1.0);
p = 4
diffusion = GalleryOperator("diffusion", p, p, mesh);
bddc = LumpedBDDC(diffusion)

# note: either syntax works
interiorinverse = LFAToolkit.getinteriorinverse(bddc);
interiorinverse = bddc.interiorinverse;

# verify
@assert interiorinverse ≈ Matrix(diffusion.elementmatrix[bddc.interiornodes, bddc.interiornodes])^-1

# output

```
"""
function getinteriorinverse(bddc::BDDC)
    # assemble if needed
    if !isdefined(bddc, :interiorinverse)
        elementmatrix = bddc.operator.elementmatrix
        interiormatrix = Matrix(elementmatrix[bddc.interiornodes, bddc.interiornodes])

        # store
        bddc.interiorinverse = interiormatrix^-1
    end

    # return
    return getfield(bddc, :interiorinverse)
end

"""
```julia
getschur(bddc)
```

Compute or retrieve the Schur complement matrix for a BDDC preconditioner

# Arguments:
- `preconditioner`:  preconditioner to compute Schur complement matrix

# Returns:
- Matrix for Schur complement

# Example:
```jldoctest
# setup
mesh = Mesh2D(1.0, 1.0);
p = 4
diffusion = GalleryOperator("diffusion", p, p, mesh);
bddc = LumpedBDDC(diffusion)

# note: either syntax works
schur = LFAToolkit.getschur(bddc);
schur = bddc.schur;

# verify
@assert schur ≈ diffusion.elementmatrix[bddc.primalnodes, bddc.primalnodes] -
    diffusion.elementmatrix[bddc.primalnodes, bddc.subassemblednodes] *
    Matrix(diffusion.elementmatrix[bddc.subassemblednodes, bddc.subassemblednodes])^-1 *
    diffusion.elementmatrix[bddc.subassemblednodes, bddc.primalnodes]

# output

```
"""
function getschur(bddc::BDDC)
    # assemble if needed
    if !isdefined(bddc, :schur)
        elementmatrix = bddc.operator.elementmatrix
        mixedmatrix = elementmatrix[bddc.primalnodes, bddc.subassemblednodes]
        schur =
            elementmatrix[bddc.primalnodes, bddc.primalnodes] -
            mixedmatrix * bddc.subassembledinverse * mixedmatrix'

        # store
        bddc.schur = schur
    end

    # return
    return getfield(bddc, :schur)
end

"""
```julia
getprimalrowmodemap(bddc)
```

Compute or retrieve the matrix mapping the rows of the primal BDDC matrix to the primal symbol matrix

# Arguments:
- `preconditioner`:  preconditioner to compute primal row mode map

# Returns:
- Matrix mapping rows of primal BDDC matrix to primal symbol matrix

# Example:
```jldoctest
# setup
mesh = Mesh2D(1.0, 1.0);
diffusion = GalleryOperator("diffusion", 3, 3, mesh);
bddc = LumpedBDDC(diffusion);

# note: either syntax works
primalmodemap = LFAToolkit.getprimalrowmodemap(bddc);
primalmodemap = bddc.primalrowmodemap;

# verify
@assert primalmodemap ≈ [1 1 1 1]

# output

```
"""
function getprimalrowmodemap(bddc::BDDC)
    # assemble if needed
    if !isdefined(bddc, :primalrowmodemap)
        primalrowmodemap = bddc.operator.rowmodemap[bddc.primalmodes, bddc.primalnodes]

        # store
        bddc.primalrowmodemap = primalrowmodemap
    end

    # return
    return getfield(bddc, :primalrowmodemap)
end

"""
```julia
getprimalcolumnmodemap(bddc)
```

Compute or retrieve the matrix mapping the columns of the primal BDDC matrix to the
    primal symbol matrix

# Arguments:
- `preconditioner`:  preconditioner to compute primal column mode map

# Returns:
- Matrix mapping columns of primal BDDC matrix to primal symbol matrix

# Example:
```jldoctest
# setup
mesh = Mesh2D(1.0, 1.0);
diffusion = GalleryOperator("diffusion", 3, 3, mesh);
bddc = LumpedBDDC(diffusion);

# note: either syntax works
primalmodemap = LFAToolkit.getprimalcolumnmodemap(bddc);
primalmodemap = bddc.primalcolumnmodemap;

# verify
@assert primalmodemap ≈ [1; 1; 1; 1]

# output

```
"""
function getprimalcolumnmodemap(bddc::BDDC)
    # assemble if needed
    if !isdefined(bddc, :primalcolumnmodemap)
        primalcolumnmodemap =
            bddc.operator.columnmodemap[bddc.primalnodes, bddc.primalmodes]

        # store
        bddc.primalcolumnmodemap = primalcolumnmodemap
    end

    # return
    return getfield(bddc, :primalcolumnmodemap)
end

"""
```julia
getsubassembledrowmodemap(bddc)
```

Compute or retrieve the matrix mapping the rows of the subassembled BDDC matrix to the
    subassembled symbol matrix

# Arguments:
- `preconditioner`:  preconditioner to compute subassembled row mode map

# Returns:
- Matrix mapping rows of subassembled BDDC matrix to subassembled symbol matrix

# Example:
```jldoctest
# setup
mesh = Mesh2D(1.0, 1.0);
diffusion = GalleryOperator("diffusion", 3, 3, mesh);
bddc = LumpedBDDC(diffusion);

# note: either syntax works
subassembledmodemap = LFAToolkit.getsubassembledrowmodemap(bddc);
subassembledmodemap = bddc.subassembledrowmodemap;

# verify
@assert subassembledmodemap ≈ bddc.operator.rowmodemap[bddc.subassembledmodes, bddc.subassemblednodes]

# output

```
"""
function getsubassembledrowmodemap(bddc::BDDC)
    # assemble if needed
    if !isdefined(bddc, :subassembledrowmodemap)
        subassembledrowmodemap =
            bddc.operator.rowmodemap[bddc.subassembledmodes, bddc.subassemblednodes]

        # store
        bddc.subassembledrowmodemap = subassembledrowmodemap
    end

    # return
    return getfield(bddc, :subassembledrowmodemap)
end

"""
```julia
getsubassembledcolumnmodemap(bddc)
```

Compute or retrieve the matrix mapping the columns of the subassembled BDDC matrix to the
    subassembled symbol matrix

# Arguments:
- `preconditioner`:  preconditioner to compute subassembled column mode map

# Returns:
- Matrix mapping columns of subassembled BDDC matrix to subassembled symbol matrix

# Example:
```jldoctest
# setup
mesh = Mesh2D(1.0, 1.0);
diffusion = GalleryOperator("diffusion", 3, 3, mesh);
bddc = LumpedBDDC(diffusion);

# note: either syntax works
subassembledmodemap = LFAToolkit.getsubassembledcolumnmodemap(bddc);
subassembledmodemap = bddc.subassembledcolumnmodemap;

# verify
@assert subassembledmodemap ≈ bddc.operator.columnmodemap[bddc.subassemblednodes, bddc.subassembledmodes]

# output

```
"""
function getsubassembledcolumnmodemap(bddc::BDDC)
    # assemble if needed
    if !isdefined(bddc, :subassembledcolumnmodemap)
        subassembledcolumnmodemap =
            bddc.operator.columnmodemap[bddc.subassemblednodes, bddc.subassembledmodes]

        # store
        bddc.subassembledcolumnmodemap = subassembledcolumnmodemap
    end

    # return
    return getfield(bddc, :subassembledcolumnmodemap)
end

"""
```julia
getmixedrowmodemap(bddc)
```

Compute or retrieve the matrix mapping the rows of the mixed BDDC matrix to the
    symbol matrix

# Arguments:
- `preconditioner`:  preconditioner to compute mixed row mode map

# Returns:
- Matrix mapping rows of mixed BDDC matrix to symbol matrix
"""
function getmixedrowmodemap(bddc::BDDC)
    # assemble if needed
    if !isdefined(bddc, :mixedrowmodemap)
        numberprimalmodes = maximum(size(bddc.primalmodes))
        numbersubassembledmodes = maximum(size(bddc.subassembledmodes))
        numbersubassemblednodes = maximum(size(bddc.subassemblednodes))
        mixedrowmodemap = [
            I(numberprimalmodes) zeros((numberprimalmodes, numbersubassemblednodes))
            zeros(numbersubassembledmodes, numberprimalmodes) bddc.subassembledrowmodemap
        ]

        # store
        bddc.mixedrowmodemap = bddc.modepermutation * mixedrowmodemap
    end

    # return
    return getfield(bddc, :mixedrowmodemap)
end

"""
```julia
getmixedcolumnmodemap(bddc)
```

Compute or retrieve the matrix mapping the columns of the mixed BDDC matrix to the
    symbol matrix

# Arguments:
- `preconditioner`:  preconditioner to compute mixed column mode map

# Returns:
- Matrix mapping columns of mixed BDDC matrix to symbol matrix
"""
function getmixedcolumnmodemap(bddc::BDDC)
    # assemble if needed
    if !isdefined(bddc, :mixedcolumnmodemap)
        numberprimalmodes = maximum(size(bddc.primalmodes))
        numbersubassembledmodes = maximum(size(bddc.subassembledmodes))
        numbersubassemblednodes = maximum(size(bddc.subassemblednodes))
        mixedcolumnmodemap = [
            I(numberprimalmodes) zeros((numberprimalmodes, numbersubassembledmodes))
            zeros(numbersubassemblednodes, numberprimalmodes) bddc.subassembledcolumnmodemap
        ]

        # store
        bddc.mixedcolumnmodemap = mixedcolumnmodemap
    end

    # return
    return getfield(bddc, :mixedcolumnmodemap)

end

"""
```julia
getmodepermutation(bddc)
```

Compute or retrieve the matrix permuting multi-component modes to standard ordering

# Arguments:
- `preconditioner`:  preconditioner to compute mode permutation matrix

# Returns:
- Matrix mapping BDDC mode ordering to standard ordering
"""
function getmodepermutation(bddc::BDDC)
    # assemble if needed
    if !isdefined(bddc, :modepermutation)
        numberprimal = maximum(size(bddc.primalmodes))
        numbersubassembled = maximum(size(bddc.subassembledmodes))
        numbermodes = numberprimal + numbersubassembled
        modepermutation = spzeros(Bool, numbermodes, numbermodes)

        # build permutation matrix
        for i = 1:numberprimal
            modepermutation[bddc.primalmodes[i], i] = true
        end
        for i = 1:numbersubassembled
            modepermutation[bddc.subassembledmodes[i], i+numberprimal] = true
        end

        # store
        bddc.modepermutation = modepermutation
    end

    # return
    return getfield(bddc, :modepermutation)
end

"""
```julia
getmixedmultiplicity(bddc)
```

Compute or retrieve the diagonal matrix of mixed interface node and primal mode
    multiplicity for the BDDC preconditioner

# Arguments:
- `preconditioner`:  preconditioner to compute mixed node/mode multiplicity

# Returns:
- Matrix of mixed multiplicity for the BDDC preconditioner
"""
function getmixedmultiplicity(bddc::BDDC)
    # assemble if needed
    if !isdefined(bddc, :mixedmultiplicity)
        numberprimalmodes = maximum(size(bddc.primalmodes))
        mixedmultiplicity = Diagonal(
            vcat(
                ones(numberprimalmodes, 1)...,
                [1 / m for m in bddc.operator.multiplicity[bddc.subassemblednodes]]...,
            ),
        )

        # store
        bddc.mixedmultiplicity = mixedmultiplicity
    end

    # return
    return getfield(bddc, :mixedmultiplicity)
end

"""
```julia
getJDT(bddc)
```

Compute or retrieve the matrix mapping the jump over subdomain interface modes

# Arguments:
- `preconditioner`:  preconditioner to compute jump mapping matrix

# Returns:
- Matrix mapping the jump over subdomain interfacemodes
"""
function getJDT(bddc::BDDC)
    # assemble if needed
    if !isdefined(bddc, :JDT)
        # setup
        numberprimalmodes = maximum(size(bddc.primalmodes))
        numbernodes, _ = size(bddc.operator.columnmodemap)
        numbermixed, _ = size(bddc.mixedcolumnmodemap)
        numberinterfacenodes = maximum(size(bddc.interfacenodes))

        # build jump map
        J_D_T_nodes = zeros(numbernodes, numbernodes)
        for i = 1:numberinterfacenodes, j = 1:numberinterfacenodes
            indxi = bddc.interfacenodes[i]
            indxj = bddc.interfacenodes[j]
            if bddc.operator.rowmodemap[:, indxi] == bddc.operator.rowmodemap[:, indxj]
                scale = -1 / bddc.operator.multiplicity[indxi]
                if i == j
                    scale = 1 + scale
                end
                J_D_T_nodes[indxi, indxj] = scale
            end
        end
        J_D_T_mixed = spzeros(numbermixed, numbermixed)
        J_D_T_mixed[numberprimalmodes+1:end, numberprimalmodes+1:end] =
            J_D_T_nodes[bddc.subassemblednodes, bddc.subassemblednodes]

        # store
        bddc.JDT = J_D_T_mixed
    end

    # return
    return getfield(bddc, :JDT)
end

# ------------------------------------------------------------------------------
# data for computing symbols
# ------------------------------------------------------------------------------

"""
```julia
computesymbolsinjection(bddc)
```

Compute or retrieve the injection operator for the BDDC symbol matrix

# Arguments:
- `preconditioner`:  preconditioner to compute injection matrix

# Returns:
- Matrix providing the injection operator of BDDC symbol matrix

"""
function computesymbolsinjection(bddc::BDDC, θ::Array)
    if (bddc.injectiontype == BDDCInjectionType.scaled)
        # lumped BDDC
        return bddc.mixedmultiplicity
    elseif (bddc.injectiontype == BDDCInjectionType.harmonic)
        # Dirichlet BDDC
        # -- validity check
        dimension = length(θ)
        if dimension != bddc.operator.inputs[1].basis.dimension
            # COV_EXCL_START
            throw(
                ArgumentError(
                    "Must provide as many values of θ as the mesh has dimensions",
                ),
            )
            # COV_EXCL_STOP
        end

        # -- setup
        numberprimalmodes = maximum(size(bddc.primalmodes))
        numbernodes, _ = size(bddc.operator.columnmodemap)
        numbermixed, _ = size(bddc.mixedcolumnmodemap)
        numberinteriornodes = maximum(size(bddc.interiornodes))
        numberinterfacenodes = maximum(size(bddc.interfacenodes))
        nodecoordinatedifferences = bddc.operator.nodecoordinatedifferences
        elementmatrix = bddc.operator.elementmatrix

        # -- harmonic operator
        A_ΓI_nodes = zeros(ComplexF64, numbernodes, numbernodes)
        for i = 1:numberinterfacenodes, j = 1:numberinteriornodes
            indxi = bddc.interfacenodes[i]
            indxj = bddc.interiornodes[j]
            A_ΓI_nodes[indxi, indxj] =
                elementmatrix[indxi, indxj] *
                ℯ^(
                    im * sum([
                        θ[k] * nodecoordinatedifferences[indxi, indxj, k] for
                        k = 1:dimension
                    ])
                )
        end
        A_II_inv = bddc.interiorinverse
        A_II_inv_nodes = zeros(ComplexF64, numbernodes, numbernodes)
        for i = 1:numberinteriornodes, j = 1:numberinteriornodes
            indxi = bddc.interiornodes[i]
            indxj = bddc.interiornodes[j]
            A_II_inv_nodes[indxi, indxj] =
                A_II_inv[i, j] *
                ℯ^(
                    im * sum([
                        θ[k] * nodecoordinatedifferences[indxi, indxj, k] for
                        k = 1:dimension
                    ])
                )
        end
        𝓗_T = A_ΓI_nodes * A_II_inv_nodes
        𝓗_T_mixed = spzeros(ComplexF64, numbermixed, numbermixed)
        𝓗_T_mixed[numberprimalmodes+1:end, numberprimalmodes+1:end] =
            𝓗_T[bddc.subassemblednodes, bddc.subassemblednodes]

        # -- jump mapping
        J_D_T_mixed = bddc.JDT

        # -- injection mapping
        return bddc.mixedmultiplicity + J_D_T_mixed * 𝓗_T_mixed
    else
        throw(ArgumentError("Injection type unknown")) # COV_EXCL_LINE
    end
end

# ------------------------------------------------------------------------------
# get/set property
# ------------------------------------------------------------------------------

function Base.getproperty(bddc::BDDC, f::Symbol)
    if f == :primalnodes
        return getprimalnodes(bddc)
    elseif f == :subassemblednodes
        return getsubassemblednodes(bddc)
    elseif f == :interfacenodes
        return getinterfacenodes(bddc)
    elseif f == :interiornodes
        return getinteriornodes(bddc)
    elseif f == :primalmodes
        return getprimalmodes(bddc)
    elseif f == :subassembledmodes
        return getsubassembledmodes(bddc)
    elseif f == :interfacemodes
        return getinterfacemodes(bddc)
    elseif f == :interiormodes
        return getinteriormodes(bddc)
    elseif f == :subassembledinverse
        return getsubassembledinverse(bddc)
    elseif f == :interiorinverse
        return getinteriorinverse(bddc)
    elseif f == :schur
        return getschur(bddc)
    elseif f == :mixedmultiplicity
        return getmixedmultiplicity(bddc)
    elseif f == :JDT
        return getJDT(bddc)
    elseif f == :primalrowmodemap
        return getprimalrowmodemap(bddc)
    elseif f == :primalcolumnmodemap
        return getprimalcolumnmodemap(bddc)
    elseif f == :subassembledrowmodemap
        return getsubassembledrowmodemap(bddc)
    elseif f == :subassembledcolumnmodemap
        return getsubassembledcolumnmodemap(bddc)
    elseif f == :mixedrowmodemap
        return getmixedrowmodemap(bddc)
    elseif f == :mixedcolumnmodemap
        return getmixedcolumnmodemap(bddc)
    elseif f == :modepermutation
        return getmodepermutation(bddc)
    else
        return getfield(bddc, f)
    end
end

function Base.setproperty!(bddc::BDDC, f::Symbol, value)
    if f == :operator || f == :injectiontype
        throw(ReadOnlyMemoryError()) # COV_EXCL_LINE
    else
        return setfield!(bddc, f, value)
    end
end

# ------------------------------------------------------------------------------
# compute symbols
# ------------------------------------------------------------------------------

"""
```julia
computesymbols(bddc, θ)
```

Compute the symbol matrix for a BDDC preconditioned operator

# Arguments:
- `bddc`:  BDDC preconditioner to compute symbol matrix for
- `θ`:     Fourier mode frequency array (one frequency per dimension)

# Returns:
- Symbol matrix for the BDDC preconditioned operator

# Lumped BDDC Example:
```jldoctest
using LinearAlgebra;

for dimension in 2:3
    # setup
    mesh = []
    if dimension == 2
        mesh = Mesh2D(1.0, 1.0);
    elseif dimension == 3
        mesh = Mesh3D(1.0, 1.0, 1.0);
    end
    diffusion = GalleryOperator("diffusion", 3, 3, mesh);
    bddc = LumpedBDDC(diffusion)

    # compute symbols
    A = computesymbols(bddc, [0.2], π*ones(dimension));

    # verify
    eigenvalues = real(eigvals(A));
    if dimension == 2
        @assert minimum(eigenvalues) ≈ 0.43999999999999995
        @assert maximum(eigenvalues) ≈ 0.8
    elseif dimension == 3
        @assert minimum(eigenvalues) ≈ -0.6319999999999972
        @assert maximum(eigenvalues) ≈ 0.8
    end
end

# output

```

# Dirichlet BDDC Example:
```jldoctest
using LinearAlgebra;

for dimension in 2:3
    # setup
    mesh = []
    if dimension == 2
        mesh = Mesh2D(1.0, 1.0);
    elseif dimension == 3
        mesh = Mesh3D(1.0, 1.0, 1.0);
    end
    diffusion = GalleryOperator("diffusion", 3, 3, mesh);
    bddc = DirichletBDDC(diffusion)

    # compute symbols
    A = computesymbols(bddc, [0.2], π*ones(dimension));

    # verify
    eigenvalues = real(eigvals(A));
    if dimension == 2
        @assert minimum(eigenvalues) ≈ 0.7999999999999998
        @assert maximum(eigenvalues) ≈ 0.8
    elseif dimension == 3
        @assert minimum(eigenvalues) ≈ 0.7801226993865031
        @assert maximum(eigenvalues) ≈ 0.8
    end
end

# output

```
"""
function computesymbols(bddc::BDDC, ω::Array, θ::Array)
    # validity check
    if length(ω) != 1
        Throw(error("exactly one parameter required for BDDC smoothing")) # COV_EXCL_LINE
    end
    dimension = length(θ)
    if dimension != bddc.operator.inputs[1].basis.dimension
        throw(ArgumentError("Must provide as many values of θ as the mesh has dimensions")) # COV_EXCL_LINE
    end

    # setup
    elementmatrix = bddc.operator.elementmatrix
    numberprimalnodes = maximum(size(bddc.primalnodes))
    numberprimalmodes = maximum(size(bddc.primalmodes))
    numbersubassemblednodes = maximum(size(bddc.subassemblednodes))
    nodecoordinatedifferences = bddc.operator.nodecoordinatedifferences

    # subdomain solver
    A_rr_inv = bddc.subassembledinverse
    A_rr_inv_nodes = zeros(ComplexF64, numbersubassemblednodes, numbersubassemblednodes)
    for i = 1:numbersubassemblednodes, j = 1:numbersubassemblednodes
        indxi = bddc.subassemblednodes[i]
        indxj = bddc.subassemblednodes[j]
        A_rr_inv_nodes[i, j] =
            A_rr_inv[i, j] *
            ℯ^(
                im * sum([
                    θ[k] * nodecoordinatedifferences[indxi, indxj, k] for k = 1:dimension
                ])
            )
    end

    # mixed subassembled primal matrices
    Â_Πr_nodes = zeros(ComplexF64, numberprimalnodes, numbersubassemblednodes)
    for i = 1:numberprimalnodes, j = 1:numbersubassemblednodes
        indxi = bddc.primalnodes[i]
        indxj = bddc.subassemblednodes[j]
        Â_Πr_nodes[i, j] =
            elementmatrix[indxi, indxj] *
            ℯ^(
                im * sum([
                    θ[k] * nodecoordinatedifferences[indxi, indxj, k] for k = 1:dimension
                ])
            )
    end
    Â_Πr_modes = bddc.primalrowmodemap * Â_Πr_nodes

    # Schur complement
    Ŝ_Π = bddc.schur
    Ŝ_Π_nodes = zeros(ComplexF64, numberprimalnodes, numberprimalnodes)
    for i = 1:numberprimalnodes, j = 1:numberprimalnodes
        indxi = bddc.primalnodes[i]
        indxj = bddc.primalnodes[j]
        Ŝ_Π_nodes[i, j] =
            Ŝ_Π[i, j] *
            ℯ^(
                im * sum([
                    θ[k] * nodecoordinatedifferences[indxi, indxj, k] for k = 1:dimension
                ])
            )
    end
    Ŝ_Π_inv_modes = (bddc.primalrowmodemap * Ŝ_Π_nodes * bddc.primalcolumnmodemap)^-1

    # subassembled nodes primal modes
    Ø = zeros((numberprimalmodes, numbersubassemblednodes))
    K_u_inv = [
        I(numberprimalmodes) Ø
        -A_rr_inv_nodes*Â_Πr_modes' I(numbersubassemblednodes)
    ]
    P_inv = [
        Ŝ_Π_inv_modes Ø
        Ø' A_rr_inv_nodes
    ]
    Â_inv_mixed = K_u_inv * P_inv * K_u_inv'

    # injection
    R_mixed = computesymbolsinjection(bddc, θ)
    R_T_Â_inv_R_modes =
        bddc.mixedrowmodemap * R_mixed' * Â_inv_mixed * R_mixed * bddc.mixedcolumnmodemap

    # return
    return I - ω[1] * R_T_Â_inv_R_modes * computesymbols(bddc.operator, θ)
end

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


text:
module TestExamples

using BenchmarkConfigSweeps
using Tables
using Test

function test_example1()
    configs = Iterators.product(
        BenchmarkConfigSweeps.nthreads.(1:2),
        BenchmarkConfigSweeps.env.(
            "BCST_A" .=> [33, 44],  # read in ./example1.jl
            "BCST_B" => "not varied",
        ),
    )
    mktempdir() do resultdir
        BenchmarkConfigSweeps.run(resultdir, joinpath(@__DIR__, "example1.jl"), configs)

        sweepresult = BenchmarkConfigSweeps.load(resultdir)

        st = BenchmarkConfigSweeps.simpletable(sweepresult)
        @test length(st.env) == 4

        ft = Tables.columns(sweepresult)
        @test ft.nthreads == [1, 1, 2, 2, 1, 1, 2, 2]
        @test ft.BCST_A == [33, 33, 33, 33, 44, 44, 44, 44]
        @test ft.k1 == [0, 33, 0, 33, 0, 44, 0, 44]
    end
end

end  # module