text:
import Base: unique, unique!, ==
using ModelingToolkit: value, operation, arguments, istree

"""
$(TYPEDEF)

A basis over the variables `u` with parameters `p` and independent variable `iv`.
It extends an `AbstractSystem` as defined in `ModelingToolkit.jl`.
`f` can either be a Julia function which is able to use ModelingToolkit variables or
a vector of `eqs`.
It can be called with the typical DiffEq signature, meaning out of place with `f(u,p,t)`
or in place with `f(du, u, p, t)`.
If `linear_independent` is set to `true`, a linear independent basis is created from all atom function in `f`.
If `simplify_eqs` is set to `true`, `simplify` is called on `f`.
Additional keyworded arguments include `name`, which can be used to name the basis, and
`observed` for defining observeables.

# Fields
$(FIELDS)

# Example

```julia
using ModelingToolkit
using DataDrivenDiffEq

@parameters w[1:2] t
@variables u[1:2](t)

Ψ = Basis([u; sin.(w.*u)], u, parameters = p, iv = t)
```

## Note

The keyword argument `eval_expression` controls the function creation
behavior. `eval_expression=true` means that `eval` is used, so normal
world-age behavior applies (i.e. the functions cannot be called from
the function that generates them). If `eval_expression=false`,
then construction via GeneralizedGenerated.jl is utilized to allow for
same world-age evaluation. However, this can cause Julia to segfault
on sufficiently large basis functions. By default eval_expression=false.

"""
mutable struct Basis <: ModelingToolkit.AbstractSystem
    """The equations of the basis"""
    eqs::Vector{Equation}
    """Dependent (state) variables"""
    states::Vector
    """Parameters"""
    ps::Vector
    observed::Vector
    """Independent variable"""
    iv::Num
    """Internal function representation of the basis"""
    f_::Function
    """Name of the basis"""
    name::Symbol
    """Internal systems"""
    systems::Vector{Basis}
end

function Basis(eqs::AbstractVector, states::AbstractVector; parameters::AbstractArray = [], iv = nothing,

    simplify = false, linear_independent = false, name = gensym(:Basis), eval_expression = false,
    observed = [],
    kwargs...)

    eqs = simplify ? ModelingToolkit.simplify.(eqs) : eqs
    eqs = linear_independent ? create_linear_independent_eqs(eqs) : eqs
    isnothing(iv) && (iv = Num(Variable(:t)))
    unique!(eqs, !simplify)

    if eval_expression
        f_oop, f_iip = eval.(build_function(eqs, value.(states), value.(parameters), [value(iv)], expression = Val{true}))
    else
        f_oop, f_iip = build_function(eqs, value.(states), value.(parameters), [value(iv)], expression = Val{false})
    end
    eqs = [Variable(:φ,i) ~ eq for (i,eq) ∈ enumerate(eqs)]
    f_(u,p,t) = f_oop(u,p,t)
    f_(du, u, p, t) = f_iip(du, u, p, t)

    return Basis(eqs, value.(states), value.(parameters), observed, value(iv), f_, name, Basis[])
end

function Basis(f::Function, states::AbstractVector; parameters::AbstractArray = [], iv = nothing, kwargs...)

    isnothing(iv) && (iv = Num(Variable(:t)))

    try
        eqs = f(states, parameters, iv)
        return Basis(eqs, states, parameters = parameters, iv = iv; kwargs...)
    catch e
        rethrow(e)
    end
end

_get_name(x::Num) = _get_name(value(x))
_get_name(x) = nameof(istree(x) ? operation(x) : x)

Base.show(io::IO, x::Basis) = print(io, "$(String.(x.name)) : $(length(x.eqs)) dimensional basis in ", "$(String.([_get_name(v) for v in x.states]))")

@inline function Base.print(io::IO, x::Basis)
    show(io, x)
    !isempty(x.ps) && println(io, "\nParameters : $(x.ps)")
    println(io, "\nIndependent variable: $(x.iv)")
    println(io, "Equations")
    for (i,eq) ∈ enumerate(x.eqs)
        if i < 5 || i == length(x.eqs)
        println(io, "$(eq.lhs) = $(eq.rhs)")
        elseif i == 5
            println(io, "...")
        else
            continue
        end
    end
end

@inline function Base.println(io::IO, x::Basis, fullview::DataType = Val{false})
    fullview == Val{false} && return print(io, x)
    show(io, x)
    !isempty(x.ps) && println(io, "\nParameters : $(x.ps)")
    println(io, "\nIndependent variable: $(x.iv)")
    println(io, "Equations")
    for (i,eq) ∈ enumerate(x.eqs)
        println(io, "$(eq.lhs) = $(eq.rhs)")
    end
end

"""

"""
variables(x::Basis) = x.states

function update!(b::Basis, eval_expression = false)

    if eval_expression
        f_oop, f_iip = eval.(ModelingToolkit.build_function([bi.rhs for bi in b.eqs], b.states, b.ps, [b.iv], expression = Val{true}))
    else
        f_oop, f_iip = ModelingToolkit.build_function([bi.rhs for bi in b.eqs], b.states, b.ps, [b.iv], expression = Val{false})
    end

    f_(u, p, t) = f_oop(u, p, t)
    f_(du, u, p, t) = f_iip(du, u, p, t)

    b.f_ = f_
    return
end

function unique(b::AbstractArray{Num}, simplify_eqs = false)
    b = simplify_eqs ? simplify.(b) : b
    returns = ones(Bool, size(b)...)
    N = maximum(eachindex(b))
    for i ∈ eachindex(b)
        returns[i] = !any([isequal(b[i], b[j]) for j in i+1:N])
    end
    return b[returns]
end

function Base.unique!(b::AbstractArray, simplify_eqs = false)
    bs = simplify_eqs ? simplify.(b) : b
    removes = zeros(Bool, size(bs)...)
    N = maximum(eachindex(bs))
    for i ∈ eachindex(bs)
        removes[i] = any([isequal(bs[i], bs[j]) for j in i+1:N])
    end
    deleteat!(b, removes)
end

function unique(b::AbstractArray{Equation}, simplify_eqs = false)
    b = simplify_eqs ? simplify.(b) : b
    returns = ones(Bool, size(b)...)
    N = maximum(eachindex(b))
    for i ∈ eachindex(b)
        returns[i] = !any([isequal(b[i].rhs, b[j].rhs) for j in i+1:N])
    end
    return b[returns]
end

function Base.unique!(b::AbstractArray{Equation}, simplify_eqs = false)
    bs = [bi.rhs for bi in b]
    bs = simplify_eqs ? simplify.(bs) : bs
    removes = zeros(Bool, size(bs)...)
    N = maximum(eachindex(bs))
    for i ∈ eachindex(bs)
        removes[i] = any([isequal(bs[i], bs[j]) for j in i+1:N])
    end
    deleteat!(b, removes)
end

function Base.unique!(b::Basis, simplify_eqs = false; eval_expression = false)
    unique!(b.eqs, simplify_eqs)
    update!(b, eval_expression)
end

function Base.unique(b::Basis; kwargs...)
    eqs = unique(map(x->x.rhs, equations(b)))
    return Basis(eqs, variables(b), parameters = parameters(b), iv = independent_variable(b), kwargs...)
end

"""
    deleteat!(basis, inds, eval_expression = false)

    Delete the entries specified by `inds` and update the `Basis` accordingly.
"""
function Base.deleteat!(b::Basis, inds; eval_expression = false)
    deleteat!(b.eqs, inds)
    update!(b, eval_expression)
    return
end

"""
    push!(basis, eq, simplify_eqs = true; eval_expression = false)

    Push the equations(s) in `eq` into the basis and update all internal fields accordingly.
    `eq` can either be a single equation or an array. If `simplify_eq` is true, the equation will be simplified.
"""
function Base.push!(b::Basis, eqs::AbstractArray, simplify_eqs = true; eval_expression = false)
    @inbounds for eq ∈ eqs
        push!(b, eq, false)
    end
    unique!(b, simplify_eqs, eval_expression = eval_expression)
    return
end

function Base.push!(b::Basis, eq, simplify_eqs = true; eval_expression = false)
    push!(b.eqs, Variable(:φ, length(b.eqs)+1)~eq)
    unique!(b, simplify_eqs, eval_expression = eval_expression)
    return
end

"""
    merge(x::Basis, y::Basis; eval_expression = false)

    Return a new `Basis`, which is defined via the union of `x` and `y` .
"""
function Base.merge(x::Basis, y::Basis; eval_expression = false)
    b =  unique(vcat([xi.rhs  for xi ∈ equations(x)], [xi.rhs  for xi ∈ equations(y)]))
    vs = unique(vcat(x.states, y.states))
    ps = unique(vcat(x.ps, y.ps))
    observed = unique(vcat(x.observed, y.observed))
    return Basis(Num.(b), vs, parameters = ps, observed = observed, eval_expression = eval_expression)
end

"""
    merge!(x::Basis, y::Basis; eval_expression = false)

    Updates `x` to include the union of both `x` and `y`.
"""
function Base.merge!(x::Basis, y::Basis; eval_expression = false)
    push!(x, map(x->x.rhs, equations(y)))
    x.states = unique(vcat(x.states, y.states))
    x.ps = unique(vcat(x.ps, y.ps))
    update!(x, eval_expression)
    return
end

Base.length(x::Basis) = length(x.eqs)
Base.size(x::Basis) = size(x.eqs)

Base.getindex(x::Basis, idx) = getindex(equations(x), idx)
Base.firstindex(x::Basis) = firstindex(equations(x))
Base.lastindex(x::Basis) = lastindex(equations(x))
Base.iterate(x::Basis) = iterate(equations(x))
Base.iterate(x::Basis, id) = iterate(equations(x), id)

function (==)(x::Basis, y::Basis)
    length(x) == length(y) || return false
    n = zeros(Bool, length(x))
    yrhs = [yi.rhs for yi in equations(y)]
    xrhs = [xi.rhs for xi in equations(x)]
    @inbounds for (i, xi) in enumerate(xrhs)
        n[i] = any(isequal.([xi], yrhs))
    end
    return all(n)
end

free_parameters(b::Basis; operations = [+]) = count_operation([xi.rhs for xi in b.eqs], operations) + length(b.eqs)

(b::Basis)(u, p::DiffEqBase.NullParameters, t) = b(u, [], t)
(b::Basis)(du, u, p::DiffEqBase.NullParameters, t) = b(du, u, [], t)
(b::Basis)(u::AbstractVector,  p::AbstractArray = [], t = nothing) = b.f_(u, isempty(p) ? parameters(b) : p, isnothing(t) ? zero(eltype(u)) : t)
(b::Basis)(du::AbstractVector, u::AbstractVector, p::AbstractArray = [], t = nothing) = b.f_(du, u, isempty(p) ? parameters(b) : p, isnothing(t) ? zero(eltype(u)) : t)

function (b::Basis)(x::AbstractMatrix, p::AbstractArray = [], t::AbstractArray = [])
    isempty(t) ? nothing : @assert size(x, 2) == length(t)

    if (isempty(p) || eltype(p) <: Num) && !isempty(parameters(b))
        pi = isempty(p) ? parameters(b) : p
        res = Array{Any}(undef,length(b), size(x)[2])
    else
        pi = p
        res = zeros(eltype(x), length(b), size(x)[2])
    end

    @inbounds for i in 1:size(x)[2]
        res[:, i] .= b.f_(x[:, i], isempty(p) ? parameters(b) : p, isempty(t) ? zero(eltype(x)) : t[i])
    end
    return res
end

function (b::Basis)(y::AbstractMatrix, x::AbstractMatrix, p::AbstractArray = [], t::AbstractArray = [])
    @assert size(x, 2) == size(y, 2) "Provide consistent arrays."
    @assert size(y, 1) == length(b) "Provide consistent arrays."
    isempty(t) ? nothing : @assert size(x, 2) == length(t)

    @inbounds for i in 1:size(x, 2)
        b.f_(view(y, :, i), view(x, :, i), isempty(p) ? parameters(b) : p, isempty(t) ? zero(eltype(x)) : t[i])
    end
end

"""
    jacobian(basis)

    Returns a function representing the jacobian matrix / gradient of the `Basis` with respect to the
    dependent variables as a function with the common signature `f(u,p,t)` for out of place and `f(du, u, p, t)` for in place computation.
"""
function jacobian(x::Basis, eval_expression = false)

    j = ModelingToolkit.jacobian([xi.rhs for xi in x.eqs], x.states)

    if eval_expression
        f_oop, f_iip = eval.(ModelingToolkit.build_function(expand_derivatives.(j), x.states, x.ps, [x.iv], expression = Val{true}))
    else
        f_oop, f_iip = ModelingToolkit.build_function(expand_derivatives.(j), x.states, x.ps, [x.iv], expression = Val{false})
    end

    f_(u, p, t) = f_oop(u, p, t)
    f_(du, u, p, t) = f_iip(du, u, p, t)

    return f_
end

"""
    dynamics(basis)

    Returns the internal function representing the dynamics of the `Basis`.
"""
function dynamics(b::Basis)
    return b.f_
end

## Symbolic Manipulation

function is_unary(f::Function)
    for m in methods(f)
        m.nargs - 1 > 1 && return false
    end
    return true
end

count_operation(x::Number, op::Function, nested::Bool = true) = 0
count_operation(x::Sym, op::Function, nested::Bool = true) = 0
count_operation(x::Num, op::Function, nested::Bool = true) = count_operation(value(x), op, nested)

function count_operation(x, op::Function, nested::Bool = true)
    if operation(x)== op
        if is_unary(op)
            # Handles sin, cos and stuff
            nested && return 1 + count_operation(arguments(x), op)
            return 1
        else
            # Handles +, *
            nested && length(arguments(x))-1 + count_operation(arguments(x), op)
            return length(arguments(x))-1
        end
    elseif nested
        return count_operation(arguments(x), op, nested)
    end
    return 0
end

function count_operation(x, ops::AbstractArray, nested::Bool = true)
    return sum([count_operation(x, op, nested) for op in ops])
end

function count_operation(x::AbstractArray, op::Function, nested::Bool = true)
    sum([count_operation(xi, op, nested) for xi in x])
end

function count_operation(x::AbstractArray, ops::AbstractArray, nested::Bool = true)
    counter = 0
    @inbounds for xi in x, op in ops
        counter += count_operation(xi, op, nested)
    end
    counter
end

function split_term!(x::AbstractArray, o, ops::AbstractArray = [+])
    if istree(o)
        n_ops = count_operation(o, ops, false)
        c_ops = 0
        @views begin
            if n_ops == 0
                x[begin]= o
            else
                counter_ = 1
                for oi in arguments(o)
                    c_ops = count_operation(oi, ops, false)
                    split_term!(x[counter_:counter_+c_ops], oi, ops)
                    counter_ += c_ops + 1
                end
            end
        end
    else
        x[begin] = o
    end
    return
end

split_term!(x::AbstractArray,o::Num, ops::AbstractArray = [+]) = split_term!(x, value(o), ops)

remove_constant_factor(x::Num) = remove_constant_factor(value(x))
remove_constant_factor(x::Number) = one(x)

function remove_constant_factor(x)
    istree(x) || return x
    n_ops = count_operation(x, *, false)+1
    ops = Array{Any}(undef, n_ops)
    @views split_term!(ops, x, [*])
    filter!(x->!isa(x, Number), ops)
    return Num(prod(ops))
end

function remove_constant_factor!(o::AbstractArray)
    for i in eachindex(o)
        o[i] = remove_constant_factor(o[i])
    end
end

function create_linear_independent_eqs(o::AbstractVector)
    o .= simplify.(o)
    remove_constant_factor!(o)
    n_ops = [count_operation(bi, +, false) for bi in o]
    n_x = sum(n_ops) + length(o)
    u_o = Array{Any}(undef, n_x)
    ind_lo, ind_up = 0, 0
    for i in eachindex(o)
        ind_lo = i > 1 ? sum(n_ops[1:i-1]) + i : 1
        ind_up = sum(n_ops[1:i]) + i
        @views split_term!(u_o[ind_lo:ind_up], o[i], [+])
    end
    remove_constant_factor!(u_o)
    unique!(u_o)
    return u_o
end


text:
################################################################################
### Helpers
################################################################################

"""
    vmuladd!(C, A, B, factor=1.0)

`C += factor * A * B` with A, B, C matrix types
"""
function vmuladd!(C::Matrix{T}, A::Matrix{T}, B::Matrix{T}, factor::T = T(1)) where {T <: Real}
    @turbo for m in 1:size(A, 1), n in 1:size(B, 2)
        Cmn = zero(eltype(C))
        for k in 1:size(A, 2)
            Cmn += A[m,k] * B[k,n]
        end
        C[m,n] += factor * Cmn
    end
end
function vmuladd!(C::Matrix{T}, A::Matrix{T}, B::Diagonal{T}, factor::T = T(1)) where {T <: Real}
    @turbo for m in 1:size(A, 1), n in 1:size(A, 2)
        C[m,n] += factor * A[m,n] * B.diag[n]
    end
end
function vmuladd!(C::Matrix{T}, A::Matrix{T}, X::Adjoint{T}, factor::T = T(1)) where {T <: Real}
    B = X.parent
    @turbo for m in 1:size(A, 1), n in 1:size(B, 2)
        Cmn = zero(eltype(C))
        for k in 1:size(A, 2)
            Cmn += A[m,k] * B[n, k]
        end
        C[m,n] += factor * Cmn
    end
end
function vmuladd!(C::Matrix{T}, X::Adjoint{T}, B::Matrix{T}, factor::T = T(1)) where {T <: Real}
    A = X.parent
    @turbo for m in 1:size(A, 1), n in 1:size(B, 2)
        Cmn = zero(eltype(C))
        for k in 1:size(A, 2)
            Cmn += A[k,m] * B[k,n]
        end
        C[m,n] += factor * Cmn
    end
end
function vmuladd!(C::Matrix{T}, X1::Adjoint{T}, X2::Adjoint{T}, factor::T = T(1)) where {T <: Real}
    A = X1.parent
    B = X2.parent
    @turbo for m in 1:size(A, 1), n in 1:size(B, 2)
        Cmn = zero(eltype(C))
        for k in 1:size(A, 2)
            Cmn += A[k,m] * B[n,k]
        end
        C[m,n] += factor * Cmn
    end
end

"""
    vmul!(C, A, B, factor)

`C = factor * A * B` with C, A, B matrices
"""
function vmul!(C::Matrix{T}, A::Matrix{T}, X::Adjoint{T}, factor::T) where {T <: Real}
    B = X.parent
    @turbo for m in 1:size(A, 1), n in 1:size(B, 2)
        Cmn = zero(eltype(C))
        for k in 1:size(A, 2)
            Cmn += A[m,k] * B[n, k]
        end
        C[m,n] = factor * Cmn
    end
end
function vmul!(C::Matrix{T}, X1::Adjoint{T}, X2::Adjoint{T}, factor::T) where {T <: Real}
    A = X1.parent
    B = X2.parent
    @turbo for m in 1:size(A, 1), n in 1:size(B, 2)
        Cmn = zero(eltype(C))
        for k in 1:size(A, 2)
            Cmn += A[k,m] * B[n,k]
        end
        C[m,n] = factor * Cmn
    end
end

################################################################################
### Overloads for complex StructArrays
################################################################################

const CMat64 = StructArray{ComplexF64, 2, NamedTuple{(:re, :im), Tuple{Matrix{Float64}, Matrix{Float64}}}, Int64}
const CVec64 = StructVector{ComplexF64, NamedTuple{(:re, :im), Tuple{Vector{Float64}, Vector{Float64}}}, Int64}

# constructors
CVec64(::UndefInitializer, n) = StructArray(Vector{ComplexF64}(undef, n))::CVec64
CMat64(::UndefInitializer, n, m) = StructArray(Matrix{ComplexF64}(undef, n, m))::CMat64
CMat64(::UniformScaling, n, m) = StructArray(Matrix{ComplexF64}(I, n, m))::CMat64

# util
function Base.copyto!(C::CMat64, ::UniformScaling)
    @turbo for i in axes(C.re, 1), j in axes(C.re, 2)
        C.re[i, j] = Float64(i == j)
    end
    @turbo for i in axes(C.re, 1), j in axes(C.re, 2)
        C.im[i, j] = 0.0
    end
    nothing
end
function Base.copyto!(C::CMat64, D::Diagonal{<: Real})
    @turbo for i in axes(C.re, 1), j in axes(C.re, 2)
        C.re[i, j] = Float64(i == j) * D.diag[i]
    end
    @turbo for i in axes(C.im, 1), j in axes(C.im, 2)
        C.im[i, j] = 0.0
    end
    nothing
end

# unoptimized - these (should) only run during initialization anyway
Base.exp(C::CMat64) = StructArray(exp(Matrix(C)))
Base.:(*)(C1::CMat64, C2::CMat64) = StructArray(Matrix(C1) * Matrix(C2))

# tested
@inline function vmul!(C::CMat64, A::CMat64, B::CMat64)
    vmul!(   C.re, A.re, B.re)       # C.re = A.re * B.re
    vmuladd!(C.re, A.im, B.im, -1.0) # C.re = C.re - A.im * B.im
    vmul!(   C.im, A.re, B.im)       # C.im = A.re * B.im
    vmuladd!(C.im, A.im, B.re)       # C.im = C.im + A.im * B.re
end

#tested
@inline function vmul!(C::CMat64, A::CMat64, B::Diagonal{<: Real})
    vmul!(C.re, A.re, B)
    vmul!(C.im, A.im, B)
end
@inline function vmul!(C::CMat64, A::CMat64, B::Diagonal{<: Real}, range)
    @views vmul!(C.re, A.re, Diagonal(B.diag[range]))
    @views vmul!(C.im, A.im, Diagonal(B.diag[range]))
end
#tested
@inline function vmul!(C::CMat64, A::CMat64, D::Diagonal{ComplexF64, <: CVec64})
    vmul!(   C.re, A.re, Diagonal(D.diag.re))
    vmuladd!(C.re, A.im, Diagonal(D.diag.im), -1.0)
    vmul!(   C.im, A.re, Diagonal(D.diag.im))
    vmuladd!(C.im, A.im, Diagonal(D.diag.re))
end
@inline function vmul!(C::CMat64, A::CMat64, D::Diagonal{ComplexF64, <: CVec64}, range)
    @views vmul!(   C.re, A.re, Diagonal(D.diag.re[range]))
    @views vmuladd!(C.re, A.im, Diagonal(D.diag.im[range]), -1.0)
    @views vmul!(   C.im, A.re, Diagonal(D.diag.im[range]))
    @views vmuladd!(C.im, A.im, Diagonal(D.diag.re[range]))
end

@inline function vmul!(C::CMat64, A::Diagonal{<: Real}, B::CMat64)
    vmul!(C.re, A, B.re)
    vmul!(C.im, A, B.im)
end
@inline function vmul!(C::CMat64, A::Diagonal{<: Real}, B::CMat64, range)
    @views vmul!(C.re, A, B.re, range)
    @views vmul!(C.im, A, B.im, range)
end
@inline function vmul!(C::CMat64, A::Diagonal{ComplexF64, <: CVec64}, B::CMat64)
    vmul!(   C.re, Diagonal(A.diag.re), B.re)
    vmuladd!(C.re, Diagonal(A.diag.im), B.im, -1.0)
    vmul!(   C.im, Diagonal(A.diag.re), B.im)
    vmuladd!(C.im, Diagonal(A.diag.im), B.re)
end
@inline function vmul!(C::CMat64, A::Diagonal{ComplexF64, <: CVec64}, B::CMat64, range)
    @views vmul!(   C.re, Diagonal(A.diag.re[range]), B.re)
    @views vmuladd!(C.re, Diagonal(A.diag.im[range]), B.im, -1.0)
    @views vmul!(   C.im, Diagonal(A.diag.re[range]), B.im)
    @views vmuladd!(C.im, Diagonal(A.diag.im[range]), B.re)
end

# tested
@inline function vmul!(C::CMat64, A::CMat64, X::Adjoint{ComplexF64, <: CMat64})
    B = X.parent
    vmul!(   C.re, A.re, adjoint(B.re))
    vmuladd!(C.re, A.im, adjoint(B.im), 1.0)
    vmul!(   C.im, A.re, adjoint(B.im), -1.0)
    vmuladd!(C.im, A.im, adjoint(B.re))
end
# tested
@inline function vmul!(C::CMat64, X::Adjoint{ComplexF64, <:CMat64}, B::CMat64)
    A = X.parent
    vmul!(   C.re, adjoint(A.re), B.re)
    vmuladd!(C.re, adjoint(A.im), B.im, 1.0)
    vmul!(   C.im, adjoint(A.re), B.im)
    vmuladd!(C.im, adjoint(A.im), B.re, -1.0)
end
# tested
@inline function vmul!(C::CMat64, X::Adjoint{ComplexF64, <: CMat64}, Y::Adjoint{ComplexF64, <: CMat64})
    A = X.parent; B = Y.parent
    vmul!(   C.re, adjoint(A.re), adjoint(B.re))
    vmuladd!(C.re, adjoint(A.im), adjoint(B.im), -1.0)
    vmul!(   C.im, adjoint(A.re), adjoint(B.im), -1.0)
    vmuladd!(C.im, adjoint(A.im), adjoint(B.re), -1.0)
end

# tested
@inline function rvmul!(A::CMat64, B::Diagonal{T}) where {T <: Real}
    rvmul!(A.re, B)
    rvmul!(A.im, B)
end
# tested
function rvmul!(A::CMat64, B::Diagonal{ComplexF64, <: CVec64})
    error("This multiplication cannot be done inplace.")
end

# tested
@inline function lvmul!(A::Diagonal{T}, B::CMat64) where {T <: Real}
    lvmul!(A, B.re)
    lvmul!(A, B.im)
end
# tested
function lvmul!(A::Diagonal{ComplexF64, <: CVec64}, B::CMat64)
    error("This multiplication cannot be done inplace.")
end

# tested
@inline rvadd!(A::CMat64, D::Diagonal{T}) where {T <: Real} = rvadd!(A.re, D)
# tested
@inline function rvadd!(A::CMat64, B::CMat64)
    rvadd!(A.re, B.re)
    rvadd!(A.im, B.im)
end

# test
@inline function vsub!(O::CMat64, A::CMat64, ::UniformScaling)
    copyto!(O.im, A.im)
    vsub!(O.re, A.re, I)
end

@bm function rdivp!(A::CMat64, T::CMat64, O::CMat64, pivot)
    # assume Diagonal is ±1!
    @inbounds begin
        N = size(A, 1)

        # Apply pivot
        for (j, p) in enumerate(pivot)
            @turbo for i in 1:N
                O.re[i, j] = A.re[i, p]
            end
            @turbo for i in 1:N
                O.im[i, j] = A.im[i, p]
            end
        end

        # do the rdiv
        # @turbo will segfault on `k in 1:0`, so pull out first loop
        # invT = conj(T[1, 1]) / abs2(T[1, 1])
        invT = 1.0 / T[1, 1]
        re = real(invT)
        im = imag(invT)

        @turbo for i in 1:N
            A.re[i, 1] = O.re[i, 1] * re
        end
        @turbo for i in 1:N
            A.re[i, 1] -= O.im[i, 1] * im
        end
        @turbo for i in 1:N
            A.im[i, 1] = O.im[i, 1] * re
        end
        @turbo for i in 1:N
            A.im[i, 1] += O.re[i, 1] * im
        end

        # TODO Is this optimal?
        for j in 2:N
            # invT = conj(T[j, j]) / abs2(T[j, j])
            invT = 1.0 / T[j, j] # TODO this might always be real?
            re = real(invT)
            im = imag(invT)

            # These loops need to be seperated. If they are not results are
            # with @turbo. I guess this is because x_re and x_im get written
            # to in parallel, i.e. values are overwritten?
            # this is worth it performance wise. 5-10x faster for 64^2 Matrix
            @turbo for i in 1:N
                x_re = O.re[i, j]
                for k in 1:j-1
                    x_re -= A.re[i, k] * T.re[k, j]
                end
                A.re[i, j] = x_re * re
                A.im[i, j] = x_re * im
            end
            @turbo for i in 1:N
                x_re = 0.0
                for k in 1:j-1
                    x_re += A.im[i, k] * T.im[k, j]
                end
                A.re[i, j] += x_re * re
                A.im[i, j] += x_re * im
            end
            @turbo for i in 1:N
                x_im = O.im[i, j]
                for k in 1:j-1
                    x_im -= A.im[i, k] * T.re[k, j]
                end
                A.re[i, j] -= x_im * im
                A.im[i, j] += x_im * re
            end
            @turbo for i in 1:N
                x_im = 0.0
                for k in 1:j-1
                    x_im -= A.re[i, k] * T.im[k, j]
                end
                A.re[i, j] -= x_im * im
                A.im[i, j] += x_im * re
            end

        end
    end
    A
end

################################################################################
### UDT
################################################################################

# tested
# annoyingly we need a temp buffer for a real and imaginary numbers along the column...
# 0 allocs in @benchmark
@inline function reflector!(x::CMat64, normu, j, n, temp::CMat64)
    # Note:
    # `temp` is a CMat64 because that will be available in udt_AVX_pivot. Could
    # also use CVec64 of length size(x, 1) or Vector{Float64}(undef, 2size(x, 1))
    @inbounds begin
        ξ1 = x[j, j]
        if iszero(normu)
            return zero(ξ1) #zero(ξ1/normu)
        end
        normu = sqrt(normu)
        ν = LinearAlgebra.copysign(normu, real(ξ1))
        ξ1 += ν
        invξ1 = 1.0 / ξ1
        x.re[j, j] = -ν
        x.im[j, j] = 0.0

        # maybe this is partial copying suboptimal? (cache)
        @turbo for i in j+1:n
            temp.re[i] = x.re[i, j]
        end
        @turbo for i in j+1:n
            temp.im[i] = x.im[i, j]
        end

        @turbo for i = j+1:n
            x.re[i, j] = temp.re[i, 1] * real(invξ1)
        end
        @turbo for i = j+1:n
            x.re[i, j] -= temp.im[i, 1] * imag(invξ1)
        end
        @turbo for i = j+1:n
            x.im[i, j] = temp.im[i, 1] * real(invξ1)
        end
        @turbo for i = j+1:n
            x.im[i, j] += temp.re[i, 1] * imag(invξ1)
        end
    end
    ξ1/ν
end

# tested, 0 allocs in @benchmark
@inline function reflectorApply!(τ::Number, A::CMat64, j, n)
    @inbounds for k = j+1:n
        # dot
        vAk_re = A.re[j, k]
        @turbo for i = j+1:n
            vAk_re += A.re[i, j] * A.re[i, k]
        end
        @turbo for i = j+1:n
            vAk_re += A.im[i, j] * A.im[i, k] # conj implied
        end

        vAk_im = A.im[j, k]
        @turbo for i = j+1:n
            vAk_im -= A.im[i, j] * A.re[i, k]
        end
        @turbo for i = j+1:n
            vAk_im += A.re[i, j] * A.im[i, k]
        end

        # conj(τ) implied
        temp = real(τ) * vAk_re + imag(τ) * vAk_im
        vAk_im = real(τ) * vAk_im - imag(τ) * vAk_re
        vAk_re = temp

        # ger
        A.re[j, k] -= vAk_re
        @turbo for i = j+1:n
            A.re[i, k] -= A.re[i, j] * vAk_re
        end
        @turbo for i = j+1:n
            A.re[i, k] += A.im[i, j] * vAk_im
        end

        A.im[j, k] -= vAk_im
        @turbo for i = j+1:n
            A.im[i, k] -= A.re[i, j] * vAk_im
        end
        @turbo for i = j+1:n
            A.im[i, k] -= A.im[i, j] * vAk_re
        end
    end

    return A
end

# tested, 0 allocs in @benchmark
@inline function indmaxcolumn(A::CMat64, j=1, n=size(A, 1))
    squared_norm = 0.0
    @turbo for k in j:n
        squared_norm += abs2(A.re[k, j])
    end
    @turbo for k in j:n
        squared_norm += abs2(A.im[k, j])
    end
    ii = j
    @inbounds for i in j+1:n
        mi = 0.0
        @turbo for k in j:n
            mi += abs2(A.re[k, i])
        end
        @turbo for k in j:n
            mi += abs2(A.im[k, i])
        end
        if abs(mi) > squared_norm
            squared_norm = mi
            ii = i
        end
    end
    return ii, squared_norm
end

@bm function udt_AVX_pivot!(
        U::CMat64,
        D::AbstractArray{Float64, 1},
        input::CMat64,
        pivot::AbstractArray{Int64, 1} = Vector(UnitRange(1:size(input, 1))),
        temp::AbstractArray{ComplexF64, 1} = Vector{ComplexF64}(undef, length(D)),
        apply_pivot::Val = Val(true)
    )
    # Assumptions:
    # - all matrices same size
    # - input can be changed (input becomes T)

    n = size(input.re, 1)
    @inbounds for i in 1:n
        pivot[i] = i
    end

    # # TODO optimize
    @inbounds for j in 1:n
        # # Find column with maximum norm in trailing submatrix
        jm, squared_max_norm = indmaxcolumn(input, j, n)

        if jm != j
            # Flip elements in pivoting vector
            tmpp = pivot[jm]
            pivot[jm] = pivot[j]
            pivot[j] = tmpp

            # Update matrix with
            @turbo for i = 1:n
                tmp = input.re[i,jm]
                input.re[i,jm] = input.re[i,j]
                input.re[i,j] = tmp
            end
            @turbo for i = 1:n
                tmp = input.im[i,jm]
                input.im[i,jm] = input.im[i,j]
                input.im[i,j] = tmp
            end
        end

        # Compute reflector of columns j
        τj = reflector!(input, squared_max_norm, j, n, U) # complex result
        temp[j] = τj

        # # Update trailing submatrix with reflector
        reflectorApply!(τj, input, j, n)
    end

    @turbo for i in axes(U.re, 1), j in axes(U.re, 2)
        U.re[i, j] = Float64(i == j)
    end
    @turbo for i in axes(U.im, 1), j in axes(U.im, 2)
        U.im[i, j] = 0.0
    end

    @inbounds begin
        U.re[n, n] -= real(temp[n])
        U.im[n, n] -= imag(temp[n])
        for k = n-1:-1:1
            for j = k:n
                vBj_re = U.re[k,j]
                @turbo for i = k+1:n
                    vBj_re += input.re[i,k] * U.re[i,j]
                end
                @turbo for i = k+1:n
                    vBj_re += input.im[i,k] * U.im[i,j]
                end

                vBj_im = U.im[k,j]
                @turbo for i = k+1:n
                    vBj_im += input.re[i,k] * U.im[i,j]
                end
                @turbo for i = k+1:n
                    vBj_im -= input.im[i,k] * U.re[i,j]
                end

                # we need the results in sync and this is not a turbo loop, so
                # temp should probably be Vector{Complex} rather than CVec64
                temp_re = real(temp[k])
                temp_im = imag(temp[k])
                re = temp_re * vBj_re - temp_im * vBj_im
                vBj_im = temp_im * vBj_re + temp_re * vBj_im
                vBj_re = re

                U.re[k,j] -= vBj_re
                U.im[k,j] -= vBj_im
                @turbo for i = k+1:n
                    U.re[i,j] -= input.re[i,k] * vBj_re
                end
                @turbo for i = k+1:n
                    U.re[i,j] += input.im[i,k] * vBj_im
                end
                @turbo for i = k+1:n
                    U.im[i,j] -= input.im[i,k] * vBj_re
                end
                @turbo for i = k+1:n
                    U.im[i,j] -= input.re[i,k] * vBj_im
                end

            end
        end
    end

    @turbo for i in 1:n
        D[i] = abs(input.re[i, i])
    end

    _apply_pivot!(input, D, pivot, apply_pivot)

    nothing
end

# function invert_pivot!(pivot::Vector{Int64})
#     from = 1
#     to = pivot[1]
#     @inbounds for _ in 1:16
#         temp = pivot[to]
#         pivot[to] = from
#         from = to
#         to = temp
#     end
#     pivot
# end

@bm function _apply_pivot!(input::CMat64, D, pivot, ::Val{true})
    n = size(input, 1)

    # without applying pivot we have a upper triangular matrix and we know
    # input.re[:, 1] = [sign(input.re[1, 1]), 0, ..., 0]
    # input.im[:, 1] = [input.im[1, 1] / D[1], 0, ..., 0]
    first_re = sign(input.re[1, 1])
    first_im = input.im[1, 1] / D[1]

    # use input.im[:, 1] as temporary storage (good cache)
    # I don't think we can do swaps like A[:, i] = A[:, j] because pivot can
    # have loops...
    @inbounds for i in 1:n
        d = 1.0 / D[i]
        for j in 1:i-1 # can be 0-length which makes @turbo crash
            input.im[pivot[j], 1] = 0.0
        end
        @turbo for j in i:n
            input.im[pivot[j], 1] = d * input.re[i, j]
        end
        @turbo for j in 1:n
            input.re[i, j] = input.im[j, 1]
        end
    end

    # Now we use the one element column from input.re as a temp and repeat
    one_column = pivot[1]
    @inbounds for i in 1:n
        d = 1.0 / D[i]
        for j in 1:i-1
            input.re[pivot[j], one_column] = 0.0
        end
        @turbo for j in i:n
            input.re[pivot[j], one_column] = d * input.im[i, j]
        end
        @turbo for j in 1:n
            input.im[i, j] = input.re[j, one_column]
        end
    end

    # Finally we restore the messed up columns
    input.re[1, one_column] = first_re
    @turbo for i in 2:n
        input.re[i, one_column] = 0.0
    end
    input.im[1, one_column] = first_im
    @turbo for i in 2:n
        input.im[i, one_column] = 0.0
    end

    input
end
@bm function _apply_pivot!(input::CMat64, D, pivot, ::Val{false})
    n = size(input, 1)
    @inbounds for i in 1:n
        d = 1.0 / D[i]
        @turbo for j in i:n
            input.re[i, j] = d * input.re[i, j]
        end
    end
    @inbounds for i in 1:n
        d = 1.0 / D[i]
        @turbo for j in i:n
            input.im[i, j] = d * input.im[i, j]
        end
    end
end

################################################################################
### GHQ
################################################################################

@inline function vmul!(C::CMat64, A::Matrix{<: Float64}, D::Diagonal{ComplexF64, <: CVec64})
    vmul!(C.re, A, Diagonal(D.diag.re))
    vmul!(C.im, A, Diagonal(D.diag.im))
end
@inline function vmul!(C::CMat64, A::Matrix{<: Float64}, D::Diagonal{ComplexF64, <: CVec64}, range)
    @views vmul!(C.re, A, Diagonal(D.diag.re[range]))
    @views vmul!(C.im, A, Diagonal(D.diag.im[range]))
end
@inline function vmul!(C::CMat64, D::Diagonal{ComplexF64, <: CVec64}, B::Matrix{<: Float64})
    vmul!(C.re, Diagonal(D.diag.re), B)
    vmul!(C.im, Diagonal(D.diag.im), B)
end
@inline function vmul!(C::CMat64, D::Diagonal{ComplexF64, <: CVec64}, B::Matrix{<: Float64}, range)
    @views vmul!(C.re, Diagonal(D.diag.re[range]), B)
    @views vmul!(C.im, Diagonal(D.diag.im[range]), B)
end

@inline function vmul!(C::CMat64, A::CMat64, B::Matrix{Float64})
    vmul!(C.re, A.re, B)
    vmul!(C.im, A.im, B)
end
@inline function vmul!(C::CMat64, A::Matrix{Float64}, B::CMat64)
    vmul!(C.re, A, B.re)
    vmul!(C.im, A, B.im)
end

ext: jl


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiNetworkingV1beta1HTTPIngressRuleValue <: SwaggerModel
    paths::Any # spec type: Union{ Nothing, Vector{IoK8sApiNetworkingV1beta1HTTPIngressPath} } # spec name: paths

    function IoK8sApiNetworkingV1beta1HTTPIngressRuleValue(;paths=nothing)
        o = new()
        validate_property(IoK8sApiNetworkingV1beta1HTTPIngressRuleValue, Symbol("paths"), paths)
        setfield!(o, Symbol("paths"), paths)
        o
    end
end # type IoK8sApiNetworkingV1beta1HTTPIngressRuleValue

const _property_map_IoK8sApiNetworkingV1beta1HTTPIngressRuleValue = Dict{Symbol,Symbol}(Symbol("paths")=>Symbol("paths"))
const _property_types_IoK8sApiNetworkingV1beta1HTTPIngressRuleValue = Dict{Symbol,String}(Symbol("paths")=>"Vector{IoK8sApiNetworkingV1beta1HTTPIngressPath}")
Base.propertynames(::Type{ IoK8sApiNetworkingV1beta1HTTPIngressRuleValue }) = collect(keys(_property_map_IoK8sApiNetworkingV1beta1HTTPIngressRuleValue))
Swagger.property_type(::Type{ IoK8sApiNetworkingV1beta1HTTPIngressRuleValue }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiNetworkingV1beta1HTTPIngressRuleValue[name]))}
Swagger.field_name(::Type{ IoK8sApiNetworkingV1beta1HTTPIngressRuleValue }, property_name::Symbol) =  _property_map_IoK8sApiNetworkingV1beta1HTTPIngressRuleValue[property_name]

function check_required(o::IoK8sApiNetworkingV1beta1HTTPIngressRuleValue)
    (getproperty(o, Symbol("paths")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiNetworkingV1beta1HTTPIngressRuleValue }, name::Symbol, val)
end


text:
@ghdef mutable struct CheckSuite
    id::Union{Integer,Nothing}
    head_sha::Union{String,Nothing}
    status::Union{String,Nothing}
    conclusion::Union{String, Nothing}
    app::Union{App,Nothing}
end
namefield(suite::CheckSuite) = suite.id

"""
    GitHub.check_suite([api,] repo::Repo, id; options...)

- https://developer.github.com/v3/checks/suites/#get-a-single-check-suite
"""
@api_default function check_suite(api::GitHubAPI, repo::Repo, id; headers = Dict(), options...)
    headers["Accept"] = "application/vnd.github.antiope-preview+json"
    result = gh_get_json(api, "/repos/$(name(repo))/check-suites/$(id)";
        headers=headers, options...)
    CheckSuite(results)
end
"""
    GitHub.check_suites([api,] repo::Repo, ref; options...)

- https://developer.github.com/v3/checks/suites/#list-check-suites-for-a-specific-ref
"""
@api_default function check_suites(api::GitHubAPI, repo::Repo, ref; headers = Dict(), options...)
    headers["Accept"] = "application/vnd.github.antiope-preview+json"
    results, page_data = gh_get_paged_json(api, "/repos/$(name(repo))/commits/$(name(ref))/check-suites";
        headers=headers, options...)
    map(CheckSuite, results["check_suites"]), page_data, results["total_count"]
end


text:
using Omega
using Test

struct MyType
  x
end

lift!(:MyType; n = 1, mod = Main)
x = MyType(uniform(0.0, 1.0))
@test x isa Omega.RandVar

struct YourType{T}
  x::T
end

lift!(:YourType, n=1, mod=Main)

x = YourType(uniform(0.0, 1.0))
rand(x)


text:
# # Tuple id type
const TupleID = NTuple{N, Int} where N

singletonid(::TupleID, i::Int) = (i,)
@inline append(a::Tuple, b::Tuple) = (a..., b...)

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
# Tests for the PWDict type.
# Author: Noah C. Benson <n@nben.net>

@testset "PWSet" begin
    els = [:b => 20.0, :c => 30.0, :a => 10.0, :d => 40.0]
    p = PWSet{Symbol,Float64}()
    for el in els
        p = push(p, el)
    end
    (lst,mst) = (first(p), pop(p))
    @test lst == :d
    (lst,mst) = (first(mst), pop(mst))
    @test lst == :c
    (lst,mst) = (first(mst), pop(mst))
    @test lst == :b
    (lst,mst) = (first(mst), pop(mst))
    @test lst == :a
    @test isempty(mst)
    # make many samples and make sure they resemble the distribution
    counts = Dict(:a => 0.0, :b => 0.0, :c => 0.0, :d => 0.0)
    for _ in 1:100000
        sym = rand(p)
        counts[sym] += 1/1000
    end
    @test  8.5 < counts[:a] < 12.5
    @test 18.5 < counts[:b] < 22.5
    @test 28.5 < counts[:c] < 32.5
    @test 38.5 < counts[:d] < 42.5
end


text:
using GtkInteract
using Base.Test

## test manipulate
@manipulate for n=1:10
    n
end

## tests widgets
opts = ["one", "two", "three"]
# write your own tests here
sl = slider(1:10)         # one from range
cb = checkbox(true, label="check")      # bool
tb = togglebutton(true, label="toggle") # bool
dd = dropdown(opts)       # 1 of n
rbs = radiobuttons(opts)  # 1 of n
sel = selectlist(opts, label="select")        # 1 of n
tbs = togglebuttons(opts) # 1 of n
bg = buttongroup(opts)    # 0,1,...,n of n
txtb = textbox("text goes here", label="textbox")

btn = button("button")
out = label()

controls = [sl, cb, tb, dd, rbs,  tbs, bg, sel, txtb]

w = window(vbox(controls...),
           hbox(halign(:end,padding(10, btn))),
           out)

## using Reactive
Reactive.foreach([x.signal for x in controls]...) do sl, cb, tb, dd, rbs, tbs, bg, sel, txtb
    push!(out, """
          $(string(sl))
          $(string(cb))
          $(string(tb))
          $(string(dd))
          $(string(rbs))
          $(string(tbs))
          $(string(bg))
          $(string(sel))
          $(string(txtb))

          """)
end

# test numeric-valued textboxes
using Gtk, Interact
w = mainwindow()
b = textbox(1)
push!(w, b)
display(w)
g = first(w.window)
widget = first(g)
setproperty!(widget, :text, "2")
@test value(signal(b)) == 1
signal_emit(widget, :activate, Void)
sleep(0.1)
@test value(signal(b)) == 2

# Images
using TestImages, Colors
w = mainwindow()
img = testimage("mandrill")
renderbuffer = convert(Matrix{RGB24}, img.data)
surf = cairoimagesurface(renderbuffer)
# Gtk.present(w.window)
push!(w, grow(surf))
display(w)
img = testimage("cameraman")
sleep(1)
push!(surf, img)

# Cleanup of signals (especially ones that run constantly!)
using Reactive
frametimer = fps(10)
w = mainwindow()
push!(w.refs, frametimer)
display(w)
@test frametimer.alive == true
destroy(w)
yield()
@test frametimer.alive == false


text:
using Convex
using Mosek, MosekTools
using Optim

using Transducers
using BenchmarkTools

function f(x, u)
    0.5 * (-x'*x + sum(dot(^)(u,2)))
end

function convex_solver(x, u_guess, ulim)
    m = length(u_guess)
    u = Convex.Variable(m)
    problem = minimize(f(x, u))
    problem.constraints += [u >= ulim[1]]
    problem.constraints += [u <= ulim[2]]
    solve!(problem, Mosek.Optimizer(); silent_solver=true)
    u.value, problem.optval
end

function ipnewton_solver(x, u_guess, ulim)
    dfc = TwiceDifferentiableConstraints(ulim[1], ulim[2])
    obj(u) = f(x, u)
    res = Optim.optimize(obj, dfc, u_guess, IPNewton())
    res.minimizer, res.minimum
end

@testset "solver" begin
    for N in [1, 10, 100]
        @show n, m = N, N
        d = 100
        xs = 1:d |> Map(i -> rand(n)) |> collect
        us = 1:d |> Map(i -> rand(m)) |> collect
        ulim = (-2*ones(m), 2*ones(m))
        # case 1: convex solver
        println("convex solver")
        @btime convex_solver(rand($n), rand($m), $ulim)
        println("non-convex solver (ipnewton)")
        @btime ipnewton_solver(rand($n), rand($m), $ulim)
    end
    nothing
end


text:
function test_convolution_layer(backend::Backend, n_group, filter_w, filter_h, pad_w, pad_h, stride_w, stride_h, n_input, freeze, T, eps)
  println("-- Testing Convolution(frozen=$freeze) on $(typeof(backend)){$T} filter=$((filter_w,filter_h))...")
  println("    > Setup")
  input_w = 16
  input_h = 10
  input_chann = 6
  input_num = 24
  n_filter = 12

  input_dims = (input_w, input_h, input_chann, input_num)
  filter_dims = (filter_w, filter_h, round(Int, input_chann/n_group), n_filter)
  bias_dims = (1, 1, n_filter, 1)

  layer = ConvolutionLayer(name="conv", kernel=(filter_w, filter_h), stride=(stride_w, stride_h),
      pad=(pad_w,pad_h), n_filter=n_filter, n_group=n_group,
      tops=Array(Symbol,n_input), bottoms=Array(Symbol,n_input))

  # convolution layer requires each input blob to be the same shape
  input = [rand(T, input_dims) for i = 1:n_input]
  inputs = Blob[make_blob(backend, x) for x in input]
  data_diffs = Blob[make_blob(backend, x) for x in input]

  state = setup(backend, layer, inputs, data_diffs)

  println("    > Forward")
  filter = rand(T, filter_dims)
  copy!(state.filter, filter)
  bias = rand(T, bias_dims)
  copy!(state.bias, bias)

  forward(backend, state, inputs)
  for i = 1:n_input
    expected_output = convolution_forward(state, filter, bias, input[i])
    @test size(expected_output) == size(state.blobs[i])
    got_output = similar(expected_output)
    copy!(got_output, state.blobs[i])
    @test all(-eps .< expected_output - got_output .< eps)
  end

  println("    > Backward")
  top_diff = [rand(T, size(state.blobs[1])) for i = 1:n_input]
  for i = 1:n_input
    copy!(state.blobs_diff[i], top_diff[i])
  end

  if freeze
    freeze!(state)
  end
  backward(backend, state, inputs, data_diffs)

  grad_filter_exp = zeros(T, filter_dims)
  grad_bias_exp = zeros(T, bias_dims)
  for i = 1:n_input
    gradients_expected = convolution_backward(state, filter, bias, input[i], top_diff[i])
    gradients_got = similar(gradients_expected[3])
    copy!(gradients_got, data_diffs[i])

    @test all(-eps .< gradients_got - gradients_expected[3] .< eps)

    grad_filter_exp += gradients_expected[1]
    grad_bias_exp += gradients_expected[2]
  end
  grad_filter_got = similar(grad_filter_exp)
  grad_bias_got = similar(grad_bias_exp)
  copy!(grad_filter_got, state.∇filter)
  copy!(grad_bias_got, state.∇bias)

  #is_grad_filter_match = all(abs(grad_filter_exp - grad_filter_got) .< eps)
  #is_grad_bias_match   = all(abs(grad_bias_exp - grad_bias_got) .< eps)

  is_grad_filter_match = all(map((x,y)->isapprox(x,y; atol=eps), grad_filter_exp, grad_filter_got))
  is_grad_bias_match   = all(map((x,y)->isapprox(x,y; atol=eps), grad_bias_exp, grad_bias_got))

  if freeze
    # when frozen, the gradients are not computed, so the got value
    # should be random un-initialized values, so should not match (with
    # very high probability)
    @test !is_grad_bias_match
    @test !is_grad_filter_match
  else
    @show maximum(abs(grad_filter_exp - grad_filter_got))
    @show eps
    @test is_grad_filter_match
    @test is_grad_bias_match
  end

  shutdown(backend, state)
end

# naive implementation of convolution forward, used to check the correctness
function convolution_forward(state, filter::Array, bias::Array, input::Array)
  width, height, channel, num = size(input)
  n_group = state.layer.n_group
  o_g = round(Int, state.layer.n_filter / n_group)
  k_g = round(Int, channel / n_group)

  output = Array(eltype(input), size(state.blobs[1]))
  output[:] = 0

  for n = 1:num
    for g = 1:n_group
      for o = 1:o_g
        for k = 1:k_g
          for y = 1:state.height_out
            for x = 1:state.width_out
              for p = 1:state.layer.kernel[2]
                for q = 1:state.layer.kernel[1]
                  in_y = (y-1) * state.layer.stride[2] - state.layer.pad[2] + p
                  in_x = (x-1) * state.layer.stride[1] - state.layer.pad[1] + q
                  if (in_y >= 1 && in_y <= height && in_x >= 1 && in_x <= width)
                    output[x, y, (g-1)*o_g+o, n] += input[in_x, in_y, (g-1)*k_g+k, n] *
                        filter[q, p, k, (g-1)*o_g+o]
                  end
                end
              end
            end
          end
        end
      end
    end
  end

  # add bias
  for n = 1:num
    for o = 1:state.layer.n_filter
      for y = 1:state.height_out
        for x = 1:state.width_out
          output[x, y, o, n] += bias[o]
        end
      end
    end
  end

  return output
end

# naive implementation of convolution backward, used to check correctness
function convolution_backward(state, filter::Array, bias::Array, input::Array, top_diff::Array)
  ∇filter = zeros(eltype(filter), size(filter))
  ∇bias   = zeros(eltype(bias), size(bias))
  ∇input  = zeros(eltype(input), size(input))

  width, height, channels, num = size(input)
  n_group = state.layer.n_group
  o_g = round(Int, state.layer.n_filter / n_group)
  k_g = round(Int, channels / n_group)

  # ∇bias
  for n = 1:num
    for o = 1:state.layer.n_filter
      for y = 1:state.height_out
        for x = 1:state.width_out
          ∇bias[o] += top_diff[x, y, o, n]
        end
      end
    end
  end

  # ∇filter and ∇input
  for n = 1:num
    for g = 1:n_group
      for o = 1:o_g
        for k = 1:k_g
          for y = 1:state.height_out
            for x = 1:state.width_out
              for p = 1:state.layer.kernel[2]
                for q = 1:state.layer.kernel[1]
                  in_y = (y-1) * state.layer.stride[2] - state.layer.pad[2] + p
                  in_x = (x-1) * state.layer.stride[1] - state.layer.pad[1] + q
                  if (in_y >= 1 && in_y <= height && in_x >= 1 && in_x <= width)
                    ∇filter[q,p,k,(g-1)*o_g+o] += top_diff[x,y,(g-1)*o_g+o,n] * input[in_x,in_y,(g-1)*k_g+k,n]
                    ∇input[in_x,in_y,(g-1)*k_g+k,n] += top_diff[x,y,(g-1)*o_g+o,n] * filter[q,p,k,(g-1)*o_g+o]
                  end
                end
              end
            end
          end
        end
      end
    end
  end

  return (∇filter, ∇bias, ∇input)
end

function test_convolution_layer(backend::Backend, n_group, filter_w, filter_h, pad_w, pad_h, stride_w, stride_h, n_input, T, eps)
  test_convolution_layer(backend, n_group, filter_w, filter_h, pad_w, pad_h, stride_w, stride_h, n_input, true, T, eps)
  test_convolution_layer(backend, n_group, filter_w, filter_h, pad_w, pad_h, stride_w, stride_h, n_input, false, T, eps)
end
function test_convolution_layer(backend::Backend, n_input, T, eps)
  test_convolution_layer(backend, 2, 3, 4, 2, 2, 1, 2, n_input, T, eps)
  #XXX: temporarily disable this test for GPU of a bug in cuDNN v2 RC2
  # see https://devtalk.nvidia.com/default/topic/806940/gpu-accelerated-libraries/cudnn-1x1-filter-convolutions
  #if !test_gpu
    test_convolution_layer(backend, 1, 1, 1, 0, 0, 1, 1, n_input, T, eps)
  #end
end

function test_convolution_layer(backend::Backend)
  test_convolution_layer(backend, 3, Float64, 1e-5)
  test_convolution_layer(backend, 3, Float32, 1e-2) # Float32 is sooo inaccurate?
end
if test_cpu
  test_convolution_layer(backend_cpu)
end
if test_gpu
  test_convolution_layer(backend_gpu)
end


text:
push!(LOAD_PATH, string(pwd(),"/src"))

using ModifiedStokesSolver
using Base.Test

@testset "NumDiff" begin

    function numdiv(F, r, h)
        hx = [h, 0]
        hy = [0, h]
        Fx = (F(r.+hx) .- F(r.-hx)) ./ (2*h)
        Fy = (F(r.+hy) .- F(r.-hy)) ./ (2*h)
        return Fx[1] + Fy[2]
    end

    function numgrad(F, r, h)
        hx = [h, 0]
        hy = [0, h]
        Fx = (F(r.+hx) - F(r.-hx)) ./ (2*h)
        Fy = (F(r.+hy) - F(r.-hy)) ./ (2*h)
        return (Fx, Fy)
    end

    function numlap(F, r, h)
        hx = [h, 0]
        hy = [0, h]
        Fxx = (F(r.+hx) .-2.*F(r) .+ F(r.-hx)) ./ (h^2)
        Fyy = (F(r.+hy) .-2.*F(r) .+ F(r.-hy)) ./ (h^2)
        return Fxx .+ Fyy
    end

    h = 1e-4
    srand(1)
    r = 1 .- 2.*rand(2)
    f = 1 .- 2.*rand(2)
    n = 1 .- 2.*rand(2)
    a = rand()
    rnorm = sqrt(sum(r.^2))
    println("a|r| = ", a*rnorm)

    p = pressure(r, f, a)
    u = stokeslet(r, f, a)
    divu = numdiv(r -> stokeslet(r, f, a), r, h)
    Lp = numlap(r -> pressure(r, f, a), r, h)
    Lu = numlap(r -> stokeslet(r, f, a), r, h)
    gradp = numgrad(r -> pressure(r, f, a), r, h)
    pde = a^2.*u .- Lu .+ gradp

    println("u = ",u)
    println("div u = ",divu)
    println("L p = ", Lp)
    println("a^2 u - Lu + gradp = ", pde)

    # Test streslet with a|r| <1 and >1
    T = stresslet(r, f, n, a)
    divT_zlt1 = numdiv(r -> stresslet(r, f, n, a), r, h)
    println("T = ", T)
    println("div T = ", divT_zlt1)

    a = 2/rnorm
    println("a|r| = ", a*rnorm)
    T = stresslet(r, f, n, a)
    divT_zgt1 = numdiv(r -> stresslet(r, f, n, a), r, h)
    println("T = ", T)
    println("div T = ", divT_zgt1)

    @test abs(divu) < 1e-7
    @test abs(Lp) < 1e-7
    @test maximum(abs.(pde)) < 1e-7
    @test abs(divT_zlt1) < 1e-7
    @test abs(divT_zgt1) < 1e-7

end


text:
using PyPlot

abstract type Render_Option end

mutable struct Render_Airspace <: Render_Option
	pause_time::Float64
	if_annotate::Bool
	function Render_Airspace(; pause_time::Float64=0.1, annotate::Bool=false)
		show()
		this = new()
		this.pause_time = pause_time
		this.if_annotate = annotate
		return this
	end
end

struct Not_Render <: Render_Option
end

function render_airspace(::Any, render_option::Not_Render, ::Any)
	nothing
end

function render_airspace(t::Int64,
				render_option::Render_Airspace,
				All_AC::Array{Aircraft, 1})

	# fig = figure("Airspace", figsize = (7, 5))

	for ac in All_AC
		if  ac.desti == false  # && ac.in_airspace == true

			if ac.advisory == :COC
				scatter(ac.dyn.x, ac.dyn.y,
					marker = "o", color = "blue", s = 8)
			else
				scatter(ac.dyn.x, ac.dyn.y,
					marker = "o", color = "red", s = 8)
			end

			th = linspace(-pi, pi, 30)
			# print(ac.sensor.sensing_range)
			if ac.advisory == :COC
				plot(ac.dyn.x + ac.sensor.sensing_range * cos.(th),
					ac.dyn.y + ac.sensor.sensing_range * sin.(th),
					linestyle = "--", color = "green",
					linewidth = 0.4)
			else
				plot(ac.dyn.x + ac.sensor.sensing_range * cos.(th),
					ac.dyn.y + ac.sensor.sensing_range * sin.(th),
					linestyle = "--", color = "red",
					linewidth = 0.4)
			end

			arrow_len = airspace_dim / 12

			arrow(
				ac.dyn.x,
				ac.dyn.y,
				cos(ac.dyn.heading) * arrow_len * 0.8, # dx
				sin(ac.dyn.heading) * arrow_len * 0.8, # dy
				head_width = airspace_dim / 70,
				width = 1,
				head_length = airspace_dim / 70,
				overhang = 0.5,
				head_starts_at_zero = "true",
				facecolor = "black",
				length_includes_head = "true")

			scatter(ac.x_dest, ac.y_dest,
					marker = ",", color = "magenta", s = 8)

			plot([ac.dyn.x; ac.x_dest], [ac.dyn.y; ac.y_dest],
				linestyle = "--", color = "black", linewidth = 1)

			if render_option.if_annotate
				annotate(string(ac.advisory),
					xy = [ac.dyn.x; ac.dyn.y],# Arrow tip
					xytext = [ac.dyn.x + arrow_len/5; ac.dyn.y + arrow_len/5 * 2],
					xycoords = "data")

				turn_rate_str = string(rad2deg(ac.dyn.turn_rate))
				if ac.dyn.turn_rate != 0
					turn_rate_str = turn_rate_str[1:4]
				end

				annotate(turn_rate_str,
					xy = [ac.dyn.x; ac.dyn.y],# Arrow tip
					xytext = [ac.dyn.x + arrow_len/5; ac.dyn.y],
					xycoords = "data")

				annotate(string(round(rad2deg(ac.dyn.heading))),
					xy = [ac.dyn.x; ac.dyn.y],# Arrow tip
					xytext = [ac.dyn.x + arrow_len/5; ac.dyn.y - arrow_len/5 * 2],
					xycoords = "data")
			end
		end
	end

	xlabel("X (ft)")
	ylabel("Y (ft)")
	title(string(t) * ": " * string(length(All_AC)))
	axis("equal")
	ax = gca()
	ax[:set_xlim]([0, airspace_dim])
	ax[:set_ylim]([0, airspace_dim])

	pause(render_option.pause_time)
	clf()
end

# function CAS_plot(t::Float64,
# 				plot_option::Plot_Evaluation_Results,
# 				stats::Stats)
# 	title = "Number of Aircraft at Each Time " * string(t)
# 	figure(title)
# 	plot(stats.t, stats.num_ac, linestyle = "-", color = "b", label = "Number of Aircraft")
# 	xlabel("Simulation Time")
# 	ylabel("Number of Aircraft at t")
# 	savefig("../Simulation_Results/eval_plots/" * title * ".pdf")

# 	title = "Cumulative Number of Aircraft Arrived" * string(t)
# 	figure(title)
# 	plot(stats.t, stats.num_ac_arrived, linestyle = "-", color = "r", label = "Number of Aircraft Arrived")
# 	xlabel("Simulation Time")
# 	ylabel("Number of Aircraft Arrived")
# 	savefig("../Simulation_Results/eval_plots/" * title * ".pdf")

# 	title = "Number of NMAC at Each Step " * string(t)
# 	figure(title)
# 	plot(stats.t, stats.num_NMAC_t, linestyle = "-", color = "b", linewidth = 0.3)
# 	xlabel("Simulation Time")
# 	ylabel("Number of NMAC at Each Step")
# 	savefig("../Simulation_Results/eval_plots/" * title * ".pdf")

# 	title = "Cumulative Number of NMAC " * string(t)
# 	figure(title)
# 	plot(stats.t, stats.num_NMAC_cum, linestyle = "-", color = "r", label = "Cumulative Number of NMAC")
# 	xlabel("Simulation Time")
# 	ylabel("Cumulative Number of NMAC")
# 	savefig("../Simulation_Results/eval_plots/" * title * ".pdf")

# 	title = "Average Flight Time " * string(t)
# 	figure(title)
# 	plot(stats.t, stats.average_flight_time, linestyle = "-", color = "b")
# 	xlabel("Simulation Time")
# 	ylabel("Average Flight Time")
# 	savefig("../Simulation_Results/eval_plots/" * title * ".pdf")

# 	show()
# end

max_stars_repo_name: parachutel/DeepCorrectionCAS


text:
"""
    AssemblyModule

Module for assemblers  of system matrices and vectors.
"""
module AssemblyModule

__precompile__(true)

using FinEtools
import FinEtools.AssemblyModule: startassembly!, assemble!, makematrix!
using SparseMatricesCSR
using  LinearAlgebra: diag

"""
    SysmatAssemblerSparseCSRSymm{T<:Number} <: AbstractSysmatAssembler

Assembler for a **symmetric square** matrix  assembled from symmetric square
matrices. The matrix is stored in the Compressed Sparse Row (CSR) format.
"""
mutable struct SysmatAssemblerSparseCSRSymm{T<:Number} <: AbstractSysmatAssembler
    # Type for assembling of a sparse global matrix from elementwise matrices.
    buffer_length:: FInt;
    matbuffer::Vector{T};
    rowbuffer::Vector{FInt};
    colbuffer::Vector{FInt};
    buffer_pointer:: FInt;
    ndofs:: FInt;
end

"""
    SysmatAssemblerSparseCSRSymm(zero::T=0.0) where {T<:Number}

Construct blank system matrix assembler for symmetric CSR matrices. The matrix
entries are of type `T`.

# Example

This is how a symmetric sparse matrix is assembled from two square dense matrices.
```
	a = SysmatAssemblerSparseCSRSymm(0.0)
	startassembly!(a, 5, 5, 3, 7, 7)
	m = [0.24406   0.599773    0.833404  0.0420141
		0.786024  0.00206713  0.995379  0.780298
		0.845816  0.198459    0.355149  0.224996]
	assemble!(a, m'*m, [5 2 1 4], [5 2 1 4])
	m = [0.146618  0.53471   0.614342    0.737833
		 0.479719  0.41354   0.00760941  0.836455
		 0.254868  0.476189  0.460794    0.00919633
		 0.159064  0.261821  0.317078    0.77646
		 0.643538  0.429817  0.59788     0.958909]
	assemble!(a, m'*m, [2 3 1 5], [2 3 1 5])
	A = makematrix!(a)
```
# See also
[]
"""
function SysmatAssemblerSparseCSRSymm(zero::T=0.0) where {T<:Number}
    return SysmatAssemblerSparseCSRSymm{T}(0,[zero],[0],[0],0,0)
end

"""
    startassembly!(self::SysmatAssemblerSparseCSRSymm{T},
      elem_mat_dim::FInt, ignore1::FInt, elem_mat_nmatrices::FInt,
      ndofs::FInt, ignore2::FInt) where {T<:Number}

Start the assembly of a symmetric square CSR global matrix.

The method makes buffers for matrix assembly. It must be called before
the first call to the method `assemble!`.
- `elem_mat_nrows`= number of rows in typical element matrix,
- `elem_mat_ncols`= number of columns in a typical element matrix,
- `elem_mat_nmatrices`= number of element matrices,
- `ndofs_row`= Total number of equations in the row direction,
- `ndofs_col`= Total number of equations in the column direction.

The values stored in the buffers are initially undefined!
"""
function startassembly!(self::SysmatAssemblerSparseCSRSymm{T}, elem_mat_dim::FInt, ignore1::FInt, elem_mat_nmatrices::FInt, ndofs::FInt, ignore2::FInt) where {T<:Number}
    self.buffer_length = elem_mat_nmatrices*elem_mat_dim^2;
    self.rowbuffer = Array{FInt, 1}(undef, self.buffer_length);
    self.colbuffer = Array{FInt, 1}(undef, self.buffer_length);
    self.matbuffer = Array{T, 1}(undef, self.buffer_length);
    self.buffer_pointer = 1;
    self.ndofs = ndofs;
    return self
end

"""
    assemble!(self::SysmatAssemblerSparseCSRSymm{T}, mat::FMat{T},  dofnums::FIntVec, ignore::FIntVec) where {T<:Number}

Assemble a square symmetric matrix.

`dofnums` are the row degree of freedom numbers, the column degree of freedom
number input is ignored (the row and column numbers are assumed to be the same).
"""
function assemble!(self::SysmatAssemblerSparseCSRSymm{T}, mat::FMat{T},  dofnums::Union{FIntVec, FIntMat}, ignore::Union{FIntVec, FIntMat}) where {T<:Number}
    # Assembly of a square symmetric matrix.
    # The method assembles the square symmetric matrix using the two vectors of
    # equation numbers for the rows and columns.
    nrows=length(dofnums); ncolumns=nrows;
    p = self.buffer_pointer
    @assert p+ncolumns*nrows <= self.buffer_length+1
    @assert size(mat) == (nrows, ncolumns)
    @inbounds for j in 1:ncolumns
        @inbounds for i in 1:nrows
            self.matbuffer[p] = mat[i,j] # serialized matrix
            self.rowbuffer[p] = dofnums[i];
            self.colbuffer[p] = dofnums[j];
            p=p+1
        end
    end
    self.buffer_pointer=p;
    return self
end

"""
    makematrix!(self::SysmatAssemblerSparseCSRSymm)

Make a sparse symmetric square matrix.
"""
function makematrix!(self::SysmatAssemblerSparseCSRSymm)
    # Make a sparse matrix.
    # The method makes a sparse matrix from the assembly buffers.
    @assert length(self.rowbuffer) >= self.buffer_pointer-1
    @assert length(self.colbuffer) >= self.buffer_pointer-1
    # Here we will go through the rows and columns, and whenever the row or
    # the column refer to indexes outside of the limits of the matrix, the
    # corresponding value will be set to 0 and assembled to row and column 1.
    @inbounds for j=1:self.buffer_pointer-1
        r = self.rowbuffer[j]
        if (r > self.ndofs) || (r <= 0)
            self.rowbuffer[j] = 1;
            self.matbuffer[j] = 0.0
        end
        c = self.colbuffer[j]
        if (c > self.ndofs) || (c <= 0)
            self.colbuffer[j] = 1;
            self.matbuffer[j] = 0.0
        end
    end
    # Use the CSR package.
    S = SparseMatricesCSR.sparsecsr(self.rowbuffer[1:self.buffer_pointer-1],
               self.colbuffer[1:self.buffer_pointer-1],
               self.matbuffer[1:self.buffer_pointer-1],
               self.ndofs, self.ndofs);
    self = SysmatAssemblerSparseCSRSymm(zero(eltype(self.matbuffer))) # get rid of the buffers
    return S
end

"""
    SysmatAssemblerVecOfMatSymm{T<:Number} <: AbstractSysmatAssembler

Assembler for a vector of matrices.

The global matrix is represented by a collection of constituent matrices, which
in this case are not assembled, but rather kept in a vector. The global matrix
is assumed to be symmetric.
"""
mutable struct SysmatAssemblerVecOfMatSymm{T<:Number} <: AbstractSysmatAssembler
    # Type for assembling of a sparse global matrix from elementwise matrices.
    elem_mat_dim::FInt
    buffer_length::FInt;
    buffer::Vector{Tuple{Vector{FInt}, Matrix{T}}};
    buffer_pointer::FInt;
    ndofs::FInt;
end

"""
    SysmatAssemblerVecOfMatSymm(zero::T=0.0) where {T<:Number}

Construct blank system matrix assembler for vector of matrices. The matrix
entries have elements of type `T`.

# Example

This is how a symmetric sparse matrix is assembled from two square dense matrices.
```
    a = SysmatAssemblerVecOfMatSymm(0.0)
    startassembly!(a, 5, 5, 3, 7, 7)
    m = [0.24406   0.599773    0.833404  0.0420141
        0.786024  0.00206713  0.995379  0.780298
        0.845816  0.198459    0.355149  0.224996]
    assemble!(a, m'*m, [5 2 1 4], [5 2 1 4])
    m = [0.146618  0.53471   0.614342    0.737833
         0.479719  0.41354   0.00760941  0.836455
         0.254868  0.476189  0.460794    0.00919633
         0.159064  0.261821  0.317078    0.77646
         0.643538  0.429817  0.59788     0.958909]
    assemble!(a, m'*m, [2 3 1 5], [2 3 1 5])
    A = makematrix!(a)
```
# See also
[]
"""
function SysmatAssemblerVecOfMatSymm(zero::T=0.0) where {T<:Number}
    return SysmatAssemblerVecOfMatSymm{T}(0,0,[],0,0)
end

"""
    startassembly!(self::SysmatAssemblerVecOfMatSymm{T}, ignore0::FInt, ignore1::FInt, elem_mat_nmatrices::FInt, ndofs::FInt, ignore2::FInt) where {T<:Number}

Start the assembly of a symmetric square CSR global matrix.

The method makes buffers for matrix assembly. It must be called before
the first call to the method `assemble!`.
- `elem_mat_dim`= dimension of the element matrix,
- `ignore1`= ignored,
- `elem_mat_nmatrices`= number of element matrices,
- `ndofs`= Total number of equations in the row and column direction.

The values stored in the buffers are initially undefined!
"""
function startassembly!(self::SysmatAssemblerVecOfMatSymm{T}, elem_mat_dim::FInt, ignore1::FInt, elem_mat_nmatrices::FInt, ndofs::FInt, ignore2::FInt) where {T<:Number}
    self.elem_mat_dim = elem_mat_dim
    self.buffer_length = elem_mat_nmatrices;
    self.buffer = Array{Tuple{Vector{FInt}, Matrix{T}}, 1}(undef, self.buffer_length);
    self.buffer_pointer = 1;
    self.ndofs = ndofs;
    return self
end

"""
    assemble!(self::SysmatAssemblerVecOfMatSymm{T}, mat::FMat{T}, dofnums::FIntVec, ignore::FIntVec) where {T<:Number}

Assemble a square symmetric matrix.

`dofnums` are the row degree of freedom numbers, the column degree of freedom
number input is ignored (the row and column numbers are assumed to be the same).
"""
function assemble!(self::SysmatAssemblerVecOfMatSymm{T}, mat::FMat{T},  dofnums::Union{FIntVec, FIntMat}, ignore::Union{FIntVec, FIntMat}) where {T<:Number}
    p = self.buffer_pointer
    @assert p <= self.buffer_length
    @assert self.elem_mat_dim == size(mat, 1)
    @assert self.elem_mat_dim == size(mat, 2)
    @assert self.elem_mat_dim == length(dofnums)
    d, m = deepcopy(vec(dofnums)), deepcopy(mat)
    fnz = 1
    @inbounds for j in 1:length(d)
        if d[j] == 0
            d[j] = fnz
            @inbounds for k in 1:length(d)
                m[j, k] = 0.0
            end
            @inbounds for k in 1:length(d)
                m[k, j] = 0.0
            end
        else
            fnz = d[j]
        end
    end
    self.buffer[p] = (d, m)
    self.buffer_pointer=p+1;
    return self
end

"""
    makematrix!(self::SysmatAssemblerVecOfMatSymm)

This function is dummy: the contents is preserved in the assembler itself.
"""
function makematrix!(self::SysmatAssemblerVecOfMatSymm)
    self.buffer_length = self.buffer_pointer-1
    return nothing
end

end


text:
using Test
using LogicCircuits
import .LogicCircuits.LoadSave:
   zoo_lc_file, zoo_vtree_file, parse_comment_line, parse_lc_header_line, parse_lc_literal_line, parse_literal_line, parse_lc_decision_line, parse_bias_line, parse_lc_file,
   CircuitFormatLine, BiasLine, DecisionLine, WeightedLiteralLine, CircuitHeaderLine, CircuitCommentLine, LCElement, CircuitFormatLines

@testset "Logistic circuit file parser tests" begin
   @test parse_comment_line("c blah blablah") isa CircuitCommentLine
   @test parse_comment_line("c") isa CircuitCommentLine
   @test parse_comment_line("c    blah blablah") isa CircuitCommentLine
   @test parse_lc_header_line("Logisitic Circuit") isa CircuitHeaderLine
   @test parse_lc_header_line("Logistic Circuit") isa CircuitHeaderLine
   @test parse_lc_literal_line("T 0 0 1 0.11139313932426485 0.5341755009918099 0.4104354811044485 0.2088029562981886 0.38317256253159404 0.21456111303752262 0.33798418436324884 0.7382343563376387 0.5769125897294547 0.13071237914862724") isa WeightedLiteralLine
   @test parse_lc_literal_line("F 1069 490 491 0.6277380017061743 -0.45260459349249044 0.34609986139917703 0.6004763090354547 0.2394524067773312 0.22081649811500942 -0.26666977618500204 0.14544044474614298 0.30372580539872435 0.2192352511676825") isa WeightedLiteralLine
   @test parse_lc_decision_line("D 1799 985 4 (472 474 0.27742886347699697 -0.0894114793745983 0.5298165134268861 0.5827938730880822 0.14116799704274996 0.3970938168763751 0.17798346381236296 0.08917988964843772 -0.05605305315306568 0.1702693902831316) (472 475 0.3833466224435187 0.8445851879217264 -0.3572571803165608 0.1793868357569113 -0.2373580813674068 0.670248227361854 -0.11119443329855791 0.13163431621813051 0.5421030929813475 0.25786192990838014) (473 474 1.0369907390437323 0.44729016983853126 -0.07892427803381961 0.38996680892303803 0.5285038536250287 0.3944289684978373 0.2762655604492141 0.556958084538147 0.2711846681681724 0.39922629776124985) (473 475 0.032883234975809694 -0.02256663542306192 0.6555725013615572 0.5140023339657676 0.11841852634121926 0.14907399101146324 -0.22404529652178906 -0.11976212824115842 -0.15206954052616856 0.0022385109727181413)") isa DecisionLine{LCElement}
   @test parse_lc_decision_line("D 10652 1001 2 (508 511 0.0008337025235152718 -6.048729079142479e-05 0.0012900540050118133 0.006382987897195768 0.00013330176570593142 -0.0034902489721742023 0.003162325487226574 -0.009619185307110537 0.043311151137203116 -0.007194862955461081) (509 511 0.023396488696149225 -0.000729066431265012 6.551173017401332e-06 0.05715185398005281 -0.008310854435718613 -0.003834142193742804 -0.005833871820252338 -0.05352747769146413 0.010573950714222884 0.03262423061844396)") isa DecisionLine{LCElement}
   @test parse_bias_line("B -0.6090213458520287 0.10061233805363132 -0.44510731039287776 -0.4536824618763301 -0.738392695523771 -0.5610245232140584 -0.4586543592164493 -0.07962059343551083 -0.2582953135054242 -0.03257926010007175") isa BiasLine
   @test parse_lc_file(zoo_lc_file("mnist-large.circuit")) isa CircuitFormatLines
end

@testset "Logistic circuit loader test" begin

   circuit = load_logic_circuit(zoo_lc_file("mnist-large.circuit"))

   @test circuit isa LogicCircuit
   @test isdecomposable(circuit)

   circuit = load_smooth_logic_circuit(zoo_lc_file("mnist-large.circuit"))

   @test circuit isa LogicCircuit
   @test isdecomposable(circuit)
   @test issmooth(circuit)

   circuit, vtree = load_struct_smooth_logic_circuit(zoo_lc_file("mnist-large.circuit"), zoo_vtree_file("balanced.vtree"))

   @test circuit isa StructLogicCircuit
   @test vtree isa PlainVtree
   @test isdecomposable(circuit)
   @test issmooth(circuit)

   # TODO: add a test and function to check whether a structured circuit respects its vtree (or another vtree)

end


text:
using Documenter, GWinstekPSP

makedocs(
  modules = [GWinstekPSP],
  clean = false,
  format = :html,
  sitename = "GWinstekPSP.jl",
  authors = "Chris Stook",
  linkcheck = !("skiplinks" in ARGS),
  pages = [
    "Home" => "index.md",
    "Instalation" => "install.md",
    "Usage" => "usage.md",
    "API" => "api.md"
    ],
  )

deploydocs(
  repo = "github.com/cstook/GWinstekPSP.jl.git",
  target = "build",
  deps = nothing,
  make = nothing,
  )


text:
# This example demonstrates one-sided communication,
# specifically activate Remote Memory Access (RMA)

using MPI

MPI.Init()
const world_sz = MPI.Comm_size(MPI.COMM_WORLD)
const rank = MPI.Comm_rank(MPI.COMM_WORLD)

# allocate memory
all_ranks = fill(-1, world_sz)
# create RMA window on all ranks
win = MPI.Win_create(all_ranks, MPI.COMM_WORLD)

#### first, let's MPI.Put on all ranks

# start the communication epoch
MPI.Win_fence(0, win)
# each rank writes to exposed windows of rank 0
# Signature: obj, target_rank, target_displacement, window
MPI.Put(rank, 0, rank, win)
# finish the communication epoch
MPI.Win_fence(0, win)
# print window content on all ranks
for j in 0:world_sz-1
    if rank == j
        println("After Put, Rank $rank:")
        @show all_ranks
    end
    MPI.Barrier(MPI.COMM_WORLD)
end
rank == 0 && println()

#### now, let's MPI.Get on all ranks

# start the communication epoch
MPI.Win_fence(0, win)
# each rank reads from exposed windows of rank 0
MPI.Get(all_ranks, 0, win)
# finish the communication epoch
MPI.Win_fence(0, win)
# print window content on all ranks
for j in 0:world_sz-1
    if rank == j
        println("After Get, Rank $rank:")
        @show all_ranks
    end
    MPI.Barrier(MPI.COMM_WORLD)
end

# free window
MPI.free(win)


text:
using Comonicon

@cast function scan_tfim(N::Int; lognsamples::Int=6, nprocs::Int=1, nthreads::Int=2)
    if nprocs > 1
        procs = addprocs(nprocs-1; exeflags=["--project=$(pkgdir(RenyiEntropyOfMagic))", "--threads=$nthreads"])
        Distributed.remotecall_eval(Main, procs, Expr(:toplevel, :(using RenyiEntropyOfMagic)))
    end
    results = scan_tfim_renyi_entropy(N; nsamples=10^lognsamples)
    serialize(pkgdir(RenyiEntropyOfMagic, "data", "$N-sites-result.dat"), results)
end

@cast function prepare_state(N::Int;nprocs::Int=1)
    if nprocs > 1
        procs = addprocs(nprocs-1; exeflags=["--project=$(pkgdir(RenyiEntropyOfMagic))"])
        Distributed.remotecall_eval(Main, procs, Expr(:toplevel, :(using RenyiEntropyOfMagic)))
    end
    prepare_states(N, 0.1:0.1:2.5)
end

@main


text:
using Symbolics, Test
import Symbolics: VarDomainPairing
import DomainSets: Interval, Ball, infimum, supremum, radius, center

domain = Interval(0, 1)
@test infimum(domain) == 0
@test supremum(domain) == 1

ball = Ball()
@test radius(ball) == 1.0
@test center(ball) == [0.0,0.0,0.0]

ball = Ball(2.5, [1,2])
@test radius(ball) == 2.5
@test center(ball) == [1,2]

@syms t x
t in domain
domains = [t ∈ Interval(0.0,1.0),
           x ∈ Interval(0.0,1.0)]
@test domains[1] isa VarDomainPairing
@test domains[2] isa VarDomainPairing

@syms y z
@test ((x,y) ∈ Ball(2.0, [0,0])) isa VarDomainPairing
@test ((x,y,z) ∈ Ball(1.5, [1,2,3])) isa VarDomainPairing

# Test tuple gets converted to Interval
var_domain_pair = t ∈ (0,1)
@test var_domain_pair isa VarDomainPairing
@test var_domain_pair.domain isa Interval

# Other types
t = Symbolics.Num(:t)
@assert (t ∈ domain) isa VarDomainPairing
t = Symbolics.variable(:t)
@assert (t ∈ domain) isa VarDomainPairing


text:
using HDF5
using Test
using Pkg

@info "libhdf5 v$(HDF5.h5_get_libversion())"

@testset "HDF5.jl" begin

include("plain.jl")
include("compound.jl")
include("custom.jl")
include("reference.jl")
include("dataspace.jl")
include("hyperslab.jl")
include("readremote.jl")
include("extend_test.jl")
include("gc.jl")
include("external.jl")
include("swmr.jl")
include("mmap.jl")
include("properties.jl")
include("table.jl")

try
    using MPI
    if HDF5.has_parallel()
        # basic MPI tests, for actual parallel tests we need to run in MPI mode
        include("mpio.jl")
    end
catch
end

end


text:
using MapImages
using Test

@testset "MapImages.jl" begin
    # Write your own tests here.
end


text:
# file for function used for node/dof numbering

# this can be generalized to 3D trivially
function populateDofNumbers(mesh::PumiMesh)
# populate the dofnums_Nptr with values calculated from
# nodenums_Nptr

#  resetAllIts2(mesh.m_ptr)
  # mesh iterator increment, retreval functions
  # TODO: checks if using mutable arrays is causing dynamic dispatch
#  iterators_inc = [incrementVertIt, incrementEdgeIt, incrementFaceIt, incrementElIt]
#  iterators_get = [getVert, getEdge, getFace, getEl]
  num_entities = mesh.numEntitiesPerType
  num_nodes_entity = mesh.numNodesPerType

  for etype = 1:(mesh.dim+1) # loop over entity types
    if (num_nodes_entity[etype] != 0)  # if no nodes on this type of entity, skip
      for entity_ptr in apf.MeshIterator(mesh.m_ptr, etype-1)  # loop over all entities of this type
	for node = 1:num_nodes_entity[etype]
          nodenum = apf.getNumberJ(mesh.nodenums_Nptr, entity_ptr, node-1, 0)
	  if nodenum != 0
	    for i=1:mesh.numDofPerNode
	      dofnum_i = (nodenum -1)*mesh.numDofPerNode + i
  	      apf.numberJ(mesh.dofnums_Nptr, entity_ptr, node-1, i-1, dofnum_i)
	    end  # end loop over dofsPerNode
	  end   # end if nodenum != 0
	end  # end loop over node

#      iterators_inc[etype]()
      end  # end loop over entitiesa
    end  # end if
  end  # end loop over entity types

#  resetAllIts2(mesh.m_ptr)

#  apf.writeVtkFiles("dofs_numbered", mesh.m_ptr)
  return nothing
end

# this can be generalized trivially
function populateNodeStatus(mesh::PumiMesh)
# populate the nodestatus_Nptr with values
# currently we set all nodes to status 3 (free)

#  resetAllIts2(mesh.m_ptr)
  # mesh iterator increment, retreval functions
  # TODO: checks if using mutable arrays is causing dynamic dispatch
#  iterators_inc = [incrementVertIt, incrementEdgeIt, incrementFaceIt, incrementElIt]
#  iterators_get = [getVert, getEdge, getFace, getEl]
  num_entities = mesh.numEntitiesPerType
  num_nodes_entity = mesh.numNodesPerType

#  mesh.numNodesPerType = num_nodes_entity

  for etype = 1:(mesh.dim + 1) # loop over entity types
    if (num_nodes_entity[etype] != 0)  # if no nodes on this type of entity, skip
      for entity_ptr in apf.MeshIterator(mesh.m_ptr, etype-1)  # loop over all entities of this type
	for node = 1:num_nodes_entity[etype]
	  apf.numberJ(mesh.nodestatus_Nptr, entity_ptr, node-1, 0, 3)
	end  # end loop over node

#      iterators_inc[etype]()
      end  # end loop over entitiesa
    end  # end if
  end  # end loop over entity types

#  resetAllIts2(mesh.m_ptr)
  return nothing
end

# this can be generalized to 3D trivially
function getDofNumbers(mesh::PumiMeshDG)
# populate array of dof numbers, in same shape as solution array u (or q)

#println("in getDofNumbers")
#println("numNodesPerElement = ", mesh.numNodesPerElement)

  mesh.dofs = Array{Int}(mesh.numDofPerNode, mesh.numNodesPerElement, mesh.shared_element_offsets[end] - 1)

  for i=1:mesh.numEl
    #TODO: use the non-allocating version
    dofnums = getGlobalNodeNumbers(mesh, i)

    for j=1:mesh.numNodesPerElement
      for k=1:mesh.numDofPerNode  # loop over dofs on the node
        mesh.dofs[k, j, i] = dofnums[k,j]
      end
    end
  end
  # get dof number of non-local elements
  # post receives first
  send_reqs = mesh.send_reqs
  recv_reqs = mesh.recv_reqs
  for i=1:mesh.npeers
    # get the pointer of the start location
    start_el = mesh.shared_element_offsets[i]
    numel = mesh.shared_element_offsets[i+1] - start_el
    ndata = mesh.numDofPerNode*mesh.numNodesPerElement*numel
    lin_idx = sub2ind(size(mesh.dofs), 1, 1, start_el)
    ptr_start_el = pointer(mesh.dofs, lin_idx)
    peer_i = mesh.peer_parts[i]
    recv_reqs[i] = MPI.Irecv!(ptr_start_el, ndata, peer_i, 1, mesh.comm)
  end

  # now send data
  dof_sendbuf = Array{Array{Int, 3}}(mesh.npeers)
  for i=1:mesh.npeers
    elnums_i = mesh.local_element_lists[i]
    numel = length(elnums_i)
    dof_sendbuf[i] = Array{Int}(mesh.numDofPerNode, mesh.numNodesPerElement, numel)
    sendbuf_i = dof_sendbuf[i]
    for j=1:length(elnums_i)
      elnum_j = elnums_i[j]
      for k=1:mesh.numNodesPerElement
        for p=1:mesh.numDofPerNode
          sendbuf_i[p, k, j] = mesh.dofs[p, k, elnum_j]
        end
      end
    end  # end loop over current list of elements

    # now do the send
    send_reqs[i] = MPI.Isend(sendbuf_i, mesh.peer_parts[i], 1, mesh.comm)
  end

  # figure out the local to global offset of the dof numbers
  # compute number of local dofs
  ndof = 0
  for i=1:(mesh.dim + 1)
    ndof += mesh.numNodesPerType[i]*mesh.numEntitiesPerType[i]*mesh.numDofPerNode
  end
  println(mesh.f, "ndof = ", ndof)
  # get all processes dof offsets
  dof_offsets = MPI.Allgather(ndof, mesh.comm)
  println(mesh.f, "dof_offsets = ", dof_offsets)
  # figure out own dof_offset
  dof_offset = 0
  for i=1:mesh.myrank  # sum all dofs up to but *not* including self
    dof_offset += dof_offsets[i]
  end
  mesh.dof_offset = dof_offset

  # figure out peer dof_offsets
  peer_dof_offsets = Array{Int}(mesh.npeers)
  for i=1:mesh.npeers
    dof_offset_i = 0
    for j=1:mesh.peer_parts[i]
      dof_offset_i += dof_offsets[j]
    end
    peer_dof_offsets[i] = dof_offset_i
  end

  println(mesh.f, "peer_parts = ", mesh.peer_parts)
  println(mesh.f, "peer_dof_offsets = ", peer_dof_offsets)
  # adjust other elements dof such that dof + dof_offset of this process = global
  # dof number
  # store local dof number + (offset of owner - offset of host) in dofs array
  dofs = mesh.dofs
  for i=1:mesh.npeers
    idx, stat =  MPI.Waitany!(recv_reqs)
    # do subtraction first to avoid overflow
    offset = peer_dof_offsets[idx] - dof_offset
    start_el = mesh.shared_element_offsets[idx]
    end_el = mesh.shared_element_offsets[idx+1] - 1
    for j=start_el:end_el
      for k=1:mesh.numNodesPerElement
        for p=1:mesh.numDofPerNode
          dofs[p, k, j] += offset
        end
      end
    end

  end  # end loop over peers

  # wait for the sends to finish before returning
  for i=1:mesh.npeers
    MPI.Wait!(send_reqs[i])
  end

  #println("mesh.dof = ", mesh.dofs)

  return nothing

end

"""
  Node numbering algorithm that loops over elements and numbers their downward
  adjacencies, generating a decent sparsity pattern.  This algorithm is
  primarily intended for CG meshes, although I think it works for DG meshes too.
  This probably only works for 2D CG meshes.

  **Inputs**

   * mesh: a PumiMesh
   * number_dofs: Bool, if true, number the dofs (mesh.dofnums_Npts), otherwise
                  number the nodes (mesh.nodenums_Nptr)
"""
function numberNodes(mesh::PumiMesh, number_dofs=false)
# assign node numbers to entire mesh, or dof numbers if number_dofs=true,
# using the correct Numbering pointer
# assumes mesh elements have already been reordered
# this works for high order elements

  @assert mesh.dim == 2

  # initally number all dofs as numDof+1 to 2*numDof
  # this allows quick check to see if somthing is labelled or not
  numberNodesElement(mesh, number_dofs=number_dofs, start_at_one=false)

  # calculate number of nodes, dofs
  numnodes = mesh.numNodes

  num_entities = mesh.numEntitiesPerType
  num_nodes_entity = mesh.numNodesPerType

#  mesh.numNodesPerType = num_nodes_entity

  if number_dofs
    numbering_ptr = mesh.dofnums_Nptr
    curr_dof = mesh.numDof + 1
    dofpernode = mesh.numDofPerNode
    numDof = mesh.numDofPerNode*numnodes
  else  # do node numbering
    numbering_ptr = mesh.nodenums_Nptr
    curr_dof = mesh.numNodes + 1
    dofpernode = 1
    numDof = numnodes
  end

  verts_i = Array{Ptr{Void}}(12)
  edges_i = Array{Ptr{Void}}(12)
#  resetAllIts2(mesh.m_ptr)
  el_i_ptr = Ptr{Void}(0)  # hold current element
# TODO: move all if statements out one for loop (check only first dof on each node)
  curr_dof = 1
  for el_i_ptr in apf.MeshIterator(mesh.m_ptr, mesh.dim)
    # get vertices, edges for this element
    numVert = apf.getDownward(mesh.m_ptr, el_i_ptr, 0, verts_i)
#    println("verts_i = ", verts_i)
#    println("mesh.verts = ", mesh.verts)
    numEdge = apf.getDownward(mesh.m_ptr, el_i_ptr, 1, edges_i)
    for j=1:3  # loop over vertices, edges  #TODO: why is this 1:3?
      vert_j = verts_i[j]
      edge_j = edges_i[j]
      for k=1:num_nodes_entity[1] # loop over vertex nodes
        for p=1:dofpernode
          dofnum_k = apf.getNumberJ(numbering_ptr, vert_j, k-1, p-1)
          if dofnum_k > numDof  # still has initial number
            # give it new (final) number
            apf.numberJ(numbering_ptr, vert_j, 0, k-1, curr_dof)
            curr_dof += 1
          end
        end
      end  # end loop over vertex dofs

      # loop over nodes on edge
      for k=1:num_nodes_entity[2]  # loop over nodes
	for p=1:dofpernode  # loop over dofs
	  dofnum_p = apf.getNumberJ(numbering_ptr, edge_j, k-1, p-1)
	  if dofnum_p > numDof  # still has initial number
	    # give it new (final) number
	    apf.numberJ(numbering_ptr, edge_j, k-1, p-1, curr_dof)
	    curr_dof += 1
	  end
	end
      end
    end  # end loop over vertices, edges
    # label face nodes
    for k=1:num_nodes_entity[3]  # loop over nodes on face
        # visit the nodes in the SBP ordering, so consider k to be the SBP
        # node number
        # this only works when there are no nodes on edges or vertices, because
        # we would have to do additional translation to take them into account

      for p=1:dofpernode  # loop over dofs
	dofnum_p = apf.getNumberJ(numbering_ptr, el_i_ptr, k-1, p-1)
	if dofnum_p > numDof
	  apf.numberJ(numbering_ptr, el_i_ptr, k-1, p-1, curr_dof)
	  curr_dof += 1
	end
      end
    end  # end loop over face nodes
  end  # end loop over elements

  if (curr_dof -1) != numDof
    error("Warning: number of dofs assigned is not equal to the expected number" *
          "number of dofs assigned = $(curr_dof-1) ,expected number =  $numDof")
  end

  return nothing

end

"""
  This function numbers the mesh entities using a simple iteration over all
  entities of each dimension.  This generates a bad sparsity pattern for CG
  meshes.  For DG meshes the sparsity pattern is ok if Pumi reordered the
  elements in its database, otherwise the sparsity pattern is bad.

  **Inputs**

   * mesh: a PumiMesh

  **Keyword Arguments**

   * number_dofs: if true, number the dofs (mesh.dofnums_Nptr), otherwise
                  number the nodes (mesh.nodenums_Nptr), default false
   * start_at_one: if true, start numbering dofs at 1, otherwise start at
                   number of dofs/nodes + 1.  This is useful as an initial
                   numbering for other algorithms, because any dof/node
                   with a number greater than the number of dofs/nodes has not
                   been relabeled yet.
"""
function numberNodesElement(mesh::PumiMesh; number_dofs=false, start_at_one=true)
# number all nodes (or dofs) by iterating over all entities, from lowest
# dimension to highest, and numbering the nodes in order

  # calculate number of nodes, dofs
  numnodes = mesh.numNodes

  num_entities = mesh.numEntitiesPerType
  num_nodes_entity = mesh.numNodesPerType

  if number_dofs
    numbering_ptr = mesh.dofnums_Nptr
    curr_dof = mesh.numDof + 1
    dofpernode = mesh.numDofPerNode
    numDof = mesh.numDofPerNode*numnodes
  else  # do node numbering
    numbering_ptr = mesh.nodenums_Nptr
    curr_dof = mesh.numNodes + 1
    dofpernode = 1
    numDof = numnodes
  end

  if (numDof > 2^30 || numDof < 0)
    error("Warning: too many dofs, renumbering will fail")
  end

#  println("num_entities = ", num_entities)
#  println("num_nodes_entity = ", num_nodes_entity)

  if start_at_one
    curr_dof = 1
  else
    curr_dof = numDof+1
  end

  if mesh.isDG
    nodemap =  mesh.nodemapSbpToPumi
  else  # fake nodemap: always 1:1, this works because we apf.iterate over the
        # entities of each dimension, so every entity is visited exactly once
        # Such an algorith generates a horrible sparsity pattern, but this
        # is the default algorithm, so whatever
    nodemap = collect(eltype(mesh.nodemapSbpToPumi), 1:mesh.numNodesPerElement)
  end

  for etype = 1:(mesh.dim + 1) # loop over entity types
    if (num_nodes_entity[etype] != 0)  # if no nodes on this type of entity, skip
      for entity_ptr in apf.MeshIterator(mesh.m_ptr, etype - 1)  # loop over all entities of this type
	for node = 1:num_nodes_entity[etype]
          pumi_node = nodemap[node]  # this appears to only work for elements
                                     # with nodes on the interior, otherwise
                                     # there would need to be different node
                                     # maps for entities of different dimension
	  for dof = 1:dofpernode
	    apf.numberJ(numbering_ptr, entity_ptr, pumi_node-1, dof-1, curr_dof)
	    curr_dof += 1
	  end  # end loop over dof
	end  # end loop over node
      end  # end loop over entitiesa
    end  # end if
  end  # end loop over entity types

  return nothing
end

function getDofNumbers(mesh::PumiMesh2)
# populate array of dof numbers, in same shape as solution array u (or q)

mesh.dofs = Array{Int32}(mesh.numDofPerNode, mesh.numNodesPerElement, mesh.numEl)

for i=1:mesh.numEl
  dofnums = getGlobalNodeNumbers(mesh, i)

  for j=1:mesh.numNodesPerElement
    for k=1:mesh.numDofPerNode  # loop over dofs on the node
      mesh.dofs[k, j, i] = dofnums[k,j]
    end
  end
end

#println("mesh.dof = ", mesh.dofs)

return nothing

end

"""
  This function numbers the dof/nodes of the mesh using and adjacency-type
  algorithm, starting from the mesh element nearest the specified coordinates.
  The dofs/nodes are numbered in an increasing order.  This only works for
  DG meshes.  This function overwrites the mesh.el_Nptr numbering with
  the new element numbering, in addition to numbering the dofs/nodes.

  Note that this function updates the mesh.verts, mesh.edges, mesh.faces,
  mesh.elements arrays so they match the numbering (ie. mesh.element[1] has
  the number 0 in mesh.el_Nptr, because the numbering is zero based).

  The numbering of dofs/nodes is 1-based, and the numbering of elements is
  zero-based.

  **Inputs**

   * mesh: a DG mesh
   * start_coords: vector of length mesh.dim specifying starting location
                   of the numbering.  The element with the closests linear
                   centroid (based on vertex locations only) is selected as
                   the starting element.
  * number_dofs: if true, numbers the dofs (mesh.dofnums_Nptr), otherwise
                 numbers the nodes (mesh.nodenums_Nptr).
"""
function numberNodesWindy(mesh::PumiMeshDG, start_coords, number_dofs=false)
# number nodes and elements starting from a particular location

  # calculate number of nodes, dofs
  numDof = mesh.numDof
  if (numDof > 2^30 || numDof < 0)
    println("Warning: too many dofs, renumbering will fail")
  end

  if number_dofs
    numbering_ptr = mesh.dofnums_Nptr
    dofpernode = mesh.numDofPerNode
    numDof = mesh.numDofPerNode*numnodes
  else  # do node numbering
    numbering_ptr = mesh.nodenums_Nptr
    dofpernode = 1
    numDof = mesh.numNodes
  end
  dim = mesh.dim
  el_Nptr = mesh.el_Nptr
  adj_els = Array{Ptr{Void}}(mesh.numFacesPerElement)
  nodemap = mesh.nodemapSbpToPumi
  numDofPerElement = dofpernode*mesh.numNodesPerElement

  # initially number all components in range (numEl+1):(2*numEl)
  curr_elnum = mesh.numEl+1
  for el_i in apf.MeshIterator(mesh.m_ptr, dim)
    apf.numberJ(el_Nptr, el_i, 0, 0, curr_elnum)
    curr_elnum += 1
#    incrementIt(dim)
  end
  @assert (curr_elnum - 1) == 2*mesh.numEl

  # do the final numbering
  start_el = getStartEl(mesh, start_coords)
  que = FIFOQueue{Ptr{Void}}(size_hint = div(mesh.numEl, 2))
  push!(que, start_el)

  curr_elnum = mesh.numEl - 1
  curr_dof = numDof
  # the dreaded while loop
  while (!isempty(que))
#    print("\n")
    curr_el = pop!(que)

    # only unlabelled entities are added to the que, and they are added
    # exactly once, so no need to check here
    apf.numberJ(el_Nptr, curr_el, 0, 0, curr_elnum)
    curr_elnum -= 1

    # now label dofs
    # We have to label the dofs on each element in ascending order, so that
    # eqn.res = eqn.res_vec
    _curr_dof = curr_dof - numDofPerElement + 1
    for i=1:mesh.numNodesPerElement
      pumi_node = nodemap[i]
      for j=1:dofpernode
        apf.numberJ(numbering_ptr, curr_el, pumi_node-1, j-1, _curr_dof)
        _curr_dof += 1
      end
    end
    curr_dof -= numDofPerElement

    # add face adjacent neighbor elements to the que
    numadj = apf.countBridgeAdjacent(mesh.m_ptr, curr_el, dim-1, dim)
    apf.getBridgeAdjacent(adj_els)

    for i=1:numadj
      el_i = adj_els[i]
      elnum = apf.getNumberJ(el_Nptr, el_i, 0, 0)

       # if not numbered and not in que
       # when adding to que, scale elnum so it is in range (2*numEl+1):3*numEl
       # as long as there are at least 3 nodes per element, node numbers will
       # overflow before this
      if (elnum > mesh.numEl) && (elnum <= 2*mesh.numEl)
        apf.numberJ(el_Nptr, el_i, 0, 0, elnum + mesh.numEl)
        push!(que, el_i)
      end
    end  # end loop over adjacent elements

    # this can only happen if either 1) this algorithm is broken, or
    # 2) there are disjoint sets of elements on this process
    # Assume it is 2 and find a new starting element
    if curr_elnum != -1 && isempty(que)
      new_startel = getStartEl(mesh, start_coords, el_Nptr, mesh.numEl - 1)
      @assert new_startel != C_NULL
      push!(que, new_startel)
    end

  end  # end while loop

  #-----------------------------------------
  # Error checking and cleanup

  if curr_dof != 0
    println("Warning: number of dofs assigned is not equal to the expected number")
    println("number of dofs remaining = ", curr_dof)
    throw(ErrorException("dof numbering failed"))
  end

  if curr_elnum != -1
    throw(ErrorException("element numbering failed"))
  end
#=
  # the element number should have been zero based (oops), so decrement
  # the element numbers
  for i=1:mesh.numEl
    el_i = mesh.elements[i]
    idx = apf.getNumberJ(mesh.el_Nptr, el_i, 0, 0)
    apf.numberJ(mesh.el_Nptr, el_i, 0, 0, idx - 1)
  end
=#
  # need to update the element pointer array because element numbering has
  # changed
  #TODO: don't reallocate the arrays
  mesh.verts, mesh.edges, mesh.faces, mesh.elements = getEntityPointers(mesh)
  return nothing
end

"""
  Finds the element with the (linear) centroid closest to the specified
  coordinates.  Can optionally skip certain elements that are already numbered.

  **Inputs**

   * mesh
   * start_coords: vector of length `mesh.dim` containing the point that
                   the element will be closest to
   * el_Nptr: an apf::Numbering* that numbers the mesh elements (optional)
   * sentinel: the sentinal value for element numbers.  Element numbers less
               than or equal to this value will be skipped.  This argument
               is required if `el_Nptr` is provided. This value should be
               zero-based

"""
function getStartEl(mesh::PumiMeshDG, start_coords, el_Nptr::Ptr{Void}=C_NULL, sentinel::Integer=-1)
  # find the element with centroid closest to start_coords

  numVertsPerElement = mesh.numTypePerElement[1]
  coords = zeros(3, numVertsPerElement)
  centroid = zeros(3)
  dim = mesh.dim  # hoist the lookup

#  resetIt(dim)

  min_norm = typemin(Float64)
  min_el = C_NULL

  for el_i  in apf.MeshIterator(mesh.m_ptr, dim)

    # skip already numbered elements
    if el_Nptr != C_NULL
      elnum = apf.getNumberJ(el_Nptr, el_i, 0, 0)
      if elnum <= sentinel
        continue
      end
    end

    getElementCoords(mesh, el_i, coords)

    # compute centroid
    fill!(centroid, 0.0)
    for j=1:dim
      for k=1:numVertsPerElement
        centroid[j] += coords[j, k]
      end
      centroid[j] = centroid[j]/numVertsPerElement
    end

    # compute norm of distance from start_coords
    nrm = 0.0
    for j=1:dim
      nrm += (centroid[j] - start_coords[j])*(centroid[j] - start_coords[j])
    end
    # skip the square root because it doesn't affect the comparison

    if nrm > min_norm
      min_norm = nrm
      min_el = el_i
    end

  end

  return min_el
end

"""
  Verifies the element numbering, specifically that every element number
  is assigned exactly once, and all numbers are in the range of 0 to
  mesh.numEl -1
"""
function checkElementNumbering(mesh, el_Nptr)

  # check all numbers are assigned, none more than once
  elnums = zeros(Int, mesh.numEl)

  for el in apf.MeshIterator(mesh.m_ptr, mesh.dim)
    elnum = apf.getNumberJ(el_Nptr, el, 0, 0) + 1
    @assert elnum >= 1
    @assert elnum <= mesh.numEl
    elnums[elnum] += 1
  end

  for i=1:mesh.numEl
    @assert elnums[i] == 1
  end

  return nothing
end

"""
  Verifies the dof numbering, specifically that every dof number is assigned
  exactly once and all numbers are in the range 1 to mesh.numDof
"""
function checkDofNumbering(mesh, dof_Nptr)

  dofnums = zeros(UInt8, mesh.numDof)

  for el in apf.MeshIterator(mesh.m_ptr, mesh.dim)
    for j=1:mesh.numNodesPerElement
      for k=1:mesh.numDofPerNode
        dofnum = apf.getNumberJ(dof_Nptr, el, j-1, k-1)
        @assert dofnum > 0
        @assert dofnum <= mesh.numDof
        dofnums[dofnum] += 1
      end
    end
  end

  for i=1:mesh.numDof
    @assert dofnums[i] == 1
  end

  return nothing
end