text:
abstract type AbstractMPS end

"""
    length(::MPS/MPO)

The number of sites of an MPS/MPO.
"""
Base.length(m::AbstractMPS) = length(m.data)

"""
    ITensors.data(::MPS/MPO)

The Vector storage of an MPS/MPO.

This is mostly for internal usage, please let us
know if there is functionality not available for
MPS/MPO you would like.
"""
data(m::AbstractMPS) = m.data

leftlim(m::AbstractMPS) = m.llim

rightlim(m::AbstractMPS) = m.rlim

function setleftlim!(m::AbstractMPS, new_ll::Int)
  m.llim = new_ll
end

function setrightlim!(m::AbstractMPS, new_rl::Int)
  m.rlim = new_rl
end

isortho(m::AbstractMPS) = leftlim(m)+1 == rightlim(m)-1

function orthocenter(m::T) where {T<:AbstractMPS}
  !isortho(m) && error("$T has no well-defined orthogonality center")
  return leftlim(m)+1
end

Base.getindex(M::AbstractMPS,
              n::Integer) = getindex(data(M),n)

function Base.setindex!(M::AbstractMPS,
                        T::ITensor,
                        n::Integer;
                        set_limits::Bool = true)
  if set_limits
    (n <= leftlim(M)) && setleftlim!(M,n-1)
    (n >= rightlim(M)) && setrightlim!(M,n+1)
  end
  data(M)[n] = T
  return M
end

Base.copy(m::AbstractMPS) = typeof(m)(copy(data(m)),
                                      leftlim(m),
                                      rightlim(m))

Base.similar(m::AbstractMPS) = typeof(m)(similar(data(m)),
                                         0,
                                         length(m))

Base.deepcopy(m::AbstractMPS) = typeof(m)(deepcopy(data(m)),
                                          leftlim(m),
                                          rightlim(m))

Base.eachindex(m::AbstractMPS) = 1:length(m)

Base.iterate(M::AbstractMPS) = iterate(data(M))

Base.iterate(M::AbstractMPS, state) = iterate(data(M), state)

"""
    unique_siteind(A::MPO, B::MPS, j::Int)
    unique_siteind(A::MPO, B::MPO, j::Int)

Get the site index of MPO `A` that is unique to `A` (not shared with MPS/MPO `B`).
"""
function unique_siteind(A::AbstractMPS, B::AbstractMPS, j::Int)
  N = length(A)
  j == 1 && return uniqueind(A[j], A[j+1], B[j])
  j == N && return uniqueind(A[j], A[j-1], B[j])
  return uniqueind(A[j], A[j-1], A[j+1], B[j])
end

"""
    unique_siteinds(A::MPO, B::MPS)
    unique_siteinds(A::MPO, B::MPO)

Get the site indices of MPO `A` that are unique to `A` (not shared with MPS/MPO `B`), as a `Vector{<:Index}`.
"""
function unique_siteinds(A::AbstractMPS, B::AbstractMPS)
  return [unique_siteind(A, B, j) for j in eachindex(A)]
end

"""
    common_siteind(A::MPO, B::MPS, j::Int)
    common_siteind(A::MPO, B::MPO, j::Int)

Get the site index of MPO `A` that is shared with MPS/MPO `B`.
"""
function common_siteind(A::AbstractMPS, B::AbstractMPS, j::Int)
  return commonind(A[j], B[j])
end

"""
    common_siteinds(A::MPO, B::MPS)
    common_siteinds(A::MPO, B::MPO)

Get the site indices of MPO `A` that are shared with MPS/MPO `B`, as a `Vector{<:Index}`.
"""
function common_siteinds(A::AbstractMPS, B::AbstractMPS)
  return [common_siteind(A, B, j) for j in eachindex(A)]
end

function Base.map!(f::Function, M::AbstractMPS)
  for i in eachindex(M)
    setindex!(M, f(M[i]), i; set_limits = false)
  end
  return M
end

Base.map(f::Function, M::AbstractMPS) = map!(f, copy(M))

for fname in (:dag,
              :prime,
              :setprime,
              :noprime,
              :addtags,
              :removetags,
              :replacetags,
              :settags)
  @eval begin
    """
        $($fname)(M::MPS, args...; kwargs...)

        $($fname)(M::MPO, args...; kwargs...)

    Apply $($fname) to all ITensors of an MPS/MPO, returning a new MPS/MPO.

    The ITensors of the MPS/MPO will be a view of the storage of the original ITensors.
    """
    $fname(M::AbstractMPS,
           args...;
           kwargs...) = map(m -> $fname(m, args...;
                                        kwargs...), M)

    """
        $($fname)!(M::MPS, args...; kwargs...)

        $($fname)!(M::MPO, args...; kwargs...)

    Apply $($fname) to all ITensors of an MPS/MPO in-place.
    """
    $(Symbol(fname, :!))(M::AbstractMPS,
                         args...;
                         kwargs...) = map!(m -> $fname(m, args...;
                                                       kwargs...), M)

  end
end

function map_linkinds!(f::Function, M::AbstractMPS)
  for i in eachindex(M)[1:end-1]
    l = linkind(M, i)
    if !isnothing(l)
      l̃ = f(l)
      setindex!(M, replaceind(M[i], l, l̃), i;
                set_limits = false)
      setindex!(M, replaceind(M[i+1], l, l̃), i+1;
                set_limits = false)
    end
  end
  return M
end

map_linkinds(f::Function, M::AbstractMPS) = map_linkinds!(f, copy(M))

function map_common_siteinds!(f::Function, M1::AbstractMPS,
                                           M2::AbstractMPS)
  length(M1) != length(M2) && error("MPOs/MPSs must be the same length")
  for i in eachindex(M1)
    s = common_siteind(M1, M2, i)
    if !isnothing(s)
      s̃ = f(s)
      setindex!(M1, replaceind(M1[i], s, s̃), i;
                set_limits = false)
      setindex!(M2, replaceind(M2[i], s, s̃), i;
                set_limits = false)
    end
  end
  return M1, M2
end

function map_common_siteinds(f::Function, M1::AbstractMPS,
                                          M2::AbstractMPS)
  return map_common_siteinds!(f, copy(M1), copy(M2))
end

function map_unique_siteinds!(f::Function, M1::AbstractMPS,
                                           M2::AbstractMPS)
  length(M1) != length(M2) && error("MPOs/MPSs must be the same length")
  for i in eachindex(M1)
    s = unique_siteind(M1, M2, i)
    if !isnothing(s)
      s̃ = f(s)
      setindex!(M1, replaceind(M1[i], s, s̃), i;
                set_limits = false)
    end
  end
  return M1
end

function map_unique_siteinds(f::Function, M1::AbstractMPS,
                                        M2::AbstractMPS)
  return map_unique_siteinds!(f, copy(M1), M2)
end

for fname in (:sim,
              :prime,
              :setprime,
              :noprime,
              :addtags,
              :removetags,
              :replacetags,
              :settags)
  fname_linkinds = Symbol(fname, :_linkinds)
  fname_linkinds_inplace = Symbol(fname_linkinds, :!)
  fname_common_siteinds = Symbol(fname, :_common_siteinds)
  fname_common_siteinds_inplace = Symbol(fname_common_siteinds, :!)
  fname_unique_siteinds = Symbol(fname, :_unique_siteinds)
  fname_unique_siteinds_inplace = Symbol(fname_unique_siteinds, :!)

  @eval begin
    """
        $($fname_linkinds)(M::MPS, args...; kwargs...)

        $($fname_linkinds)(M::MPO, args...; kwargs...)

    Apply $($fname) to all link indices of an MPS/MPO, returning a new MPS/MPO.

    The ITensors of the MPS/MPO will be a view of the storage of the original ITensors.
    """
    $fname_linkinds(M::AbstractMPS,
                    args...;
                    kwargs...) = map_linkinds(i -> $fname(i, args...;
                                                         kwargs...), M)

    """
        $($fname_linkinds)!(M::MPS, args...; kwargs...)

        $($fname_linkinds)!(M::MPO, args...; kwargs...)

    Apply $($fname) to all link indices of the ITensors of an MPS/MPO in-place.
    """
    function $fname_linkinds_inplace(M::AbstractMPS,
                                     args...;
                                     kwargs...)
      return map_linkinds!(i -> $fname(i, args...;
                                      kwargs...), M)
    end

    """
        $($fname_common_siteinds)(M1::MPO, M2::MPS, args...; kwargs...)

        $($fname_common_siteinds)(M1::MPO, M2::MPO, args...; kwargs...)

    Apply $($fname) to the site indices that are shared by `M1` and `M2`.

    Returns new MPSs/MPOs. The ITensors of the MPSs/MPOs will be a view of the storage of the original ITensors.
    """
    function $fname_common_siteinds(M1::AbstractMPS,
                                    M2::AbstractMPS,
                                    args...;
                                    kwargs...)
      return map_common_siteinds(i -> $fname(i, args...;
                                             kwargs...), M1, M2)
    end

    """
        $($fname_common_siteinds)!(M1::MPO, M2::MPS, args...; kwargs...)

        $($fname_common_siteinds)!(M1::MPO, M2::MPO, args...; kwargs...)

    Apply $($fname) to the site indices that are shared by `M1` and `M2`. Returns new MPSs/MPOs.

    Modifies the input MPSs/MPOs in-place.
    """
    function $fname_common_siteinds_inplace(M1::AbstractMPS,
                                            M2::AbstractMPS,
                                            args...;
                                            kwargs...)
      return map_common_siteinds!(i -> $fname(i, args...;
                                              kwargs...), M1, M2)
    end

    """
        $($fname_unique_siteinds)(M1::MPO, M2::MPS, args...; kwargs...)

    Apply $($fname) to the site indices of `M1` that are not shared with `M2`. Returns new MPSs/MPOs.

    The ITensors of the MPSs/MPOs will be a view of the storage of the original ITensors.
    """
    function $fname_unique_siteinds(M1::AbstractMPS,
                                    M2::AbstractMPS,
                                    args...;
                                    kwargs...)
      return map_unique_siteinds(i -> $fname(i, args...;
                                             kwargs...), M1, M2)
    end

    """
        $($fname_unique_siteinds)!(M1::MPO, M2::MPS, args...; kwargs...)

    Apply $($fname) to the site indices of `M1` that are not shared with `M2`. Modifies the input MPSs/MPOs in-place.
    """
    function $fname_unique_siteinds_inplace(M1::AbstractMPS,
                                            M2::AbstractMPS,
                                            args...;
                                            kwargs...)
      return map_unique_siteinds!(i -> $fname(i, args...;
                                              kwargs...), M1, M2)
    end
  end
end

"""
    maxlinkdim(M::MPS)

    maxlinkdim(M::MPO)

Get the maximum link dimension of the MPS or MPO.
"""
function maxlinkdim(M::AbstractMPS)
  md = 0
  for b ∈ eachindex(M)[1:end-1]
    l = linkind(M, b)
    linkdim = isnothing(l) ? 0 : dim(l)
    md = max(md, linkdim)
  end
  md
end

"""
    linkind(M::MPS, j::Int)

    linkind(M::MPO, j::Int)

Get the link or bond Index connecting the
MPS or MPO tensor on site j to site j+1
"""
function linkind(M::AbstractMPS, j::Int)
  N = length(M)
  j ≥ length(M) && return nothing
  return commonind(M[j], M[j+1])
end

function _log_or_not_dot(M1::MPST,
                         M2::MPST,
                         loginner::Bool;
                         make_inds_match::Bool = true)::Number where {MPST <: AbstractMPS}
  N = length(M1)
  if length(M2) != N
    throw(DimensionMismatch("inner: mismatched lengths $N and $(length(M2))"))
  end
  M1dag = dag(M1)
  sim_linkinds!(M1dag)
  if make_inds_match
    replace_siteinds!(M1dag, siteinds(M2))
  end
  O = M1dag[1] * M2[1]

  if loginner
    normO = norm(O)
    log_inner_tot = log(normO)
    O ./= normO
  end

  for j in eachindex(M1)[2:end]
    O = (O*M1dag[j])*M2[j]

    if loginner
      normO = norm(O)
      log_inner_tot += log(normO)
      O ./= normO
    end

  end

  if loginner
    return log_inner_tot
  end

  return O[]
end

"""
    dot(A::MPS, B::MPS; make_inds_match = true)
    inner(A::MPS, B::MPS; make_inds_match = true)

    dot(A::MPO, B::MPO)
    inner(A::MPO, B::MPO)

Compute the inner product `<A|B>`. If `A` and `B` are MPOs, computes the Frobenius inner product.

If `make_inds_match = true`, the function attempts to make
the site indices match before contracting (so for example, the
inputs can have different site indices, as long as they
have the same dimensions or QN blocks).

For now, `make_inds_match` is only supported for MPSs.

See also `logdot`/`loginner`.
"""
function LinearAlgebra.dot(M1::MPST,
                           M2::MPST;
                           kwargs...) where {MPST <: AbstractMPS}
  return _log_or_not_dot(M1, M2, false; kwargs...)
end

"""
    logdot(A::MPS, B::MPS; make_inds_match = true)
    loginner(A::MPS, B::MPS; make_inds_match = true)

    logdot(A::MPO, B::MPO)
    loginner(A::MPO, B::MPO)

Compute the logarithm of the inner product `<A|B>`. If `A` and `B` are MPOs, computes the logarithm of the Frobenius inner product.

This is useful for larger MPS/MPO, where in the limit of large numbers of sites the inner product can diverge or approach zero.

If `make_inds_match = true`, the function attempts to make
the site indices match before contracting (so for example, the
inputs can have different site indices, as long as they
have the same dimensions or QN blocks).

For now, `make_inds_match` is only supported for MPSs.
"""
function logdot(M1::MPST,
                M2::MPST;
                kwargs...) where {MPST <: AbstractMPS}
  return _log_or_not_dot(M1, M2, true; kwargs...)
end

inner(M1::MPST,
      M2::MPST;
      kwargs...) where {MPST <: AbstractMPS} = dot(M1, M2; kwargs...)

loginner(M1::MPST,
         M2::MPST;
         kwargs...) where {MPST <: AbstractMPS} = logdot(M1, M2; kwargs...)

"""
    norm(A::MPS)

    norm(A::MPO)

Compute the norm of the MPS or MPO.

See also `lognorm`.
"""
function LinearAlgebra.norm(M::AbstractMPS)
  return sqrt(dot(M, M))
end

"""
    lognorm(A::MPS)

    lognorm(A::MPO)

Compute the logarithm of the norm of the MPS or MPO.

This is useful for larger MPS/MPO that are not gauged, where in the limit of large numbers of sites the norm can diverge or approach zero.

See also `norm` and `loginner`/`logdot`.
"""
function lognorm(M::AbstractMPS)
  return 0.5 * logdot(M, M)
end

function plussers(::Type{T},
                  left_ind::Index,
                  right_ind::Index,
                  sum_ind::Index) where {T<:Array}
  total_dim    = dim(left_ind) + dim(right_ind)
  total_dim    = max(total_dim, 1)
  # TODO: I am not sure if we should be using delta
  # tensors for this purpose? I think we should consider
  # not allowing them to be made with different index sizes
  #left_tensor  = δ(left_ind, sum_ind)
  left_tensor  = diagITensor(1.0,left_ind, sum_ind)
  right_tensor = ITensor(right_ind, sum_ind)
  for i in 1:dim(right_ind)
    right_tensor[right_ind(i), sum_ind(dim(left_ind) + i)] = 1
  end
  return left_tensor, right_tensor
end

"""
    add(A::MPS, B::MPS; kwargs...)
    +(A::MPS, B::MPS; kwargs...)

    add(A::MPO, B::MPO; kwargs...)
    +(A::MPO, B::MPO; kwargs...)

Add two MPS/MPO with each other, with some optional
truncation.
"""
function Base.:+(A::T, B::T; kwargs...) where {T <: AbstractMPS}
  N = length(A)
  length(B) != N && throw(DimensionMismatch("lengths of MPOs A ($N) and B ($(length(B))) do not match"))
  orthogonalize!(A, 1; kwargs...)
  orthogonalize!(B, 1; kwargs...)
  C = similar(A)
  rand_plev = 13124
  lAs = [linkind(A, i) for i in 1:N-1]
  prime!(A, rand_plev, "Link")

  first  = Vector{ITensor{2}}(undef,N-1)
  second = Vector{ITensor{2}}(undef,N-1)
  for i in 1:N-1
    lA = linkind(A, i)
    lB = linkind(B, i)
    r  = Index(dim(lA) + dim(lB), tags(lA))
    f, s = plussers(typeof(data(A[1])), lA, lB, r)
    first[i]  = f
    second[i] = s
  end
  C[1] = A[1] * first[1] + B[1] * second[1]
  for i in 2:N-1
      C[i] = dag(first[i-1]) * A[i] * first[i] + dag(second[i-1]) * B[i] * second[i]
  end
  C[N] = dag(first[N-1]) * A[N] + dag(second[N-1]) * B[N]
  prime!(C, -rand_plev, "Link")
  truncate!(C; kwargs...)
  return C
end

add(A::T, B::T;
    kwargs...) where {T <: AbstractMPS} = +(A, B; kwargs...)

"""
    sum(A::Vector{MPS}; kwargs...)

    sum(A::Vector{MPO}; kwargs...)

Add multiple MPS/MPO with each other, with some optional
truncation.
"""
function Base.sum(A::Vector{T};
                  kwargs...) where {T <: AbstractMPS}
  length(A) == 0 && return T()
  length(A) == 1 && return A[1]
  length(A) == 2 && return +(A[1], A[2]; kwargs...)
  nsize = isodd(length(A)) ? (div(length(A) - 1, 2) + 1) : div(length(A), 2)
  newterms = Vector{T}(undef, nsize)
  np = 1
  for n in 1:2:length(A) - 1
    newterms[np] = +(A[n], A[n+1]; kwargs...)
    np += 1
  end
  if isodd(length(A))
    newterms[nsize] = A[end]
  end
  return sum(newterms; kwargs...)
end

"""
    orthogonalize!(M::MPS, j::Int; kwargs...)

    orthogonalize!(M::MPO, j::Int; kwargs...)

Move the orthogonality center of the MPS
to site `j`. No observable property of the
MPS will be changed, and no truncation of the
bond indices is performed. Afterward, tensors
`1,2,...,j-1` will be left-orthogonal and tensors
`j+1,j+2,...,N` will be right-orthogonal.
"""
function orthogonalize!(M::AbstractMPS,
                        j::Int;
                        kwargs...)
  while leftlim(M) < (j-1)
    (leftlim(M) < 0) && setleftlim!(M, 0)
    b = leftlim(M)+1
    linds = uniqueinds(M[b],M[b+1])
    L,R = factorize(M[b], linds)
    M[b] = L
    M[b+1] *= R

    setleftlim!(M,b)
    if rightlim(M) < leftlim(M)+2
      setrightlim!(M, leftlim(M)+2)
    end
  end

  N = length(M)

  while rightlim(M) > (j+1)
    (rightlim(M) > (N+1)) && setrightlim!(M,N+1)
    b = rightlim(M)-2
    rinds = uniqueinds(M[b+1],M[b])
    L,R = factorize(M[b+1], rinds)
    M[b+1] = L
    M[b] *= R

    setrightlim!(M,b+1)
    if leftlim(M) > rightlim(M)-2
      setleftlim!(M, rightlim(M)-2)
    end
  end
end

"""
    truncate!(M::MPS; kwargs...)

    truncate!(M::MPO; kwargs...)

Perform a truncation of all bonds of an MPS/MPO,
using the truncation parameters (cutoff,maxdim, etc.)
provided as keyword arguments.
"""
function NDTensors.truncate!(M::AbstractMPS;
                             kwargs...)
  N = length(M)

  # Left-orthogonalize all tensors to make
  # truncations controlled
  orthogonalize!(M,N)

  # Perform truncations in a right-to-left sweep
  for j in reverse(2:N)
    rinds = uniqueinds(M[j],M[j-1])
    U,S,V = svd(M[j],rinds;kwargs...)
    M[j] = U
    M[j-1] *= (S*V)
    setrightlim!(M,j)
  end

end

NDTensors.contract(A::AbstractMPS,
                   B::AbstractMPS;
                   kwargs...) = *(A, B; kwargs...)

"""
    *(x::Number, M::MPS)

    *(x::Number, M::MPO)

Scale the MPS or MPO by the provided number.

Note: right now this just naively scales the
middle tensor. In the future, the plan would be
to scale the tensors between the left limit
and right limit by `x^(1/N)` where `N` is the distance
between the left limit and right limit.
"""
function Base.:*(x::Number, M::AbstractMPS)
  N = deepcopy(M)
  c = div(length(N), 2)
  N[c] .*= x
  return N
end

Base.:-(M::AbstractMPS) = Base.:*(-1,M)

"""
    hasqns(M::MPS)

    hasqns(M::MPO)

Return true if the MPS or MPO has
tensors which carry quantum numbers.
"""
hasqns(M::AbstractMPS) = hasqns(M[1])

"""
    flux(M::MPS)

    flux(M::MPO)

    totalqn(M::MPS)

    totalqn(M::MPO)

For an MPS or MPO which conserves quantum
numbers, compute the total QN flux. For
a tensor network such as an MPS or MPO,
the flux is the sum of fluxes of each of
the tensors in the network. The name
`totalqn` is an alias for `flux`.
"""
function flux(M::AbstractMPS)
  hasqns(M) || return nothing
  q = QN()
  for j=M.llim+1:M.rlim-1
    q += flux(M[j])
  end
  return q
end

totalqn(M::AbstractMPS) = flux(M)

function checkflux(M::AbstractMPS)
  for m in M
    checkflux(m)
  end
  return nothing
end

function Base.show(io::IO, M::AbstractMPS)
  print(io,"$(typeof(M))")
  (length(M) > 0) && print(io,"\n")
  for (i, A) ∈ enumerate(data(M))
    if order(A) != 0
      println(io,"[$i] $(inds(A))")
    else
      println(io,"[$i] ITensor()")
    end
  end
end

@deprecate orthoCenter(args...;
                       kwargs...) orthocenter(args...; kwargs...)

import .NDTensors.store

@deprecate store(m::AbstractMPS) data(m)

@deprecate replacesites!(args...;
                         kwargs...) ITensors.replace_siteinds!(args...; kwargs...)

@deprecate applyMPO(args...; kwargs...) contract(args...; kwargs...)

@deprecate applympo(args...; kwargs...) contract(args...; kwargs...)

@deprecate errorMPOprod(args...;
                        kwargs...) error_contract(args...;
                                                  kwargs...)

@deprecate error_mpoprod(args...;
                         kwargs...) error_contract(args...;
                                                   kwargs...)

@deprecate error_mul(args...;
                     kwargs...) error_contract(args...;
                                               kwargs...)

@deprecate multMPO(args...; kwargs...) contract(args...; kwargs...)

import Base.sum

@deprecate sum(A::AbstractMPS,
               B::AbstractMPS; kwargs...) add(A, B; kwargs...)

@deprecate multmpo(args...;
                   kwargs...) contract(args...; kwargs...)

@deprecate set_leftlim!(args...;
                        kwargs...) ITensors.setleftlim!(args...;
                                                        kwargs...)

@deprecate set_rightlim!(args...;
                         kwargs...) ITensors.setrightlim!(args...;
                                                          kwargs...)

@deprecate tensors(args...;
                   kwargs...) ITensors.data(args...; kwargs...)

@deprecate primelinks!(args...;
                       kwargs...) ITensors.prime_linkinds!(args...;
                                                          kwargs...)

@deprecate simlinks!(args...;
                     kwargs...) ITensors.sim_linkinds!(args...;
                                                      kwargs...)

@deprecate mul(A::AbstractMPS,
               B::AbstractMPS;
               kwargs...) contract(A, B; kwargs...)


text:
"""
Represents an aberration. This uses Zernike polynomials.

See also: [`Piston`](@ref), [`Tilt`](@ref), [`Tip`](@ref), [`Defocus`](@ref), [`ObliqueAstigmatism`](@ref), [`VerticalAstigmatism`](@ref), [`HorizontalComa`](@ref), [`VerticalComa`](@ref), [`ObliqueTrefoil`](@ref), [`VerticalTrefoil`](@ref), [`Spherical`](@ref), [`AddAberration`](@ref), [`Astigmatism`](@ref), [`Coma`](@ref), [`Trefoil `](@ref)
"""
abstract type Aberration end

"""
Piston aberration.
"""
struct Piston <: Aberration
    piston::Real
end
Piston() = Piston(1)

"""
Tilt aberration.
"""
struct Tilt <: Aberration
    tilt::Real
end
Tilt() = Tilt(1)

"""
Tip aberration.
"""
struct Tip <: Aberration
    tip::Real
end
Tip() = Tip(1)

"""
Defocus aberration.
"""
struct Defocus <: Aberration
    defocus::Real
end
Defocus() = Defocus(1)

"""
ObliqueAstigmatism aberration.
"""
struct ObliqueAstigmatism <: Aberration
    oblique::Real
end
ObliqueAstigmatism() = ObliqueAstigmatism(1)

"""
VerticalAstigmatism aberration.
"""
struct VerticalAstigmatism <: Aberration
    vertical::Real
end
VerticalAstigmatism() = VerticalAstigmatism(1)

"""
HorizontalComa aberration.
"""
struct HorizontalComa <: Aberration
    horizontal::Real
end
HorizontalComa() = HorizontalComa(1)

"""
VerticalComa aberration.
"""
struct VerticalComa <: Aberration
    vertical::Real
end
VerticalComa() = VerticalComa(1)

"""
ObliqueTrefoil aberration.
"""
struct ObliqueTrefoil <: Aberration
    oblique::Real
end
ObliqueTrefoil() = ObliqueTrefoil(1)

"""
VerticalTrefoil aberration.
"""
struct VerticalTrefoil <: Aberration
    vertical::Real
end
VerticalTrefoil() = VerticalTrefoil(1)

"""
Spherical aberration.
"""
struct Spherical <: Aberration
    spherical::Real
end
Spherical() = Spherical(1)
abstract type AbstractAddAberration <: Aberration end

"""
The addition of two aberrations. It is created when `+` operator is used on two
[`Aberration`](@ref).
"""
struct AddAberration <: AbstractAddAberration
    left::Aberration
    right::Aberration
end

"""
Astigmatism aberration.
"""
struct Astigmatism <: AbstractAddAberration
    left::ObliqueAstigmatism
    right::VerticalAstigmatism
end
Astigmatism(oblique::Real, vertical::Real) = Astigmatism(ObliqueAstigmatism(oblique), VerticalAstigmatism(vertical))
Astigmatism() = Astigmatism(ObliqueAstigmatism(), VerticalAstigmatism())

"""
Coma aberration.
"""
struct Coma <: AbstractAddAberration
    left::HorizontalComa
    right::VerticalComa
end
Coma(horizontal::Real, vertical::Real) = Coma(HorizontalComa(horizontal), VerticalComa(vertical))
Coma() = Coma(HorizontalComa(), VerticalComa())

"""
Trefoil aberration.
"""
struct Trefoil <: AbstractAddAberration
    left::ObliqueTrefoil
    right::VerticalTrefoil
end
Trefoil(oblique::Real, vertical::Real) = Trefoil(ObliqueTrefoil(oblique), VerticalTrefoil(vertical))
Trefoil() = Trefoil(ObliqueTrefoil(), VerticalTrefoil())

"""
    phase(aberration, coordinates=:cartesian)

Create a function mapping the phase of the `aberration` in the given `coordinates`
system (`:cartesian` or `:polar`).

See also: [`intensity`](@ref), [`amplitude`](@ref)
"""
phase(ab::Aberration, coordinates::Symbol=:cartesian) = phase(convert(Coordinates, coordinates), ab)
phase(t::Coordinates, ab::Aberration) = phase(typeof(t), ab::Aberration)
phase(::Type{Cartesian}, ab::Aberration) = begin
    ϕ = phase(Polar(), ab)
    (x,y) -> begin
        r,θ = cartesian_to_polar(x,y)
        ϕ(r,θ)
    end
end
phase(::Type{Polar}, ab::Piston) = (r,θ) -> ab.piston
phase(::Type{Cartesian}, ab::Piston) = (x,y) -> ab.piston
phase(::Type{Polar}, ab::Tilt) = (r,θ) -> ab.tilt*2r*sin(θ)
phase(::Type{Polar}, ab::Tip) = (r,θ) -> ab.tip*2r*cos(θ)
phase(::Type{Polar}, ab::Defocus) = (r,θ) -> ab.defocus * √3(2r^2-1)
phase(::Type{Polar}, ab::ObliqueAstigmatism) = (r,θ) -> ab.oblique * √6 * r^2 * sin(2θ)
phase(::Type{Polar}, ab::VerticalAstigmatism) = (r,θ) -> ab.vertical * √6 * r^2 * cos(2θ)
phase(::Type{Polar}, ab::HorizontalComa) = (r,θ) -> ab.horizontal * 2*√2*(3r^3-2r)cos(θ)
phase(::Type{Polar}, ab::VerticalComa) = (r,θ) -> ab.vertical * 2*√2*(3r^3-2r)sin(θ)
phase(::Type{Polar}, ab::ObliqueTrefoil) = (r,θ) -> ab.oblique * 2*√2*r^3 * cos(3θ)
phase(::Type{Polar}, ab::VerticalTrefoil) = (r,θ) -> ab.vertical * 2*√2*r^3 * sin(3θ)
phase(::Type{Polar}, ab::Spherical) = (r,θ) -> ab.spherical * √5(6r^4-6r^2+1)
phase(c::Type{Polar}, ab::AbstractAddAberration) = begin
    ϕ1 = phase(c, ab.left)
    ϕ2 = phase(c, ab.right)
    (r,θ) -> ϕ1(r,θ) + ϕ2(r,θ)
end

Base.:+(a::T, b::U) where {T<:Aberration, U<:Aberration} = AddAberration(a,b)

Base.:*(a::T, b::Piston) where {T<:Real} = Piston(a * b.piston)
Base.:*(a::T, b::Tilt) where {T<:Real} = Tilt(a * b.tilt)
Base.:*(a::T, b::Tip) where {T<:Real} = Tip(a * b.tip)
Base.:*(a::T, b::Defocus) where {T<:Real} = Defocus(a * b.defocus)
Base.:*(a::T, b::ObliqueAstigmatism) where {T<:Real} = ObliqueAstigmatism(a * b.oblique)
Base.:*(a::T, b::VerticalAstigmatism) where {T<:Real} = VerticalAstigmatism(a * b.vertical)
Base.:*(a::T, b::HorizontalComa) where {T<:Real} = HorizontalComa(a * b.horizontal)
Base.:*(a::T, b::VerticalComa) where {T<:Real} = VerticalComa(a * b.vertical)
Base.:*(a::T, b::ObliqueTrefoil) where {T<:Real} = ObliqueTrefoil(a * b.oblique)
Base.:*(a::T, b::VerticalTrefoil) where {T<:Real} = VerticalTrefoil(a * b.vertical)
Base.:*(a::T, b::Spherical) where {T<:Real} = Spherical(a * b.spherical)
Base.:*(a::T, b::U) where {T<:Real, U<:AbstractAddAberration} = U(a * b.left, a * b.right)

"""
    intensity(aberration, coordinates=:cartesian; mask=(x,y)->1)

Creates a function mapping the intensity in the given `coordinates` system
(`:cartesian` or `:polar`). This is equivalent of building two `Wavefront`, one
with no aberration and the other with the given `aberration`.

See also: [`phase`](@ref), [`amplitude`](@ref)
"""
intensity(ab::Aberration=Piston(), coordinates::Symbol=:cartesian; mask::Function=((a,b)->1)) = begin
    c = convert(Coordinates, coordinates)
    ϕ = phase(c, ab)
    (a,b) -> abs(1 + mask(a,b)*exp(1im*ϕ(a,b)))^2
end

"""
    project(aberration, data; mask=(x,y)->1)

Project the `data` onto the Zernike polynomial corresponding to the `aberration`
in order to get the fitting [`Aberration`](@ref). A `mask` can be used to fit only
in certain parts of the data (typically usefull if your interference pattern is
only relevant in a circular region). The coordinates are assumed to be normalized
such that (x,y) ∈ [-1;1]².

See also: [`correct`](@ref)
"""
project(ab::Type{<:Aberration}, data::Array{<:Real, 2}; mask::Function=((x,y)->1)) = project(ab(), data, mask=mask)
function project(ab::Aberration, data::Array{<:Real, 2}; mask::Function=((x,y)->1))
    ϕ = phase(Cartesian, ab)
    x,y = range(-1,1,length=size(data)[1]), range(-1,1,length=size(data)[2])
    X = repeat(reshape(x, 1, :), length(y), 1)
    Y = repeat(y, 1, length(x))
    mask_values = mask.(X,Y)
    ab_values = ϕ.(X,Y) .* mask_values
    sum(data .* ab_values) / sum(mask_values) * ab
end
function project(ab::AbstractAddAberration, data::Array{<:Real, 2}; mask::Function=((x,y)->1))
    typeof(ab)(project(ab.left, data, mask=mask),project(ab.right, data, mask=mask))
end

"""
    correct(aberration, data; mask=(x,y)->1)

Correct the given aberration in the data by [`project`](@ref)ing the data on the
aberration.

See also: [`project`](@ref)
"""
function correct(ab::Type{<:Aberration}, data::Array{<:Real, 2}; mask::Function=((x,y)->1))
    aberration = project(ab,data,mask=mask)
    ϕ = phase(Cartesian, aberration)
    x,y = range(-1,1,length=size(data)[1]), range(-1,1,length=size(data)[2])
    X = repeat(reshape(x, 1, :), length(y), 1)
    Y = repeat(y, 1, length(x))
    mask_values = mask.(X,Y)
    ab_values = ϕ.(X,Y) .* mask_values
    data - ab_values
end

Base.show(a::AddAberration) = begin
    Base.show(a.left)
    Base.show(" + ")
    Base.show(a.right)
end
Base.show(a::Astigmatism) = Base.show("Astigmatism($(a.left.oblique), $(a.right.vertical))")
Base.show(a::Trefoil) = Base.show("Astigmatism($(a.left.oblique), $(a.right.vertical))")
Base.show(a::Coma) = Base.show("Coma($(a.left.horizontal), $(a.right.vertical))")

ext: jl


text:
using Graphics
using Compat
using Compat.Test
using Compat.LinearAlgebra

@testset "Geometry" begin
    ## Point-vector identity (typealias)
    @test Point(2, 1) == Vec2(2, 1)

    ## Vector operators
    @test Vec2(1, 2) + Vec2(3, 4) == Vec2(4, 6)
    @test Vec2(1, 2) - Vec2(3, 4) == Vec2(-2, -2)
    @test Vec2(1, 2) * 1 == Vec2(1, 2) # identity
    @test Vec2(1, 2) * 0 == Vec2(0, 0) # null
    @test Vec2(1, 2) * 2 == Vec2(2, 4) # scalar multiplication
    @test Vec2(1, 2) / 2 == Vec2(0.5, 1)

    ## rotation: rotate()
    @test rotate(Vec2(1, 2), π, Vec2(1, 1)).x ≈ 1
    @test rotate(Vec2(1, 2), π, Vec2(1, 1)).y ≈ 0
    @test rotate(Vec2(1, 2), 2π).x ≈ 1
    @test rotate(Vec2(1, 2), 2π).y ≈ 2

    ## Euclidean norm/magnitude: norm()
    @test norm(Vec2(1, 2)) == sqrt(5)
    @test norm(Vec2(0, 0)) == 0

    ## Bounding boxes: BoundingBox
    BBT_point_1 = Vec2(1, 1)
    BBT_point_2 = Vec2(2, 3)
    BBT_point_3 = Vec2(4, 5)
    BBT = BoundingBox(BBT_point_1, BBT_point_2, BBT_point_3)

    ### BoundingBox attributes
    @test BBT == BoundingBox(1, 4, 1, 5)
    @test height(BBT) == 4
    @test width(BBT) == 3
    @test diagonal(BBT) == 5
    @test aspect_ratio(BBT) ≈ 1.333333333333
    @test xmin(BBT) == 1
    @test xmax(BBT) == 4
    @test center(BBT) == Vec2(2.5, 3)
    @test xrange(BBT) == (1, 4)
    @test yrange(BBT) == (1, 5)

    ### BoundingBox operations
    BBT_1 = BoundingBox(2, 3, 4, 5)
    BBT_2 = BoundingBox(6, 7, 8, 9)
    BBT_3 = BoundingBox(NaN, NaN, NaN, NaN)
    BBT_4 = BoundingBox(-Inf, -Inf, -Inf, -Inf)

    #### BoundingBox (+)
    @test BBT_1 + BBT_2 == BoundingBox(2, 7, 4, 9)
    @test BBT_3 + BBT_4 == BoundingBox(-Inf, -Inf, -Inf, -Inf)

    #### BoundingBox (&)
    @test BBT_1 & BBT_2 == BoundingBox(6, 3, 8, 5)
    @test BBT_3 & BBT_4 == BoundingBox(-Inf, -Inf, -Inf, -Inf)

    #### deform()
    @test deform(BBT_1, -1, 2, -3, 4) == BoundingBox(1, 5, 1, 9)

    #### shift()
    @test shift(BBT_1, -1, 2) == BoundingBox(1, 2, 6, 7)

    #### scale()
    @test BBT_1 * 3 == BoundingBox(1, 4, 3, 6)

    #### rotate()
    BBT_3 = rotate(BBT_1, π, Point(0, 1))
    @test BBT_3.xmin ≈ -3
    @test BBT_3.ymin ≈ -3
    @test BBT_3.xmax ≈ -2
    @test BBT_3.ymax ≈ -2

    #### with_aspect_ratio()
    @test with_aspect_ratio(BBT_1, 2) == BoundingBox(2.25, 2.75, 4, 5)
    @test with_aspect_ratio(BBT_1, 0.5) == BoundingBox(2, 3, 4.25, 4.75)

    #### isinside()
    @test isinside(BBT_1, Point(2.5, 4.5))
    @test isinside(BBT_1, Point(2, 4))
    @test isinside(BBT_1, Point(1, 3)) == false
end

# using Gtk.ShortNames, Cairo

# @testset "Canvas" begin
#     win = Window() |> (cnvs = Canvas(300, 280))
#     draw(cnvs) do c
#         set_coords(getgc(cnvs), BoundingBox(0, 1, 0, 1))
#     end
#     showall(win)
#     sleep(0.5)
#     mtrx = Cairo.get_matrix(getgc(cnvs))
#     @test mtrx.xx == 300
#     @test mtrx.yy == 280
#     @test mtrx.xy == mtrx.yx == mtrx.x0 == mtrx.y0 == 0
# end


text:
import Term.Measures: Measure, width, height
import Term: remove_markup, Panel

@testset "\e[34mMeasure - str" begin
    @test size(Measure("a"^10)) == (10, 1)

    strings = (
        "asadasda"^2,
        "~~±±||__ASdaSDvxcvxc\nasdfsiaudfhsdiufhndskjv",
        "dasdadsa\nsdfsd"^7,
        "asdasdsa|dfvcxvashfusdn\nfidsuhfsdf",
    )
    for string in strings
        m = Measure(string)
        @test m.w == lw(string)
        @test m.h == nlines(string)
    end

    strings = (
        "[red]is my color"^2,
        "~~±±||__{red on_green}aasdas{bold}asdsad{/bold}asdas{/red on_green}xc\nasdfsiaudfhsdiufhndskjv",
        "dasda{#ffffff}dsa{#ffffff}\nsdfsd"^7,
        "asdasdsa|dfvcxvashf{bold}usdn\nfid[/bold]suhfsdf",
    )
    for string in strings
        m = Measure(string)
        @test m.w == lw(remove_markup(string))
        @test m.h == nlines(remove_markup(string))
    end
end

# TODO test with renderables

@testset "\r34mMeasure funcs" begin
    @test width("test") == 4
    @test height("test") == 1

    w = h = 5
    p = Panel(width = w, height = h)
    @test width(p) == w
    @test height(p) == h

    s = first(p.segments)
    @test width(s) == w
    @test height(s) == 1
end


text:
###########################
# An n-dimensional simplex
###########################

abstract type Simplex{T,C} <: Domain{T} end

isclosedset(::Simplex{T,:closed}) where {T} = true
isclosedset(::Simplex{T,:open}) where {T} = false

isopenset(d::Simplex) = !isclosedset(d)

abstract type UnitSimplex{T,C} <: Simplex{T,C} end
const ClosedUnitSimplex{T} = UnitSimplex{T,:closed}
const OpenUnitSimplex{T} = UnitSimplex{T,:open}

UnitSimplex(n::Int) = DynamicUnitSimplex(n)
UnitSimplex(::Val{N}) where {N} = EuclideanUnitSimplex{N}()

UnitSimplex{T}(n::Int) where {T <: StaticTypes} = StaticUnitSimplex{T}(n)
UnitSimplex{T}(::Val{N}) where {N,T} = StaticUnitSimplex{T}(Val(N))
UnitSimplex{T}() where {T <: StaticTypes} = StaticUnitSimplex{T}()
UnitSimplex{T}(n::Int) where {T} = DynamicUnitSimplex{T}(n)

UnitSimplex{T,C}(n::Int) where {T <: StaticTypes,C} = StaticUnitSimplex{T,C}(n)
UnitSimplex{T,C}(::Val{N}) where {N,T,C} = StaticUnitSimplex{T,C}(Val(N))
UnitSimplex{T,C}() where {T <: StaticTypes,C} = StaticUnitSimplex{T,C}()
UnitSimplex{T,C}(n::Int) where {T,C} = DynamicUnitSimplex{T,C}(n)

insimplex_closed(x) = mapreduce( t-> t >= 0, &, x) && norm(x,1) <= 1
insimplex_open(x) = mapreduce( t-> t > 0, &, x) && norm(x,1) < 1
insimplex_closed(x, tol) = mapreduce( t-> t >= -tol, &, x) && norm(x,1) <= 1+tol
insimplex_open(x, tol) = mapreduce( t-> t > -tol, &, x) && norm(x,1) < 1+tol

indomain(x, ::ClosedUnitSimplex) = insimplex_closed(x)
indomain(x, ::OpenUnitSimplex) = insimplex_open(x)
approx_indomain(x, ::ClosedUnitSimplex, tolerance) = insimplex_closed(x, tolerance)
approx_indomain(x, ::OpenUnitSimplex, tolerance) = insimplex_open(x, tolerance)

isempty(::UnitSimplex) = false

==(d1::UnitSimplex, d2::UnitSimplex) =
    isclosedset(d1)==isclosedset(d2) && dimension(d1)==dimension(d2)

boundingbox(d::UnitSimplex{T}) where {T} = UnitCube{T}(dimension(d))

struct StaticUnitSimplex{T,C} <: UnitSimplex{T,C}
end

StaticUnitSimplex(::Val{N}) where {N} = StaticUnitSimplex{SVector{N,Float64}}()

StaticUnitSimplex{T}() where {T} = StaticUnitSimplex{T,:closed}()

StaticUnitSimplex{T}(n::Int) where {T} =
    (@assert n == euclideandimension(T); StaticUnitSimplex{T}())
StaticUnitSimplex{T}(::Val{N}) where {N,T} =
    (@assert N == euclideandimension(T); StaticUnitSimplex{T}())

StaticUnitSimplex{T,C}(n::Int) where {T,C} =
    (@assert n == euclideandimension(T); StaticUnitSimplex{T,C}())
StaticUnitSimplex{T,C}(::Val{N}) where {N,T,C} =
    (@assert N == euclideandimension(T); StaticUnitSimplex{T,C}())

const EuclideanUnitSimplex{N,T,C} = StaticUnitSimplex{SVector{N,T},C}

EuclideanUnitSimplex{N}() where {N} = EuclideanUnitSimplex{N,Float64}()

struct DynamicUnitSimplex{T,C} <: UnitSimplex{T,C}
    dimension   ::  Int

    DynamicUnitSimplex{T,C}(n::Int) where {T,C} = new(n)
    DynamicUnitSimplex{T,C}(n::Int) where {T<:StaticTypes,C} =
        (@assert n == euclideandimension(T); new(n))
end

DynamicUnitSimplex(n::Int) = DynamicUnitSimplex{Vector{Float64}}(n)
DynamicUnitSimplex{T}(n::Int) where {T} = DynamicUnitSimplex{T,:closed}(n)

dimension(d::DynamicUnitSimplex) = d.dimension

const VectorUnitSimplex{T,C} = DynamicUnitSimplex{Vector{T},C}

VectorUnitSimplex(dimension) = VectorUnitSimplex{Float64}(dimension)

show(io::IO, d::EuclideanUnitSimplex{N,Float64,:closed}) where {N} = print(io, "UnitSimplex(Val($(N)))")
show(io::IO, d::VectorUnitSimplex{Float64,:closed}) = print(io, "UnitSimplex($(dimension(d)))")

center(d::EuclideanUnitSimplex{N,T}) where {N,T} = ones(SVector{N,T})/N
center(d::VectorUnitSimplex{T}) where {T} = ones(T, dimension(d))/dimension(d)

corners(::EuclideanUnitSimplex{N,T}) where {N,T} =
    vcat([zero(SVector{N,T})], [ SVector{N,T}(ntuple( i -> i==j, N)) for j in 1:N])
corners(d::VectorUnitSimplex{T}) where {T} =
    vcat([zeros(T,dimension(d))], [ (z = zeros(T, dimension(d)); z[j]=1; z) for j in 1:dimension(d)])

interior(d::EuclideanUnitSimplex{N,T}) where {N,T} = EuclideanUnitSimplex{N,T,:open}()
closure(d::EuclideanUnitSimplex{N,T}) where {N,T} = EuclideanUnitSimplex{N,T,:closed}()
interior(d::VectorUnitSimplex{T}) where {T} = VectorUnitSimplex{T,:open}(dimension(d))
closure(d::VectorUnitSimplex{T}) where {T} = VectorUnitSimplex{T,:closed}(dimension(d))

# We pick the center point, because it belongs to the domain regardless of
# whether it is open or closed.
point_in_domain(d::UnitSimplex) = center(d)

similardomain(d::StaticUnitSimplex{S,C}, ::Type{T}) where {S,T,C} =
    StaticUnitSimplex{T,C}()
similardomain(d::DynamicUnitSimplex{S,C}, ::Type{T}) where {S,T,C} =
    DynamicUnitSimplex{T,C}(d.dimension)


text:
# Default implementations of generative model functions

# This file is long, messy and fragile. One day a more robust paradigm should be used.

function implemented(f::typeof(generate_s), TT::Type)
    if !hasmethod(f, TT)
        return false
    end
    m = which(f, TT)
    if m.module == POMDPs && !implemented(transition, Tuple{TT.parameters[1:end-1]...})
        return false
    else # a more specific implementation exists
        return true
    end
end

@generated function generate_s(p::Union{MDP,POMDP}, s, a, rng::AbstractRNG)
    impl = quote
        td = transition(p, s, a)
        return rand(rng, td)
    end
    if implemented(transition, Tuple{p, s, a})
        return impl
    else
        treq = @req transition(::p,::s,::a)
        reqs = [(implemented(treq...), treq...)]
        this = @req(generate_s(::p, ::s, ::a, ::rng))
        return quote
            try
                $impl # trick the compiler to put the right backedges in
            catch
                failed_synth_warning($this, $reqs)
                throw(MethodError(generate_s, (p,s,a,rng)))
            end
        end
    end
end

function implemented(f::typeof(generate_sr), TT::Type)
    if !hasmethod(f, TT)
        return false
    end
    m = which(f, TT)
    reqs_met = implemented(generate_s, TT) && implemented(reward, Tuple{TT.parameters[1:end-1]..., TT.parameters[2]})
    if m.module == POMDPs && !reqs_met
        return false
    else # a more specific implementation exists
        return true
    end
end

@generated function generate_sr(p::Union{POMDP,MDP}, s, a, rng::AbstractRNG)
    impl = quote
        sp = generate_s(p, s, a, rng)
        return sp, reward(p, s, a, sp)
    end

    if implemented(generate_s, Tuple{p, s, a, rng}) && implemented(reward, Tuple{p, s, a, s})
        return impl
    else
        reqs = [@req(transition(::p,::s,::a)), @req(reward(::p, ::s, ::a, ::s))]
        cl = [(implemented(r...), r...) for r in reqs]
        greqs = [@req(generate_s(::p, ::s, ::a, ::AbstractRNG)), @req(reward(::p, ::s, ::a, ::s))]
        gcl = [(implemented(r...), r...) for r in greqs]
        this = @req(generate_sr(::p, ::s, ::a, ::rng))
        return quote
            try
                $impl # trick to get the compiler to put the right backedges in
            catch
                failed_synth_warning($this, $cl, $gcl)
                throw(MethodError(generate_sr, (p,s,a,rng)))
            end
        end
    end
end

function implemented(f::typeof(generate_o), TT::Type)
    if !hasmethod(f, TT)
        return false
    end
    m = which(f, TT)
    if m.module == POMDPs && !implemented(observation, Tuple{TT.parameters[1:end-1]...})
        return false
    else # a more specific implementation exists
        return true
    end
end

@generated function generate_o(p::POMDP, s, a, sp, rng::AbstractRNG)
    impl = quote
        od = observation(p, s, a, sp)
        return rand(rng, od)
    end

    if implemented(observation, Tuple{p, s, a, sp})
        return impl
    else
        oreq = @req observation(::p, ::s, ::a, ::sp)
        reqs = [(implemented(oreq...), oreq...)]
        this = @req(generate_o(::p, ::s, ::a, ::sp, ::rng))
        return quote
            try
                $impl # trick to get the compiler to put the right backedges in
            catch
                failed_synth_warning($this, $reqs)
                throw(MethodError(generate_o, (p, s, a, sp, rng)))
            end
        end
    end
end

function implemented(f::typeof(generate_so), TT::Type)
    if !hasmethod(f, TT)
        return false
    end
    m = which(f, TT)
    reqs_met = implemented(generate_s, TT) && implemented(generate_o, Tuple{TT.parameters[1:end-1]..., TT.parameters[2], TT.parameters[end]})
    if m.module == POMDPs && !reqs_met
        return false
    else # a more specific implementation exists
        return true
    end
end

@generated function generate_so(p::POMDP, s, a, rng::AbstractRNG)
    impl = quote
        sp = generate_s(p, s, a, rng)
        return sp, generate_o(p, s, a, sp, rng)
    end

    if implemented(generate_s, Tuple{p, s, a, rng}) && implemented(generate_o, Tuple{p, s, a, s, rng})
        return impl
    else
        reqs = [@req(transition(::p,::s,::a)), @req(observation(::p, ::s, ::a, ::s))]
        cl = [(implemented(r...), r...) for r in reqs]
        greqs = [@req(generate_s(::p,::s,::a,::AbstractRNG)), @req(generate_o(::p, ::s, ::a, ::s, ::AbstractRNG))]
        gcl = [(implemented(r...), r...) for r in greqs]
        this = @req(generate_so(::p, ::s, ::a, ::rng))
        return quote
            try
                $impl # trick to get the compiler to put the right backedges in
            catch
                failed_synth_warning($this, $cl, $gcl)
                throw(MethodError(generate_so, (p,s,a,rng)))
            end
        end
    end
end

function implemented(f::typeof(generate_sor), TT::Type)
    if !hasmethod(f, TT)
        return false
    end
    m = which(f, TT)
    so_reqs_met = implemented(generate_so, TT) && implemented(reward, Tuple{TT.parameters[1:end-1]..., TT.parameters[2]})
    sr_reqs_met = implemented(generate_sr, TT) && implemented(generate_o, Tuple{TT.parameters[1:end-1]..., TT.parameters[2], TT.parameters[end]})
    if m.module == POMDPs && !so_reqs_met && !sr_reqs_met
        return false
    else # a more specific implementation exists
        return true
    end
end

@generated function generate_sor(p::POMDP, s, a, rng::AbstractRNG)
    so_impl = quote
        sp, o = generate_so(p, s, a, rng)
        sp, o, reward(p, s, a, sp)
    end

    sr_impl = quote
        sp, r = generate_sr(p, s, a, rng)
        o = generate_o(p, s, a, sp, rng)
        sp, o, r
    end

    if implemented(generate_so, Tuple{p, s, a, rng}) && implemented(reward, Tuple{p, s, a, s})
        return so_impl
    elseif implemented(generate_sr, Tuple{p, s, a, rng}) && implemented(generate_o, Tuple{p, s, a, s, rng})
        return sr_impl
    else
        reqs = [@req(transition(::p,::s,::a)), @req(observation(::p, ::s, ::a, ::s)), @req(reward(::p,::s,::a,::s))]
        cl = [(implemented(r...), r...) for r in reqs]
        greqs = [@req(generate_sr(::p,::s,::a,::AbstractRNG)), @req(generate_o(::p, ::s, ::a, ::s, ::AbstractRNG))]
        gcl = [(implemented(r...), r...) for r in greqs]
        this = @req(generate_sor(::p, ::s, ::a, ::rng))
        return quote
            try # dirty trick to get the compiler to insert the right backedges
                $so_impl
                $sr_impl
            catch
                failed_synth_warning($this, $cl, $gcl)
                throw(MethodError(generate_sor, (p,s,a,rng)))
            end
        end
    end
end

function implemented(f::typeof(generate_or), TT::Type)
    if !hasmethod(f, TT)
        return false
    end
    m = which(f, TT)
    reqs_met = implemented(generate_o, TT) && implemented(reward, Tuple{TT.parameters[1:end-1]..., TT.parameters[2]})
    if m.module == POMDPs && !reqs_met
        return false
    else # a more specific implementation exists
        return true
    end
end

@generated function generate_or(p::POMDP, s, a, sp, rng::AbstractRNG)
    impl = quote
        o = generate_o(p, s, a, sp, rng)
        return o, reward(p, s, a, sp)
    end

    if implemented(generate_o, Tuple{p, s, a, sp, rng}) && implemented(reward, Tuple{p, s, a, sp})
        return impl
    else
        reqs = [@req(observation(::p, ::s, ::a, ::sp)), @req(reward(::p,::s,::a,::sp))]
        cl = [(implemented(r...), r...) for r in reqs]
        greqs = [@req(generate_o(::p, ::s, ::a, ::s, ::AbstractRNG)), @req(reward(::p,::s,::a,::sp))]
        gcl = [(implemented(r...), r...) for r in greqs]
        this = @req(generate_or(::p, ::s, ::a, ::rng))
        return quote
            try
                $impl # trick to get the compiler to insert the right backedges
            catch
                failed_synth_warning($this, $cl, $gcl)
                throw(MethodError(generate_or, (p,s,a,sp,rng)))
            end
        end
    end
end

function implemented(f::typeof(initialstate), TT::Type)
    if !hasmethod(f, TT)
        return false
    end
    m = which(f, TT)
    if m.module == POMDPs && !implemented(initialstate_distribution, Tuple{TT.parameters[1]})
        return false
    else
        return true
    end
end

@generated function initialstate(p::Union{POMDP,MDP}, rng::AbstractRNG)
    impl = quote
        d = initialstate_distribution(p)
        return rand(rng, d)
    end

    if implemented(initialstate_distribution, Tuple{p})
        return impl
    else
        req = @req initialstate_distribution(::p)
        reqs = [(implemented(req...), req...)]
        this = @req(initialstate(::p, ::rng))
        return quote
            try
                $impl # trick to get the compiler to insert the right backedges
            catch
                failed_synth_warning($this, $reqs)
                throw(MethodError(initialstate, (p, rng)))
            end
        end
    end
end

function failed_synth_warning(gen::Tuple, reqs::Vector, greqs::Vector=[])
    io = IOBuffer()
    ioc = IOContext(io, stderr)
    print(ioc, "Hint: Either implement ")
    printstyled(ioc, format_method(gen...), color=:blue)
    println(ioc, " directly, or, to automatically synthesize it, implement the following methods from the explicit interface:\n")
    show_checked_list(ioc, reqs)
    if !isempty(greqs)
        println(ioc, """
    \nOR implement the following methods from the generative interface:
    """)
        show_checked_list(ioc, greqs)
    end
    println(ioc, "\n([✔] = already implemented correctly; [X] = missing)")
    @warn("""
          POMDPs.jl: Could not find or synthesize $(format_method(gen...)). Note that this is only a warning to help diagnose problems; consider fixing errors first.

          $(String(take!(io)))
          """)
end


text:
################################################################################
# pwdict.jl
#
# The Persistent weighted collection PWDict which is built on top of the
# persistent heap type in PHeap.jl# transient heap type.
#
# @author Noah C. Benson
#
# MIT License
# Copyright (c) 2020 Noah C. Benson

# ==============================================================================
# AbstractPWDict
"""
    AbstractPWDict

AbstractPWDict is a subtype of AbstractPDict that is overloaded by classes that
represent weighted persistent dictionaries.

See also: [`AbstractPDict`](@ref), [`PWDict`](@ref), [`AbstractPWSet`](@ref).
"""
abstract type AbstractPWDict{K,V,W<:Number} <: AbstractPDict{K,V} end

# ==============================================================================
# PWDict
_to_pwtup(p::Tuple{K,V,W}) where {K,V,W<:Number} = p
_to_pwtup(p::Tuple{K,Tuple{V,W}}) where {K,V,W<:Number} = begin
    (k,(v,w)) = p
    return (k,v,w)
end
_to_pwtup(p::Pair{K,Tuple{V,W}}) where {K,V,W<:Number} = begin
    (k,(v,w)) = p
    return (k,v,w)
end
_to_pwtup(p::Pair{K,Pair{V,W}}) where {K,V,W<:Number} = begin
    (k,(v,w)) = p
    return (k,v,w)
end
_to_pwtup(p::Tuple{Tuple{K,V},W}) where {K,V,W<:Number} = begin
    ((k,v),w) = p
    return (k,v,w)
end
_to_pwtup(p::Tuple{Pair{K,V},W}) where {K,V,W<:Number} = begin
    ((k,v),w) = p
    return (k,v,w)
end
_to_pwtup(p::Pair{Pair{K,V},W}) where {K,V,W<:Number} = begin
    ((k,v),w) = p
    return (k,v,w)
end
const PWDictHeap{K,W,D} = PHeap{K,W,typeof(>),D} where {K,W<:Number,D<:AbstractPDict{K,Int}}
macro _pwdict_code(name::Symbol, dicttype::Symbol)
    eq = gensym()
    h = gensym()
    return quote
        struct $name{K,V,W<:Number} <: AbstractPWDict{K,V,W}
            heap::PWDictHeap{K,W,$dicttype{K,Int}}
            dict::$dicttype{K,V}
        end
        function $name{K,V,W}() where {K,V,W<:Number}
            return $name{K,V,W}(PWDictHeap{K,W,$dicttype{K,Int}}(>),
                                $dicttype{K,V}())
        end
        function $name{K,V,W}(d::$name{K,V,W}) where {K,V,W}
            return d
        end
        function $name{K,V,W}(ps::Union{Tuple,Pair}...) where {K,V,W<:Number}
            return reduce(push, ps, init=$name{K,V,W}())
        end
        function $name{K,V}(::Tuple{}) where {K,V}
            return $name{K,V,Float64}(PWDictHeap{K,Float64,$dicttype{K,Int}}(>),
                                      $dicttype{K,V}())
        end
        function $name{K,V}() where {K,V}
            return $name{K,V,Float64}(PWDictHeap{K,Float64,$dicttype{K,Int}}(>),
                                      $dicttype{K,V}())
        end
        function $name{K,V}(d::$name{K,V,W}) where {K,V,W<:Number}
            return d
        end
        function $name{K,V}(ps::Union{Tuple,Pair}...) where {K,V}
            return $name{K,V,Float64}(ps...)
        end
        function $name(d::$name{K,V,W}) where {K,V,W<:Number}
            return d
        end
        function $name(ps::Union{Tuple,Pair}...)
            return $name{Any,Any,Float64}(ps...)
        end
        function $name(::Tuple{})
            return $name{Any,Any,Float64}()
        end
        function $name()
            return $name{Any,Any,Float64}(PWDictHeap{Any,Float64,$dicttype{Any,Int}}(),
                                          $dicttype{Any,Any}())
        end
        # Document the equal/hash types.
        equalfn(::Type{T}) where {T <: $name} = equalfn($dicttype)
        hashfn(::Type{T}) where {T <: $name} = hashfn($dicttype)
        equalfn(::$name) = equalfn($dicttype)
        hashfn(::$name) = hashfn($dicttype)
        # Generic construction.
        # Base methods.
        #Base.empty(s::$name{K,V,W}, ::Type{J}=K, ::Type{U}=V, ::Type{X}=W) where {K,V,W,J,U,X} = $name{J,U,X}()
        Base.length(s::$name) = length(s.heap)
        Base.iterate(u::$name{K,V,W}) where {K,V,W} = begin
            (length(u) == 0) && return nothing
            return iterate(u, u.heap)
        end
        Base.iterate(u::$name{K,V,W}, x) where {K,V,W} = begin
            nxt = iterate(u.heap, x)
            (nxt === nothing) && return nothing
            (k,b) = nxt
            return (k => u.dict[k], b)
        end
        Base.get(u::$name{K,V,W}, k, df) where {K,V,W} = get(u.dict, k, df)
        Base.in(kv::Pair, u::$name{K,V,W}, eqfn::Function) where {K,V,W} = in(kv, u.dict, eqfn)
        Base.in(kv::Pair, u::$name{K,V,W}) where {K,V,W} = in(kv, u.dict)
        # And the Air methods.
        Air.push(u::$name{K,V,W}, x::Pair{J,U}) where {K,V,W,J,U} = begin
            (k,v,w) = _to_pwtup(x)
            heap = push(u.heap, (k,w))
            dict = push(u.dict, k => v)
            (heap === u.heap && dict === u.dict) && return u
            return $name{K,V,W}(heap, dict)
        end
        Air.push(u::$name{K,V,W}, x::Tuple) where {K,V,W,J,U} = begin
            (k,v,w) = _to_pwtup(x)
            heap = push(u.heap, (k,w))
            dict = push(u.dict, k => v)
            (heap === u.heap && dict === u.dict) && return u
            return $name{K,V,W}(heap, dict)
        end
        Air.pop(u::$name{K,V,W}) where {K,V,W} = begin
            (length(u) == 0) && throw(ArgumentError("PWDict must be non-empty"))
            x = first(u.heap)
            dict = delete(u.dict, x)
            return $name{K,V,W}(pop(u.heap), dict)
        end
        Air.first(u::$name{K,V,W}) where {K,V,W} = begin
            (length(u) == 0) && throw(ArgumentError("PWDict must be non-empty"))
            k = first(u.heap)
            return k => u.dict[k]
        end
        Air.setindex(u::$name{K,V,W}, v, k::J) where {K,V,W,J} = begin
            (get(u.heap._index, k, 0) == 0) && throw(
                ArgumentError("when inserting into PWDict, value and weight are required"))
            d = setindex(u.dict, v, k)
            (d === u.dict) && return u
            return $name{K,V,W}(u.heap, d)
        end
        Air.setindex(u::$name{K,V,W}, v::Tuple{VV,WW}, k::J) where {K,V,W,VV,WW<:Number,J} = begin
            return push(u, k => v)
        end
        Air.setindex(u::$name{K,V,W}, v::Pair{VV,WW}, k::J) where {K,V,W,VV,WW<:Number,J} = begin
            return push(u, k => v)
        end
        Air.delete(u::$name{K,V,W}, k::J) where {K,V,W,J} = begin
            heap = delete(u.heap, k)
            (heap === u.heap) && return u
            dict = delete(u.dict, k)
            return $name{K,V,W}(heap, dict)
        end
        Air.getweight(u::$name{K,V,W}, x::J) where {K,V,W,J} = getweight(u.heap, x)
        Air.setweight(u::$name{K,V,W}, x::J, w::X) where {K,V,W,J,X<:Number} = begin
            heap = setweight(u.heap, x, w)
            (heap === u.heap) && return u
            return $name{K,V,W}(heap, u.dict)
        end
        Random.rand(u::$name{K,V,W}) where {K,V,W} = begin
            (length(u) > 0) || throw(ArgumentError("PWDict must be non-empty"))
            k = rand(u.heap)
            return k => u.dict[k]
        end
    end |> esc
end
# Declare the types.
@_pwdict_code PWDict PDict
@_pwdict_code PWIdDict PIdDict
# And document them.
@doc """
    PWDict{K,V}

A persistent dictionary with weighted pairs. As such, a `PWDict` supports the
typical operations of a `PDict` as well as the following:
* The `first` function yields the key-value pair with the highest weight.
* The `pop` function yields a copy of the dictionary without the pair that has
  the highest weight.
* Iteration occurs in the order of greatest to least weight.
* The weights can be changed with the `getweight` and `setweight` functions;
  `setweight` yields a duplicate dictionary with updated weights.

See also: [`PDict`](@ref), [`getweight`](@ref), [`setweight`](@ref).

# Examples

```@meta
DocTestSetup = quote
    using Air
end
```

```jldoctest; filter=r"PWDict{Symbol, ?Int64, ?Float64} with 3 entries:"
julia> PWDict{Symbol,Int}(:a => (1,0.1), :b => (2, 0.2), :c => (3, 0.3))
PWDict{Symbol,Int64,Float64} with 3 entries:
  :c => 3
  :b => 2
  :a => 1
```
""" PWDict
@doc """
    PWIdDict{K,V}

A persistent dictionary with weighted pairs. As such, a `PWIdDict` supports the
typical operations of a `PIdDict` as well as the following:
* The `first` function yields the key-value pair with the highest weight.
* The `pop` function yields a copy of the dictionary without the pair that has
  the highest weight.
* Iteration occurs in the order of greatest to least weight.
* The weights can be changed with the `getweight` and `setweight` functions;
  `setweight` yields a duplicate dictionary with updated weights.

See also: [`PIdDict`](@ref), [`PWDict`](@ref), [`getweight`](@ref),
[`setweight`](@ref).

# Examples

```@meta
DocTestSetup = quote
    using Air
end
```

```jldoctest; filter=r"PWIdDict{Symbol, ?Int64, ?Float64} with 3 entries:"
julia> PWIdDict{Symbol,Int}(:a => (1,0.1), :b => (2, 0.2), :c => (3, 0.3))
PWIdDict{Symbol,Int64,Float64} with 3 entries:
  :c => 3
  :b => 2
  :a => 1
```
""" PWIdDict
# For some reason, these need to be out of the macro above.
Base.empty(s::PWDict{K,V,W}, ::Type{J}=K, ::Type{U}=V, ::Type{X}=W) where {K,V,W,J,U,X} = PWDict{J,U,X}()
Base.empty(s::PWIdDict{K,V,W}, ::Type{J}=K, ::Type{U}=V, ::Type{X}=W) where {K,V,W,J,U,X} = PWIdDict{J,U,X}()
# Equality functions also.
Base.isequal(s::DS, t) where {DS <: AbstractPWDict} = false
Base.isequal(t, s::DS) where {DS <: AbstractPWDict} = false
Base.isequal(s::DS, t::DT) where {
    KT,VT,KS,VS,
    DS <: AbstractPWDict{KS,VS},
    DT <: AbstractPWDict{KT,VT}} = begin
    (length(t) == length(s)) || return false
    for (k,v) in s
        tt = get(t, k, t)
        (tt === t) && return false
        isequal(tt, v) || return false
        (getweight(s, k) == getweight(t, k)) || return false
    end
    (equalfn(DS) === equalfn(DT)) && return true
    for (k,v) in t
        ss = get(s, k, s)
        (ss === s) && return false
        isequal(ss, v) || return false
    end
    return true
end
Base.hash(u::PWD) where {K,V,W,PWD<:AbstractPWDict{K,V,W}} = begin
    h = length(u)
    for (k,v) in u.dict
        w = getweight(u.heap, k)
        h += hash(k) ⊻ hash(v) ⊻ hash(w)
    end
    return h
end

# Export the relevant symbols.
export PWDict, PWIdDict


text:
function int_p_pprime(mm,nn,q,z,N::Integer=matsize4([mm;nn],q))
    P = per_coeffs([mm;nn],q,N)
    map(product(1:length(mm),1:length(nn))) do x
        m = mm[x[1]]
        n = nn[x[2]]
        Pm = P[:,x[1]]
        Pn = P[:,length(mm)+x[2]]
        res = zero(z)
        if iseven(m)&&iseven(n)
            for k=0:N-1, l=0:N-1
                if k==l
                    @. res += Pm[k+1]*Pn[l+1]*cos(4k*z)/4
                else
                    @. res += Pm[k+1]*Pn[l+1]*0.5l*( cos(2(k+l)*z)/(k+l) - cos(2(k-l)*z)/(k-l) )
                end
            end
        elseif iseven(m)&&isodd(n)
            for k=0:N-1, l=0:N-1
                if k==l+1
                    @. res += Pm[k+1]*Pn[l+1]*( k*z + sin(4k*z)/4 )
                else
                    @. res += Pm[k+1]*Pn[l+1]*0.5(l+1)*( sin(2(k-l-1)*z)/(k-l-1) + sin(2(k+l+1)*z)/(k+l+1) )
                end
            end
        elseif isodd(m)&&iseven(n)
            for k=0:N-1, l=0:N-1
                if k+1==l
                    @. res += Pm[k+1]*Pn[l+1]*( sin(4(k+1)*z)/4 - (k+1)*z )
                else
                    @. res += Pm[k+1]*Pn[l+1]*0.5l*( sin(2(k+l+1)*z)/(k+l+1) - sin(2(k-l+1)*z)/(k-l+1) )
                end
            end
        else
            for k=0:N-1, l=0:N-1
                if k==l
                    @. res -= Pm[k+1]*Pn[l+1]*cos(4(k+1)*z)/4
                else
                    @. res -= Pm[k+1]*Pn[l+1]*0.5(l+1)*( cos(2(k-l)*z)/(k-l) + cos(2(k+l+2)*z)/(k+l+2) )
                end
            end
        end
        return (res[2]-res[1])
    end
end


text:
"""
# Method p51_skyline

Plane or axisymmetric strain analysis of an elastic solid (Plane
structural element) using 3-, 6-, 10- or 15-node right-angled triangles
(Triangle finite elements) or 4-, 8- or 9-node rectangular quadrilaterals
(Quadrilateral finite elements). Mesh numbered in x(r)- or y(z)- direction.

### Constructors
```julia
p51(data)
```
### Arguments
```julia
* `data::Dict{Symbol, Any}` : Dictionary containing all input data
```

### Required data dictionary keys
```julia
* struc_el::StructuralElement                          : Structural element
* support::Array{Tuple{Int,Array{Int,1}},1}        : Fixed-displacements vector
* loaded_nodes::Array{Tuple{Int,Array{Float64,1}},1} : Node load vector
* properties::Vector{Float64}                          : Material properties
* x_coords::FloatRange{Floalt64}                       : x-coordinate vector
* y_coords::FloatRange{Floalt64}                       : y-coordinate vector
* thickness:: Float64                                  : Thickness of plate
```

### Optional additional data dictionary keys
```julia
* penalty = 1e20               : Penalty used for fixed degrees of freedoms
* etype::Vector{Int}           : Element material vector if np_types > 1
```

### Return values
```julia
* fem                          : Fem object
```

### Related help
```julia
?StructuralElement             : List of available structural element types
?Plane                         : Help on a Plane structural element
?FiniteElement                 : List finite element types
?Quadrilateral                 : Help on Quadrilateral finite element
```
"""
function p51_skyline(data::Dict{Symbol, Any})

  # Setup basic dimensions of arrays

  # Parse & check FEdict data

  if :struc_el in keys(data)
    struc_el = data[:struc_el]
  else
    println("No fin_el type specified.")
    return
  end

  ndim = struc_el.ndim
  nst = struc_el.nst

  # Add radial stress
  struc_el.axisymmetric && (nst = 4)

  fin_el = struc_el.fin_el
  @assert typeof(fin_el) <: FiniteElement

  if typeof(fin_el) == Line
    (nels, nn) = mesh_size(fin_el, struc_el.nxe)
  elseif typeof(fin_el) == Triangle || typeof(fin_el) == Quadrilateral
    (nels, nn) = mesh_size(fin_el, struc_el.nxe, struc_el.nye)
  elseif typeof(fin_el) == Hexahedron
    (nels, nn) = mesh_size(fin_el, struc_el.nxe, struc_el.nye, struc_el.nze)
  else
    println("$(typeof(fin_el)) is not a known finite element.")
    return
  end

  nodof = fin_el.nodof           # Degrees of freedom per node
  ndof = fin_el.nod * nodof      # Degrees of freedom per fin_el

  # Update penalty if specified in FEdict

  penalty = 1e20
  if :penalty in keys(data)
    penalty = data[:penalty]
  end

  # Allocate all arrays

  # Start with arrays to be initialized from FEdict

  if :properties in keys(data)
    prop = zeros(size(data[:properties], 1), size(data[:properties], 2))
    for i in 1:size(data[:properties], 1)
      prop[i, :] = data[:properties][i, :]
    end
  else
    println("No :properties key found in FEdict")
  end

  nf = ones(Int, nodof, nn)
  if :support in keys(data)
    for i in 1:size(data[:support], 1)
      nf[:, data[:support][i][1]] = data[:support][i][2]
    end
  end

  x_coords = zeros(nn)
  if :x_coords in keys(data)
    x_coords = data[:x_coords]
  end

  y_coords = zeros(nn)
  if :y_coords in keys(data)
    y_coords = data[:y_coords]
  end

  z_coords = zeros(nn)
  if :z_coords in keys(data)
    z_coords = data[:z_coords]
  end

  etype = ones(Int, nels)
  if :etype in keys(data)
    etype = data[:etype]
  end

  # All other arrays

  points = zeros(struc_el.nip, ndim)
  g = zeros(Int, ndof)
  g_coord = zeros(ndim,nn)
  fun = zeros(fin_el.nod)
  coord = zeros(fin_el.nod, ndim)
  gamma = zeros(nels)
  jac = zeros(ndim, ndim)
  g_num = zeros(Int, fin_el.nod, nels)
  der = zeros(ndim, fin_el.nod)
  deriv = zeros(ndim, fin_el.nod)
  bee = zeros(nst,ndof)
  km = zeros(ndof, ndof)
  mm = zeros(ndof, ndof)
  gm = zeros(ndof, ndof)
  kg = zeros(ndof, ndof)
  eld = zeros(ndof)
  weights = zeros(struc_el.nip)
  g_g = zeros(Int, ndof, nels)
  num = zeros(Int, fin_el.nod)
  actions = zeros(ndof, nels)
  displacements = zeros(size(nf, 1), ndim)
  gc = ones(ndim)
  dee = zeros(nst,nst)
  sigma = zeros(nst)
  axial = zeros(nels)

  formnf!(nodof, nn, nf)
  neq = maximum(nf)
  kdiag = zeros(Int, neq)

  # Find global array sizes

  for iel in 1:nels
    geom_rect!(fin_el, iel, x_coords, y_coords, coord, num, struc_el.direction)
    num_to_g!(num, nf, g)
    g_num[:, iel] = num
    g_coord[:, num] = coord'
    g_g[:, iel] = g
    fkdiag!(kdiag, g)
  end

  for i in 2:neq
    kdiag[i] = kdiag[i] + kdiag[i-1]
  end

  kv = zeros(kdiag[neq])
  gv = zeros(kdiag[neq])

  println("There are $(neq) equations and the skyline storage is $(kdiag[neq]).")

  sample!(fin_el, points, weights)

  for iel in 1:nels
    deemat!(dee, prop[etype[iel], 1], prop[etype[iel], 2])
    num = g_num[:, iel]
    coord = g_coord[:, num]'              # Transpose
    g = g_g[:, iel]
    km = zeros(ndof, ndof)
    for i in 1:struc_el.nip
      shape_fun!(fun, points, i)
      shape_der!(der, points, i)
      jac = der*coord
      detm = det(jac)
      jac = inv(jac)
      deriv = jac*der
      beemat!(bee, deriv)
      if struc_el.axisymmetric
        gc = fun'*coord
        bee[4, 1:2:(ndof-1)] = fun[:]/gc[1]
      end
      km += (bee')*dee*bee*detm*weights[i]*gc[1]
    end
    fsparv(kv, km, g, kdiag)
  end

  loads = zeros(neq + 1)
  if :loaded_nodes in keys(data)
    for i in 1:size(data[:loaded_nodes], 1)
      loads[nf[:, data[:loaded_nodes][i][1]]+1] = data[:loaded_nodes][i][2]
    end
  end

  fixed_freedoms = 0
  if :fixed_freedoms in keys(data)
    fixed_freedoms = size(data[:fixed_freedoms], 1)
  end
  no = zeros(Int, fixed_freedoms)
  node = zeros(Int, fixed_freedoms)
  sense = zeros(Int, fixed_freedoms)
  value = zeros(Float64, fixed_freedoms)
  if :fixed_freedoms in keys(data) && fixed_freedoms > 0
    for i in 1:fixed_freedoms
      no[i] = nf[data[:fixed_freedoms][i][2], data[:fixed_freedoms][i][1]]
      value[i] = data[:fixed_freedoms][i][3]
    end
    kv[kdiag[no]] = kv[kdiag[no]] + penalty
    loads[no+1] = kv[kdiag[no]] .* value
  end

  sparin!(kv, kdiag)
  loads[2:end] = spabac!(kv, loads[2:end], kdiag)
  loads[1] = 0.0
  nf1 = deepcopy(nf) + 1

  if struc_el.axisymmetric
    println("\nNode     r-disp          z-disp")
  else
    println("\nNode     x-disp          y-disp")
  end

  tmp = []
  for i in 1:nn
    tmp = vcat(tmp, loads[nf1[:,i]])
    xstr = @sprintf("%+.4e", loads[nf1[1,i]])
    ystr = @sprintf("%+.4e", loads[nf1[2,i]])
    println("  $(i)    $(xstr)     $(ystr)")
  end

  struc_el.nip = 1
  points = zeros(struc_el.nip, ndim)
  weights = zeros(struc_el.nip)
  sample!(fin_el, points, weights)
  println("\nThe integration point (nip = $(struc_el.nip)) stresses are:")
  if struc_el.axisymmetric
    println("\nElement  r-coord   z-coord     sig_r        sig_z        tau_rz        sig_t")
  else
    println("\nElement  x-coord   y-coord     sig_x        sig_y        tau_xy")
  end
  for iel in 1:nels
    deemat!(dee, prop[etype[iel], 1], prop[etype[iel], 2])
    num = g_num[:, iel]
    coord = g_coord[:, num]'
    g = g_g[:, iel]
    eld = loads[g+1]
    for i in 1:struc_el.nip
      shape_fun!(fun, points, i)
      shape_der!(der, points, i)
      gc = fun'*coord
      jac = inv(der*coord)
      deriv = jac*der
      beemat!(bee, deriv)
      if struc_el.axisymmetric
        gc = fun'*coord
        bee[4, 1:2:(ndof-1)] = fun[:]/gc[1]
      end
      sigma = dee*(bee*eld)
      gc1 = @sprintf("%+.4f", gc[1])
      gc2 = @sprintf("%+.4f", gc[2])
      s1 = @sprintf("%+.4e", sigma[1])
      s2 = @sprintf("%+.4e", sigma[2])
      s3 = @sprintf("%+.4e", sigma[3])
      println("   $(iel)     $(gc1)   $(gc2)   $(s1)  $(s2)  $(s3)")
    end
  end
  println()

  FEM(struc_el, fin_el, ndim, nels, nst, ndof, nn, nodof, neq, penalty,
    etype, g, g_g, g_num, kdiag, nf, no, node, num, sense, actions,
    bee, coord, gamma, dee, der, deriv, displacements, eld, fun, gc,
    g_coord, jac, km, mm, gm, kv, gv, loads, points, prop, sigma, value,
    weights, x_coords, y_coords, z_coords, axial)

end

max_stars_repo_path: src/5 Elastic Solids/p51_skyline.jl

max_issues_repo_path: src/5 Elastic Solids/p51_skyline.jl

max_forks_repo_path: src/5 Elastic Solids/p51_skyline.jl


text:
using Tikz

colors = Dict( "blueblue" => [50.,60.,180.] )

style = tikzStyle( "style1", lw = 2.0 )

obj = tikzObj( "example1.tex", standalone = true, colors = colors, styles = [style] )

plot( obj, [0., 5., 2.], [1.,6.,-1], smooth=false, linestyle="--", color="blueblue", style=style )

node( obj, 1.0, 3.0, color = "blueblue", extracommands="circle,fill=gray,inner sep = 0pt,minimum size = 5pt" )

closeObj( obj )

max_stars_repo_name: mopg/julia-tikz

max_issues_repo_name: mopg/julia-tikz

max_forks_repo_name: mopg/julia-tikz


text:
@static if Sys.isapple()
    struct Cpasswd
        pw_name::Cstring
        pw_passwd::Cstring
        pw_uid::Cint
        pw_gid::Cint
        pw_change::Cint
        pw_class::Cstring
        pw_gecos::Cstring
        pw_dir::Cstring
        pw_shell::Cstring
        pw_expire::Cint
        pw_fields::Cint
    end
elseif Sys.islinux()
    struct Cpasswd
       pw_name::Cstring
       pw_passwd::Cstring
       pw_uid::Cint
       pw_gid::Cint
       pw_gecos::Cstring
       pw_dir::Cstring
       pw_shell::Cstring
    end
else
    struct Cpasswd
        pw_name::Cstring
        pw_uid::Cint
        pw_gid::Cint
        pw_dir::Cstring
        pw_shell::Cstring
    end

    Cpasswd() = Cpasswd(pointer("NA"), 0, 0, pointer("NA"), pointer("NA"))
end

struct Cgroup
    gr_name::Cstring
    gr_passwd::Cstring
    gr_gid::Cint
end

Cgroup() = Cgroup(pointer("NA"), pointer("NA"), 0)

struct User
    name::String
    uid::UInt64
    gid::UInt64
    dir::String
    shell::String
end

function User(ps::Cpasswd)
    User(
        unsafe_string(ps.pw_name),
        UInt64(ps.pw_uid),
        UInt64(ps.pw_gid),
        unsafe_string(ps.pw_dir),
        unsafe_string(ps.pw_shell)
    )
end

User(passwd::Ptr{Cpasswd}) = User(unsafe_load(passwd))
Base.show(io::IO, user::User) = print(io, "$(user.uid) ($(user.name))")

@static if Sys.isunix()
    function User(name::String)
        Libc.errno(0)
        ps = ccall(:getpwnam, Ptr{Cpasswd}, (Ptr{UInt8},), name)
        ret = Libc.errno()

        systemerror(:getpwnam, !iszero(ret))
        ps == C_NULL && throw(ArgumentError("User $name not found."))

        return User(ps)
    end

    function User(uid::UInt)
        Libc.errno(0)
        ps = ccall(:getpwuid, Ptr{Cpasswd}, (UInt64,), uid)
        ret = Libc.errno()

        systemerror(:getpwuid, !iszero(ret))
        ps == C_NULL && throw(ArgumentError("User $uid not found."))

        return User(ps)
    end

    User() = User(UInt(ccall(:geteuid, Cint, ())))
else
    User(name::String) = User(Cpasswd())
    User(uid::UInt) = User(Cpasswd())
    User() = User(UInt64(0))
end

struct Group
    name::String
    gid::UInt64
end

Group(gr::Cgroup) = Group(unsafe_string(gr.gr_name), UInt64(gr.gr_gid))
Group(group::Ptr{Cgroup}) = Group(unsafe_load(group))
Base.show(io::IO, group::Group) = print(io, "$(group.gid) ($(group.name))")

@static if Sys.isunix()
    function Group(name::String)
        Libc.errno(0)
        gr = ccall(:getgrnam, Ptr{Cgroup}, (Ptr{UInt8},), name)
        ret = Libc.errno()

        systemerror(:getgrnam, !iszero(ret))
        gr == C_NULL && throw(ArgumentError("Group $name not found."))
        return Group(gr)
    end

    function Group(gid::UInt)
        Libc.errno(0)
        gr = ccall(:getgrgid, Ptr{Cgroup}, (UInt64,), gid)
        ret = Libc.errno()

        systemerror(:getgrgid, !iszero(ret))
        gr == C_NULL && throw(ArgumentError("Group $gid not found."))
        return Group(gr)
    end

    Group() = Group(UInt64(ccall(:getegid, Cint, ())))
else
    Group(name::String) = Group(Cgroup())
    Group(uid::UInt) = Group(Cgroup())
    Group() = Group(UInt64(0))
end


text:
workspace()
using Gaston

include("BayesianBlocks.jl")

sig = [randn(100); randn(100)+10]
edges = bayesian_blocks(sig)
# chps = extract_changepoints(last)
plot(sig,sig+5)

max_stars_repo_name: av-maslov/Recurrency

max_issues_repo_name: av-maslov/Recurrency

max_forks_repo_name: av-maslov/Recurrency


text:
module SyncBarriersBenchmarks

using BenchmarkTools: Benchmark, BenchmarkGroup

include("bench_uniform_loops.jl")
include("bench_fuzzy.jl")
include("bench_cml.jl")
include("bench_gcm.jl")

function setup(; uniform_loops = (), fuzzy = (), cml = (), gcm = ())
    suite = BenchmarkGroup()
    suite["uniform_loops"] = BenchUniformLoops.setup(; uniform_loops...)
    suite["fuzzy"] = BenchFuzzy.setup(; fuzzy...)
    suite["cml"] = BenchCML.setup(; cml...)
    suite["gcm"] = BenchGCM.setup(; gcm...)
    return suite
end

function set_smoke_params!(bench)
    bench.params.seconds = 0.001
    bench.params.evals = 1
    bench.params.samples = 1
    bench.params.gctrial = false
    bench.params.gcsample = false
    return bench
end

foreach_benchmark(f!, bench::Benchmark) = f!(bench)
function foreach_benchmark(f!, group::BenchmarkGroup)
    for x in values(group)
        foreach_benchmark(f!, x)
    end
end

function setup_smoke(; ntasks = Threads.nthreads())
    suite = setup(
        fuzzy = (ntasks = ntasks, n = ntasks, m = ntasks, niters = 10),
        cml = (ntasks = ntasks, nsites = ntasks + 2, nsteps = 10),
        gcm = (ntasks = ntasks, nsites = ntasks + 2, nsteps = 10),
    )
    foreach_benchmark(set_smoke_params!, suite)
    return suite
end

function clear()
    BenchUniformLoops.clear()
    BenchFuzzy.clear()
    BenchCML.clear()
    BenchGCM.clear()
end

end  # module SyncBarriersBenchmarks


text:
using BinaryBuilder, SHA

version_dict = Dict()
for f in readdir("products/")
    if !endswith(f, ".tar.gz")
        continue
    end
    @info("Hashing $f....")
    name, version, platkey = extract_name_version_platform_key(f)

    if !(version in keys(version_dict))
        version_dict[version] = Dict()
    end
    hash = open(joinpath("products", f)) do f
        bytes2hex(sha256(f))
    end
    version_dict[version][triplet(platkey)] = (f => hash)
end

name = "GCC"
tag_name = BinaryBuilder.get_tag_name()
bin_path = "https://github.com/staticfloat/GCCBuilder/releases/download/$(tag_name)"
products = [ExecutableProduct(Prefix("."), "gcc", :gcc)]
for version in keys(version_dict)
    print_buildjl(pwd(), name, version, products, version_dict[version], bin_path)
end

@info("Uploading to GitHub....")
run(`ghr -replace $(tag_name) products/`)

max_stars_repo_name: staticfloat/GCCBuilder

max_issues_repo_name: staticfloat/GCCBuilder

max_forks_repo_name: staticfloat/GCCBuilder


text:
using Highbury
using Test

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


text:
module DebugPrintln

export dprintln, debug, nodebug

DebugFlag = true

function set_debug_flag(value)
	global DebugFlag
	DebugFlag = value
end
nodebug() = set_debug_flag(false)
debug() = set_debug_flag(true)

function dprintln(args...)
	if DebugPrintln.DebugFlag == true
		println(args...)
	end
end

end


text:
module SOtoZulip
using Base64
using HTTP
using JSON3
using Underscores
using CodecZlib
using Gumbo
using Cascadia
using Cascadia: matchFirst
using SQLite
using MD5
using Dates
using StructTypes

# Zulip related things
export ZulipClient, sendMessage, global_zulip!
export SOClient, searchtag, getquestions, getqanswers, getallqanswers
export show_query, getdb
export process_question, process_answers
export invalidate_answer, invalidate_question
export get_questions, clean, get_qids

include("zulipclient.jl")
include("soclient.jl")
include("utils.jl")
include("dbwrapper.jl")

end # module


text:
#
# Breadth-first constructive approach
#

function breadth_first_constructive(indsT::Vector)
  ntensors = length(indsT)
  if ntensors ≤ 16
    return breadth_first_constructive(UInt16, DimT, indsT)
  elseif ntensors ≤ 32
    return breadth_first_constructive(UInt32, DimT, indsT)
  elseif ntensors ≤ 64
    return breadth_first_constructive(UInt64, DimT, indsT)
  elseif ntensors ≤ 128
    return breadth_first_constructive(UInt128, DimT, indsT)
  else
    return breadth_first_constructive(BitSet, DimT, indsT)
  end
end

breadth_first_constructive(T::Tuple) = breadth_first_constructive(collect(T))

function breadth_first_constructive(::Type{TensorSetT}, ::Type{DimT},
                                    T::Vector{IndexSetT}) where {IndexSetT, TensorSetT, DimT}
  labels, alldims = label_dims(DimT, T)
  nlabels = length(alldims)
  if nlabels ≤ 16
    return breadth_first_constructive(TensorSetT, UInt16, labels, alldims)
  elseif nlabels ≤ 32
    return breadth_first_constructive(TensorSetT, UInt32, labels, alldims)
  elseif nlabels ≤ 64
    return breadth_first_constructive(TensorSetT, UInt64, labels, alldims)
  elseif nlabels ≤ 128
    return breadth_first_constructive(TensorSetT, UInt128, labels, alldims)
  else
    return breadth_first_constructive(TensorSetT, BitSet, labels, alldims)
  end
end

function breadth_first_constructive(::Type{TensorSetT}, ::Type{LabelSetT}, labels::Vector, alldims::Vector) where {TensorSetT, LabelSetT}
  return breadth_first_constructive(TensorSetT, map(label -> bitset(LabelSetT, label), labels), alldims)
end

# TODO: delete?
#function breadth_first_constructive(::Type{TensorSetT}, ::Type{LabelSetT}, ::Type{DimT},
#                                    T::Vector{<: ITensor}) where {TensorSetT, LabelSetT, DimT}
#  indsT = [inds(Tₙ) for Tₙ in T]
#  return breadth_first_constructive(TensorSetT, LabelSetT, DimT, indsT)
#end

function breadth_first_constructive(::Type{TensorSetT}, ::Type{LabelSetT}, ::Type{DimT},
                                    T::Vector{IndexSetT}) where {IndexSetT, TensorSetT, LabelSetT, DimT}
  labels, alldims = label_dims(DimT, T)
  return breadth_first_constructive(TensorSetT, map(label -> bitset(LabelSetT, label), labels), alldims)
end

# A type storing information about subnetworks
const SubNetwork{LabelSetT, DimT} = NamedTuple{(:inds, :cost, :sequence), Tuple{LabelSetT, DimT, Vector{Any}}}

function breadth_first_constructive(::Type{TensorSetT},
                                    T::Vector{LabelSetT},
                                    alldims::Vector{DimT}) where {TensorSetT, LabelSetT, DimT}
  components = connectedcomponents(T, alldims)
  N = length(components)
  if N == 1
    return _breadth_first_constructive(TensorSetT, collect(1:length(T)), T, alldims)
  end
  sequences = Vector{Any}(undef, N)
  for n in 1:N
    componentsₙ = components[n]
    if length(componentsₙ) == 1
      sequences[n] = only(componentsₙ)
      continue
    elseif length(componentsₙ) == 2
      sequences[n] = componentsₙ
      continue
    end
    sequences[n] = _breadth_first_constructive(TensorSetT, componentsₙ, T[componentsₙ], alldims)
  end
  return sequences
end

# Apply breadth_first_constructive to a single disconnected subnetwork
# Based on: https://arxiv.org/abs/1304.6112 and https://github.com/Jutho/TensorOperations.jl/blob/v3.1.0/src/indexnotation/optimaltree.jl
function _breadth_first_constructive(::Type{TensorSetT},
                                     Tlabels::Vector,
                                     T::Vector{LabelSetT},
                                     alldims::Vector{DimT}) where {TensorSetT, LabelSetT, DimT}
  n = length(T)

  # `cache[c]` is the set of all objects made up by
  # contracting `c` unique tensors from the original tensors `1:n`.
  cache = Vector{Dict{TensorSetT, SubNetwork{LabelSetT, DimT}}}(undef, n)
  for c in 1:n
    # Initialized to empty
    cache[c] = eltype(cache)()
  end
  # Fill the first cache with trivial data
  for i in 1:n
    cache[1][bitset(TensorSetT, [Tlabels[i]])] = (inds = T[i], cost = 0, sequence = Any[])
  end

  # TODO: pick a reasonable maxcost, the product of all dimensions
  # of tensors in the network. Could also be `typemax(DimT)`?
  maxcost = try
    # This may overflow, so we catch the error and return typemax(DimT)
    Base.Checked.checked_mul(reduce(Base.Checked.checked_mul, alldims; init = one(DimT)), maximum(alldims))
  catch
    typemax(DimT)
  end

  # TODO: pick a reasonable initialcost
  # Maybe use the cost of the trivial contraction [4, [3, [2, 1]]]?
  tensordims = Vector{DimT}(undef, n)
  for k in 1:n
    tensordims[k] = dim(T[k], alldims)
  end
  _initialcost, overflow = Base.Checked.mul_with_overflow(maximum(tensordims), minimum(tensordims))
  _initialcost = overflow ? typemax(DimT) : _initialcost
  initialcost = min(maxcost, _initialcost)

  # Factor to increase the cost cap by each iteration
  costfac = maximum(alldims)

  currentcost = initialcost
  previouscost = zero(initialcost)

  while isempty(cache[n])
    nextcost = maxcost

    # c is the total number of tensors being contracted
    # in the current sequence
    for c in 2:n
      # For each pair of sets Sᵈ, Sᶜ⁻ᵈ, 1 ≤ d ≤ ⌊c/2⌋
      for d in 1:c÷2
        for a in keys(cache[d]), b in keys(cache[c-d])
          if d == c-d && _isless(b, a)
            # When d == c-d (the subset sizes are equal), check that
            # b > a so that that case (a,b) and (b,a) are not repeated
            continue
          end

          if !_isemptyset(_intersect(a, b))
            # Check that each element of S¹ appears
            # at most once in (TᵃTᵇ).
            continue
          end

          # Use previously computed cost of contracting network `ab` and compare against the previouscost
          ab = _union(a, b)
          cache_c = @inbounds cache[c]
          cache_ab = get(cache_c, ab, nothing)
          currentcost_ab = isnothing(cache_ab) ? currentcost : cache_ab.cost
          if currentcost_ab ≤ previouscost
            continue
          end

          # Determine the cost μ of contracting Tᵃ, Tᵇ
          # These dictionary calls and `contraction_cost` take
          # up most of the time.
          cache_a = cache[d][a]
          cache_b = cache[c-d][b]

          if dim(_intersect(cache_a.inds, cache_b.inds), alldims) < 2
            # XXX: For now, ignore outer products contractions.
            # In the future, handle this in a more sophisticated way.
            continue
          end

          cost, inds_ab = contraction_cost(cache_a.inds, cache_b.inds, alldims)
          if iszero(cost)
            # If the cost is zero, that means the multiplication overflowed
            continue
          end

          if d > 1
            # Add to cost of contracting the subnetwork `a`
            cost, overflow = Base.Checked.add_with_overflow(cost, cache_a.cost)
            overflow && continue
          end
          if c-d > 1
            # Add to cost of contracting the subnetwork `b`
            cost, overflow = Base.Checked.add_with_overflow(cost, cache_b.cost)
            overflow && continue
          end

          if cost ≤ currentcost_ab
            cost_ab = cost
            if d == 1
              sequence_a = _only(a)
            else
              sequence_a = cache_a.sequence
            end
            if c-d == 1
              sequence_b = _only(b)
            else
              sequence_b = cache_b.sequence
            end
            sequence_ab = Any[sequence_a, sequence_b]

            # XXX: this call is pretty slow (maybe takes 1/3 of total time in large n limit)
            cache_c[ab] = (inds = inds_ab, cost = cost_ab, sequence = sequence_ab)
          end # if cost ≤ currentcost_ab
        end # for a in S[d], b in S[c-d]
      end # for d in 1:c÷2
    end # for c in 2:n
    previouscost = currentcost
    currentcost = min(maxcost, nextcost*costfac)

    # Reset all tensors to old
    for i in 1:n
      for a in eachindex(cache[i])
        cache_a = cache[i][a]
        cache[i][a] = (inds = cache_a.inds, cost = cache_a.cost, sequence = cache_a.sequence)
      end
    end
  end # while isempty(S[n])
  Sⁿ = bitset(TensorSetT, Tlabels)
  return cache[n][Sⁿ].sequence
end


text:
export proxPositive!, normPositive

"""
    proxPositive!(x::Vector{T},λ::Float64=1.0;kargs...) where T

    enforce positivity and realness of solution `x`.
"""
function proxPositive!(x::Vector{T},λ::Float64=1.0;kargs...) where T
  enfReal!(x)
  enfPos!(x)
end

"""
    returns the value of the characteristic function of real, positive numbers.
    normPositive(x) = (isreal(x)&&x>0) ? 0 : Inf
"""
function normPositive(x::Vector{T},λ::Float64=1.0;kargs...) where T
  y = copy(x)
  proxPositive!(y)
  if y != x
    return Inf
  end
  return 0.
end


text:
module TestExamples3dP4est

using Test
using Trixi

include("test_trixi.jl")

# pathof(Trixi) returns /path/to/Trixi/src/Trixi.jl, dirname gives the parent directory
EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "p4est_3d_dgsem")

@testset "P4estMesh" begin
  @trixi_testset "elixir_advection_basic.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_basic.jl"),
      l2   = [0.00013446460962856976],
      linf = [0.0012577781391462928])
  end

  @trixi_testset "elixir_advection_unstructured_curved.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_unstructured_curved.jl"),
      l2   = [0.0006244885699399409],
      linf = [0.04076651402041587])
  end

  @trixi_testset "elixir_advection_nonconforming.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_nonconforming.jl"),
      l2   = [0.0024774648310858928],
      linf = [0.021727876954353964])
  end

  @trixi_testset "elixir_advection_amr.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_amr.jl"),
      l2   = [9.773852895157622e-6],
      linf = [0.0005853874124926162])
  end

  @trixi_testset "elixir_advection_amr_unstructured_curved.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_amr_unstructured_curved.jl"),
      l2   = [0.00014665036779554962],
      linf = [0.00845936405684372])
  end

  @trixi_testset "elixir_advection_cubed_sphere.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_cubed_sphere.jl"),
      l2   = [0.0077828979483271195],
      linf = [0.08759188779479488])
  end

  @trixi_testset "elixir_advection_restart.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_restart.jl"),
      l2   = [0.0281388160824776],
      linf = [0.08740635193023694])
  end

  @trixi_testset "elixir_euler_source_terms_nonperiodic.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms_nonperiodic.jl"),
      l2   = [0.0018268326813744103, 0.001745601029521995, 0.001745601029521962, 0.0017456010295218891, 0.003239834454817457],
      linf = [0.014660503198892005, 0.01506958815284798, 0.01506958815283821, 0.015069588152864632, 0.02700205515651044])
  end

  @trixi_testset "elixir_euler_free_stream.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_free_stream.jl"),
      l2   = [7.839921025049303e-15, 4.3320953097202227e-14, 4.61318381170666e-14, 5.466388734250446e-14, 9.381101254280883e-14],
      linf = [1.0269562977782698e-12, 4.9040771443742415e-12, 6.0748628349927e-12, 7.261524714863299e-12, 1.1615597372838238e-11],
      tspan = (0.0, 0.1))
  end

  @trixi_testset "elixir_euler_free_stream_extruded.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_free_stream_extruded.jl"),
      l2   = [1.6667578354338437e-14, 8.351207753024527e-15, 1.805859250264852e-14, 1.90966219766217e-14, 1.90966219766217e-14],
      linf = [1.0413891970983968e-13, 4.988787161153141e-13, 2.0694557179012918e-13, 6.221689829999377e-13, 1.0587086762825493e-12])
  end

  @trixi_testset "elixir_euler_ec.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_ec.jl"),
      l2   = [0.009411352784942132, 0.007257089010979736, 0.007478550902224631, 0.007472675228591207, 0.023851056842379016],
      linf = [0.19307851386933228, 0.27472353046951714, 0.3039307725248127, 0.3063060931738112, 0.5144856763103725])
  end
end

end # module


text:
using Word2Vec
using LinearAlgebra

function loadNonStopWords(fileName::String)::Tuple{Vector{String}, Dict{String, Int}}
    file = open(fileName)
    nonStopWords = Vector{String}()
    nonStopWordIDs = Dict{String, Int}()
    for phrase in readlines(file)
        if phrase != ""
            push!(nonStopWords, phrase)
            nonStopWordIDs[phrase] = length(nonStopWords)
        end
    end
    close(file)
    (nonStopWords, nonStopWordIDs)
end

function loadPhraseIDGroups(fileName::String, nonStopWordIDs::Dict{String,Int})::Vector{Vector{Int}}
    file = open(fileName)
    phraseIDGroups = Vector{Vector{Int}}()
    for line in readlines(file)
        if length(line) == 0
            continue
        end
        phrases = split(line)
        idGroup = Vector{Int}()
        for phrase in phrases
            if haskey(nonStopWordIDs, phrase)
                push!(idGroup, nonStopWordIDs[phrase])
            end
        end
        push!(phraseIDGroups, idGroup)
    end
    close(file)
    phraseIDGroups
end

function normalizeVectors(wordVectors::WordVectors, nonStopWords::Vector{String})::Array{Float64,2}
    numDims, numWords = size(wordVectors)
    numVectors = length(nonStopWords)
    result = Array{Float64,2}(undef, (numDims, numVectors))
    for i in 1:numVectors
        vector = get_vector(wordVectors, nonStopWords[i])
        w = 1.0 / norm(vector)
        result[:,i] = vector * w
    end
    result
end

function getSimilarityMatrix(vectors::Array{Float64,2})::Array{Float64,2}
    0.5 .+ 0.5 .* (vectors' * vectors)
end

function getPhrasePairFreqs(nonStopWords::Vector{String}, nonStopWordIDs::Dict{String,Int}, phraseIDGroups::Vector{Vector{Int}})::Dict{Tuple{Int,Int},Int}
    numPhrases = length(nonStopWords)
    idMap = Dict{Int,Int}()
    for (newID, phrase) in enumerate(nonStopWords)
        idMap[nonStopWordIDs[phrase]] = newID
    end

    pairFreqs = Dict{Tuple{Int,Int},Int}()
    for idGroup in phraseIDGroups
        numIDsInGroup = length(idGroup)
        convertedIDs = [idMap[id] for id in idGroup]
        for (j1, id1) in enumerate(convertedIDs)
            for j2 in j1 + 1:numIDsInGroup
                id2 = convertedIDs[j2]
                pair = if id1 < id2
                    (id1, id2)
                else
                    (id2, id1)
                end
                pairFreqs[pair] = get(pairFreqs, pair, 0) + 1
            end
        end
    end

    pairFreqs
end

function loadAspects(fileName::String, nonStopWordIDs::Dict{String,Int})::Tuple{Vector{Set{Int}}, Vector{Int}}
    file = open(fileName)
    aspects = Vector{Set{Int}}()
    numPhrases = length(nonStopWordIDs)
    aspectIDs = Vector{Int}(undef, numPhrases)
    for line in readlines(file)
        fields = split(line)
        if length(fields) != 3
            continue
        end
        aspectID = parse(Int, fields[1])
        card = parse(Int, fields[2])
        phrase = fields[3]
        phraseID = nonStopWordIDs[phrase]
        if aspectID > length(aspects)
            @assert aspectID == length(aspects) + 1
            push!(aspects, Set{Int}())
        end
        push!(aspects[aspectID], phraseID)
        aspectIDs[phraseID] = aspectID
    end
    close(file)
    (aspects, aspectIDs)
end

function findDualAspectPairFreqs(pairFreqs::Dict{Tuple{Int,Int},Int}, aspectIDs::Vector{Int})::Vector{Tuple{Tuple{Int,Int},Int}}
    dualAspectPairFreqs = Vector{Tuple{Tuple{Int,Int},Int}}()
    for (pair, freq) in pairFreqs
        aspectID1 = aspectIDs[pair[1]]
        aspectID2 = aspectIDs[pair[2]]
        if aspectID1 != aspectID2
            push!(dualAspectPairFreqs, (pair, freq))
        end
    end
    dualAspectPairFreqs
end

function getDualAspectConn(simMat::Array{Float64}, dualAspectPairFreqs::Vector{Tuple{Tuple{Int,Int},Int}}, pairID1::Int, pairID2::Int)::Float64
    if pairID1 != pairID2
        pairFreq1 = dualAspectPairFreqs[pairID1]
        pair1 = pairFreq1[1]
        freq1 = pairFreq1[2]
        pairFreq2 = dualAspectPairFreqs[pairID2]
        pair2 = pairFreq2[1]
        freq2 = pairFreq2[2]
        sim1 = simMat[pair1[1], pair2[1]] * simMat[pair1[2], pair2[2]]
        sim2 = simMat[pair1[1], pair2[2]] * simMat[pair1[2], pair2[1]]
        return min(sim1, sim2) * freq1 * freq2
    else
        return 0.0
    end
end

function sparsifyDualAspectConn(simMat::Array{Float64}, dualAspectPairFreqs::Vector{Tuple{Tuple{Int,Int},Int}}, numSingleSideNeighbors::Int)::Vector{Dict{Int,Float64}}
    numPoints = length(dualAspectPairFreqs)
    singleSideConns = Vector{Dict{Int,Float64}}(undef, numPoints)
    numCPUs = length(Sys.cpu_info())
    jobs = Channel{Tuple{Int,Int}}(numCPUs)
    Segments = Vector{Tuple{Int,Int}}
    jobOutputs = Channel{Tuple{Segments,Vector{Vector{Dict{Int,Float64}}}}}(numCPUs)

    function makeJobs(batchSize::Int)
        for i in 1:batchSize:numPoints
            put!(jobs, (i, min(i + batchSize - 1, numPoints)))
        end
    end

    function runJob(iCPU::Int)
        mySegments = Segments()
        myOutputs = Vector{Vector{Dict{Int,Float64}}}()
        conn1s = Vector{Float64}(undef, numPoints)
        for job in jobs
            push!(mySegments, job)
            segOutput = Vector{Dict{Int,Float64}}(undef, job[2] - job[1] + 1)
            # println("cpu $iCPU: computing sumConn for $job")
            for pairID1 in job[1]:job[2]
                sumConn1 = 0.0
                for pairID2 in 1:numPoints
                    conn1s[pairID2] = getDualAspectConn(simMat, dualAspectPairFreqs, pairID1, pairID2)
                end
                perm = sortperm(conn1s, rev=true)
                neighbors = Dict{Int, Float64}()
                for i in 1:numSingleSideNeighbors
                    neighbors[perm[i]] = conn1s[perm[i]]
                end
                segOutput[pairID1 - job[1] + 1] = neighbors
                if pairID1 % 1000 == 0
                    println("cpu $iCPU: $pairID1 of $numPoints have been processed")
                end
            end
            push!(myOutputs, segOutput)
            # println("cpu $iCPU: finish computing sumConn for $job")
        end # job
        put!(jobOutputs, (mySegments, myOutputs))
    end # runJob

    bind(jobs, @async makeJobs(100))
    for iCPU in 1:numCPUs
        Threads.@spawn runJob(iCPU)
    end

    for iCPU in 1:numCPUs
        segmentedOutputs = take!(jobOutputs)
        for (segID, segment) in enumerate(segmentedOutputs[1])
            singleSideConns[segment[1]:segment[2]] = segmentedOutputs[2][segID]
        end
    end

    doubleSideConns = [Dict{Int,Float64}() for i in 1:numPoints]
    for i in 1: numPoints
        for (j, conn) in singleSideConns[i]
            doubleSideConns[i][j] = conn
            doubleSideConns[j][i] = conn
        end
    end
    doubleSideConns
end

function saveDualAspectPairFreqs(fileName::String, dualAspectPairFreqs::Vector{Tuple{Tuple{Int,Int},Int}})
    file = open(fileName, "w")
    for pairFreq in dualAspectPairFreqs
        pair = pairFreq[1]
        freq = pairFreq[2]
        i = pair[1]
        j = pair[2]
        println(file, "$i $j $freq")
    end
    close(file)
end

function saveConns(fileName::String, sparseConns::Vector{Dict{Int,Float64}})
    file = open(fileName, "w")
    for (i, connsI) in enumerate(sparseConns)
        for (j, conn) in connsI
            println(file, "$i $j $conn")
        end
    end
    close(file)
end

numSingleSideNeighbors = 100
wordVectors = wordvectors("DualAspectTopicRepresentation/workspace/NIPS.awe.model.txt")
nonStopWords, nonStopWordIDs = loadNonStopWords("DualAspectTopicRepresentation/workspace/NIPS.nonStopWords.txt")
phraseIDGroups = loadPhraseIDGroups("DualAspectTopicRepresentation/workspace/NIPS.awe.training.txt", nonStopWordIDs)
pairFreqs = getPhrasePairFreqs(nonStopWords, nonStopWordIDs, phraseIDGroups)
vectors = normalizeVectors(wordVectors, nonStopWords)
simMat = getSimilarityMatrix(vectors)
aspects, aspectIDs = loadAspects("DualAspectTopicRepresentation/workspace/NIPS.aspects.1.txt", nonStopWordIDs)
dualAspectPairFreqs = findDualAspectPairFreqs(pairFreqs, aspectIDs)
saveDualAspectPairFreqs("DualAspectTopicRepresentation/workspace/NIPS.DAPairFreq.txt", dualAspectPairFreqs)
println("sparsifying connections ...")
sparseConns = sparsifyDualAspectConn(simMat, dualAspectPairFreqs, numSingleSideNeighbors)
saveConns("DualAspectTopicRepresentation/workspace/NIPS.DAPairConns.txt", sparseConns)

max_stars_repo_path: src/04 sparsify_connections.jl

max_issues_repo_path: src/04 sparsify_connections.jl

max_forks_repo_path: src/04 sparsify_connections.jl