text:
export pseudo_basis, basis_pmat

################################################################################
#
#  Basic field access
#
################################################################################

doc"""
***
      nf(O::NfRelOrd) -> RelativeExtension

> Returns the ambient number field of $\mathcal O$.
"""
nf(O::NfRelOrd) = O.nf

doc"""
    parent(O::NfRelOrd) -> NfRelOrdSet

Returns the parent of $\mathcal O$, that is, the set of orders of the ambient
number field.
"""
parent(O::NfRelOrd) = O.parent

doc"""
    isequation_order(O::NfRelOrd) -> Bool

> Returns whether $\mathcal O$ is the equation order of the ambient number
> field.
"""
isequation_order(O::NfRelOrd) = O.isequation_order

ismaximal_known(O::NfRelOrd) = O.ismaximal != 0

ismaximal(O::NfRelOrd) = O.ismaximal == 1

issimple(O::NfRelOrd) = issimple(nf(O))

elem_type(::Type{NfRelOrd{T, S}}) where {T, S} = NfRelOrdElem{T}

################################################################################
#
#  "Assure" functions for fields
#
################################################################################

function assure_has_basis_pmat(O::NfRelOrd{T, S}) where {T, S}
  if isdefined(O, :basis_pmat)
    return nothing
  end
  if !isdefined(O, :pseudo_basis)
    error("No pseudo_basis and no basis_pmat defined.")
  end
  pb = pseudo_basis(O, Val{false})
  L = nf(O)
  M = zero_matrix(base_ring(L), degree(O), degree(O))
  C = Vector{S}()
  for i = 1:degree(O)
    elem_to_mat_row!(M, i, pb[i][1])
    push!(C, deepcopy(pb[i][2]))
  end
  O.basis_pmat = PseudoMatrix(M, C)
  return nothing
end

function assure_has_pseudo_basis(O::NfRelOrd{T, S}) where {T, S}
  if isdefined(O, :pseudo_basis)
    return nothing
  end
  if !isdefined(O, :basis_pmat)
    error("No pseudo_basis and no basis_pmat defined.")
  end
  P = basis_pmat(O, Val{false})
  L = nf(O)
  pseudo_basis = Vector{Tuple{elem_type(L), S}}()
  for i = 1:degree(O)
    a = elem_from_mat_row(L, P.matrix, i)
    push!(pseudo_basis, (a, deepcopy(P.coeffs[i])))
  end
  O.pseudo_basis = pseudo_basis
  return nothing
end

function assure_has_basis_nf(O::NfRelOrd{T, S}) where {T, S}
  if isdefined(O, :basis_nf)
    return nothing
  end
  L = nf(O)
  pb = pseudo_basis(O)
  basis_nf = Vector{elem_type(L)}()
  for i = 1:degree(O)
    push!(basis_nf, pb[i][1])
  end
  O.basis_nf = basis_nf
  return nothing
end

function assure_has_basis_mat(O::NfRelOrd)
  if isdefined(O, :basis_mat)
    return nothing
  end
  O.basis_mat = basis_pmat(O).matrix
  return nothing
end

function assure_has_basis_mat_inv(O::NfRelOrd)
  if isdefined(O, :basis_mat_inv)
    return nothing
  end
  O.basis_mat_inv = inv(basis_mat(O, Val{false}))
  return nothing
end

function assure_has_inv_coeff_ideals(O::NfRelOrd)
  if isdefined(O, :inv_coeff_ideals)
    return nothing
  end
  pb = pseudo_basis(O, Val{false})
  O.inv_coeff_ideals = [ inv(pb[i][2]) for i in 1:degree(O) ]
  return nothing
end

################################################################################
#
#  Pseudo basis / basis pseudo-matrix
#
################################################################################

doc"""
***
      pseudo_basis(O::NfRelOrd{T, S}) -> Vector{Tuple{RelativeElement{T}{T}, S}}

> Returns the pseudo-basis of $\mathcal O$.
"""
function pseudo_basis(O::NfRelOrd, copy::Type{Val{T}} = Val{true}) where T
  assure_has_pseudo_basis(O)
  if copy == Val{true}
    return deepcopy(O.pseudo_basis)
  else
    return O.pseudo_basis
  end
end

doc"""
***
      basis_pmat(O::NfRelOrd) -> PMat

> Returns the basis pseudo-matrix of $\mathcal O$ with respect to the power basis
> of the ambient number field.
"""
function basis_pmat(O::NfRelOrd, copy::Type{Val{T}} = Val{true}) where T
  assure_has_basis_pmat(O)
  if copy == Val{true}
    return deepcopy(O.basis_pmat)
  else
    return O.basis_pmat
  end
end

doc"""
***
      inv_coeff_ideals(O::NfRelOrd{T, S}) -> Vector{S}

> Returns the inverses of the coefficient ideals of the pseudo basis of $O$.
"""
function inv_coeff_ideals(O::NfRelOrd, copy::Type{Val{T}} = Val{true}) where T
  assure_has_inv_coeff_ideals(O)
  if copy == Val{true}
    return deepcopy(O.inv_coeff_ideals)
  else
    return O.inv_coeff_ideals
  end
end

################################################################################
#
#  Basis / (inverse) basis matrix
#
################################################################################

doc"""
***
      basis_nf(O::NfRelOrd) -> Array{RelativeElement, 1}

> Returns the elements of the pseudo-basis of $\mathcal O$ as elements of the
> ambient number field.
"""
function basis_nf(O::NfRelOrd, copy::Type{Val{T}} = Val{true}) where T
  assure_has_basis_nf(O)
  if copy == Val{true}
    return deepcopy(O.basis_nf)
  else
    return O.basis_nf
  end
end

doc"""
***
      basis_mat(O::NfRelOrd{T, S}) -> Generic.Mat{T}

> Returns the basis matrix of $\mathcal O$ with respect to the power basis
> of the ambient number field.
"""
function basis_mat(O::NfRelOrd, copy::Type{Val{T}} = Val{true}) where T
  assure_has_basis_mat(O)
  if copy == Val{true}
    return deepcopy(O.basis_mat)
  else
    return O.basis_mat
  end
end

doc"""
***
      basis_mat_inv(O::NfRelOrd{T, S}) -> Generic.Mat{T}

> Returns the inverse of the basis matrix of $\mathcal O$.
"""
function basis_mat_inv(O::NfRelOrd, copy::Type{Val{T}} = Val{true}) where T
  assure_has_basis_mat_inv(O)
  if copy == Val{true}
    return deepcopy(O.basis_mat_inv)
  else
    return O.basis_mat_inv
  end
end

################################################################################
#
#  String I/O
#
################################################################################

function show(io::IO, S::NfRelOrdSet)
  print(io, "Set of orders of the number field ")
  print(io, S.nf)
end

function show(io::IO, O::NfRelOrd)
  compact = get(io, :compact, false)
  if compact
    if ismaximal_known(O) && ismaximal(O)
      print(io, "Relative maximal order with pseudo-basis ")
    else
      print(io, "Relative order with pseudo-basis ")
    end
    pb = pseudo_basis(O, Val{false})
    for i = 1:degree(O)
      print(io, "($(pb[i][1])) * ")
      showcompact(io, pb[i][2])
      if i != degree(O)
        print(io, ", ")
      end
    end
  else
    if ismaximal_known(O) && ismaximal(O)
      print(io, "Relative maximal order of ")
    else
      print(io, "Relative order of ")
    end
    println(io, nf(O))
    print(io, "with pseudo-basis ")
    pb = pseudo_basis(O, Val{false})
    for i = 1:degree(O)
      print(io, "\n(")
      print(io, pb[i][1])
      print(io, ", ")
      showcompact(io, pb[i][2])
      print(io, ")")
    end
  end
end

################################################################################
#
#  Discriminant
#
################################################################################

function assure_has_discriminant(O::NfRelOrd{T, S}) where {T, S}
  if isdefined(O, :disc_abs) || isdefined(O, :disc_rel)
    return nothing
  end
  d = det(trace_matrix(O))
  pb = pseudo_basis(O, Val{false})
  a = pb[1][2]^2
  for i = 2:degree(O)
    a *= pb[i][2]^2
  end
  disc = d*a
  simplify(disc)
  if T == nf_elem
    O.disc_abs = numerator(disc)
  else
    O.disc_rel = numerator(disc)
  end
  return nothing
end

function discriminant(O::NfRelOrd{nf_elem, S}) where S
  assure_has_discriminant(O)
  return deepcopy(O.disc_abs)
end

function discriminant(O::NfRelOrd{T, S}) where {T <: RelativeElement{U} where U, S}
  assure_has_discriminant(O)
  return deepcopy(O.disc_rel)
end

################################################################################
#
#  Degree
#
################################################################################

doc"""
***
      degree(O::NfRelOrd) -> Int

> Returns the degree of $\mathcal O$.
"""
degree(O::NfRelOrd) = degree(nf(O))

################################################################################
#
#  Deepcopy
#
################################################################################

doc"""
***
      deepcopy(O::NfRelOrd) -> NfRelOrd

> Makes a copy of $\mathcal O$.
"""
function Base.deepcopy_internal(O::NfRelOrd{T, S}, dict::ObjectIdDict) where {T, S}
  z = NfRelOrd{T, S}(O.nf)
  for x in fieldnames(O)
    if x != :nf && x != :parent && isdefined(O, x)
      setfield!(z, x, Base.deepcopy_internal(getfield(O, x), dict))
    end
  end
  z.nf = O.nf
  z.parent = O.parent
  return z
end

################################################################################
#
#  Inclusion of number field elements
#
################################################################################

function _check_elem_in_order(a::RelativeElement{T}, O::NfRelOrd{T, S}, short::Type{Val{V}} = Val{false}) where {T, S, V}
  b_pmat = basis_pmat(O, Val{false})
  t = zero_matrix(base_ring(nf(O)), 1, degree(O))
  elem_to_mat_row!(t, 1, a)
  t = t*basis_mat_inv(O, Val{false})
  if short == Val{true}
    for i = 1:degree(O)
      if !(t[1, i] in b_pmat.coeffs[i])
        return false
      end
    end
    return true
  else
    for i = 1:degree(O)
      if !(t[1, i] in b_pmat.coeffs[i])
        return false, Vector{T}()
      end
    end
    v = Vector{T}(degree(O))
    for i in 1:degree(O)
      v[i] = deepcopy(t[1, i])
    end
    return true, v
  end
end

doc"""
***
      in(a::RelativeElement, O::NfRelOrd) -> Bool

> Checks whether $a$ lies in $\mathcal O$.
"""
function in(a::RelativeElement{T}, O::NfRelOrd{T, S}) where {T, S}
  return _check_elem_in_order(a, O, Val{true})
end

################################################################################
#
#  Construction
#
################################################################################

doc"""
***
      Order(K::RelativeExtension{T}, M::Generic.Mat{T}) -> NfRelOrd

> Returns the order which has basis matrix $M$ with respect to the power basis
> of $K$.
"""
function Order(L::RelativeExtension{nf_elem}, M::Generic.Mat{nf_elem})
  # checks
  return NfRelOrd{nf_elem, NfOrdFracIdl}(L, deepcopy(M))
end

function Order(L::RelativeExtension{S}, M::Generic.Mat{S}) where S <: RelativeElement{T} where T
  # checks
  return NfRelOrd{elem_type(base_ring(L)), NfRelOrdFracIdl{T}}(L, deepcopy(M))
end

doc"""
***
      Order(K::RelativeExtension, M::PMat) -> NfRelOrd

> Returns the order which has basis pseudo-matrix $M$ with respect to the power basis
> of $K$.
"""
function Order(L::RelativeExtension{T}, M::PMat{T, S}) where {T, S}
  # checks
  return NfRelOrd{T, S}(L, deepcopy(M))
end

doc"""
***
      EquationOrder(L::RelativeExtension) -> NfRelOrd

> Returns the equation order of the number field $L$.
"""
function EquationOrder(L::RelativeExtension)
  M = identity_matrix(base_ring(L), degree(L))
  PM = PseudoMatrix(M)
  O = Order(L, PM)
  O.basis_mat_inv = M
  O.isequation_order = true
  return O
end

doc"""
***
      maximal_order(L::RelativeExtension) -> NfRelOrd

> Returns the maximal order of $L$.
"""
function maximal_order(L::RelativeExtension)
  try
    O = _get_maximal_order_of_nf_rel(L)
    return O
  catch e
    if !isa(e, AccessorNotSetError)
      rethrow(e)
    end
    O = MaximalOrder(L)
    _set_maximal_order_of_nf_rel(L, O)
    return O
  end
end

MaximalOrder(L::RelativeExtension) = MaximalOrder(EquationOrder(L))

function maximal_order_via_absolute(L::NfRel)
  Labs, lToLabs, kToLabs = absolute_field(L)
  Oabs = maximal_order(Labs)
  return relative_order(Oabs, lToLabs)
end

function maximal_order_via_absolute(m::NfRelToNf)
  Oabs = maximal_order(codomain(m))
  return relative_order(Oabs, m)
end

function maximal_order_via_simple(L::NfRel_ns)
  Ls, m = simple_extension(L)
  Os = maximal_order(Ls)
  return non_simple_order(Os, m)
end

function maximal_order_via_simple(m::NfRelToNfRel_nsMor)
  Os = maximal_order(domain(m))
  return non_simple_order(Os, m)
end

################################################################################
#
#  Equality
#
################################################################################

function ==(R::NfRelOrd, S::NfRelOrd)
  nf(R) != nf(S) && return false
  return basis_pmat(R, Val{false}) == basis_pmat(S, Val{false})
end

################################################################################
#
#  Trace matrix
#
################################################################################

doc"""
***
      trace_matrix(O::NfRelOrd{T, S}) -> Generic.Mat{T}

> Returns the trace matrix of $\mathcal O$.
"""
function trace_matrix(O::NfRelOrd)
  if isdefined(O, :trace_mat)
    return deepcopy(O.trace_mat)
  end
  L = nf(O)
  K = base_ring(L)
  b = basis_nf(O, Val{false})
  d = degree(L)
  g = zero_matrix(K, d, d)
  for i = 1:d
    t = trace(b[i]*b[i])
    g[i, i] = t
    for j = (i + 1):d
      t = trace(b[i]*b[j])
      g[i, j] = t
      g[j, i] = t
    end
  end
  O.trace_mat = g
  return deepcopy(g)
end

################################################################################
#
#  Dedekind criterion
#
################################################################################

function nf_elem_poly_to_fq_poly(R::FqPolyRing, m::Union{NfToFqMor, NfRelToFqMor}, f::Generic.Poly{T}) where {T <: Union{nf_elem, NfRelElem}}
  @assert codomain(m) == base_ring(R)
  @assert domain(m) == base_ring(parent(f))

  g = zero(R)
  for i = 0:degree(f)
    setcoeff!(g, i, m(coeff(f, i)))
  end
  return g
end

function nf_elem_poly_to_fq_nmod_poly(R::FqNmodPolyRing, m::NfToFqNmodMor, f::Generic.Poly{nf_elem})
  @assert codomain(m) == base_ring(R)
  @assert domain(m) == base_ring(parent(f))

  g = zero(R)
  for i = 0:degree(f)
    setcoeff!(g, i, m(coeff(f, i)))
  end
  return g
end

function fq_nmod_poly_to_nf_elem_poly(R::Generic.PolyRing{nf_elem}, m::InverseMap, f::fq_nmod_poly)
  @assert codomain(m) == base_ring(R)
  @assert domain(m) == base_ring(parent(f))

  g = zero(R)
  for i = 0:degree(f)
    setcoeff!(g, i, m(coeff(f, i)))
  end
  return g
end

function fq_poly_to_nf_elem_poly(R::Generic.PolyRing{T}, m::InverseMap, f::fq_poly) where {T <: Union{nf_elem, NfRelElem}}
  @assert codomain(m) == base_ring(R)
  @assert domain(m) == base_ring(parent(f))

  g = zero(R)
  for i = 0:degree(f)
    setcoeff!(g, i, m(coeff(f, i)))
  end
  return g
end

# Algorithm IV.6. in "Berechnung relativer Ganzheitsbasen mit dem
# Round-2-Algorithmus" by C. Friedrichs.
function dedekind_test(O::NfRelOrd, p::Union{NfOrdIdl, NfRelOrdIdl}, compute_order::Type{Val{S}} = Val{true}) where S
  !isequation_order(O) && error("Order must be an equation order")
  !issimple(O) && error("Not implemented for non-simple extensions")

  L = nf(O)
  K = base_ring(L)
  T = L.pol
  Kx = parent(T)
  OK = maximal_order(K)
  F, mF = ResidueField(OK, p)
  mmF = extend(mF, K)
  immF = inv(mmF)
  Fy, y = PolynomialRing(F,"y", cached=false)

  Tmodp = nf_elem_poly_to_fq_poly(Fy, mmF, T)
  fac = factor(Tmodp)
  g = Kx(1)
  for (t, e) in fac
    mul!(g, g, fq_poly_to_nf_elem_poly(Kx, immF, t))
  end
  gmodp = nf_elem_poly_to_fq_poly(Fy, mmF, g)
  hmodp = divexact(Tmodp, gmodp)
  h = fq_poly_to_nf_elem_poly(Kx, immF, hmodp)
  a = anti_uniformizer(p)
  f = a*(g*h - T)
  fmodp = nf_elem_poly_to_fq_poly(Fy, mmF, f)

  d = gcd(fmodp, gcd(gmodp, hmodp))

  if compute_order == Val{false}
    return isone(d)
  else
    if isone(d)
      return true, O
    end

    Umodp = divexact(Tmodp, d)
    U = fq_poly_to_nf_elem_poly(Kx, immF, Umodp)
    PM = PseudoMatrix(representation_matrix(a*U(gen(L))), [ K(1)*OK for i = 1:degree(O) ])
    PN = vcat(basis_pmat(O), PM)
    PN = sub(pseudo_hnf(PN, :lowerleft, true), degree(O) + 1:2*degree(O), 1:degree(O))
    OO = Order(L, PN)
    OO.isequation_order = false
    return false, OO
  end
end

dedekind_ispmaximal(O::NfRelOrd, p::Union{NfOrdIdl, NfRelOrdIdl}) = dedekind_test(O, p, Val{false})

dedekind_poverorder(O::NfRelOrd, p::Union{NfOrdIdl, NfRelOrdIdl}) = dedekind_test(O, p)[2]

################################################################################
#
#  p-overorder
#
################################################################################

doc"""
***
      poverorder(O::NfRelOrd, p::Union{NfOrdIdl, NfRelOrdIdl}) -> NfRelOrd

> This function tries to find an order that is locally larger than $\mathcal O$
> at the prime $p$.
"""
function poverorder(O::NfRelOrd, p::Union{NfOrdIdl, NfRelOrdIdl})
  if isequation_order(O) && issimple(O)
    return dedekind_poverorder(O, p)
  else
    return ring_of_multipliers(pradical(O, p))
  end
end

################################################################################
#
#  p-maximal overorder
#
################################################################################

doc"""
***
      pmaximal_overorder(O::NfRelOrd, p::Union{NfOrdIdl, NfRelOrdIdl}) -> NfRelOrd

> This function finds a $p$-maximal order $R$ containing $\mathcal O$.
"""
function pmaximal_overorder(O::NfRelOrd, p::Union{NfOrdIdl, NfRelOrdIdl})
  d = discriminant(O)
  OO = poverorder(O, p)
  dd = discriminant(OO)
  while d != dd
    d = dd
    OO = poverorder(OO, p)
    dd = discriminant(OO)
  end
  return OO
end

function MaximalOrder(O::NfRelOrd)
  disc = discriminant(O)
  fac = factor(disc)
  OO = deepcopy(O)
  for (p, e) in fac
    if e == 1
      continue
    end
    OO += pmaximal_overorder(O, p)
  end
  OO.ismaximal = 1
  return OO
end

################################################################################
#
#  Addition of orders
#
################################################################################

doc"""
***
      +(R::NfRelOrd, S::NfRelOrd) -> NfRelOrd

> Given two orders $R$, $S$ of $K$, this function returns the smallest order
> containing both $R$ and $S$.
"""
function +(a::NfRelOrd{T, S}, b::NfRelOrd{T, S}) where {T, S}
  # checks
  @assert nf(a) == nf(b)
  aB = basis_pmat(a)
  bB = basis_pmat(b)
  d = degree(a)
  PM = sub(pseudo_hnf(vcat(aB, bB), :lowerleft, true), d + 1:2*d, 1:d)
  return NfRelOrd{T, S}(nf(a), PM)
end

################################################################################
#
#  Absolute to relative
#
################################################################################

function relative_order(O::NfOrd, m::NfRelToNf)
  L = domain(m)
  Labs = codomain(m)
  @assert nf(O) == Labs
  K = base_ring(L)
  OK = maximal_order(K)
  mm = inv(m)
  B = basis(O, Val{false})
  d = degree(L)
  dabs = degree(Labs)
  M = zero_matrix(K, dabs, d)
  for i = 1:dabs
    elem_to_mat_row!(M, i, mm(Labs(B[i])))
  end
  PM = sub(pseudo_hnf(PseudoMatrix(M), :lowerleft, true), (dabs - d + 1):dabs, 1:d)
  return NfRelOrd{typeof(PM.matrix[1, 1]), typeof(PM.coeffs[1])}(L, PM)
end

################################################################################
#
#  Simple to non-simple
#
################################################################################

function non_simple_order(O::NfRelOrd, m::NfRelToNfRel_nsMor)
  L = domain(m)
  L_ns = codomain(m)
  @assert nf(O) == L
  K = base_ring(L)
  B = basis_nf(O, Val{false})
  d = degree(L)
  M = zero_matrix(K, d, d)
  for i = 1:d
    elem_to_mat_row!(M, i, m(L(B[i])))
  end
  PM = pseudo_hnf(PseudoMatrix(M, Hecke.basis_pmat(O).coeffs), :lowerleft, true)
  return NfRelOrd{typeof(PM.matrix[1, 1]), typeof(PM.coeffs[1])}(L_ns, PM)
end

################################################################################
#
#  Denominator in an order
#
################################################################################

doc"""
    denominator(a::RelativeElement, O::NfRelOrd) -> fmpz

Returns the smallest positive integer $k$ such that $k \cdot a$ is contained in
$\mathcal O$.
"""
function denominator(a::RelativeElement, O::NfRelOrd)
  t = zero_matrix(base_ring(nf(O)), 1, degree(O))
  elem_to_mat_row!(t, 1, a)
  t = t*basis_mat_inv(O, Val{false})
  d = fmpz(1)
  icv = inv_coeff_ideals(O, Val{false})
  for i = 1:degree(O)
    tt = icv[i]*t[1, i]
    tt = simplify(tt)
    d = lcm(d, denominator(tt))
  end
  return d
end

################################################################################
#
#  Random elements
#
################################################################################

function rand(O::NfRelOrd, B::Int)
  pb = pseudo_basis(O, Val{false})
  z = nf(O)()
  for i = 1:degree(O)
    t = rand(pb[i][2], B)
    z += t*pb[i][1]
  end
  return O(z)
end

###############################################################################
#
#  From the order of a relative extension
#  to the absolute maximal order
#
###############################################################################

function _maximal_absolute_order_from_relative(L::NfRel_ns)

  #We compute the absolute extension and the maps
  S,mS=simple_extension(L)
  K,mK=absolute_field(S, false)

  #we compute the relative maximal order of L and of the base field
  OL=maximal_order(L)
  O=maximal_order(L.base_ring)

  #take the basis
  basisL=[OL.pseudo_basis[i][1]*denominator(OL.pseudo_basis[i][2]) for i=1:degree(L.pol)]#basis(S)
  basisO=[L(O.basis_nf[i]) for i=1:degree(O)]

  #and bring them to K
  nbasisL=[mK(mS\(el)) for el in basis(L)]#[mK(el) for el in basisL]
  nbasisO=[mK(mS\(el)) for el in basisO]

  #construct the order generated by the products
  cbasis=[x*y for x in nbasisL for y in nbasisO]

  Ord=Order(K, cbasis)
  #println("The elements give an order")
  return MaximalOrder(Ord)

end


text:
module ContinuumArrays
using IntervalSets, LinearAlgebra, LazyArrays, FillArrays, BandedMatrices, QuasiArrays, InfiniteArrays
import Base: @_inline_meta, @_propagate_inbounds_meta, axes, getindex, convert, prod, *, /, \, +, -, ==, ^,
                IndexStyle, IndexLinear, ==, OneTo, tail, similar, copyto!, copy, diff,
                first, last, show, isempty, findfirst, findlast, findall, Slice, union, minimum, maximum, sum, _sum,
                getproperty, isone, iszero, zero, abs, <, ≤, >, ≥, string
import Base.Broadcast: materialize, BroadcastStyle, broadcasted
import LazyArrays: MemoryLayout, Applied, ApplyStyle, flatten, _flatten, colsupport, most, combine_mul_styles, AbstractArrayApplyStyle,
                        adjointlayout, arguments, _mul_arguments, call, broadcastlayout, layout_getindex, UnknownLayout,
                        sublayout, sub_materialize, ApplyLayout, BroadcastLayout, combine_mul_styles, applylayout,
                        simplifiable, _simplify
import LinearAlgebra: pinv, dot, norm2
import BandedMatrices: AbstractBandedLayout, _BandedMatrix
import FillArrays: AbstractFill, getindex_value, SquareEye
import ArrayLayouts: mul
import QuasiArrays: cardinality, checkindex, QuasiAdjoint, QuasiTranspose, Inclusion, SubQuasiArray,
                    QuasiDiagonal, MulQuasiArray, MulQuasiMatrix, MulQuasiVector, QuasiMatMulMat,
                    ApplyQuasiArray, ApplyQuasiMatrix, LazyQuasiArrayApplyStyle, AbstractQuasiArrayApplyStyle, AbstractQuasiLazyLayout,
                    LazyQuasiArray, LazyQuasiVector, LazyQuasiMatrix, LazyLayout, LazyQuasiArrayStyle, _factorize
import InfiniteArrays: Infinity

export Spline, LinearSpline, HeavisideSpline, DiracDelta, Derivative, ℵ₁, Inclusion, Basis, WeightedBasis, grid, transform, affine

####
# Interval indexing support
####
struct AlephInfinity{N} <: Integer end

isone(::AlephInfinity) = false
iszero(::AlephInfinity) = false

==(::AlephInfinity, ::Int) = false
==(::Int, ::AlephInfinity) = false

*(::AlephInfinity{N}, ::AlephInfinity{N}) where N = AlephInfinity{N}()
*(::AlephInfinity{N}, ::Infinity) where N = AlephInfinity{N}()
*(::Infinity, ::AlephInfinity{N}) where N = AlephInfinity{N}()
function *(a::Integer, b::AlephInfinity)
    a > 0 || throw(ArgumentError("$a is negative"))
    b
end

*(a::AlephInfinity, b::Integer) = b*a

abs(a::AlephInfinity) = a
zero(::AlephInfinity) = 0

for OP in (:<, :≤)
    @eval begin
        $OP(::Real, ::AlephInfinity) = true
        $OP(::AlephInfinity, ::Real) = false
    end
end

for OP in (:>, :≥)
    @eval begin
        $OP(::Real, ::AlephInfinity) = false
        $OP(::AlephInfinity, ::Real) = true
    end
end

const ℵ₁ = AlephInfinity{1}()

string(::AlephInfinity{1}) = "ℵ₁"

show(io::IO, F::AlephInfinity{1}) where N =
    print(io, "ℵ₁")

const QMul2{A,B} = Mul{<:Any, <:Any, <:A,<:B}
const QMul3{A,B,C} = Mul{<:AbstractQuasiArrayApplyStyle, <:Tuple{A,B,C}}

cardinality(::AbstractInterval) = ℵ₁
*(ℵ::AlephInfinity) = ℵ

Inclusion(d::AbstractInterval{T}) where T = Inclusion{float(T)}(d)
first(S::Inclusion{<:Any,<:AbstractInterval}) = leftendpoint(S.domain)
last(S::Inclusion{<:Any,<:AbstractInterval}) = rightendpoint(S.domain)

norm2(x::Inclusion{T,<:AbstractInterval}) where T = sqrt(dot(x,x))

function dot(x::Inclusion{T,<:AbstractInterval}, y::Inclusion{V,<:AbstractInterval}) where {T,V}
    x == y || throw(DimensionMismatch("first quasivector has axis $(x) which does not match the axis of the second, $(y)."))
    TV = promote_type(T,V)
    isempty(x) && return zero(TV)
    a,b = endpoints(x.domain)
    convert(TV, b^3 - a^3)/3
end

for find in (:findfirst, :findlast)
    @eval $find(f::Base.Fix2{typeof(isequal)}, d::Inclusion) = f.x in d.domain ? f.x : nothing
end

function findall(f::Base.Fix2{typeof(isequal)}, d::Inclusion)
    r = findfirst(f,d)
    r === nothing ? eltype(d)[] : [r]
end

function checkindex(::Type{Bool}, inds::Inclusion{<:Any,<:AbstractInterval}, r::Inclusion{<:Any,<:AbstractInterval})
    @_propagate_inbounds_meta
    isempty(r) | (checkindex(Bool, inds, first(r)) & checkindex(Bool, inds, last(r)))
end

###
# Maps
###

# Affine map represents A*x .+ b
abstract type AbstractAffineQuasiVector{T,AA,X,B} <: AbstractQuasiVector{T} end

show(io::IO, ::MIME"text/plain", a::AbstractAffineQuasiVector) = print(io, "$(a.A) * $(a.x) .+ ($(a.b))")

struct AffineQuasiVector{T,AA,X,B} <: AbstractAffineQuasiVector{T,AA,X,B}
    A::AA
    x::X
    b::B
end

AffineQuasiVector(A::AA, x::X, b::B) where {AA,X,B} =
    AffineQuasiVector{promote_type(eltype(AA), eltype(X), eltype(B)),AA,X,B}(A,x,b)

AffineQuasiVector(A, x) = AffineQuasiVector(A, x, zero(promote_type(eltype(A),eltype(x))))
AffineQuasiVector(x) = AffineQuasiVector(one(eltype(x)), x)

AffineQuasiVector(A, x::AffineQuasiVector, b) = AffineQuasiVector(A*x.A, x.x, A*x.b .+ b)

axes(A::AbstractAffineQuasiVector) = axes(A.x)
affine_getindex(A, k) = A.A*A.x[k] .+ A.b
getindex(A::AbstractAffineQuasiVector, k::Number) = affine_getindex(A, k)
function getindex(A::AbstractAffineQuasiVector, k::Inclusion)
    @boundscheck A.x[k] # throws bounds error if k ≠ x
    A
end

getindex(A::AbstractAffineQuasiVector, ::Colon) = copy(A)

copy(A::AbstractAffineQuasiVector) = A

inbounds_getindex(A::AbstractAffineQuasiVector{<:Any,<:Any,<:Inclusion}, k::Number) = A.A*k .+ A.b
isempty(A::AbstractAffineQuasiVector) = isempty(A.x)
==(a::AbstractAffineQuasiVector, b::AbstractAffineQuasiVector) = a.A == b.A && a.x == b.x && a.b == b.b

BroadcastStyle(::Type{<:AbstractAffineQuasiVector}) = LazyQuasiArrayStyle{1}()

for op in(:*, :\, :+, :-)
    @eval broadcasted(::LazyQuasiArrayStyle{1}, ::typeof($op), a::Number, x::Inclusion) = broadcast($op, a, AffineQuasiVector(x))
end
for op in(:/, :+, :-)
    @eval broadcasted(::LazyQuasiArrayStyle{1}, ::typeof($op), x::Inclusion, a::Number) = broadcast($op, AffineQuasiVector(x), a)
end

broadcasted(::LazyQuasiArrayStyle{1}, ::typeof(*), a::Number, x::AbstractAffineQuasiVector) = AffineQuasiVector(a, x)
broadcasted(::LazyQuasiArrayStyle{1}, ::typeof(\), a::Number, x::AbstractAffineQuasiVector) = AffineQuasiVector(inv(a), x)
broadcasted(::LazyQuasiArrayStyle{1}, ::typeof(/), x::AbstractAffineQuasiVector, a::Number) = AffineQuasiVector(inv(a), x)
broadcasted(::LazyQuasiArrayStyle{1}, ::typeof(+), a::Number, x::AbstractAffineQuasiVector) = AffineQuasiVector(one(eltype(x)), x, a)
broadcasted(::LazyQuasiArrayStyle{1}, ::typeof(+), x::AbstractAffineQuasiVector, b::Number) = AffineQuasiVector(one(eltype(x)), x, b)
broadcasted(::LazyQuasiArrayStyle{1}, ::typeof(-), a::Number, x::AbstractAffineQuasiVector) = AffineQuasiVector(-one(eltype(x)), x, a)
broadcasted(::LazyQuasiArrayStyle{1}, ::typeof(-), x::AbstractAffineQuasiVector, b::Number) = AffineQuasiVector(one(eltype(x)), x, -b)

function checkindex(::Type{Bool}, inds::Inclusion{<:Any,<:AbstractInterval}, r::AbstractAffineQuasiVector)
    @_propagate_inbounds_meta
    isempty(r) | (checkindex(Bool, inds, first(r)) & checkindex(Bool, inds, last(r)))
end

affine_igetindex(d, x) = d.A \ (x .- d.b)
igetindex(d::AbstractAffineQuasiVector, x) = affine_igetindex(d, x)

for find in (:findfirst, :findlast, :findall)
    @eval $find(f::Base.Fix2{typeof(isequal)}, d::AbstractAffineQuasiVector) = $find(isequal(igetindex(d, f.x)), d.x)
end

minimum(d::AbstractAffineQuasiVector) = signbit(d.A) ? last(d) : first(d)
maximum(d::AbstractAffineQuasiVector) = signbit(d.A) ? first(d) : last(d)

union(d::AbstractAffineQuasiVector) = Inclusion(minimum(d)..maximum(d))

struct AffineMap{T,D,R} <: AbstractAffineQuasiVector{T,T,D,T}
    domain::D
    range::R
end

AffineMap(domain::AbstractQuasiVector{T}, range::AbstractQuasiVector{V}) where {T,V} =
    AffineMap{promote_type(T,V), typeof(domain),typeof(range)}(domain,range)

measure(x::Inclusion) = last(x)-first(x)

function getproperty(A::AffineMap, d::Symbol)
    d == :x && return A.domain
    d == :A && return measure(A.range)/measure(A.domain)
    d == :b && return (last(A.domain)*first(A.range) - first(A.domain)*last(A.range))/measure(A.domain)
    getfield(A, d)
end

function getindex(A::AffineMap, k::Number)
    # ensure we exactly hit range
    k == first(A.domain) && return first(A.range)
    k == last(A.domain) && return last(A.range)
    affine_getindex(A, k)
end

function igetindex(A::AffineMap, k::Number)
    # ensure we exactly hit range
    k == first(A.range) && return first(A.domain)
    k == last(A.range) && return last(A.domain)
    affine_igetindex(A, k)
end

first(A::AffineMap) = first(A.range)
last(A::AffineMap) = last(A.range)

affine(a::AbstractQuasiVector, b::AbstractQuasiVector) = AffineMap(a, b)
affine(a, b::AbstractQuasiVector) = affine(Inclusion(a), b)
affine(a::AbstractQuasiVector, b) = affine(a, Inclusion(b))
affine(a, b) = affine(Inclusion(a), Inclusion(b))

# mapped vectors
const AffineMappedQuasiVector = SubQuasiArray{<:Any, 1, <:Any, <:Tuple{AbstractAffineQuasiVector}}
const AffineMappedQuasiMatrix = SubQuasiArray{<:Any, 2, <:Any, <:Tuple{AbstractAffineQuasiVector,Slice}}

==(a::AffineMappedQuasiVector, b::AffineMappedQuasiVector) = parentindices(a) == parentindices(b) && parent(a) == parent(b)

_sum(V::AffineMappedQuasiVector, ::Colon) = parentindices(V)[1].A \ sum(parent(V))

# pretty print for bases
show(io::IO, P::AffineMappedQuasiMatrix) = print(io, "$(parent(P)) affine mapped to $(parentindices(P)[1].x.domain)")
show(io::IO, P::AffineMappedQuasiVector) = print(io, "$(parent(P)) affine mapped to $(parentindices(P)[1].x.domain)")
show(io::IO, ::MIME"text/plain", P::AffineMappedQuasiMatrix) = show(io, P)
show(io::IO, ::MIME"text/plain", P::AffineMappedQuasiVector) = show(io, P)

const QInfAxes = Union{Inclusion,AbstractAffineQuasiVector}

sub_materialize(_, V::AbstractQuasiArray, ::Tuple{QInfAxes}) = V
sub_materialize(_, V::AbstractQuasiArray, ::Tuple{QInfAxes,QInfAxes}) = V
sub_materialize(_, V::AbstractQuasiArray, ::Tuple{Any,QInfAxes}) = V
sub_materialize(_, V::AbstractQuasiArray, ::Tuple{QInfAxes,Any}) = V

include("operators.jl")
include("bases/bases.jl")

end

ext: jl


text:
using CSTParser, Test, CSTParser.Tokenize
using CSTParser: parse, remlineinfo!, span, headof, kindof, valof

@testset "CSTParser" begin
    include("spec.jl")
    include("parser.jl")
    include("interface.jl")
    include("display.jl")
    include("iterate.jl")
    if VERSION >= v"1.6"
        include("check_base.jl")
        include("errparse.jl")
    end
end


text:
# julia1
using Plots
:)mozda


text:
"""
```julia
@def_sim sim_name default_type_params begin
    obs_typedef
end
```

Given a type definition for a single observation in a simulation
(`obs_typedef`), evaluate that type definition as is, but also creates a second
type named `sim_name` as well as various methods on the new type.

The fields of `sim_name` will have the same name as the fields of
`obs_typedef`, but will be arrays of whatever the type of the corresponding
`obs_typedef` field was. The intention is for `sim_name` to be a struct of
arrays (see https://en.wikipedia.org/wiki/AOS_and_SOA). If you want an array
of structs, just simply collect an array of instances of the type defined in
`obs_typedef`. The struct of arrays storage format has better cache efficiency
and data locality if you want to operate on all values of a particular field at
once, rather than all the fields of a particular value.

In addition to the new type `sim_name`, the following methods will be defined:

- `sim_name(sz::NTuple{N,Int})`. This is a constructor for `sim_name` that
  allocates arrays of size `sz` for each field. If `obs_typedef` inlcuded any
  type parameters, then the default values (specified in `default_type_params`)
  will be used.
- `Base.endof(::sim_name)`: equal to the length of any of its fields
- `Base.length(::sim_name)`: equal to the length of any of its fields
- The iterator protocol for `sim_name`. The type of each element of the
  iterator is the type defined in `obs_typedef`. This amounts tho defining the
  following methods
    - `Base.start(::sim_name)::Int`
    - `Base.next(::sim_name, ::Int)::Tuple{Observation,Int}`
    - `Base.done(::sim_name, ::Int)::Bool`
- `Base.getindex(sim::sim_name, ix::Int)`. This implements _linear indexing_
  for `sim_name` and will return an instance of the type defined in
  `obs_typedef`

## Example

NOTE: the `using MacroTools`  and call to `MacroTools.prettify` is not
necessary and is only used here to clean up the output so it is easier to read

```
julia> using MacroTools

julia> macroexpand(:(@def_sim Simulation (T => Float64,) struct Observation{T<:Number}
           c::T
           k::T
           i_z::Int
       end
       )) |> MacroTools.prettify
quote
    struct Simulation{prairiedog, T <: Number}
        c::Array{T, prairiedog}
        k::Array{T, prairiedog}
        i_z::Array{Int, prairiedog}
    end
    function Simulation{prairiedog}(sz::NTuple{prairiedog, Int})
        c = Array{Float64, prairiedog}(sz)
        k = Array{Float64, prairiedog}(sz)
        i_z = Array{Int, prairiedog}(sz)
        Simulation(c, k, i_z)
    end
    struct Observation{T <: Number}
        c::T
        k::T
        i_z::Int
    end
    Base.endof(sim::Simulation) = length(sim.c)
    Base.length(sim::Simulation) = endof(sim)
    Base.start(sim::Simulation) = 1
    Base.next(sim::Simulation, ix::Int) = (sim[ix], ix + 1)
    Base.done(sim::Simulation, ix::Int) = ix >= length(sim)
    function Base.getindex(sim::Simulation, ix::Int)
        \$(Expr(:boundscheck, true))
        if ix > length(sim)
            throw(BoundsError("\$(length(sim))-element Simulation at index \$(ix)"))
        end
        \$(Expr(:boundscheck, :pop))
        \$(Expr(:inbounds, true))
        out = Observation(sim.c[ix], sim.k[ix], sim.i_z[ix])
        \$(Expr(:inbounds, :pop))
        return out
    end
end
```
"""
macro def_sim(sim_name, default_type_params, obs_typedef)

    N_sym = gensym("N")
    # construct default_type_param map
    tp_map = Dict{Symbol,Symbol}()
    if default_type_params.head != :tuple
        m = "`default_type_params` must be a tuple of `Pair`s."
        m = string(m, " If you ony have one pair, use the notation `(a => b,)`")
        error(m)
    end
    for pair in default_type_params.args
        if pair.head != :call || pair.args[1] != :(=>)
            error("Expected tuple of the form (t1 => t2, x1 => x2)")
        end
        tp_map[pair.args[2]] = pair.args[3]
    end

    ex = obs_typedef   # simplify name
    obs_name = ex.args[2].args[1]

    obs_fields = ex.args[3].args
    sim_fields = Expr(:block)
    for field in obs_fields
        isa(field, LineNumberNode) && continue

        if field.head == :(::)
            name = field.args[1]
            typ = field.args[2]
            push!(sim_fields.args, :($(name)::Array{$typ,$N_sym}))
        end
    end

    sim_typename = Expr(:curly, sim_name, N_sym)
    if ex.args[2].head == :curly
        type_params = ex.args[2].args[2:end]
        append!(sim_typename.args, type_params)
    else
        type_params = []
    end
    sim_type = Expr(:struct, ex.args[1], sim_typename, sim_fields)
    # sim_type = quote
    #     struct $(sim_name){N,$(type_params...)}
    #         $sim_fields
    #     end
    # end

    body = Expr(:block, map(sim_fields.args) do expr
        tp_name = expr.args[2].args[2]
        _eltype = get(tp_map, tp_name, tp_name)
        arr_type = :(Array{$_eltype,$N_sym})
        Expr(:(=), expr.args[1], Expr(:call, arr_type, :undef, :sz))
    end..., :($(sim_name)($([i.args[1] for i in sim_fields.args]...))))

    sim_constructor = :(function $(sim_name)(sz::NTuple{$N_sym,Int}) where $N_sym
    $body
    end)

    getindex_out_expr = Expr(:call, obs_name, map(sim_fields.args) do expr
        Expr(:ref, :(sim.$(expr.args[1])), :ix)
    end...)

    others = quote
        function Base.getindex(sim::$(sim_name), ix::Int)
            @boundscheck begin
                if ix > length(sim)
                    throw(BoundsError("$(length(sim))-element Simulation at index $ix"))
                end
            end
            @inbounds out = $getindex_out_expr
            return out
        end
        function Base.iterate(sim::$(sim_name), ix::Int=1)
            ix >= length(sim) && return nothing
            (sim[ix], ix+1)
        end
        Base.length(sim::$(sim_name)) = length(sim.$(sim_fields.args[1].args[1]))
    end
    out = esc(Expr(:block, sim_type, sim_constructor, obs_typedef, others))
end


text:
import Pkg; Pkg.activate("../repo/TuringBnpBenchmarks")

using Turing
using Distributions
using DelimitedFiles
import Random

idx, Y = let
  data = readdlm("patients.csv", ',')
  Int.(data[:, 1]), data[:, 2:end]
end

# FAM without missing data.
@model fam(y, idx, K, L0, L1, phi=0) = begin
  I = length(unique(idx))
  Nsum = length(idx)
  J = size(y, 2)

  alpha ~ Gamma(1, 1)
  v ~ filldist(Beta(alpha/K, 1), K)
  Z ~ arraydist([Bernoulli(v[k]) for j in 1:J, k in 1:K])
  lambda = tzeros(Int, Nsum)

  sigma ~ filldist(LogNormal(0, 1), I)
  delta0 ~ filldist(LogNormal(0, 0.1), L0)
  delta1 ~ filldist(LogNormal(0, 0.1), L1)
  mu0 = -cumsum(delta0)
  mu1 = cumsum(delta1)

  W ~ filldist(Dirichlet(K, 1/K), I)

  eta0 = Array{Float64, 3}(undef, I, J, L0)
  eta1 = Array{Float64, 3}(undef, I, J, L1)
  for i in 1:I
    for j in 1:J
      eta0[i, j, :] ~ Dirichlet(L0, 1)
      eta1[i, j, :] ~ Dirichlet(L1, 1)
    end
  end

  for j in 1:J
    for n in 1:Nsum
      i = idx[n]

      lambda[n] ~ Categorical(W[:, i])
      z = Z[j, lambda[n]]
      Lz = (z == 0) ? L0 : L1
      muz = (z == 0) ? mu0 : mu1
      etaz = (z == 0) ? eta0 : eta1

      y[n, j] ~ UnivariateGMM(muz, fill(sigma[i], Lz), Categorical(etaz[i, j, :]))
    end
  end
end

Random.seed!(0)
@time chain = begin
    K = 25
    L0 = 5
    L1 = 5
    nburn = 10
    nsamps = 10
    iterations = nburn + nsamps
    sample(fam(Y, idx, K, L0, L1), Gibbs(PG(100)), iterations)
end;

nothing

#=
group(chain, :Z).value.data
=#


text:
using Recruiters
using PeaceCypher
using DemeNet: Signer, DemeSpec, Deme, ID, AbstractID, Certificate, Intent, Profile

demespec = DemeSpec("PeaceDeme",:default,:PeaceCypher,:default,:PeaceCypher,:Recruiters)
deme = Deme(demespec)

maintainer = Signer(deme,"maintainer")
server = Signer(deme,"server")

MAINTAINER_ID = maintainer.id
SERVER_ID = server.id
TOOKEN_PORT = 2006
CERTIFIER_PORT = 2007

config = CertifierConfig(MAINTAINER_ID,SERVER_ID,TOOKEN_PORT,CERTIFIER_PORT)
certifier = Certifier(config,deme,server)

sleep(1)

tooken = 123333

addtooken(config,deme,tooken,maintainer)

# Now the maintainer shares the tooken, demespec and ledger port with new member

profile = Profile(Dict("uuid"=>1234431))

member = Signer(deme,"member")
profilecert = Certificate(profile,member)

certify(CERTIFIER_PORT,SERVER_ID,deme,profilecert,tooken)

# The recruiter then can get registration
@show take!(certifier.tookencertifier.tickets)

# ### Let's now test the ticket API

tooken = 123456

# ### The maintainer API

addtooken(config,deme,tooken,maintainer)
invite = ticket(demespec,CERTIFIER_PORT,SERVER_ID,tooken)

sleep(1)

# ### Now from users point of view
profile = Profile(Dict("uuid"=>11223344))
register(invite,profile)
@show take!(certifier.tookencertifier.tickets)

# And now lastly the sending of the ticket over email

# sendinvite(config,deme,"akels14@gmail.com",SMTPConfig(),maintainer)


text:
# This hack will go away once CLBLAS doesn't
# rely on the clunky cl.CLArray type in OpenCL anymore. Sadly, we need cl.Buffer for that, which  overwrites the default
# constructor registering a finalizer, which we don't want here, so we need to
# get the real constructor.
@generated function clbuffer(ptr::OwnedPtr{T}, len::Integer) where T
    Expr(:new, cl.Buffer{T}, true, :(pointer(ptr)), :(len), :(ptr.mapped), Ptr{T}(C_NULL))
end

function clbuffer(A::CLArray)
    ptr = pointer(A)
    clbuffer(ptr, length(A))
end

import CLFFT

# figure out a gc safe way to store plans.
# weak refs won't work, since the caching should keep them alive.
# But at the end, we need to free all of these, otherwise CLFFT will crash
# at closing time.
# An atexit hook here, which will empty the dictionary seems to introduce racing
# conditions.
#const plan_dict = Dict()
import Base: *, plan_ifft!, plan_fft!, plan_fft, plan_ifft, size, plan_bfft, plan_bfft!

immutable CLFFTPlan{Direction, Inplace, T, N} <: Base.FFTW.FFTWPlan{T, Direction, Inplace}
    plan::CLFFT.Plan{T}
    function (::Type{CLFFTPlan{Direction, Inplace}})(A::CLArray{T, N}) where {T, N, Direction, Inplace}
        ctx = global_context(A)
        p = CLFFT.Plan(T, ctx, size(A))
        CLFFT.set_layout!(p, :interleaved, :interleaved)
        if Inplace
            CLFFT.set_result!(p, :inplace)
        else
            CLFFT.set_result!(p, :outofplace)
        end
        CLFFT.set_scaling_factor!(p, Direction, 1f0)
        CLFFT.bake!(p, global_queue(A))
        new{Direction, Inplace, T, N}(p)
    end
end
size(x::CLFFTPlan) = (CLFFT.lengths(x.plan)...,)

# ignore flags, but have them to make it base compatible.
# TODO can we actually implement the flags?
function plan_fft(A::CLArray; flags = nothing, timelimit = Inf)
    CLFFTPlan{:forward, false}(A)
end
function plan_fft!(A::CLArray; flags = nothing, timelimit = Inf)
    CLFFTPlan{:forward, true}(A)
end
function plan_bfft(A::CLArray, region; flags = nothing, timelimit = Inf)
    CLFFTPlan{:backward, false}(A)
end
function plan_bfft!(A::CLArray, region; flags = nothing, timelimit = Inf)
    CLFFTPlan{:backward, true}(A)
end

const _queue_ref = Vector{cl.CmdQueue}(1)
function *(plan::CLFFTPlan{Direction, true, T, N}, A::CLArray{T, N}) where {T, N, Direction}
    _queue_ref[] = global_queue(A)
    CLFFT.enqueue_transform(plan.plan, Direction, _queue_ref, clbuffer(A), nothing)
    A
end
function *(plan::CLFFTPlan{Direction, false, T, N}, A::CLArray{T, N}) where {T, N, Direction}
    _queue_ref[] = global_queue(A)
    y = typeof(A)(size(plan))
    CLFFT.enqueue_transform(plan.plan, Direction, _queue_ref, clbuffer(A), clbuffer(y))
    y
end

# Enable BLAS support via GPUArrays
import CLBLAS
import GPUArrays: blas_module, blasbuffer
blas_module(::CLArray) = CLBLAS
blasbuffer(A::CLArray) = cl.CLArray(clbuffer(A), global_queue(A), size(A))


text:
@testset "CacheImage" begin
    @test (CacheImage <: AugmentorBase.AffineOperation) == false
    @test (CacheImage <: AugmentorBase.ImageOperation) == true
    @test typeof(@inferred(CacheImage())) <: CacheImage
    @test str_show(CacheImage()) == "AugmentorBase.CacheImage()"
    @test str_showconst(CacheImage()) == "CacheImage()"
    @test str_showcompact(CacheImage()) == "Cache into temporary buffer"

    @test @inferred(AugmentorBase.applyeager(CacheImage(),square)) === square
    @test @inferred(AugmentorBase.applyeager(CacheImage(),(square,))) === (square,)
    @test @inferred(AugmentorBase.applyeager(CacheImage(),(square,square2))) === (square,square2)

    # SubArray
    v = view(square, :, :)
    img = @inferred AugmentorBase.applyeager(CacheImage(), v)
    @test typeof(img) <: Array
    @test eltype(img) == eltype(v)
    @test img !== square
    @test img == square
    # Identidy ranges
    v = view(square, IdentityRange(1:3), IdentityRange(1:3))
    img = @inferred AugmentorBase.applyeager(CacheImage(), v)
    @test typeof(img) <: OffsetArray
    @test eltype(img) == eltype(v)
    @test img == OffsetArray(square, 0, 0)
    # Affine
    v = AugmentorBase.prepareaffine(square)
    img = @inferred AugmentorBase.applyeager(CacheImage(), v)
    @test typeof(img) <: OffsetArray
    @test eltype(img) == eltype(v)
    @test img == OffsetArray(square, 0, 0)
    # Array and SubArray
    v = view(square, :, :)
    tmp,img = @inferred AugmentorBase.applyeager(CacheImage(), (square,v))
    @test typeof(img) <: Array
    @test eltype(img) == eltype(v)
    @test tmp === square
    @test img !== square
    @test img == square
    # OffsetArray
    o = OffsetArray(square, (-1,2))
    @test @inferred(AugmentorBase.applyeager(CacheImage(),o)) === o
    @test @inferred(AugmentorBase.applyeager(CacheImage(),(o,square))) === (o,square)

    @test AugmentorBase.supports_eager(CacheImage) === true
    @test AugmentorBase.supports_lazy(CacheImage) === false
    @test AugmentorBase.supports_view(CacheImage) === false
    @test AugmentorBase.supports_stepview(CacheImage) === false
    @test AugmentorBase.supports_permute(CacheImage) === false
    @test AugmentorBase.supports_affine(CacheImage) === false
    @test AugmentorBase.supports_affineview(CacheImage) === false

    @test_throws MethodError AugmentorBase.applylazy(CacheImage(), v)
    @test_throws MethodError AugmentorBase.applylazy(CacheImage(), (v,v))
    @test_throws MethodError AugmentorBase.applyview(CacheImage(), v)
    @test_throws MethodError AugmentorBase.applystepview(CacheImage(), v)
    @test_throws MethodError AugmentorBase.applypermute(CacheImage(), v)
    @test_throws MethodError AugmentorBase.applyaffine(CacheImage(), v)
    @test_throws MethodError AugmentorBase.applyaffineview(CacheImage(), v)
end

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

@testset "CacheImageInto" begin
    @test_throws UndefVarError CacheImageInto
    @test (AugmentorBase.CacheImageInto <: AugmentorBase.AffineOperation) == false
    @test (AugmentorBase.CacheImageInto <: AugmentorBase.ImageOperation) == true
    @test_throws MethodError AugmentorBase.CacheImageInto()

    @testset "single image" begin
        buf = copy(rect)
        @test typeof(@inferred(CacheImage(buf))) <: AugmentorBase.CacheImageInto
        op = @inferred CacheImage(buf)
        @test AugmentorBase.CacheImageInto(buf) === op
        @test str_show(op) == "AugmentorBase.CacheImageInto(::Array{Gray{N0f8},2})"
        @test str_showconst(op) == "CacheImage(Array{Gray{N0f8}}(2, 3))"
        op2 = @inferred CacheImage(Array{Gray{N0f8}}(undef, 2, 3))
        @test typeof(op) == typeof(op2)
        @test typeof(op.buffer) == typeof(op2.buffer)
        @test size(op.buffer) == size(op2.buffer)
        @test str_showcompact(op) == "Cache into preallocated 2×3 Array{Gray{N0f8},2} with eltype Gray{Normed{UInt8,8}}"

        v = AugmentorBase.applylazy(Resize(2,3), camera)
        res = @inferred AugmentorBase.applyeager(op, v)
        @test res == v
        @test typeof(res) <: OffsetArray
        @test parent(res) === op.buffer

        res = @inferred AugmentorBase.applyeager(op, rect)
        @test res == rect
        @test res === op.buffer

        res = @inferred AugmentorBase.applylazy(op, v)
        @test res == v
        @test typeof(res) <: OffsetArray
        @test parent(res) === op.buffer

        res = @inferred AugmentorBase.applylazy(op, rect)
        @test res == rect
        @test res === op.buffer

        @test_throws ArgumentError AugmentorBase.applyeager(op, camera)
        @test_throws MethodError AugmentorBase.applyview(CacheImage(buf), v)
        @test_throws MethodError AugmentorBase.applystepview(CacheImage(buf), v)
        @test_throws MethodError AugmentorBase.applypermute(CacheImage(buf), v)
        @test_throws MethodError AugmentorBase.applyaffine(CacheImage(buf), v)
        @test_throws MethodError AugmentorBase.applyaffineview(CacheImage(buf), v)
    end

    @testset "multiple images" begin
        buf1 = copy(square)
        buf2 = copy(rgb_rect)
        @test typeof(@inferred(CacheImage(buf1,buf2))) <: AugmentorBase.CacheImageInto
        op = @inferred CacheImage(buf1,buf2)
        @test op === @inferred CacheImage((buf1,buf2))
        @test AugmentorBase.CacheImageInto((buf1,buf2)) === op
        @test str_show(op) == "AugmentorBase.CacheImageInto((::Array{Gray{N0f8},2}, ::Array{RGB{N0f8},2}))"
        @test str_showconst(op) == "CacheImage(Array{Gray{N0f8}}(3, 3), Array{RGB{N0f8}}(2, 3))"
        op2 = @inferred CacheImage(Array{Gray{N0f8}}(undef, 3, 3), Array{RGB{N0f8}}(undef, 2, 3))
        @test typeof(op) == typeof(op2)
        @test typeof(op.buffer) == typeof(op2.buffer)
        @test size.(op.buffer) === size.(op2.buffer)
        @test str_showcompact(op) == "Cache into preallocated (3×3 Array{Gray{N0f8},2} with eltype Gray{Normed{UInt8,8}}, 2×3 Array{RGB{N0f8},2} with eltype RGB{Normed{UInt8,8}})"
        @test buf1 == square
        @test buf2 == rgb_rect
        v1 = AugmentorBase.applylazy(Resize(3,3), camera)
        v2 = AugmentorBase.applylazy(Resize(2,3), RGB.(camera))
        res = @inferred AugmentorBase.applyeager(op, (v1,v2))
        @test buf1 != square
        @test buf2 != rgb_rect
        @test res == (v1, v2)
        @test typeof(res) <: NTuple{2,OffsetArray}
        @test parent.(res) === (op.buffer[1], op.buffer[2])

        @test_throws ArgumentError AugmentorBase.applyeager(op, (camera,buf1)) #1
        @test_throws MethodError AugmentorBase.applylazy(op, v1)
        @test_throws ArgumentError AugmentorBase.applylazy(op, (buf2,buf1)) #3
        @test_throws BoundsError AugmentorBase.applylazy(op, (buf1,))
        # ?
        @test_throws DimensionMismatch AugmentorBase.applylazy(op, (v1,v1))
    end

    @test AugmentorBase.supports_eager(AugmentorBase.CacheImageInto) === true
    @test AugmentorBase.supports_lazy(AugmentorBase.CacheImageInto) === true
    @test AugmentorBase.supports_view(AugmentorBase.CacheImageInto) === false
    @test AugmentorBase.supports_stepview(AugmentorBase.CacheImageInto) === false
    @test AugmentorBase.supports_permute(AugmentorBase.CacheImageInto) === false
    @test AugmentorBase.supports_affine(AugmentorBase.CacheImageInto) === false
    @test AugmentorBase.supports_affineview(AugmentorBase.CacheImageInto) === false
end


text:
using AxisArrays
using FFTW

export rates, scales, nrates, nscales, default_rates, default_scales,
  cortical, cycoct, co, scalefilter, ratefilter

# re-express the spectral and cortical dimensions to have meta data specific to
# an axis and then make it possible to add an axis to an existing array (maybe
# have a flag to allow multiple axes of the same type)

@dimension Sc "Sc" Scale
@refunit cycoct "cyc/oct" CyclesPerOct Sc false

# NOTE: the `low` and `high` fields are used to determine which filters should
# be low- and high-pass, rather than band-pass
struct ScaleAxis
  low::typeof(1.0cycoct)
  high::typeof(1.0cycoct)
end

struct RateAxis
  low::typeof(1.0Hz)
  high::typeof(1.0Hz)
end

rates(x::MetaUnion{AxisArray}) =
  axisvalues(AxisArrays.axes(x,Axis{:rate}))[1]
nrates(x) = length(rates(x))

scales(x::MetaUnion{AxisArray}) =
  axisvalues(AxisArrays.axes(x,Axis{:scale}))[1]
nscales(x) = length(scales(x))

const default_rates = sort([-2 .^ (1:0.5:5); 2 .^ (1:0.5:5)]).*Hz
const default_scales = (2 .^ (-2:0.5:3)).*cycoct
const spect_rate = 24

struct DefaultValue
end
const default = DefaultValue()
withdefault(x,default) = x
withdefault(::DefaultValue,default) = default
function checkmatch(pos,kwd,name)
  if pos !== kwd
    error("Cannot specify both a position and keyword value for $name.")
  end
end
checkmatch(::DefaultValue,x,str) = nothing

ascycoct(x::Number) = x*cycoct
ascycoct(x::Quantity) = uconvert(cycoct,x)
ascycoct(x::Quantity{<:Any,Sc}) = x

abstract type CorticalFilter
end
abstract type CorticalFilterInv
end

function DSP.filt(filter::CorticalFilter,y::MetaAxisArray)
  @assert all(∈(axisnames(y)),(:time,:freq))
  if axisname(filter) ∈ y
    ax = axisname(filter)
    error("Input already has an axis named `$ax`. If you intended ",
          "to add a new dimension, you will have to change the name of the ",
          "axis. When you define the filter you can specify the axis name ",
          "using the `axis` keyword argument.")
  end

  fir = FIRFiltering(y,Axis{fromaxis(filter)})
  cr = initfilter(y,filter)
  for (I,H) in list_filters(fir,cr,filter)
    cr[I] = view(apply(fir,H),Base.axes(y)...)
  end

  cr
end

function DSP.filt(cinv::CorticalFilterInv,cr::MetaAxisArray)
  z_cum = FFTCum(cr,axisnames(cinv))

  filters = list_filters(z_cum,cr,cinv)

  inner_dims = size(cr)[2:end-1]
  if length(inner_dims) != ndims(cinv)
    error("When computing the inverse you must invert all scale/rate "*
          "dimensions; partial inverses are not supported.")
  end

  for (I,filter) in filters
    addfft!(z_cum,view(cr,:,I,:),filter)
  end

  removed = removeaxes(getmeta(cr),axisnames(cinv)...)
  norm = normalize!(z_cum,cr,cinv.norm)
  MetaAxisArray(removed,AxisArray(norm,AxisArrays.axes(cr)[[1,end]]))
end

# rate filters
struct TimeRateFilter <: CorticalFilter
  data::Vector{typeof(1.0Hz)}
  bandonly::Bool
  axis::Symbol
end
axisname(x::TimeRateFilter) = x.axis
fromaxis(x::TimeRateFilter) = :time
AxisArrays.axisnames(x::TimeRateFilter) = (x.axis,)
Base.length(x::TimeRateFilter) = length(x.data)

"""
    ratefilter(rates,bandonly=false,axis=:rate)

Define a temporal-rate filter bank at the given rates (assumed to be in
Hz by default, but can be specified in a different unit via `Unitful`).
The rates can be either a positional or keyword argument. To compute
use inverse of the filter-bank call `inv`. Uses default_rates
if the rates are unspecified.

Apply the filter bank using `filt`, like so:

```julia
using FileIO
using CorticalSpectralTemporalResponses

x = load("testsound.wav")
X = filt(audiospect,x)
ratef = ratefilter()
R_cr = filt(ratef,X)
```

"""
function ratefilter(r=default;rates=withdefault(r,default_rates),
  bandonly=false,axis=:rate)
  checkmatch(r,rates,"rates")
  rates = asHz.(rates)

  if axis != :rate && !occursin("rate",string(axis))
    error("Rate axis name `$axis` must contain the word 'rate'.")
  end
  TimeRateFilter(rates,bandonly,axis)
end
list_filters(fir,cr,filter::TimeRateFilter) =
  ((Axis{axisname(filter)}(i), HR)
   for (i,HR) in enumerate(rate_filters(fir,cr,axisname(filter))))

# inverse of rate filters
struct TimeRateFilterInv <: CorticalFilterInv
  rates::TimeRateFilter
  norm::Float64
end
Base.ndims(x::TimeRateFilterInv) = 1
axisname(x::TimeRateFilterInv) = axisname(x.rates)
AxisArrays.axisnames(x::TimeRateFilterInv) = axisnames(x.rates)
Base.inv(rates::TimeRateFilter;norm=0.9) = TimeRateFilterInv(rates,norm)
list_filters(z_cum,cr,rateinv::TimeRateFilterInv) =
  ((i, HR)
   for (i,HR) in enumerate(rate_filters(z_cum,cr,axisname(rateinv),use_conj=true)))

# scale filters
struct FreqScaleFilter <: CorticalFilter
  data::Vector{typeof(1.0cycoct)}
  bandonly::Bool
  axis::Symbol
end
axisname(x::FreqScaleFilter) = x.axis
fromaxis(x::FreqScaleFilter) = :freq
AxisArrays.axisnames(x::FreqScaleFilter) = (x.axis,)
list_filters(fir,cs,scales::FreqScaleFilter) =
  ((Axis{axisname(scales)}(i), [HS; zero(HS)]')
   for (i,HS) in enumerate(scale_filters(fir,cs,axisname(scales))))

"""
    scalefilter(scales,bandonly=false,axis=:rate)

Define a spectral-scale filter bank at the given scales (assumed to be in
cycles/octave). The scales can be either a positional or keyword argument. To
compute use inverse of the filter-bank call `inv`. Uses default_scales if the
scales are unspecified.

Apply the filter bank using `filt`, like so:

```julia
using FileIO
using CorticalSpectralTemporalResponses

x = load("testsound.wav")
X = filt(audiospect,x)
scalef = scalefilter()
S_cr = filt(scalef,X)
```

"""
function scalefilter(s=default;scales=withdefault(s,default_scales),
    bandonly=false,axis=:scale)
  checkmatch(s,scales,"scales")
  scales = ascycoct.(scales)

  if axis != :scale && !occursin("scale",string(axis))
    error("Scale axis name `$axis` must contain the word 'scale'.")
  end
  FreqScaleFilter(scales,bandonly,axis)
end

# inverse of scale filters
struct FreqScaleFilterInv <: CorticalFilterInv
  scales::FreqScaleFilter
  norm::Float64
end
Base.ndims(x::FreqScaleFilterInv) = 1
axisname(x::FreqScaleFilterInv) = axisname(x.scales)
AxisArrays.axisnames(x::FreqScaleFilterInv) = (axisname(x.scales),)
Base.inv(scales::FreqScaleFilter;norm=0.9) = FreqScaleFilterInv(scales,norm)
list_filters(z_cum,cr,scaleinv::FreqScaleFilterInv) =
  ((i, [HS; zero(HS)]')
   for (i,HS) in enumerate(scale_filters(z_cum,cr,axisname(scaleinv))))

# combination of both scales and rates
struct ScaleRateFilter <: CorticalFilter
  scales::FreqScaleFilter
  rates::TimeRateFilter
end

"""
    cortical(scales,rates,bandonly=false,axis=:rate)

Define a joint scale-rate filter bank at the given scales and rates (these
are assumed to be in cycles per octave and Hz respectively, but can be
specified with different units via `Unitful`). The scales and rates can be
either both as positional arguments or both as keyword arguments. To compute
the inverse of the filter-bank call `inv`.

Apply the filter bank using `filt`, like so:

```julia
using FileIO
using CorticalSpectralTemporalResponses

x = load("testsound.wav")
X = filt(audiospect,x)
cortf = cortical()
cr = filt(cortf,X)
```
"""
function cortical(s=default,r=default;
    scales=withdefault(s,default_scales),
    rates=withdefault(r,default_rates),
    bandonly=false,axes=(:scale,:rate))

    checkmatch(s,scales,"scales")
    checkmatch(r,rates,"rates")
    scales = ascycoct.(scales)
    rates = asHz.(rates)

    ScaleRateFilter(scalefilter(scales,bandonly=bandonly,axis=axes[1]),
                    ratefilter(rates,bandonly=bandonly,axis=axes[2]))
end
cortical(scales::FreqScaleFilter,rates::TimeRateFilter) =
  ScaleRateFilter(scales,rates)

function DSP.filt(cort::ScaleRateFilter,cr::MetaAxisArray,progresbar=true)
  filt(cort.rates,filt(cort.scales,cr))
end

# inverse of scale-rate filters
struct ScaleRateFilterInv <: CorticalFilterInv
  scales::FreqScaleFilterInv
  rates::TimeRateFilterInv
  norm::Float64
end
Base.ndims(x::ScaleRateFilterInv) = 2
Base.inv(cf::ScaleRateFilter;norm=0.9) =
  ScaleRateFilterInv(inv(cf.scales),inv(cf.rates),norm)
AxisArrays.axisnames(x::ScaleRateFilterInv) =
  (axisname(x.scales),axisname(x.rates))

function list_filters(z_cum,cr,cf::ScaleRateFilterInv)
  ((CartesianIndex(i,j), HR.*HS)
   for (i,HR) in list_filters(z_cum,cr,cf.rates)
   for (j,HS) in list_filters(z_cum,cr,cf.scales))
end

################################################################################
# private helper functions

function find_fft_dims(y)
  @assert axisdim(y,Axis{:freq}) == ndims(y)
  @assert axisdim(y,Axis{:time}) == 1
  find_fft_dims(size(y))
end
find_fft_dims(y::NTuple{N,Int}) where {N} =
  (nextprod([2,3,5],y[1]),y[2:end-1]...,nextprod([2,3,5],y[end]))

struct FIRFiltering{T,N,P}
  Y::Array{T,N}
  plan::P
end

function FIRFiltering(y,axis)
  dims = map(AxisArrays.axes(y)) do ax
    if AxisArrays.axes(y,axis) == ax
      2nextprod([2,3,5],length(ax))
    else
      length(ax)
    end
  end

  along = axisdim(y,axis)
  Y = fft(pad(y,dims),along)
  FIRFiltering(Y,plan_ifft(Y,along))
end
apply(fir::FIRFiltering,H) = fir.plan * (fir.Y .* H)
Base.size(x::FIRFiltering,i...) = size(x.Y,i...)
Base.ndims(x::FIRFiltering) = ndims(x.Y)

function initfilter(y,ratefilter::TimeRateFilter)
  bandonly = ratefilter.bandonly
  rates = sort(ratefilter.data)
  r = Axis{axisname(ratefilter)}(rates)
  ax = AxisArrays.axes(y)
  newax = ax[1],r,ax[2:end]...

  axar = AxisArray(zeros(complex(eltype(y)),length.(newax)...),newax...)
  arates = sort!(unique!(abs.(rates)))
  rate_axis = bandonly ?
     RateAxis(-Inf*Hz,Inf*Hz) : RateAxis(first(arates),last(arates))
  axis_meta = addaxes(getmeta(y);Dict(axisname(ratefilter) => rate_axis)...)
  MetaAxisArray(axis_meta,axar)
end

function initfilter(y,scalefilter::FreqScaleFilter)
  bandonly = scalefilter.bandonly
  scales = sort(scalefilter.data)
  s = Axis{axisname(scalefilter)}(scales)
  ax = AxisArrays.axes(y)
  newax = ax[1:end-1]...,s,ax[end]

  axar = AxisArray(zeros(complex(eltype(y)),length.(newax)...),newax...)
  scale_axis = bandonly ?
    ScaleAxis(-Inf*cycoct,Inf*cycoct) :
    ScaleAxis(first(scales),last(scales))
  axis_meta = addaxes(getmeta(y);Dict(axisname(scalefilter) => scale_axis)...)
  MetaAxisArray(axis_meta,axar)
end

# TODO: do this for rates as well
reshape_for(v::Array{T,3},cr::AxisArray{T,3}) where T = v
reshape_for(v::Array{T,4},cr::AxisArray{T,4}) where T = v
reshape_for(v::Array{T,3},cr::AxisArray{T,4}) where T =
    reshape(v,ntimes(cr),1,nfrequencies(cr))

# keeps track of cumulative sum of FIR filters
# in frequency-space so we can readily normalize the result.
struct FFTCum{T,P,N,M}
  z::Array{Complex{T},N}
  z_cum::Array{Complex{T},M}
  h_cum::Array{T,M}
  nfrequencies::Int
  ntimes::Int
  plan::P
end

# TODO: working on generalizing FFTCum to working
# for any number of scale and rate axes
function withoutdim(dims,without)
  dims[1:without-1]...,dims[without+1:end]...
end

function FFTCum(cr::MetaAxisArray,withoutaxes)
  dims = find_fft_dims((size(cr,1),size(cr,ndims(cr))))
  mult = fill(ndims(cr),1)
  mult[1] += any(ax -> startswith(string(ax),"scale"),withoutaxes)
  mult[end] += any(ax -> startswith(string(ax),"rate"),withoutaxes)
  z = zeros(eltype(cr),(dims .* mult)...)

  cumsize = (size(z,1),size(z,2))
  z_cum = zeros(eltype(z),cumsize)
  h_cum = zeros(real(eltype(z)),cumsize)
  plan = plan_fft(z)
  FFTCum(z,z_cum,h_cum,nfrequencies(cr),ntimes(cr),plan)
end

Base.size(x::FFTCum,i...) = size(x.z_cum,i...)
Base.ndims(x::FFTCum) = ndims(x.z_cum)

function addfft!(x::FFTCum,cr,h)
  @assert x.ntimes == size(cr,1)
  @assert x.nfrequencies == size(cr,2)

  x.z[1:x.ntimes,1:x.nfrequencies] = cr
  Z = x.plan * x.z
  x.h_cum .+= abs2.(h)
  x.z_cum .+= h .* Z

  x
end

function normalize!(x::FFTCum,cr,norm)
  x.h_cum[:,1] .*= 2
  old_sum = sum(x.h_cum[:,nfrequencies(cr)])
  x.h_cum .= norm.*x.h_cum .+ (1 .- norm).*maximum(x.h_cum)
  x.h_cum .*= old_sum ./ sum(view(x.h_cum,:,nfrequencies(cr)))
  x.z_cum ./= x.h_cum

  spectc = view((x.plan \ x.z_cum),1:ntimes(cr),1:nfrequencies(cr))
  max.(real.(2 .* spectc),0)
end

pad(x,lens) = pad(x,lens...)
function pad(x,lens::T...) where T <: Number
  @assert all(size(x) .<= lens)
  y = zeros(eltype(x),lens)
  y[Base.axes(x)...] = x
  y
end

# transforms a bandpass frequency response into either a high or low pass
# response (or leaves it untouched)
function askind(H,len,maxi,kind,nonorm)
  if kind == :band
    H
  else
    old_sum = sum(H)
    if kind == :low
      H[1:maxi-1] .= 1
    elseif kind == :high
      H[maxi+1:len] .= 1
    else
      error("Unexpected filter kind '$kind'.")
    end
    if !nonorm
      H .= H ./ sum(H) .* old_sum
    end

    H
  end
end

function scale_filters(Y,x,scaleax)
  N_f = size(Y,ndims(Y)) >> 1
  scaleparam = getproperty(x,scaleax)
  map(scales(x)) do scale
	  scale_filter(ustrip(uconvert(cycoct,scale)), N_f, spect_rate,
                 scale <= scaleparam.low ? :low :
                 scale < scaleparam.high ? :band : :high)
  end
end

# create the frequency-scale filter (filter along spectral axis)
function scale_filter(scale,len,ts,kind)
  f2 = ((0:len-1)./len.*ts ./ 2 ./ abs(scale)).^2
  H = f2 .* exp.(1 .- f2)

  askind(H,len,argmax(H),kind,false)
end

function rate_filters(Y,x,rateax;use_conj=false)
  N_t = size(Y,1) >> 1
  rateparam = getproperty(x,rateax)

  map(rates(x)) do rate
    rate_filter(ustrip(uconvert(Hz,rate)), N_t, x.time.Δ,
                abs(rate) <= rateparam.low ? :low :
                abs(rate) < rateparam.high ? :band : :high,use_conj)
  end
end

# create the temporal-rate filter (filter along temporal axis)
function rate_filter(rate,len,Δt,kind,use_conj=false,return_partial=false)
  t = (0:len-1)*ustrip(uconvert(s,Δt))*abs(rate)
  h = @. sin(2π*t) * t^2 * exp(-3.5t)
  h .-= mean(h)

  H0 = view(fft(pad(h,2len)),1:len)
  A = angle.(H0)
  H = abs.(H0)

  maxH,maxi = findmax(H)
  H ./= maxH
  HR = askind(H,len,maxi,kind,true) .* exp.(A*im)

  if use_conj
    HR = conj.(HR)
  end

  if rate >= 0
    HR = pad(HR,2length(HR))
	else
    HR = pad(HR,2length(HR))
		HR[2:end] .= conj.(reverse(HR[2:end]))
		HR[len+1] = abs(HR[len+2])
	end

  if return_partial
    HR,h
  else
    HR
  end
end


text:
function complement(f::CSGNode, coords)
    ev = evaluate(f, coords)
    ev.signint = -1*ev.signint
    ev.value = -1*ev.value
    return ev
end

function intersection(f::CSGNode, g::CSGNode, coords)
    return max(evaluate(f, coords), evaluate(g, coords))
end

function subtraction(f::CSGNode, g::CSGNode, coords)
    return max(evaluate(f, coords), complement(g, coords))
end

function union(f::CSGNode, g::CSGNode, coords)
    return min(evaluate(f, coords), evaluate(g, coords))
end

const CSGOperations = [complement, intersection, subtraction, union]
const CSGopsyms = [:complement, :intersection, :subtraction, :union]
const opDict = Dict(:complement=>complement,
                :intersection=>intersection,
                :subtraction=>subtraction,
                :union=>union)

function complement(f::CachedCSGNode, coords, ind)
    return flipsign(evaluate(f, coords, ind))
end

function intersection(f::CachedCSGNode, g::CachedCSGNode, coords, ind)
    return max(evaluate(f, coords, ind), evaluate(g, coords, ind))
end

function subtraction(f::CachedCSGNode, g::CachedCSGNode, coords, ind)
    return max(evaluate(f, coords, ind), complement(g, coords, ind))
end

function union(f::CachedCSGNode, g::CachedCSGNode, coords, ind)
    return min(evaluate(f, coords, ind), evaluate(g, coords, ind))
end


text:
using Unicode

abstract type CallInfo; end

# Call could be resolved to a singular MI
struct MICallInfo <: CallInfo
    mi::MethodInstance
    rt
end
get_mi(ci::MICallInfo) = ci.mi

# Failed
struct FailedCallInfo <: CallInfo
    sig
    rt
end

function get_mi(ci::FailedCallInfo)
    @error "MethodInstance extraction failed" ci.sig ci.rt
    return nothing
end

# Generated
struct GeneratedCallInfo <: CallInfo
    sig
    rt
end
function get_mi(genci::GeneratedCallInfo)
    @error "Can't extract MethodInstance from call to generated functions" genci.sig genci.rt
    return nothing
end

struct MultiCallInfo <: CallInfo
    sig
    rt
    callinfos::Vector{CallInfo}
end
# actual code-error
get_mi(ci::MultiCallInfo) = error("Can't extract MethodInstance from multiple call informations")

struct DeoptimizedCallInfo <: CallInfo
    accurate::CallInfo
    deoptimized::CallInfo
end
get_mi(ci::DeoptimizedCallInfo) = get_mi(ci.accurate)

struct TaskCallInfo <: CallInfo
    ci::CallInfo
end
get_mi(tci::TaskCallInfo) = get_mi(tci.ci)

# Special handling for ReturnTypeCall
struct ReturnTypeCallInfo <: CallInfo
    called_mi::CallInfo
end
get_mi(rtci::ReturnTypeCallInfo) = get_mi(rtci.called_mi)

# CUDA callsite
struct CuCallInfo <: CallInfo
     cumi::MICallInfo
end
get_mi(gci::CuCallInfo) = get_mi(gci.cumi)

struct Callsite
    id::Int # ssa-id
    info::CallInfo
end
get_mi(c::Callsite) = get_mi(c.info)

# Callsite printing
mutable struct TextWidthLimiter
    io::IO
    width::Int
    limit::Int
end
TextWidthLimiter(io::IO, limit) = TextWidthLimiter(io, 0, limit)
has_space(limiter::TextWidthLimiter, width::Int) = limiter.width + width < limiter.limit - 1
has_space(limiter::TextWidthLimiter, s) = has_space(limiter, textwidth(string(s)))
function Base.print(io::TextWidthLimiter, s::String)
    io.width == io.limit && return 0
    width = textwidth(s::String)
    if has_space(io, width)
        print(io.io, s)
        io.width += width
        return
    else
        for c in graphemes(s)
            cwidth = textwidth(c)
            if has_space(io, cwidth)
                print(io, c)
                io.width += cwidth
            else
                break
            end
        end
        print(io, '…')
        io.width += 1
    end
end

function headstring(@nospecialize(T))
    if T isa Union || T === Union{}
        return string(T)
    elseif T isa UnionAll
        return headstring(Base.unwrap_unionall(T))
    else
        return string(T.name)
    end
end

function __show_limited(limiter, name, tt, rt)
    if !has_space(limiter, name)
        print(limiter, '…')
        return
    end
    print(limiter, string(name))
    pstrings = map(string, tt)
    headstrings = map(headstring, tt)
    print(limiter, "(")
    if length(pstrings) != 0
        # See if we have space to print all the parameters fully
        if has_space(limiter, sum(textwidth, pstrings) + 3*length(pstrings))
            print(limiter, join(map(T->string("::", T), pstrings), ","))
        # Alright, see if we at least have enough space for each head
        elseif has_space(limiter, sum(textwidth, headstrings) + 6*length(pstrings))
            print(limiter, join(map(T->string("::", T, "{…}"), headstrings), ","))
        # Fine, what about just indicating the number of arguments
        elseif has_space(limiter, 2*(length(tt)))
            print(limiter, join(map(T->"…", pstrings), ","))
        else
            print(limiter, "…")
        end
    end
    print(limiter, ")")

    # If we have space for the return type, print it
    rts = string(rt)
    if has_space(limiter, textwidth(rts)+2)
        print(limiter, string("::", rts))
    end
end

function show_callinfo(limiter, mici::MICallInfo)
    mi = mici.mi
    tt = (Base.unwrap_unionall(mi.specTypes)::DataType).parameters[2:end]
    name = mi.def.name
    rt = mici.rt
    __show_limited(limiter, name, tt, rt)
end

function show_callinfo(limiter, ci::Union{MultiCallInfo, FailedCallInfo, GeneratedCallInfo})
    types = (ci.sig::DataType).parameters
    f = types[1]
    if f isa Union
        name = string(f)
    else
        name = nameof(f)
    end
    tt = types[2:end]
    rt = ci.rt
    __show_limited(limiter, name, tt, rt)
end

function Base.show(io::IO, c::Callsite)
    limit = get(io, :limit, false)
    cols = limit ? displaysize(io)[2] : typemax(Int)
    limiter = TextWidthLimiter(io, cols)
    print(limiter, string("%", c.id, " "))
    if isa(c.info, MICallInfo)
        print(limiter, " = invoke ")
        show_callinfo(limiter, c.info)
    elseif c.info isa MultiCallInfo
        print(limiter, " = call ")
        show_callinfo(limiter, c.info)
    elseif isa(c.info, DeoptimizedCallInfo)
        deoptstr = sprint(context = IOContext(io, :color => true)) do tmpio
            print(tmpio, " = ")
            printstyled(tmpio, "deoptimized"; color = :red)
            print(tmpio, " ")
        end
        print(limiter, deoptstr)
        show_callinfo(limiter, c.info.accurate)
    elseif c.info isa FailedCallInfo ||
           c.info isa GeneratedCallInfo
        print(limiter, " = call ")
        show_callinfo(limiter, c.info)
    elseif c.info isa TaskCallInfo
        print(limiter, " = task < ")
        show_callinfo(limiter, c.info.ci)
        print(limiter, " >")
    elseif isa(c.info, ReturnTypeCallInfo)
        print(limiter, " = return_type < ")
        show_callinfo(limiter, c.info.called_mi)
        print(limiter, " >")
    elseif isa(c.info, CuCallInfo)
        print(limiter, " = cucall < ")
        show_callinfo(limiter, c.info.cumi)
        print(limiter, " >")
    end
end

is_callsite(cs::Callsite, mi::MethodInstance) = is_callsite(cs.info, mi)
is_callsite(info::MICallInfo, mi::MethodInstance) = info.mi == mi
is_callsite(info::DeoptimizedCallInfo, mi::MethodInstance) = is_callsite(info.accurate, mi)
is_callsite(info::TaskCallInfo, mi::MethodInstance) = is_callsite(info.ci, mi)
is_callsite(info::ReturnTypeCallInfo, mi::MethodInstance) = is_callsite(info.called_mi, mi)
is_callsite(info::CuCallInfo, mi::MethodInstance) = is_callsite(info.cumi, mi)
function is_callsite(info::MultiCallInfo, mi::MethodInstance)
    for csi in info.callinfos
        is_callsite(csi, mi) && return true
    end
    return false
end
is_callsite(::CallInfo, mi::MethodInstance) = false


text:
##################################################################################
# This file is part of ModelBaseEcon.jl
# BSD 3-Clause License
# Copyright (c) 2020-2022, Bank of Canada
# All rights reserved.
##################################################################################

abstract type AbstractEquation end

# equations must have these fields: expr, vinds, vsyms, eval_resid, eval_RJ
for fn in (:expr, :vinds, :vsyms, :eval_resid, :eval_RJ)
    local qnfn = QuoteNode(fn)
    eval(quote
        @inline $fn(eqn::AbstractEquation) = getfield(eqn, $qnfn)
    end)
end

# equations might have these fields. If not, we provide defaults
@inline flags(eqn::AbstractEquation) = hasfield(typeof(eqn), :flags) ? getfield(eqn, :flags) : nothing
@inline flag(eqn::AbstractEquation, f::Symbol) = (flgs = flags(eqn); hasfield(typeof(flgs), f) ? getfield(flgs, f) : false)
@inline doc(eqn::AbstractEquation) = :doc in fieldnames(typeof(eqn)) ? getfield(eqn, :doc) : ""

# @inline eval_resid(eqn::AbstractEquation, x) = eval_resid(eqn)(x)
# @inline eval_RJ(eqn::AbstractEquation, x) = eval_RJ(eqn)(x)

#
function Base.show(io::IO, eqn::AbstractEquation)
    flagstr = ""
    flgs = flags(eqn)
    for f in fieldnames(typeof(flgs))
        if getfield(flgs, f)
            flagstr *= "@$(f) "
        end
    end
    docstr = ""
    if !isempty(doc(eqn)) && !get(io, :compact, false)
        docstr = "\"$(doc(eqn))\" "
    end
    print(io, docstr, flagstr, expr(eqn))
end

Base.:(==)(e1::AbstractEquation, e2::AbstractEquation) = flags(e1) == flags(e2) && expr(e1) == expr(e2)
Base.hash(e::AbstractEquation, h::UInt) = hash((flags(e), expr(e)), h)

abstract type AbstractModel end

# a subtype of AbstractModel is expected to have a number of fields.
# If it doesn't, the creater of the new model type must define the
# access methods that follow.

@inline variables(m::AM) where {AM<:AbstractModel} = getfield(m, :variables)
@inline nvariables(m::AM) where {AM<:AbstractModel} = length(variables(m))

@inline shocks(m::AM) where {AM<:AbstractModel} = getfield(m, :shocks)
@inline nshocks(m::AM) where {AM<:AbstractModel} = length(shocks(m))

@inline allvars(m::AM) where {AM<:AbstractModel} = vcat(variables(m), shocks(m))
@inline nallvars(m::AM) where {AM<:AbstractModel} = length(variables(m)) + length(shocks(m))

@inline sstate(m::AM) where {AM<:AbstractModel} = getfield(m, :sstate)

@inline parameters(m::AM) where {AM<:AbstractModel} = getfield(m, :parameters)

@inline equations(m::AM) where {AM<:AbstractModel} = getfield(m, :equations)
@inline nequations(m::AM) where {AM<:AbstractModel} = length(equations(m))

@inline alleqns(m::AM) where {AM<:AbstractModel} = getfield(m, :equations)
@inline nalleqns(m::AM) where {AM<:AbstractModel} = length(equations(m))

export parameters
export variables, nvariables
export shocks, nshocks
export equations, nequations
export sstate
#######

# @inline moduleof(f::Function) = parentmodule(f)
@inline moduleof(e::AbstractEquation) = parentmodule(eval_resid(e))
function moduleof(m::AbstractModel)
    eqns = equations(m)
    if isempty(eqns)
        error("Unable to determine the module containing the given model. Try adding equations to it and call `@initialize`.")
    end
    return moduleof(first(eqns))
end


text:
"""
    drifters_ElipotEtAl16(pth,lst;chnk=Inf,rng=(-Inf,Inf))

Read near-surface [drifter data](https://doi.org/10.1002/2016JC011716) from the
[Global Drifter Program](https://doi.org/10.25921/7ntx-z961) into a DataFrame.

```
pth,lst=drifters_ElipotEtAl16()
df=drifters_ElipotEtAl16( pth*lst[end], rng=(2014.1,2014.2) )
```
"""
function drifters_ElipotEtAl16(fil::String;chnk=1000,rng=(-Inf,Inf))
    t=ncread(fil,"TIME")
    t_u=ncgetatt(fil,"TIME","units")
    lo=ncread(fil,"LON")
    la=ncread(fil,"LAT")

    ##

    ii=findall(isfinite.(lo.*la.*t))

    t=t[ii]
    lo=lo[ii]
    la=la[ii]
    ID=ncread(fil,"ID")[ii]
    DROGUE=ncread(fil,"DROGUE")[ii]
    U=ncread(fil,"U")[ii]
    V=ncread(fil,"V")[ii]

    ##

    t=timedecode(t, t_u)
    tmp=dayofyear.(t)+(hour.(t) + minute.(t)/60 ) /24
    t=year.(t)+tmp./daysinyear.(t)

    ii=findall( (t.>rng[1]).&(t.<=rng[2]).*(DROGUE.==1) )

    t=t[ii]
    lo=lo[ii]
    la=la[ii]
    ID=ID[ii]
    DROGUE=DROGUE[ii]
    U=U[ii]
    V=V[ii]

    ##

    df = DataFrame(ID=Int[], lon=Float64[], lat=Float64[],
                u=Float64[], v=Float64[], t=Float64[])
    !isinf(chnk) ? nn=Int(ceil(length(ii)/chnk)) : nn=1
    for jj=1:nn
       #println([jj nn])
       !isinf(chnk) ? i=(jj-1)*chnk.+(1:chnk) : i=(1:length(ii))
       i=i[findall(i.<length(ii))]
       append!(df,DataFrame(lon=lo[i], lat=la[i], t=t[i], u=U[i], v=V[i], ID=Int.(ID[i])))
    end

    return df
end

"""
    drifters_ElipotEtAl16()

Path name and file list for near-surface [drifter data](https://doi.org/10.1002/2016JC011716)
from the [Global Drifter Program](https://doi.org/10.25921/7ntx-z961)

```
pth,lst=drifters_ElipotEtAl16()
```
"""
function drifters_ElipotEtAl16()
    pth="Drifter_hourly_v013/"
    lst=["driftertrajGPS_1.03.nc","driftertrajWMLE_1.02_block1.nc","driftertrajWMLE_1.02_block2.nc",
       "driftertrajWMLE_1.02_block3.nc","driftertrajWMLE_1.02_block4.nc","driftertrajWMLE_1.02_block5.nc",
       "driftertrajWMLE_1.02_block6.nc","driftertrajWMLE_1.03_block7.nc"]
    return pth,lst
 end

"""
    drifters_ElipotEtAl16(t0::Number,t1::Number)

Loop over all files and call drifters_ElipotEtAl16 with rng=(t0,t1)

```
@everywhere using OceanRobots, CSV
@distributed for y in 2005:2020
    df=drifters_ElipotEtAl16(y+0.0,y+1.0)
    fil="Drifter_hourly_v013/csv/drifters_"*string(y)*".csv"
    CSV.write(fil, df)
end
```
"""
function drifters_ElipotEtAl16( t0::Number,t1::Number )
    pth,lst=drifters_ElipotEtAl16()
    df = DataFrame([fill(Int, 1) ; fill(Float64, 5)], [:ID, :lon, :lat, :u, :v, :t])
    for fil in lst
       println(fil)
       append!(df,drifters_ElipotEtAl16( pth*fil,chnk=10000,rng=(t0,t1) ))
    end
    return df
 end


text:
@testset "reifiedinset.jl" begin
    @testset "ReifiedInSet()" begin
        trailer = SeaPearl.Trailer()
        a = SeaPearl.IntSetVar(2, 5, "a", trailer)
        x = SeaPearl.IntVar(0, 3, "x", trailer)
        b = SeaPearl.BoolVar("b", trailer)

        constraint = SeaPearl.ReifiedInSet(x, a, b, trailer)

        @test constraint in a.onDomainChange
        @test constraint in x.onDomainChange
        @test constraint in b.onDomainChange
        @test constraint.active.value
    end

    @testset "propagate!(::ReifiedInSet)" begin
        trailer = SeaPearl.Trailer()
        a = SeaPearl.IntSetVar(2, 5, "a", trailer)
        x = SeaPearl.IntVar(0, 3, "x", trailer)
        b = SeaPearl.BoolVar("b", trailer)
        constraint = SeaPearl.ReifiedInSet(x, a, b, trailer)

        toPropagate = Set{SeaPearl.Constraint}()
        prunedDomains = SeaPearl.CPModification()

        @test SeaPearl.propagate!(constraint, toPropagate, prunedDomains)

        @test length(x.domain) == 4
        @test 0 in x.domain
        @test constraint.active.value
        @test prunedDomains == SeaPearl.CPModification()
        @test constraint.active.value

        ### Filtering b
        # b -> true
        SeaPearl.assign!(x, 3)
        SeaPearl.require!(a.domain, 3)

        @test SeaPearl.propagate!(constraint, toPropagate, prunedDomains)

        @test SeaPearl.isbound(b)
        @test SeaPearl.assignedValue(b)
        @test !constraint.active.value

        # b -> false
        trailer = SeaPearl.Trailer()
        a = SeaPearl.IntSetVar(2, 5, "a", trailer)
        x = SeaPearl.IntVar(0, 1, "x", trailer)
        b = SeaPearl.BoolVar("b", trailer)
        constraint = SeaPearl.ReifiedInSet(x, a, b, trailer)

        toPropagate = Set{SeaPearl.Constraint}()
        prunedDomains = SeaPearl.CPModification()
        @test SeaPearl.propagate!(constraint, toPropagate, prunedDomains)

        @test SeaPearl.isbound(b)
        @test !SeaPearl.assignedValue(b)
        @test !constraint.active.value

        ### Filtering a or x
        # b == true
        trailer = SeaPearl.Trailer()
        a = SeaPearl.IntSetVar(2, 5, "a", trailer)
        x = SeaPearl.IntVar(0, 3, "x", trailer)
        b = SeaPearl.BoolVar("b", trailer)
        constraint = SeaPearl.ReifiedInSet(x, a, b, trailer)

        SeaPearl.assign!(b, true)

        toPropagate = Set{SeaPearl.Constraint}()
        prunedDomains = SeaPearl.CPModification()

        @test SeaPearl.propagate!(constraint, toPropagate, prunedDomains)

        @test length(x.domain) == 2
        @test !(0 in x.domain)
        @test !(1 in x.domain)
        @test 2 in x.domain
        @test constraint.active.value
        @test prunedDomains == SeaPearl.CPModification("x" => [0, 1])

        SeaPearl.assign!(x, 3)

        @test SeaPearl.propagate!(constraint, toPropagate, prunedDomains)

        @test SeaPearl.required_values(a.domain) == Set{Int}([3])
        @test SeaPearl.possible_not_required_values(a.domain) == Set{Int}([2, 4, 5])
        @test !constraint.active.value

        # b == false
        trailer = SeaPearl.Trailer()
        a = SeaPearl.IntSetVar(2, 5, "a", trailer)
        x = SeaPearl.IntVar(0, 4, "x", trailer)
        b = SeaPearl.BoolVar("b", trailer)
        constraint = SeaPearl.ReifiedInSet(x, a, b, trailer)

        SeaPearl.assign!(b, false)
        SeaPearl.require!(a.domain, 2)
        SeaPearl.require!(a.domain, 3)

        @test SeaPearl.propagate!(constraint, toPropagate, prunedDomains)
        @test length(x.domain) == 3
        @test 0 in x.domain
        @test 1 in x.domain
        @test 4 in x.domain
        @test !(2 in x.domain)
        @test constraint.active.value

        SeaPearl.require!(a.domain, 4)
        @test SeaPearl.propagate!(constraint, toPropagate, prunedDomains)
        @test !constraint.active.value

        ### Infeasible
        trailer = SeaPearl.Trailer()
        a = SeaPearl.IntSetVar(2, 5, "a", trailer)
        x = SeaPearl.IntVar(2, 3, "x", trailer)
        b = SeaPearl.BoolVar("b", trailer)
        constraint = SeaPearl.ReifiedInSet(x, a, b, trailer)

        SeaPearl.assign!(b, false)
        SeaPearl.require!(a.domain, 2)
        SeaPearl.require!(a.domain, 3)

        @test !SeaPearl.propagate!(constraint, toPropagate, prunedDomains)

        trailer = SeaPearl.Trailer()
        a = SeaPearl.IntSetVar(2, 5, "a", trailer)
        x = SeaPearl.IntVar(0, 1, "x", trailer)
        b = SeaPearl.BoolVar("b", trailer)
        constraint = SeaPearl.ReifiedInSet(x, a, b, trailer)

        SeaPearl.assign!(b, true)

        @test !SeaPearl.propagate!(constraint, toPropagate, prunedDomains)
    end
end


text:
import Base: +, -, abs, abs2, acos, acosh, asin, asinh, atan, atanh, ceil, cos, cosh, exp, exp10, exp2, expm1, floor, log, log10, log1p, log2, one, round, sign, sin, sinh, sqrt, tan, tanh, trunc, zero
import Base.Math: cbrt, cospi, sinpi
import SpecialFunctions: besselj, besselj0, besselj1, bessely, bessely0, bessely1, digamma, erf, erfc, erfcinv, erfcx, erfinv, gamma, lgamma, trigamma
import Base.Broadcast: broadcasted
using Knet.LibKnet8: @knet8, unary_ops, unary_ops_with_int_degree
# include("broadcast.jl") ## Bcasted

# unary.jl: Unary Array->Array operations.

function unary_op(f, j=f, o...)
    J=Symbol(j)
    M = which(@__MODULE__, J)
    for S in (32,64)
        T = Symbol("Float$S")
        F = "$(f)_$S"
        @eval begin
            function broadcasted(::typeof($J),x::KnetArray{$T})
                y = similar(x)
                @knet8($F,(Cint,Ptr{$T},Ptr{$T}),length(y),x,y)
                return y
            end
            # Bcasted methods
            ($M).$J(x::Bcasted{<:KnetArray{$T}}) = broadcasted($J, x.value) |> Bcasted
            broadcasted(::typeof($J),x::Bcasted{<:KnetArray{$T}}) = broadcasted($J, x.value) |> Bcasted
        end
    end
    @eval begin # so we do not trigger some default Base implementation
        ($M).$J(x::Bcasted) = throw(MethodError($J,(x,)))
        broadcasted(::typeof($J),x::Bcasted) = throw(MethodError($J,(x,)))
    end
end

for f in unary_ops
    if !isa(f,Tuple); f=(f,); end
    j = length(f) >= 2 ? f[2] : f[1]
    ## Process non-Base functions in the appropriate submodule:
    if !isdefined(@__MODULE__, Symbol(j))
        # @warn "$j not defined in $(@__MODULE__), skipping."
        # invx, elu, gelu, relu, selu, sigm
        continue
    end
    unary_op(f...)
end

function unary_op_with_int_degree(f, j=f, o...)
    J=Symbol(j)
    M = which(@__MODULE__, J)
    for S in (32,64)
        T = Symbol("Float$S")
        F = "$(f)_$S"
        @eval begin
            function broadcasted(::typeof($J),d::Integer,x::KnetArray{$T})
                y = similar(x)
                @knet8($F,(Cint,Cint,Ptr{$T},Ptr{$T}),length(y),d,x,y)
                return y
            end
            # Bcasted methods
            ($M).$J(d,x::Bcasted{<:KnetArray{$T}}) = broadcasted($J, d, x.value) |> Bcasted
            broadcasted(::typeof($J),d,x::Bcasted{<:KnetArray{$T}}) = broadcasted($J, d, x.value) |> Bcasted
        end
    end
    @eval begin # so we do not trigger some default Base implementation
        ($M).$J(d,x::Bcasted) = throw(MethodError($J,(x,)))
        broadcasted(::typeof($J),d,x::Bcasted) = throw(MethodError($J,(x,)))
    end
end

for f in unary_ops_with_int_degree
    if !isa(f,Tuple); f=(f,); end
    unary_op_with_int_degree(f...)
end

# Unary plus and minus

broadcasted(::typeof(+), a::KnetArray)=a
+(a::KnetArray)=a
-(a::KnetArray)=broadcasted(-,a)

# Identity: Issue #335
broadcasted(::typeof(identity), a::KnetArray)=a

# TODO: SpecialFunctions loggamma update #486.
# SpecialFunctions 0.8: lgamma(x::Real)` is deprecated, use `(logabsgamma(x))[1]` instead.
# Other alternative is loggamma, throws a DomainError if gamma(x) is negative.
# However neither are equivalent to cuda lgamma.
# commenting this out until we implement a better solution.
# if !isdefined(SpecialFunctions, :loggamma) && isdefined(SpecialFunctions, :lgamma)
#     loggamma(x) = lgamma(x)
# end
# if isdefined(SpecialFunctions, :loggamma) && !isdefined(SpecialFunctions, :lgamma)
#     lgamma(x) = loggamma(x)
# end
# if isempty(methods(loggamma,(AutoGrad.Value,))) # Should be fixed after AutoGrad 1.1.4
#     @primitive loggamma(x),dy,y (dy.*(digamma.(x)))
# end

# Unbroadcasted zero works on arrays: this moved to karray.jl
# import Base: zero
# zero(x::KnetArray)=zero.(x)


text:
@recipe function f(solution::MovingWindowSolution; vars = nothing, label = "")
    fontfamily --> "Computer Modern"
    framestyle --> :box
    gridalpha --> 0.2
    linewidth --> 1.5
    minorgrid --> 0.1
    minorgridstyle --> :dash
    seriestype --> :path
    tick_direction --> :out
    xwiden --> false
    M = length(solution)
    for m in 1:M
        vars --> vars
        @series begin
            seriescolor --> m
            solution[m]
        end
        @series begin
            seriescolor --> m
            seriestype := :vline
            linestyle --> :dash
            [solution[m][end][end].t[end]]
        end
    end
end


text:
import Pkg; Pkg.activate("."); Pkg.instantiate()
using Documenter

#include("../src/Pitches.jl")
push!(LOAD_PATH, "../src/")
using Pitches

makedocs(
    sitename="Pitches.jl",
    pages = [
        "index.md",
        "tutorial.md",
        "interface.md",
        "Pitch and Interval Types" => [
            "spelled.md",
            "midi.md",
            "logfreq.md",
        ]
    ]
)


text:
using BinaryProvider # requires BinaryProvider 0.3.0 or later

# Parse some basic command-line arguments
const verbose = "--verbose" in ARGS
const prefix = Prefix(get([a for a in ARGS if a != "--verbose"], 1, joinpath(@__DIR__, "usr")))
products = [
    LibraryProduct(joinpath(@__DIR__, "usr"), ["liblsl"], :liblsl),
]

# Download binaries from hosted location
bin_prefix = "https://github.com/samuelpowell/LSLBuilder/releases/download/v1.13.0-b13"

# Listing of files generated by BinaryBuilder:
download_info = Dict(
    Linux(:aarch64, libc=:glibc) => ("$bin_prefix/liblsl.v1.13.0-b13.aarch64-linux-gnu.tar.gz", "48c55cb92e300cfa26c08fe507fc2cc20159f04513a542678eb794040c49a452"),
    Linux(:aarch64, libc=:musl) => ("$bin_prefix/liblsl.v1.13.0-b13.aarch64-linux-musl.tar.gz", "172955b82e88ad0fb9862fa30382d7f3811896fad31f1850f99f210445373ddf"),
    Linux(:armv7l, libc=:glibc, call_abi=:eabihf) => ("$bin_prefix/liblsl.v1.13.0-b13.arm-linux-gnueabihf.tar.gz", "3f63e740da7ac0977daa8469ee621c69853da5657992781d78c10e1fa63daf54"),
    Linux(:armv7l, libc=:musl, call_abi=:eabihf) => ("$bin_prefix/liblsl.v1.13.0-b13.arm-linux-musleabihf.tar.gz", "33358a73b71c34f0c8fe21b5c701169de1bbf4305357efb97270c0eea4f12d10"),
    Linux(:i686, libc=:glibc) => ("$bin_prefix/liblsl.v1.13.0-b13.i686-linux-gnu.tar.gz", "5019a298d6f3fd1f36eb9a200a721ad29cc83b9697df707061e820dcba224460"),
    Linux(:i686, libc=:musl) => ("$bin_prefix/liblsl.v1.13.0-b13.i686-linux-musl.tar.gz", "4960c3a3ec498fc3e1dc855687f78dcb7330284b0db3c2329b9b921e1360c369"),
    Windows(:i686) => ("$bin_prefix/liblsl.v1.13.0-b13.i686-w64-mingw32.tar.gz", "b2323a2fe16b2f061a4330573e41394f7cf1d09d54c0ce4f689f93f32085f41e"),
    Linux(:powerpc64le, libc=:glibc) => ("$bin_prefix/liblsl.v1.13.0-b13.powerpc64le-linux-gnu.tar.gz", "279f3460c90b01d0db68058af7c69e59c623b57136da14229dfcc237be443b92"),
    MacOS(:x86_64) => ("$bin_prefix/liblsl.v1.13.0-b13.x86_64-apple-darwin14.tar.gz", "ef3f353ee80856d47178b5a9aff579cae66bcb81381300cda3d59b79c1e4e3d8"),
    Linux(:x86_64, libc=:glibc) => ("$bin_prefix/liblsl.v1.13.0-b13.x86_64-linux-gnu.tar.gz", "65e5d90fdc40024f404f553fcb49a285cceb27297c6b93482719a016460b6375"),
    Linux(:x86_64, libc=:musl) => ("$bin_prefix/liblsl.v1.13.0-b13.x86_64-linux-musl.tar.gz", "40ccd453480aa09c5dfca5df73efda9f4a23468ddd4702597004c427eb6e491b"),
    FreeBSD(:x86_64) => ("$bin_prefix/liblsl.v1.13.0-b13.x86_64-unknown-freebsd11.1.tar.gz", "2a133452203fc808c2c62413851a116815a98eb7d40fdfc181806a837ed20054"),
    Windows(:x86_64) => ("$bin_prefix/liblsl.v1.13.0-b13.x86_64-w64-mingw32.tar.gz", "eac664995173453adfb2145f80a1923a28be5a7f8123c675629be297d73ad234"),
)

# Install unsatisfied or updated dependencies:
unsatisfied = any(!satisfied(p; verbose=verbose) for p in products)
dl_info = choose_download(download_info, platform_key_abi())
if dl_info === nothing && unsatisfied
    # If we don't have a compatible .tar.gz to download, complain.
    # Alternatively, you could attempt to install from a separate provider,
    # build from source or something even more ambitious here.
    error("Your platform (\"$(Sys.MACHINE)\", parsed as \"$(triplet(platform_key_abi()))\") is not supported by this package!")
end

# If we have a download, and we are unsatisfied (or the version we're
# trying to install is not itself installed) then load it up!
if unsatisfied || !isinstalled(dl_info...; prefix=prefix)
    # Download and install binaries
    install(dl_info...; prefix=prefix, force=true, verbose=verbose)
end

# Cludge official binary library to avoid problems with BinaryBuilder
@static if Sys.iswindows()

    dllprefix = joinpath(@__DIR__, "usr")
    exe7z = joinpath(Sys.BINDIR, "7z.exe")
    dlpath = joinpath(@__DIR__, "usr", "downloads", "liblsl-1.13.0-official.7z")

    if Sys.WORD_SIZE == 64
        dlurl = "https://github.com/sccn/liblsl/releases/download/1.13.0-b13/liblsl-1.13.0-Win64.7z"
        dllfn = "liblsl64.dll"
    else
        dlurl = "https://github.com/sccn/liblsl/releases/download/1.13.0-b13/liblsl-1.13.0-Win32.7z"
        dllfn = "liblsl32.dll"
    end

    download(dlurl, dlpath)

    run(`$exe7z e $dlpath -o$(dllprefix) $dllfn -r -aoa`)
    mv(joinpath(dllprefix, dllfn), joinpath(dllprefix, "liblsl.dll"), force=true)
end

# Write out a deps.jl file that will contain mappings for our products
write_deps_file(joinpath(@__DIR__, "deps.jl"), products, verbose=verbose)


text:
using LightGraphs

export dim, pdim

struct IndexManager
    v_data::OrderedDict{Int, UnitRange{Int}}
    e_data::OrderedDict{Int, UnitRange{Int}}
    v_para::OrderedDict{Int, UnitRange{Int}}
    e_para::OrderedDict{Int, UnitRange{Int}}
    function IndexManager()
        dict() = OrderedDict{Int, UnitRange{Int}}()
        new((dict() for i in 1:4)...)
    end
end

const EXECUTION_STYLES = (:seq, :threaded)

struct Network{EX,NL,VTup}
    "vertex batches of same function"
    vertexbatches::VTup
    "network layer"
    nl::NL
    "index manager"
    im::IndexManager
end

dim(nw::Network) = full_data_range(nw.im)[end]
pdim(nw::Network) = full_para_range(nw.im)[end]

struct NetworkLayer{GT,CTup,AF}
    "graph/toplogy of layer"
    g::GT
    "batches of edges with same color"
    colorbatches::CTup
    "accumulation function"
    accumulator::AF
    "depth/dim of accumulation"
    accdim::Int
    "lazy cache pool"
    cachepool::CachePool
end

accumulator_cache(nl::NetworkLayer, T) = getcache(nl.cachepool, T, nl.accdim * nv(nl.g))

struct ColorBatch{ETup}
    "edge indices (as in LG edge iterator) contained in batch"
    edges::Vector{Int}
    "edge batches of same function"
    edgebatches::ETup
end

struct EdgeBatch{F}
    "edge indices (as in LG edge iterator) contained in batch"
    edges::Vector{Int}
    "edge function"
    fun::F
    "edge dimension"
    dim::Int
    "first index of first batch element in flat state vector (-1 for static edge)"
    firstidx::Int
    "parameter dimensions"
    pdim::Int
    "first index of first batch element in flat parameter Vector"
    pfirstidx::Int
    "src: idx in g, idx in data, dim ; dst: idx in g, idx in data, dim"
    vertex_indices::Vector{Int}
end

component_idxs(eb::EdgeBatch) = eb.edges

@inline function vertex_indices(batch::EdgeBatch, i)
    v = batch.vertex_indices
    idx = (i-1)*6
    @inbounds begin
        (v[idx + 1], v[idx + 2], v[idx + 3], v[idx + 4], v[idx + 5], v[idx + 6])
    end
end

@inline function src_dst_ranges(layer::NetworkLayer, batch::EdgeBatch, i)
    (src, src_dat, src_dim, dst, dst_dat, dst_dim) = vertex_indices(batch, i)
    # ranges of src and data idx in data array
    src_r = src_dat : src_dat + src_dim - 1
    dst_r = dst_dat : dst_dat + dst_dim - 1

    # range of src vertex in accumulator array
    src_acc_first = 1 + (src-1) * layer.accdim
    src_acc = src_acc_first : src_acc_first + layer.accdim - 1

    # range of dst vertex in accumulator array
    dst_acc_first = 1 + (dst-1) * layer.accdim
    dst_acc = dst_acc_first : dst_acc_first + layer.accdim - 1

    # range of edge in flat state array
    edge_first = batch.firstidx + (i-1) * batch.dim
    edge_r = edge_first : edge_first + batch.dim - 1

    return (src_r, dst_r, src_acc, dst_acc, edge_r)
end

struct VertexBatch{F}
    "vertex indices contained in batch"
    vertices::Vector{Int}
    "vertex function"
    fun::F
    "vertex dimension"
    dim::Int
    "first index of first batch element in flat state Vector"
    firstidx::Int
    "parameter dimensions"
    pdim::Int
    "first index of first batch element in flat parameter Vector"
    pfirstidx::Int
end

component_idxs(vb::VertexBatch) = vb.vertices

@inline function vertex_ranges(layer::NetworkLayer, batch::VertexBatch, i)
    # range of the vertex in data array
    vstart = batch.firstidx + (i-1) * batch.dim
    v = vstart : vstart + batch.dim - 1

    # range in accumulator array
    v_idx = batch.vertices[i]
    astart = 1 + (v_idx-1) * layer.accdim
    acc = astart : astart + layer.accdim - 1

    return (v, acc)
end

@inline function parameter_range(batch::Union{EdgeBatch,VertexBatch}, i)
    pstart = batch.pfirstidx + (i-1) * batch.pdim
    p = pstart : pstart + batch.pdim - 1
end

Base.length(cb::ColorBatch) = length(cb.edges)
Base.length(eb::EdgeBatch) = length(eb.edges)
Base.length(vb::VertexBatch) = length(vb.vertices)

LightGraphs.nv(nw::Network) = sum(nv.(nw.vertexbatches))
LightGraphs.nv(vb::VertexBatch) = length(vb.vertices)

LightGraphs.ne(nw::Network) = size(ne.(nl))
LightGraphs.ne(nl::NetworkLayer) = ne(nl.g)

_lastinrange(d::OrderedDict{Int, UnitRange{Int}}) = isempty(d) ? 0 : d[d.keys[end]][end]

_lastindex_data(im::IndexManager) = max(_lastinrange(im.v_data), _lastinrange(im.e_data))
_lastindex_para(im::IndexManager) = max(_lastinrange(im.v_para), _lastinrange(im.e_para))

function _register_comp!(im::IndexManager, data, para, idxs, dim, pdim)
    firstidx = datapos = _lastindex_data(im) + 1
    pfirstidx = parapos = _lastindex_para(im) + 1
    for i in idxs
        if haskey(data, i) || haskey(para, i)
            error("Index $i allready in $data or $para")
        end
        if dim > 0
            data[i] = datapos : datapos + dim - 1
            datapos += dim
        end
        if pdim > 0
            para[i] = parapos : parapos + pdim - 1
            parapos += pdim
        end
    end
    return firstidx, pfirstidx
end

register_edges!(im::IndexManager, idxs, dim, pdim) = _register_comp!(im, im.e_data, im.e_para, idxs, dim, pdim)
register_vertices!(im::IndexManager, idxs, dim, pdim) = _register_comp!(im, im.v_data, im.v_para, idxs, dim, pdim)

vertex_data_range(im::IndexManager, i) = im.v_data[i]
vertex_para_range(im::IndexManager, i) = im.v_para[i]
edge_data_range(im::IndexManager, i) = im.e_data[i]
edge_para_range(im::IndexManager, i) = im.e_para[i]

function _isdense(last, dicts...)
    a = Vector{Int}(undef, last)
    empty!(a)
    for d in dicts
        for v in values(d)
            append!(a, v)
        end
    end
    sort!(a)
    return a == 1:last
end

isdense(im::IndexManager) = (_isdense(_lastindex_data(im), im.v_data, im.e_data) &&
    _isdense(_lastindex_para(im), im.v_para, im.e_para))

full_data_range(im::IndexManager) = 1:_lastindex_data(im)
full_para_range(im::IndexManager) = 1:_lastindex_para(im)


text:
################################################
# Unconstrained
################################################
function fitnessToMass(fitness::Vector{Float64})
    m = minimum(fitness)

    if m < 0
        fitness = 2 * abs(m) .+ fitness
    end

    fitness = 2 * maximum(fitness) .- fitness

    return fitness
end

"""
    getMass(U, max_fitness, epsilon)
    get the mass in a vector `m[i] = f[i] + 2*max_fitness*sum_vio[i] `, tolerance
    in equality constraints are given by `epsilon`
"""
function getMass(U::Array{xf_indiv,1})
    n, d = length(U), length(U[1].x)

    fitness = zeros(Float64, n)

    for i = 1:n
        fitness[i] = U[i].f
    end

    return fitnessToMass(fitness)
end