text:
export genus, representative, rank, det, uniformizer, det_representative,
       gram_matrix, representative, genus, genera_hermitian,
       local_genera_hermitian, rank

################################################################################
#
#  Local genus symbol
#
################################################################################

add_verbose_scope(:Lattice)
add_assert_scope(:Lattice)

# Need to make this type stable once we have settled on a design
mutable struct LocalGenusHerm{S, T}
  E::S                                # Field
  p::T                                # prime of base_field(E)
  data::Vector{Tuple{Int, Int, Int}}  # data
  norm_val::Vector{Int}               # additional norm valuation
                                      # (for the dyadic case)
  isdyadic::Bool                      # 2 in p
  isramified::Bool                    # p ramified in E
  non_norm_rep                        # u in K*\N(E*)
  ni::Vector{Int}                     # ni for the ramified, dyadic case

  function LocalGenusHerm{S, T}() where {S, T}
    z = new{S, T}()
    return z
  end
end

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

function Base.show(io::IO, ::MIME"text/plain", G::LocalGenusHerm)
  compact = get(io, :compact, false)
  if !compact
    print(io, "Local genus symbol at ")
    print(IOContext(io, :compact => true), prime(G), ":")
    print(io, "\n")
  end
  if isdyadic(G) && isramified(G)
    for i in 1:length(G)
      print(io, "(", scale(G, i), ", ", rank(G, i), ", ",
            det(G, i) == 1 ? "+" : "-", ", ", norm(G, i), ")")
    end
  else
    for i in 1:length(G)
      print(io, "(", scale(G, i), ", ", rank(G, i), ", ",
            det(G, i) == 1 ? "+" : "-",  ")")
    end
  end
end

function Base.show(io::IO, G::LocalGenusHerm)
  if isdyadic(G) && isramified(G)
    for i in 1:length(G)
      print(io, "(", scale(G, i), ", ", rank(G, i), ", ",
            det(G, i) == 1 ? "+" : "-", ", ", norm(G, i), ")")
      if i < length(G)
        print(io, " ")
      end
    end
  else
    for i in 1:length(G)
      print(io, "(", scale(G, i), ", ", rank(G, i),
            ", ", det(G, i) == 1 ? "+" : "-",  ")")
      if i < length(G)
        print(io, " ")
      end
    end
  end
end

################################################################################
#
#  Basic properties
#
################################################################################

@doc Markdown.doc"""
    scale(G::LocalGenusHerm, i::Int) -> Int

Given a genus symbol for Hermitian lattices over $E/K$, return the $\mathfrak
P$-valuation of the scale of the $i$th Jordan block of $G$, where $\mathfrak
P$ is a prime ideal lying over $\mathfrak p$.
"""
scale(G::LocalGenusHerm, i::Int) = G.data[i][1]

@doc Markdown.doc"""
    scales(G::LocalGenusHerm, i::Int) -> Vector{Int}

Given a genus symbol for Hermitian lattices over $E/K$, return the $\mathfrak
P$-valuation of the scales of the Jordan blocks of $G$, where $\mathfrak
P$ is a prime ideal lying over $\mathfrak p$.
"""
scales(G::LocalGenusHerm) = map(i -> scale(G, i), 1:length(G))::Vector{Int}

@doc Markdown.doc"""
    rank(G::LocalGenusHerm, i::Int)

Given a genus symbol for Hermitian lattices over $E/K$, return the rank of the
$i$th Jordan block of $G$.
"""
rank(G::LocalGenusHerm, i::Int) = G.data[i][2]

@doc Markdown.doc"""
    rank(G::LocalGenusHerm) -> Int

Given a genus symbol $G$, return the rank of a lattice in $G$.
"""
function rank(G::LocalGenusHerm)
  return sum(rank(G, i) for i in 1:length(G))
end

@doc Markdown.doc"""
    ranks(G::LocalGenusHerm)

Given a genus symbol for Hermitian lattices over $E/K$, return the ranks of the
Jordan block of $G$.
"""
ranks(G::LocalGenusHerm) = map(i -> rank(G, i), 1:length(G))::Vector{Int}

@doc Markdown.doc"""
    det(G::LocalGenusHerm, i::Int) -> Int

Given a genus symbol for Hermitian lattices over $E/K$, return the determinant of the
$i$th Jordan block of $G$. This will be `1` or `-1` depending on whether the
determinant is local norm or not.
"""
det(G::LocalGenusHerm, i::Int) = G.data[i][3]

@doc Markdown.doc"""
    det(G::LocalGenusHerm) -> Int

Given a genus symbol $G$, return the determinant of a lattice in $G$. This will be
`1` or `-1` depending on whether the determinant is a local norm or not.
"""
function det(G::LocalGenusHerm)
  return prod(det(G, i) for i in 1:length(G))
end

@doc Markdown.doc"""
    dets(G::LocalGenusHerm) -> Vector{Int}

Given a genus symbol for Hermitian lattices over $E/K$, return the determinants
of the Jordan blocks of $G$. These will be `1` or `-1` depending on whether the
determinant is local norm or not.
"""
dets(G::LocalGenusHerm) = map(i -> det(G, i), 1:length(G))::Vector{Int}

@doc Markdown.doc"""
    norm(G::LocalGenusHerm, i::Int) -> Int

Given a dyadic genus symbol for Hermitian lattices over $E/K$ at a prime
$\mathfrak p$, return the $\mathfrak p$-valuation of the norm of the $i$ Jordan
block of $G$.
"""
norm(G::LocalGenusHerm, i::Int) = begin @assert isdyadic(G); G.norm_val[i] end # this only works if it is dyadic

@doc Markdown.doc"""
    norms(G::LocalGenusHerm) -> Vector{Int}

Given a dyadic genus symbol for Hermitian lattices over $E/K$ at a prime
$\mathfrak p$, return the $\mathfrak p$-valuations of the norms of the Jordan
blocks of $G$.
"""

@doc Markdown.doc"""
    isramified(G::LocalGenus) -> Bool

Given a genus symbol for Hermitian lattices at a prime $\mathfrak p$, return
whether $\mathfrak p$ is ramified.
"""
isramified(G::LocalGenusHerm) = G.isramified

@doc Markdown.doc"""
    isdyadic(G::LocalGenus) -> Bool

Given a genus symbol for Hermitian lattices at a prime $\mathfrak p$, return
whether $\mathfrak p$ is dyadic.
"""
isdyadic(G::LocalGenusHerm) = G.isdyadic

#data(G::LocalGenusHerm) = G.data

@doc Markdown.doc"""
    length(G::LocalGenusHerm) -> Int

Return the number of Jordan blocks in a Jordan decomposition of $G$.
"""
length(G::LocalGenusHerm) = length(G.data)

@doc Markdown.doc"""
    base_field(G::LocalGenusHerm) -> NumField

Given a local genus symbol of Hermitian lattices over $E/K$, return $E$.
"""
base_field(G::LocalGenusHerm) = G.E

@doc Markdown.doc"""
    prime(G::LocalGenus) -> NfOrdIdl

Given a local genus symbol of Hermitian lattices at a prime $\mathfrak p$,
return $\mathfrak p$.
"""
prime(G::LocalGenusHerm) = G.p

################################################################################
#
#  A local non-norm
#
################################################################################

function _non_norm_rep(G::LocalGenusHerm)
  if isdefined(G, :non_norm_rep)
    return G.non_norm_rep
  end

  z = _non_norm_rep(base_field(G), base_field(base_field(G)), prime(G))
  G.non_norm_rep = z
  return z
end

################################################################################
#
#  A local unifomizer
#
################################################################################

@doc Markdown.doc"""
    uniformizer(G::LocalGenusHerm) -> NumFieldElem

Given a local genus symbol of Hermitian lattices over $E/K$ at a prime $\mathfrak p$,
return a generator for the largest invariant ideal of $E$.
"""
function uniformizer(G::LocalGenusHerm)
  E = base_field(G)
  K = base_field(E)
  if isramified(G)
    lP = prime_decomposition(maximal_order(E), prime(G))
    @assert length(lP) == 1 && lP[1][2] == 2
    Q = lP[1][1]
    pi = p_uniformizer(Q)
    A = automorphisms(E)
    uni = A[1](elem_in_nf(pi)) * A[2](elem_in_nf(pi))
    @assert iszero(coeff(uni, 1))
    @assert islocal_norm(E, coeff(uni , 0), prime(G))
    return coeff(uni, 0)
  else
    return uniformizer(prime(G))
  end
end

################################################################################
#
#  Additional dyadic information
#
################################################################################

# Get the "ni" for the ramified dyadic case
function _get_ni_from_genus(G::LocalGenusHerm)
  @assert isramified(G)
  @assert isdyadic(G)
  t = length(G)
  z = Vector{Int}(undef, t)
  for i in 1:t
    ni = minimum(2 * max(0, scale(G, i) - scale(G, j)) + 2 * norm(G, j) for j in 1:t)
    z[i] = ni
  end
  return z
end

################################################################################
#
#  Determinant representative
#
################################################################################

# TODO: I don't understand what this is doing. I know that d tells me if it is
#       a norm or not. Why do I have to multiply with a uniformizer
@doc Markdown.doc"""
    det_representative(G::LocalGenusHerm) -> NumFieldElem

Return a representative for the norm class of the determinant of $G$.
"""
function det_representative(G::LocalGenusHerm)
  z = G.data
  d = det(G)
  v = sum(scale(G, i) * rank(G, i) for i in 1:length(G))
  if isramified(G)
    v = div(v, 2)
  end
  if d == 1
    u = one(base_field(base_field(G)))
  else
    @assert isramified(G)
    u = _non_norm_rep(G)
  end
  return u * uniformizer(G)^v
end

################################################################################
#
#  Gram matrix
#
################################################################################

@doc Markdown.doc"""
    gram_matrix(G::LocalGenusHerm)

Return a matrix $M$, such that a lattice with Gram matrix $M$ is an element of
the given genus.
"""
function gram_matrix(G::LocalGenusHerm)
  return diagonal_matrix(dense_matrix_type(base_field(G))[gram_matrix(G, i) for i in 1:length(G)])
end

function gram_matrix(G::LocalGenusHerm, l::Int)
  i = scale(G, l)
  m = rank(G, l)
  d = det(G, l)
  E = base_field(G)
  K = base_field(E)
  p = elem_in_nf(p_uniformizer(prime(G)))
  A = automorphisms(E)
  _a = gen(E)
  conj = A[1](_a) == _a ? A[2] : A[1]
  if d == 1
    u = one(K)
  else
    u = _non_norm_rep(G)
  end

  if !isramified(G)
    return diagonal_matrix([E(p)^i for j in 1:m])
  end

  # ramified

  lQ = prime_decomposition(maximal_order(E), prime(G))
  @assert length(lQ) == 1 && lQ[1][2] == 2
  Q = lQ[1][1]
  pi = elem_in_nf(p_uniformizer(Q))
  H = matrix(E, 2, 2, [0, pi^i, conj(pi)^i, 0])

  if !isdyadic(G)
    # non-dyadic
    if iseven(i)
      return diagonal_matrix(push!([E(p)^div(i, 2) for j in 1:(m - 1)], u * E(p)^div(i, 2)))
    else
      @assert iseven(m)
      return diagonal_matrix([H for j in 1:div(m, 2)])
    end
  end

  # dyadic

  k = norm(G, l)

  # I should cache this e
  e = valuation(different(maximal_order(E)), Q)

  if isodd(m)
    # odd rank
    @assert 2*k == i
    r = div(m - 1, 2)
    nn = mod(m, 4)
    if d == 1
      discG = K(1)
    else
      discG = elem_in_nf(_non_norm_rep(G))
    end
    if nn == 0 || nn == 1
      discG = discG
    else
      discG = -discG
    end
    @assert iseven(i)
    if islocal_norm(E, discG*p^(-m * div(i, 2)), prime(G))
      u = K(1)
    else
      u = _non_norm_rep(G)
    end

    U = matrix(E, 1, 1, [u * p^(div(i, 2))])
    return diagonal_matrix(append!(typeof(U)[U], [H for j in 1:r]))
  else
    # even rank
    r = div(m, 2) - 1
    if islocal_norm(E, K((-1)^div(m, 2)), prime(G)) == (d == 1)
      @assert i + e >= 2 * k
      @assert 2 * k >= i
      U = matrix(E, 2, 2, [p^k, pi^i, conj(pi)^i, 0])
      return diagonal_matrix(append!(typeof(U)[U], [H for j in 1:r]))
    else
      @assert i + e > 2 * k
      @assert 2 * k >= i
      u = _non_norm_rep(G)
      u0 = u - 1
      U = matrix(E, 2, 2, [p^k, pi^i, conj(pi)^i, -p^(i -k) * u0])
      return diagonal_matrix(append!(typeof(U)[U], [H for j in 1:r]))
    end
  end
end

################################################################################
#
#  Representative
#
################################################################################

@doc Markdown.doc"""
    representative(G::LocalGenusHerm) -> HermLat

Given a local genus, return a Hermitian lattice contained in this genus.
"""
function representative(G::LocalGenusHerm)
  E = G.E
  return lattice(hermitian_space(E, gram_matrix(G)))
end

################################################################################
#
#  Constructor
#
################################################################################

function genus(::Type{HermLat}, E::S, p::T, data::Vector{Tuple{Int, Int, Int}}) where {S <: NumField, T}
  z = LocalGenusHerm{S, T}()
  z.E = E
  z.p = p
  z.isdyadic = isdyadic(p)
  z.isramified = isramified(maximal_order(E), p)
  @assert !(isramified(z) && isdyadic(z))
  z.data = data
  return z
end

function genus(::Type{HermLat}, E::S, p::T, data::Vector{Tuple{Int, Int, Int, Int}}) where {S <: NumField, T}
  z = LocalGenusHerm{S, T}()
  z.E = E
  z.p = p
  z.isdyadic = isdyadic(p)
  z.isramified = isramified(maximal_order(E), p)
  if isramified(z) && isdyadic(z)
    z.data = Tuple{Int, Int, Int}[Base.front(v) for v in data]
    z.norm_val = Int[v[end] for v in data]
    z.ni = _get_ni_from_genus(z)
  else
    z.data = Tuple{Int, Int, Int}[Base.front(v) for v in data]
  end
  return z
end

################################################################################
#
#  Genus symbol of a lattice
#
################################################################################

# TODO: caching
# TODO: better documentation

@doc Markdown.doc"""
    genus(L::HermLat, p::NfOrdIdl) -> LocalGenusHerm

Returns the genus of $L$ at the prime ideal $\mathfrak p$.

See [Kir16, Definition 8.3.1].
"""
function genus(L::HermLat, q)
  if order(q) != base_ring(base_ring(L))
    p = minimum(q)
  else
    p = q
  end

  sym = _genus_symbol(L, p)
  G = genus(HermLat, nf(base_ring(L)), p, sym)
  # Just for debugging
  @hassert :Lattice 1 begin
    if isdyadic(G) && isramified(G)
      GG = _genus_symbol_kirschmer(L, p)
      all(i -> GG[i][4] == G.ni[i], 1:length(G))
    else
      true
    end
  end
  return G
end

################################################################################
#
#  Test if lattice is contained in local genus
#
################################################################################

@doc Markdown.doc"""
    in(L::HermLat, G::LocalGenusHerm) -> Bool

Test if the lattice $L$ is contained in the local genus $G$.
"""
Base.in(L::HermLat, G::LocalGenusHerm) = genus(L, prime(G)) == G

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

@doc Markdown.doc"""
    ==(G1::LocalGenusHerm, G2::LocalGenusHerm)

Test if two local genera are equal.
"""
function ==(G1::LocalGenusHerm, G2::LocalGenusHerm)
  if base_field(G1) != base_field(G2)
    return false
  end

  if prime(G1) != prime(G2)
    return false
  end

  if length(G1) != length(G2)
    return false
  end

  t = length(G1)

  p = prime(G1)

  # We now check the Jordan type

  if any(i -> scale(G1, i) != scale(G2, i), 1:t)
    return false
  end

  if any(i -> (rank(G1, i) != rank(G2, i)), 1:t)
    return false
  end

  if det(G1) != det(G2) # rational spans must be isometric
    return false
  end

  if any(i -> (rank(G1, i) != rank(G2, i)), 1:t)
    return false
  end

  if !isramified(G1) # split or unramified
    return true
    # Everything is normal and the Jordan decomposition types agree
    #return all(det(G1, i) == det(G2, i) for i in 1:t)
  end

  if isramified(G1) && !isdyadic(G1) # ramified, but not dyadic
    # If s(L_i) is odd, then L_i = H(s(L_i))^(rk(L_i)/2) = H(s(L_i'))^(rk(L_i')/2) = L_i'
    # So all L_i, L_i' are locally isometric, in particular L_i is normal if and only if L_i' is normal
    # If s(L_i) = s(L_i') is even, then both L_i and L_i' have orthgonal bases, so they are
    # in particular normal.

    # Thus we only have to check Theorem 3.3.6 4.
    return all(i -> det(G1, i) == det(G2, i), 1:t)
    # TODO: If think I only have to do something if the scale is even. Check this!
  end

  # Dyadic ramified case

  # First check if L_i is normal if and only if L_i' is normal, i.e.,
  # that the Jordan decompositions agree
  for i in 1:t
    if (scale(G1, i) == 2 * norm(G1, i)) != (scale(G2, i) == 2 * norm(G2, i))
      return false
    end
  end

  if any(i -> G1.ni[i] != G2.ni[i], 1:t)
    return false
  end

  E = base_field(G1)
  lQ = prime_decomposition(maximal_order(E), p)
  @assert length(lQ) == 1 && lQ[1][2] == 2
  Q = lQ[1][1]

  e = valuation(different(maximal_order(E)), Q)

  for i in 1:(t - 1)
    dL1prod = prod(det(G1, j) for j in 1:i)
    dL2prod = prod(det(G2, j) for j in 1:i)

    d = dL1prod * dL2prod

    if d != 1
      if 2 * (e - 1) < G1.ni[i] + G1.ni[i + 1] - 2 * scale(G1, i)
        return false
      end
    end
  end

  return true
end

function _genus_symbol(L::HermLat, q)
  if order(q) != base_ring(base_ring(L))
    p = minimum(q)
  else
    p = q
  end
  @assert order(p) == base_ring(base_ring(L))
  B, G, S = jordan_decomposition(L, p)
  R = base_ring(L)
  E = nf(R)
  K = base_field(E)
  if !isdyadic(p) || !isramified(R, p)
    # The last entry is a dummy to make the compiler happier
    sym = Tuple{Int, Int, Int, Int}[ (S[i], nrows(B[i]), islocal_norm(E, coeff(det(G[i]), 0), p) ? 1 : -1, 0) for i in 1:length(B)]
  else
    P = prime_decomposition(R, p)[1][1]
    pi = E(K(uniformizer(p)))
    sym = Tuple{Int, Int, Int, Int}[]
    for i in 1:length(B)
      normal = _get_norm_valuation_from_gram_matrix(G[i], P) == S[i]
      GG = diagonal_matrix([pi^(max(0, S[i] - S[j])) * G[j] for j in 1:length(B)])
      v = _get_norm_valuation_from_gram_matrix(GG, P)
      @assert v == valuation(R * norm(lattice(hermitian_space(E, GG), identity_matrix(E, nrows(GG)))), P)
      r = nrows(B[i]) # rank
      s = S[i] # P-valuation of scale(L_i)
      det_class = islocal_norm(E, coeff(det(G[i]), 0), p) ? 1 : -1  # Norm class of determinant
      normi = _get_norm_valuation_from_gram_matrix(G[i], P) # P-valuation of norm(L_i)
      @assert mod(normi, 2) == 0 # I only want p-valuation
      push!(sym, (s, r, det_class, div(normi, 2)))
    end
  end
  return sym
end

################################################################################
#
#  Global genus
#
################################################################################

mutable struct GenusHerm
  E
  primes::Vector
  LGS::Vector
  rank::Int
  signatures

  function GenusHerm(E, r, LGS::Vector, signatures)
    primes = Vector(undef, length(LGS))

    for i in 1:length(LGS)
      primes[i] = prime(LGS[i])
      @assert r == rank(LGS[i])
    end
    z = new(E, primes, LGS, r, signatures)
    return z
  end
end

@doc Markdown.doc"""
    ==(G1::GenusHerm, G2::GenusHerm)

Test if two genera are equal.
"""
function ==(G1::GenusHerm, G2::GenusHerm)
  if G1.E != G2.E
    return false
  end

  if length(G1.primes) != length(G2.primes)
    return false
  end

  if G1.signatures != G2.signatures
    return false
  end

  for p in G1.primes
    if !(p in G2.primes)
      return false
    end
    if G1[p] != G2[p]
      return false
    end
  end

  return true
end

function Base.show(io::IO, G::GenusHerm)
  print(io, "Genus symbol over")
  print(io, G.E)
  print(io, "\n", "and local genera",)
  for g in G.LGS
    print(io, "\n")
    print(IOContext(io, :compact => true), prime(g), " => ", g)
  end
  print(io, "\n and signature")
  for (pl, v) in G.signatures
    print(io, "\n")
    _print_short(io, pl.r)
    print(io, " => ")
    print(io, v)
  end
end

function _print_short(io::IO, a::arb)
  r = BigFloat(a)
  s = string(r)
  if length(s) >= 10
    ss = s[1:9] * "…"
  else
    ss = s
  end
  print(io, ss)
end

function _print_short(io::IO, a::acb)
  _print_short(io, real(a))
  if !iszero(imag(a))
    print(io, " + ")
    _print_short(io, imag(a))
    print(io, " * i")
  end
end

################################################################################
#
#  Test if lattice is contained in genus
#
################################################################################

@doc Markdown.doc"""
    in(L::HermLat, G::LocalGenusHerm) -> Bool

Test if the lattice $L$ is contained in the local genus $G$.
"""
Base.in(L::HermLat, G::GenusHerm) = genus(L) == G

# This could be made more efficient

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

function Base.show(io::IO, ::MIME"text/plain", G::GenusHerm)
  print(io, "Global genus symbol\n")
  for i in 1:length(G.primes)
    print(IOContext(io, :compact => true), G.primes[i], " => ", G.LGS[i],)
    if i < length(G.primes)
      print(io, "\n")
    end
  end
  for (pl, v) in G.signatures
    print(io, "\n")
    _print_short(io, pl.r)
    print(io, " => ")
    print(io, v)
  end
end

################################################################################
#
#  Constructor
#
################################################################################

@doc Markdown.doc"""
    genus(L::Vector{LocalGenusHerm},
          signatures::Vector{Tuple{InfPlc, Int}}) -> GenusHerm

Given a set of local genera and signatures, return the corresponding global
genus.
"""
function genus(L::Vector, signatures::Dict{InfPlc, Int})
  @assert !isempty(L)
  @assert all(N >= 0 for (_, N) in signatures)
  if !_check_global_genus(L, signatures)
    throw(error("Invariants violate the product formula."))
  end
  r = rank(first(L))
  E = base_field(first(L))
  return GenusHerm(E, r, L, signatures)
end

function genus(L::Vector, signatures::Vector{Tuple{InfPlc, Int}})
  return genus(L, Dict(signatures))
end

@doc Markdown.doc"""
    genus(L::HermLat) -> GenusHerm

Given a Hermitian lattice, return the genus it belongs to.
"""
function genus(L::HermLat)
  bad = bad_primes(L)
  for p in support(discriminant(base_ring(L)))
    if isdyadic(p) && !(p in bad)
      push!(bad, p)
    end
  end
  S = real_places(base_field(base_field(L)))
  D = diagonal(rational_span(L))
  signatures = Dict{InfPlc, Int}(s => count(d -> isnegative(d, s), D) for s in S)
  return genus([genus(L, p) for p in bad], signatures)
end

################################################################################
#
#  Basic access
#
################################################################################

@doc Markdown.doc"""
    primes(G::GenusHerm) -> Vector{NfOrdIdl}

Return the primes of a global genus symbol.
"""
primes(G::GenusHerm) = G.primes

@doc Markdown.doc"""
    signatures(G::GenusHerm) -> Dict{InfPlc, Int}

Return the signatures at the infinite places, which for each real place is
given by the negative index of inertia.
"""
signatures(G::GenusHerm) = G.signatures

################################################################################
#
#  Consistency check
#
################################################################################

function _check_global_genus(LGS, signatures)
  _non_norm = _non_norm_primes(LGS)
  P = length(_non_norm)
  I = length([(s, N) for (s, N) in signatures if mod(N, 2) == 1])
  if mod(P + I, 2) == 1
    return false
  end
  return true
end

################################################################################
#
#  Primes at which the determinant is not a local norm
#
################################################################################

function _non_norm_primes(LGS::Vector)
  z = []
  for g in LGS
    p = prime(g)
    d = det(g)
    if d != 1
      push!(z, p)
    end
  end
  return z
end

################################################################################
#
#  Convenient functions
#
################################################################################

function Base.getindex(G::GenusHerm, P)
  i = findfirst(isequal(P), G.primes)
  if i === nothing
    throw(error("No local genus symbol at $P"))
  end
  return G.LGS[i]
end

################################################################################
#
#  Compute representatives of genera
#
################################################################################

function _hermitian_form_with_invariants(E, dim, P, N)
  #@show E, dim, N, [minimum(p) for p in P]
  K = base_field(E)
  R = maximal_order(K)
#  require forall{n: n in N | n in {0..dim}}: "Number of negative entries is impossible";
  infinite_pl = [ p for p in real_places(K) if _decomposition_number(E, p) == 1 ]
  length(N) != length(infinite_pl) && error("Wrong number of real places")
  S = maximal_order(E)
  prim = [ p for p in P if length(prime_decomposition(S, p)) == 1 ] # only take non-split primes
  #@show prim
  I = [ p for p in keys(N) if isodd(N[p]) ]
  !iseven(length(I) + length(P)) && error("Invariants do not satisfy the product formula")
  e = gen(E)
  x = 2 * e - trace(e)
  b = coeff(x^2, 0) # b = K(x^2)
  a = _find_quaternion_algebra(b, prim, I)
  #@show a
  D = elem_type(E)[]
  for i in 1:(dim - 1)
    if length(I) == 0
      push!(D, one(E))
    else
      push!(D, E(_weak_approximation(infinite_pl, [N[p] >= i ? -1 : 1 for p in infinite_pl])))
    end
  end
  push!(D, a * prod(D))
  Dmat = diagonal_matrix(D)
  dim0, P0, N0 = _hermitian_form_invariants(Dmat)
  @assert dim == dim0
  @assert Set(prim) == Set(P0)
  @assert N == N0
  return Dmat
end

function _hermitian_form_invariants(M)
  E = base_ring(M)
  K = base_field(E)
  @assert degree(E) == 2
  A = automorphisms(E)
  a = gen(E)
  v = A[1](a) == a ? A[2] : A[1]

  @assert M == transpose(_map(M, v))
  d = coeff(det(M), 0) # K(det(M))
  P = Dict()
  for p in keys(factor(d * maximal_order(K)))
    if islocal_norm(E, d, p)
      continue
    end
    P[p] = true
  end
  for p in keys(factor(discriminant(maximal_order(E))))
    if islocal_norm(E, d, p)
      continue
    end
    P[p] = true
  end
  D = diagonal(_gram_schmidt(M, v)[1])
  I = Dict([ p=>length([coeff(d, 0) for d in D if isnegative(coeff(d, 0), p)]) for p in real_places(K) if _decomposition_number(E, p) == 1])
  return ncols(M), collect(keys(P)), I
end

base_field(G::GenusHerm) = G.E

rank(G::GenusHerm) = G.rank

function representative(G::GenusHerm)
  P = _non_norm_primes(G.LGS)
  E = base_field(G)
  V = hermitian_space(E, _hermitian_form_with_invariants(base_field(G), rank(G), P, G.signatures))
  @vprint :Lattice 1 "Finding maximal integral lattice ..."
  M = maximal_integral_lattice(V)
  @vprint :Lattice 1 "done\n"
  for g in G.LGS
    p = prime(g)
    @vprint :Lattice 1 "Finding representative for $g ..."
    L = representative(g)
    @vprint :Lattice 1 "done\n"
    @hassert :Lattice 1 genus(L, p) == g
    #@show coefficient_ideals(pseudo_matrix(L))
    #@show matrix(pseudo_matrix(L))
    @vprint :Lattice 1 "Finding sublattice ..."
    M = find_lattice(M, L, p)
    @vprint :Lattice 1 "done\n"
  end
  return M
end

################################################################################
#
#  Enumeration of local genera
#
################################################################################

@doc Markdown.doc"""
    local_genera_hermitian(E::NumField, p::NfOrdIdl, rank::Int,
                 det_val::Int, max_scale::Int) -> Vector{LocalGenusHerm}

Return all local genera of Hermitian lattices over $E$ at $\mathfrak p$ with
rank `rank`, scale valuation bounded by `max_scale` and determinant valuation
bounded by `det_val`.
"""
function local_genera_hermitian(E, p, rank, det_val, max_scale, is_ramified = isramified(maximal_order(E), p))
  #@show E, p, rank, det_val, max_scale, is_ramified
  is_inert = !is_ramified && length(prime_decomposition(maximal_order(E), p)) == 1
  if is_ramified
    # the valuation is with respect to p
    # but the scale is with respect to P
    # in the ramified case p = P**2 and thus
    # the determinant of a block is
    # P^(scale*rank) = p^(scale*rank/2)
    det_val *= 2
  end

  K = number_field(order(p))

  scales_rks = Vector{Tuple{Int, Int}}[] # possible scales and ranks

  for rkseq in _integer_lists(rank, max_scale + 1)
    d = 0
    pgensymbol = Tuple{Int, Int}[]
    for i in 0:(max_scale + 1) - 1
      d += i * rkseq[i + 1]
      if rkseq[i + 1] != 0
        push!(pgensymbol, (i, rkseq[i + 1]))
      end
    end
    if d == det_val
        push!(scales_rks, pgensymbol)
    end
  end

  if !is_ramified
    # I add the 0 to make the compiler happy
    symbols = Vector{LocalGenusHerm{typeof(E), typeof(p)}}(undef, length(scales_rks))
    for i in 1:length(scales_rks)
      g = scales_rks[i]
      z = Tuple{Int, Int, Int, Int}[]
      for b in g
        # We have to be careful.
        # If p is inert, then the norm is not surjective.
        if !is_inert || iseven(b[1] * b[2])
          push!(z, (b[1], b[2], 1, 0))
        else
          push!(z, (b[1], b[2], -1, 0))
        end
      end
      symbols[i] = genus(HermLat, E, p, z)
    end
    return symbols
  end

  scales_rks = Vector{Tuple{Int, Int}}[g for g in scales_rks if all((mod(b[1]*b[2], 2) == 0) for b in g)]

  symbols = LocalGenusHerm{typeof(E), typeof(p)}[]
  hyperbolic_det = hilbert_symbol(K(-1), gen(K)^2//4 - 1, p)
  if !isdyadic(p) # non-dyadic
    for g in scales_rks
      n = length(g)
      dets = Vector{Int}[]
      for b in g
        if mod(b[1], 2) == 0
          push!(dets, [1, -1])
        end
        if mod(b[1], 2) == 1
          push!(dets, [hyperbolic_det^(div(b[2], 2))])
        end
      end

      for d in Iterators.product(dets...)
        g2 = Vector{Tuple{Int, Int, Int, Int}}(undef, length(g))
        for k in 1:n
          # Again the 0 for dummy purposes
          g2[k] = (g[k]..., d[k], 0)
        end
        push!(symbols, genus(HermLat, E, p, g2))
      end
    end
    return symbols
  end

  # Ramified case
  lp = prime_decomposition(maximal_order(E), p)
  @assert length(lp) == 1 && lp[1][2] == 2
  P = lp[1][1]

  e = valuation(different(maximal_order(E)), P)
  # only for debugging
  scales_rks = reverse(scales_rks)

  for g in scales_rks
    n = length(g)
    det_norms = Vector{Vector{Int}}[]
    for b in g
      if mod(b[2], 2) == 1
        @assert iseven(b[1])
        push!(det_norms, [[1, div(b[1], 2)], [-1, div(b[1], 2)]])
      end
      if mod(b[2], 2) == 0
        dn = Vector{Int}[]
        i = b[1]
        # (i + e) // 2 => k >= i/2
        for k in (ceil(Int, Int(i)//2)):(div(Int(i + e), 2) - 1)
          push!(dn, Int[1, k])
          push!(dn, Int[-1, k])
        end
        push!(dn, Int[hyperbolic_det^(div(b[2], 2)), div(i + e, 2)])
        if mod(i + e, 2) == 1
          push!(dn, Int[-hyperbolic_det^(div(b[2], 2)), div(i + e, 2)])
        end
        push!(det_norms, dn)
      end
    end
    for dn in Iterators.product(det_norms...)
      g2 = Vector{Tuple{Int, Int, Int, Int}}(undef, length(g))
      for k in 1:n
        g2[k] = (g[k]..., dn[k]...)
      end
      h = genus(HermLat, E, p, g2)
      if !(h in symbols)
        push!(symbols, h)
      end
    end
  end
  return symbols
end

@doc Markdown.doc"""
    local_genera_hermitian(E::NumField, p::NfOrdIdl, rank::Int,
                 det_val::Int, max_scale::Int) -> Vector{LocalGenusHerm}

Return all local genera of Hermitian lattices over $E$ at $\mathfrak p$ with
rank `rank`, scale valuation bounded by `max_scale` and determinant valuation
bounded by `det_val`.
"""

function genera_hermitian(E, rank, signatures, determinant; max_scale = nothing)
  K = base_field(E)
  OE = maximal_order(E)
  if max_scale === nothing
    _max_scale = determinant
  else
    _max_scale = max_scale
  end

  primes = support(discriminant(OE))
  for p in support(norm(determinant))
    if !(p in primes)
      push!(primes, p)
    end
  end

  local_symbols = Vector{LocalGenusHerm{typeof(E), typeof(primes[1])}}[]

  ms = norm(_max_scale)
  ds = norm(determinant)
  for p in primes
    det_val = valuation(ds, p)
    mscale_val = valuation(ms, p)
    det_val = div(det_val, 2)
    if !isramified(OE, p)
      mscale_val = div(mscale_val, 2)
    end
    push!(local_symbols, local_genera_hermitian(E, p, rank, det_val, mscale_val))
  end

  res = []
  it = Iterators.product(local_symbols...)
  for gs in it
    c = collect(gs)
    b = _check_global_genus(c, signatures)
    if b
      push!(res, GenusHerm(E, rank, c, signatures))
    end
  end

  return res
end

################################################################################
#
#  Genus representatives
#
################################################################################

# Return b, p, bad
# b = isdefinite(L)
# p = prime ideal of base_ring(L) which can be used for the neighbor method
# bad = bad primes of L, where L,p is not modular or p is dyadic and dividing disc(S)
function smallest_neighbour_prime(L)
  S = base_ring(L)
  R = base_ring(S)
  lp = bad_primes(L)
  bad = [p for p in lp if !ismodular(L, p)[1] ]
  for (p,_) in factor(discriminant(S))
    if isdyadic(p) && !(p in bad)
      push!(bad, p)
    end
  end

  if !isdefinite(L)
    return false, 1*S, bad
  end

  # TODO: This does not find the prime ideal with smallest norm,
  # but with smallest minimum ...

  m = rank(L)
  p = 1
  P = ideal_type(R)[]
  while length(P) == 0
    p = next_prime(p)
    lp = [ p[1] for p in prime_decomposition(R, p)]
    P = setdiff(lp, bad)
    if m == 2
      P = filter(p -> isisotropic(L, p), P)
    end
  end
  Q = prime_decomposition(S, P[1])[1][1]

  if isempty(bad)
    I = 1 * S
  else
    I = prod(bad) * S
  end
  n = absolute_norm(Q)
  if n >= 1000
    PP = prime_ideals_up_to(S, 1000)
    for QQ in PP
      if !iscoprime(QQ, I)
        continue
      end

      if isisotropic(L, QQ)
        return true, QQ, bad
      end
    end
  end
  PP = prime_ideals_up_to(S, n)
  for QQ in PP
    if !iscoprime(QQ, I)
      continue
    end
    if isisotropic(L, QQ)
      return true, QQ, bad
    end
  end
  throw(error("Impossible"))
end

function genus_generators(L)
  R = base_ring(L)
  E = nf(R)
  D = different(R)
  b, P0, bad = smallest_neighbour_prime(L)
  if isempty(bad)
    bad_prod = 1 * base_ring(R)
  else
    bad_prod = prod(bad)
  end

  # First the ideals coming from the C/C0 quotient
  Eabs, EabstoE, _ = absolute_field(E)
  Rabs = maximal_order(Eabs)
  C, h = class_group(Rabs)
  RR = base_ring(R)
  C0 = support(D)
  CC, hh = class_group(RR)
  for p in find_gens(pseudo_inv(h), PrimesSet(2, -1))[1]
    @show p
    if !(p in C0)
      push!(C0, p * R)
    end
  end
  Q0, q0 = quo(C, [ h\ideal(Rabs, [Rabs(EabstoE\b) for b in absolute_basis(i)]) for i in C0])
  q00 = pseudo_inv(q0) * h
  PP = []
  if iseven(rank(L))
    for (P, e) in factor(D)
      G = genus(L, P)
      if any(i -> isodd(rank(G, i)), 1:length(G))
        continue
      elseif !isdyadic(P)
        if any(i -> iseven(scale(G, i)), 1:length(G))
          continue
        end
      else
        if any(i -> isodd(e  + scale(G, i)) || (e + scale(G, i) != G.ni[i]), 1:length(G))
          continue
        end
      end
      push!(PP, P)
    end

    if !isempty(PP)
      U, f = unit_group_fac_elem(Rabs)
      UU, ff = unit_group_fac_elem(RR)
      nnorm = hom(U, UU, [ff\FacElem(nf(RR)(norm(f(U[i])))) for i in 1:ngens(U)])
      l = length(PP)
      VD = Int[ valuation(D, P) for P in PP ]
      K, k = kernel(nnorm)
      F = GF(2, cached = false)
      V = VectorSpace(F, length(PP))
      S = elem_type(V)[]
      for u in gens(K)
        z = elem_type(F)[]
        for i in 1:length(PP)
          zz = R(EabstoE(evaluate(f(k(u))))) - 1
          if iszero(zz) || (valuation(zz, PP[i]) >= VD[i])
            push!(z, F(0))
          else
            push!(z, F(1))
          end
        end
        push!(S, V(z))
      end
      #S = [ V(elem_type(F)[valuation(R(EabstoE(evaluate(f(k(u)))) - 1), PP[i]) >= VD[i] ? F(0) : F(1) for i in 1:length(PP)]) for u in gens(K)]
      _T, _ = sub(V, S)
      W, w = quo(V, _T)
      @show dim(W)
      if dim(W) == 0
        PP = []
      end
    end
  end

  Gens = []
  if isempty(PP)
    S = GrpAbFinGenElem[]
    Q, q = quo(Q0, S)
    Work = isdefinite(L) ? [ P0 ] : []
    p = 2
    while order(Q) > 1
      while isempty(Work)
        p = next_prime(p)
        Work = [ QQ for QQ in support(p * R) if issplit(QQ) && valuation(bad, QQ) == 0 ]
      end
      P = popfirst!(Work)
      c = q00\P
      o = order(q(c))
      if o != 1
        push!(S, c)
        Q, q = quo(Q0, S)
        push!(Gens, (P, o))
      end
    end
  else
    ll = Int(order(Q0))
    cocycle = Matrix(undef, ll, ll)
    for i in 1:ll
      for j in 1:ll
        cocycle = zero(W)
      end
    end
    C = collect(Q0)
    ideals = [ q00(C[i]) for i in 1:length(C) ]
    for i in 1:ll
      for j in 1:ll
        ij = findfirst(isequal(C[i] * C[j]), C)
        I = ideals[i] * ideals[j] * inv(ideals[ij])
        J = I * inv(a(I))
        ok, x = isprincipal(J)
        u = f(nnorm(-(ff\FacElem(nf(RR)(norm(x))))))
        x = x * u
        @assert norm(x) == 1
        y = w(V([ valuation(x - 1, PP[i]) >= VD[i] ? F(0) : F(1) for i in 1:length(PP)]))
        cocycle[i, j] = y
        cocycle[j, i] = y
      end
    end

    S = [(id(Q0), zero(W))]
    Work = isdefinite(L) ? [ P0 ] : []
    p = 2
    while length(S) != order(Q0) * length(W)
      while isempty(Work)
        p = next_prime(p)
        Work = [ QQ for QQ in support(p * R) if issplit(QQ) && valuation(bad, QQ) == 0 ]
      end
      P = popfirst!(Work)
      c = q00\P
      i = findfirst(isequal(c), C)
      I = P * inv(ideals[i])
      J = I * inv(a(I))
      ok, x = isprincipal(J)
      u = f(nnorm(-(ff\FacElem(nf(RR)(norm(x))))))
      x = x * u
      @assert norm(x) == 1
      y = V([ valuation(x - 1, PP[i]) >= VD[i] ? F(0) : F(1) for i in 1:length(PP)])
      idx = findfirst(isequal(P), PP)
      if idx !== nothing
        y = V([i == idx ? y[i] : y[i] + 1] for i in 1:dim(V)) #w(V([y[idx] = y[idx] + 1
      end
      elt = (c, w(y))
      elt1 = elt
      o = 1
      siz = length(S)
      while !(elt1 in S)
        j = findfirst(isequal(elt1[1]), C)
        @assert !isnothing(j)
        for l in 1:siz
          elt2 = S[l]
          k = findfirst(isequal(elt2[1]), C)
          @assert !isnothing(k)
          prod = (elt1[1] * elt2[1], elt1[2] + elt2[2] + cycycle[j, k])
          if !(prod in S)
            push!(S, prod)
          end
        end
        elt1 = (elt[1] * elt1[1], elt[2] + elt[1] + cocycle[i, j])
        o = o + 1
      end
      @assert length(S) == siz * o
      if o != 1
        push!(Gens, (P, o))
      end
    end
  end

  if isdefinite(L)
    return Gens, P0
  else
    return Gens, false
  end
end

# TODO: add max keyword
function genus_representatives(L)
  rank(L) < 2 && error("Rank of the lattice must be a least 2")
  definite = isdefinite(L)
  gens, P0 = genus_generators(L)
  a = involution(L)
  LL = [ L ]
  for g in gens
    if definite && g[1] == P0
      continue
    end
    I = g[1]^(g[2] - 1)
    J = a(I)
    N = neighbours_with_ppower(L, g[1], g[2] - 1)
    inter= []
    for i in 2:length(LL)
      M = pseudo_matrix(LL[i])
      IM = I * M
      JM = J * M
      inter = append!(inter, lattice(ambient_space(L), meet(IM + pseudo_matrix(x), JM) for x in N))
    end
    LL = vcat(LL, N, inter)
  end
  @assert length(LL) == prod(Int[g[2] for g in gens if !definite || g[1] != P0])
  @assert all(X -> genus(X) == genus(L), LL)

  if definite
    result = []
    for L in LL
      # Should never happen!
      @assert all(X -> !isisometric(X, L), result)
      result = append!(result, iterated_neighbours(L, P0))# : AutoOrbits:= AutoOrbits, Max:= Max);
      #max = max - length(result)
    end
    for i in 1:length(result)
      for j in 1:i-1
        @assert !isometric(result[i], result[j])
      end
    end
  else
    result = LL
  end
  return result
end

function _neighbours(L, P, result, max, callback = stdcallback)
#//  "Entering _Neighbours, #Result=", #Result, "Max=", Max;
  ok, scale = ismodular(L, P);
  !ok && error("Non-modular lattice!")
  R = base_ring(L)
  K = nf(R)
  a = involution(L)
  C = a(P)
  B = local_basis_matrix(L, P, type = :submodule)
  k, h = ResidueField(R, C)
  form = gram_matrix(ambient_space(L))
  special = false
  if scale != 0
    if isramified(R, minimum(P))
      special = isodd(scale)
      scale = div(scale + 1, 2)
    end
    form = E(elem_in_nf(uniformizer(minimum(P))))^(-scale) * form
  end
  n = rank(L)
  W = vector_space(k, n)

  LO = enumerate_lines(k, n)

  keep = true
  cont = true

  if P != C
    _G = T * form * _map(T, a)
    G = _map(_G, h)
    pi = p_uniformizer(C)
    pih = h(pi)
    for i in 1:length(LO)
      w = LO[i]
      x = [ sum(B[i, j] * (h\w[i]) for i in 1:n) for j in 1:ncols(B)]
      LL = neighbour(L, B, phi * w * G, pi * x, h, P, C, true)
      keep, cont = callback(result, LL)
      if keep
        push!(result, LL)
      end
      if !cont || lenght(result) >= max
        break
      end
    end
  elseif special
    pi = uniformizer(P)
    _G = T * form * _map(T, a)
    G = _map(_G, h)
    for i in 1:length(LO)
      w = LO[i]
      Gw = G * matrix(k, length(w), 1, w)
      ok = 0
      for d in 1:n
        if !iszero(Gw[d])
          ok = d
          break
        end
      end
      @assert ok != 0
      x = [ sum(B[i, j] * (h\w[i]) for i in 1:n) for j in 1:ncols(B)]
      nrm = _inner_product(form, x, x, a)
      b = h\(-h(nrm) / (2*Gw[d, 1]))
      x = x + b * pi * B[d]
      nrm = _inner_product(form, x, x, a)
      LL = neighbour(L, B, w * G, x, h, P, P, false)
      keep, cont = callback(result, LL)
      if keep
        push!(result, LL)
      end
      if !cont || lenght(result) >= max
        break
      end
    end
  else
    _G = T * form * _map(T, a)
    G = _map(_G, h)
    ram = isramified(R, minimum(P))
    if ram
      pi = uniformizer(P)
      S = [ h\x for x in k ]
    else
      p = minimum(P)
      pi = uniformizer(p)
      kp, hp = ResidueField(order(p), p)
      alpha = h\(degree(k) == 1 ? one(k) : gen(k))
      T = matrix(kp, 2, 1, [2, hp(tr(alpha))])
    end
    for i in 1:length(LO)
      w = LO[i]
      x = [ sum(B[i, j] * (h\w[i]) for i in 1:n) for j in 1:ncols(B)]
      Gw = G * matrix(k, length(w), 1, w)
      ok = 0
      for j in 1:n
        if !iszero(Gw[j])
          ok = j
          break
        end
      end
      @assert ok != 0
      NL = []
      if !ram
        s, V = solution(T, vector(kp, [hp(-nrm//pi)]))
        l = a(h\(1//wG[j]))
        S = [ l * (h\((s + v)[1]) + (h\(s + v)[2])*alpha) for v in v ]
      end
      for s in S
        LL = neighbour(L, B, wG, x+pi*s*B[j], h, P, P, false);
      end
      keep, cont = callback(result, LL)
      if keep
        push!(result, LL)
      end
      if !cont || lenght(result) >= max
        break
      end
    end
  end
  return result
end

function neighbours(L, P, max = inf)
#{The immediate P-neighbours of L}
#  require Order(P) eq BaseRing(L): "Arguments are incompatible";
#  require IsPrime(P): "Second argument must be prime";
#  require not IsRamified(P) or Minimum(P) ne 2: "Second argument cannot be a ramified prime over 2";
#  require IsModular(L, P) : "The lattice must be locally modular";
#  require Rank(L) ge 2: "The rank of the lattice must be at least 2";
#  require IsIsotropic(L, P): "The lattice must be locally isotropic";
#
  return _neighbours(L, P, [], max)
end

function stdcallback(list, L)
  keep = all(LL -> isisometric(LL, L), list)
#//  keep, #List;
  return keep, true;
end

function iterated_neighbours(L::HermLat, P; max = inf)# AutoOrbits, CallBack:= false, Max:= Infinity()) -> []
  #require Order(P) eq BaseRing(L): "Arguments are incompatible";
  #require IsPrime(P): "Second argument must be prime";
  #require not IsRamified(P) or Minimum(P) ne 2: "Second argument cannot be a ramified prime over 2";
  #require IsModular(L, P) : "The lattice must be locally modular";
  #require Rank(L) ge 2: "The rank of the lattice must be at least 2";
  #require IsIsotropic(L, P): "The lattice must be locally isotropic";

  result = [ L ]
  i = 1
  while length(result) < max && i <= length(result)
    # _Neighbours and the callback only add new lattices if not isometric to known ones and stop if Max is reached.
    # So we have nothing to at all.
    result = _neighbours(result[i], P, result, max)# : CallBack:= CallBack);
    i = i + 1
  end
  return result
end

################################################################################
#
#  First attempt, which mirrors Markus' Magma code
#
#  This is only used for debugging purposes
#
################################################################################

mutable struct LocalGenusSymbol{S}
  P
  data
  x
  iseven::Bool
  E
  isramified
  non_norm
end

prime(G::LocalGenusSymbol) = G.P

uniformizer(G::LocalGenusSymbol{QuadLat}) = G.x

data(G::LocalGenusSymbol) = G.data

base_field(G::LocalGenusSymbol) = G.E

function Base.show(io::IO, G::LocalGenusSymbol)
  print(io, "Local genus symbol at\n")
  print(IOContext(io, :compact => true), G.P)
  compact = get(io, :compact, false)
  if !compact
    print(io, "\nwith base field\n")
    print(io, base_field(G))
  end
  println(io, "\nWith data ", data(G))
  !G.iseven ? println(io, "and unifomizer\n", G.x) : nothing
end

# TODO: I have to redo this
function _genus_symbol_kirschmer(L::QuadLat, p::NfOrdIdl; uniformizer = zero(order(p)))
  O = order(p)
  nf(O) != base_field(L) && throw(error("Prime ideal must be an ideal of the base field of the lattice"))
  # If you pull from cache, you might have to adjust the symbol according
  # to the uniformizer flag

  J, G, E = jordan_decomposition(L, p)
  if !iszero(uniformizer)
    unif = uniformizer
    if valuation(unif, p) != 1
      throw(error("Wrong uniformizer"))
    end
  else
    unif = elem_in_nf(Hecke.uniformizer(p))
  end

  if minimum(p) != 2
    _, _h = ResidueField(O, p)
    h = extend(_h, nf(O))
    Gs = [ h(prod(diagonal(G[i]))//unif^(E[i] * nrows(J[i]))) for i in 1:length(J)]
    @assert !(0 in Gs)
    x  = [ (nrows(J[i]), E[i], issquare(Gs[i])[1] ? 1 : -1) for i in 1:length(J)]
    return LocalGenusSymbol{QuadLat}(p, x, unif, false, base_field(L), nothing, nothing)
  else
    t = length(G)
    sL = [ minimum(iszero(g[i, j]) ? inf : valuation(g[i, j], p) for j in 1:ncols(g) for i in 1:j) for g in G]
    @assert sL == E
    e = ramification_index(p)
    aL = []
    uL = []
    wL = []
    for i in 1:t
      GG = diagonal_matrix([ j < i ? unif^(2*(sL[i] - sL[j])) * G[j] : G[j] for j in 1:t])
      D = diagonal(GG)
      m, pos = findmin([valuation(d, p) for d in D])
      if e + sL[i] <= m
        push!(aL, unif^(e + sL[i]))
      else
        push!(aL, D[pos])
      end
      push!(uL, valuation(aL[i], p))
      push!(wL, min(e + sL[i], minimum(uL[i] + quadratic_defect(d//aL[i], p) for d in D)))
    end
    fL = []
    for k in 1:(t - 1)
      exp = uL[k] + uL[k + 1]
      push!(fL, (isodd(exp) ? exp : min(quadratic_defect(aL[k] * aL[k + 1], p), uL[k] + wL[k + 1], uL[k + 1], wL[k], e + div(exp, 2) + sL[k])) - 2*sL[k])
    end
    return LocalGenusSymbol{QuadLat}(p, ([nrows(G[k]) for k in 1:t], sL, wL, aL, fL, G), unif, true, base_field(L), nothing, nothing)
  end
end

# This is the "Magma" Genus symbol
function _genus_symbol_kirschmer(L::HermLat, p; uniformizer = zero(order(p)))
  @assert order(p) == base_ring(base_ring(L))

  B, G, S = jordan_decomposition(L, p)
  R = base_ring(L)
  E = nf(R)
  K = base_field(E)
  if !isdyadic(p) || !isramified(R, p)
    sym = [ (nrows(B[i]), S[i], islocal_norm(E, coeff(det(G[i]), 0), p)) for i in 1:length(B)]
  else
    P = prime_decomposition(R, p)[1][1]
    pi = E(K(Hecke.uniformizer(p)))
    sym = []
    for i in 1:length(B)
      normal = _get_norm_valuation_from_gram_matrix(G[i], P) == S[i]
      GG = diagonal_matrix([pi^(max(0, S[i] - S[j])) * G[j] for j in 1:length(B)])
      v = _get_norm_valuation_from_gram_matrix(GG, P)
      @assert v == valuation(R * norm(lattice(hermitian_space(E, GG), identity_matrix(E, nrows(GG)))), P)
      s = (nrows(B[i]), S[i], normal, v, coeff(det(diagonal_matrix([G[j] for j in 1:i])), 0))
      push!(sym, s)
    end
  end
  return sym
end


text:
#
# PBE correlation (without LDA part) - spin-polarized
# iflag = 1: J.P.Perdew, K.Burke, M.Ernzerhof, PRL 77, 3865 (1996).
#
function XC_c_pbe_spin( rho, zeta, grho )

    ga = 0.031091
    be = 0.06672455060314922
    third = 1.0/3.0
    pi34 = 0.6203504908994

    xkf = 1.919158292677513
    xks = 1.128379167095513

    rs = pi34/rho^third

    ec, vcup, vcdw = XC_c_pw_spin( rho, zeta )

    kf = xkf / rs
    ks = xks * sqrt(kf)
    fz = 0.5 * ( (1.0 + zeta)^(2.0 / 3.0) + (1.0 - zeta)^(2.0 / 3.0) )

    fz2 = fz * fz
    fz3 = fz2 * fz
    fz4 = fz3 * fz
    dfz = ( (1.0 + zeta)^(-1.0/3.0) - (1.0 - zeta)^( -1.0 / 3.0) ) / 3.0

    t = sqrt(grho) / (2.0 * fz * ks * rho)
    expe = exp( - ec / (fz3 * ga) )
    af = be / ga * (1.0 / (expe-1.0) )

    bfup = expe * (vcup - ec) / fz3
    bfdw = expe * (vcdw - ec) / fz3

    y = af * t * t

    xy = (1.0 + y) / (1.0 + y + y * y)
    qy = y * y * (2.0 + y) / (1.0 + y + y * y)^2

    s1 = 1.0 + be / ga * t * t * xy
    h0 = fz3 * ga * log(s1)

    dh0up = be * t * t * fz3 / s1 * ( -7.0 / 3.0 * xy - qy * (af * bfup / be - 7.0 / 3.0) )

    dh0dw = be * t * t * fz3 / s1 * ( -7.0 / 3.0 * xy - qy * (af * bfdw / be - 7.0 / 3.0) )

    dh0zup = (3.0 * h0 / fz - be * t * t * fz2 / s1 * (2.0 * xy -
             qy * (3.0 * af * expe * ec / fz3 / be + 2.0) ) ) * dfz * (1.0 - zeta)

    dh0zdw = -(3.0 * h0 / fz - be * t * t * fz2 / s1 * (2.0 * xy -
             qy * (3.0 * af * expe * ec / fz3 / be + 2.0) ) ) * dfz * (1.0 + zeta)

    ddh0 = be * fz / (2.0 * ks * ks * rho) * (xy - qy) / s1

    sc = rho * h0
    v1cup = h0 + dh0up + dh0zup
    v1cdw = h0 + dh0dw + dh0zdw
    v2c = ddh0

    return sc, v1cup, v1cdw, v2c

end

ext: jl


text:
import Printf: @sprintf

"""Returns a human readable string representing seconds in terms of days, hours, minutes or seconds."""
function readable_secs(secs::Real)
    days = Int(floor(secs/3600/24))
    hours = Int(floor((secs/3600) % 24))
    minutes = Int(floor((secs/60) % 60))
    seconds = Int(floor(secs%3600%60))
    secs1f = @sprintf "%.1fs" secs%3600%60
    secs2f = @sprintf "%.2fs" secs%3600%60

    if days > 0
        return "$(days)d, $(hours)h"
    elseif hours > 0
        return "$(hours)h, $(minutes)min"
    elseif minutes > 0
        return "$(minutes)min, $(seconds)s"
    elseif seconds > 10
        return secs1f
    else
        return secs2f
    end
end

max_stars_repo_name: milankl/BernoulliMap

max_issues_repo_name: milankl/BernoulliMap


text:
oddturns(w,d) = Int(floor(round(w/d;digits=6)))
eventurns(w,d) = Int(floor(round((w-d/2)/d;digits=6)))

oddwidth(n::Integer,d::Unitful.Length) = n*d
evenwidth(n::Integer,d::Unitful.Length) = (n+0.5)*d
function layerheight(n::Integer,d::Unitful.Length)
    @assert n>0
    y = d
    if n>1
        y += (n-1)*(d/2)*√3
    end
    y
end

"""
Calculate the orthocyclic coil winding for a give wire diameter, returning a `WindingLayers`.

ideal_fill(w, h, d)

    w => Horizontal width of coil for layers
    h => Vertical height to fill with layers
    d => Diameter of Wire

    ←           w             →
    ###########################
    #           ↑             #
    #                         #
    #           h             #
    #                         #
    #           ↓             #
    ###########################
                ↑
              (id)
"""
function ideal_fill(w,h,d)

    # layers=length(d/2:d*sqrt(3)/2:h-d/2)
    #TODO Improve this to the rounding isn't necessary
    layers = ideal_layers(h,d)
    odd_layer= oddturns(w,d)
    even_layer= eventurns(w,d)

    odd_turns = length(1:2:layers)*odd_layer
    even_turns = length(2:2:layers)*even_layer

    WindingLayers(layers, odd_layer, even_layer, odd_turns, even_turns, odd_turns+even_turns)
end

"""
ideal_fill(coil::CoilGeometry, awg::AWG, [fill::Real=1.0])

    coil => CoilGeometry
    awg => AWG of the wire
    fillfactor => fillfactor of the wire
"""
function ideal_fill(coil::CoilGeometry, awg::AWG, fillfactor::T=1.0) where {T<:Real}
    zero(T) < fillfactor <= one(T) || throw(DomainError(fillfactor))
    w=coil.len
    h=(coil.od-coil.id)/2
    ideal_fill(w,h,awg.d_insul/sqrt(fillfactor))
end

"""
Calculate the number of layers of wire.  Given the `h` height of the layter, and the `d` diameter of the wire.

    ideal_layers(h, d)

    ###########################
    #           ↑             #
    #                         #
    #           h             #
    #                         #
    #           ↓             #
    ###########################
                ↑
              (id)
"""
ideal_layers(h,d) = h≳d ? 1+floor(round( (h-d)/(d*√3/2) ;digits=6)) : 0
ideal_layers(coil::CoilGeometry,d::Unitful.Length) = ideal_layers((coil.od-coil.id)/2,d)

"""
Construct the minimum size coil to fit the specified winding

    enclosewinding(layers, oddlayer, evenlayer, id::Unitful.Length, d::Unitful.Length)
    enclosewinding(layers, oddlayer, evenlayer, id::Unitful.Length, gauge::AWG)
"""
function enclosewinding(layers, oddlayer, evenlayer, id::Unitful.Length, d::Unitful.Length)
    @assert layers>0
    @assert evenlayer == oddlayer || evenlayer == oddlayer-1
    # od = id + 2d + (layers-1)*d*√3
    # len = oddlayer > evenlayer ? d*oddlayer : d/2 + d*evenlayer
    w=oddlayer>evenlayer ? oddwidth(oddlayer,d) : evenwidth(evenlayer,d)
    od=2layerheight(layers,d) + id
    f = unit(id)
    CoilGeometry(id,f(od),f(w))
end
enclosewinding(layers, oddlayer, evenlayer, id::Unitful.Length, gauge::AWG) =
    enclosewinding(layers, oddlayer, evenlayer, id, gauge.d_insul)

"""
Construct the minimum size `coil` to fit the `winding` within the limits of the existing coil.

enclosewinding(coil::CoilGeometry, winding::WindingLayers, d::Unitful.Length)
enclosewinding(coil::CoilGeometry, winding::WindingLayers, gauge::AWG)
"""
function enclosewinding(coil::CoilGeometry, winding::WindingLayers, d::Unitful.Length)
#TODO Set check to make sure the calculated winding is within the limits of the existing coil
    new_coil = enclosewinding(winding.layers, winding.odd_layers, winding.even_layers, coil.od, d)
    @assert new_coil.id >= coil.id && new_coil.od <= coil.od && new_coil.len <= coil.len
    @assert isvalidcoil(new_coil)
    new_coil
end
enclosewinding(coil::CoilGeometry, winding::WindingLayers, gauge::AWG) = enclosewinding(coil,winding,gauge.d_insul)


text:
function wordcount(sentence)

    function isnon(c)
        return !isletter(c) & !isdigit(c)
    end
    tmp = lowercase(sentence)
    #tmp = replace(tmp, "' " => "")
    #tmp = replace(tmp, " '" => "")
    tmp = replace(tmp, "'" => "T")
    words = split(tmp,[',',' ','\t','\n'])
    res = Dict{String,Int32}()
    for w in words
        tst = replace(w, isnon => "")
        if last(tst,1)=="T"
            tst=chop(tst)
        end
        if first(tst,1) =="T"
            tst=reverse(chop(reverse(tst)))
        end

        tst = replace(tst, "T" => "'")

        if tst != "T" && length(tst)>0
            #println("out $tst")
            tmp = get!(res, tst, 0)
            push!(res, tst => tmp + 1)
        end
    end
    return res
end
function wordcounts(sentence)
    prev = ""
    prevclass = 0
    res = Dict{String,Int32}()

    for c in lowercase(sentence)

        if isspace(c) || c == ','
            # separate word
            if length(prev) > 0
                tmp = get!(res, prev, 0)
                push!(res, prev => tmp + 1)
                prev = ""
            end

            # elseif isnumeric(c)
            #     print(c * " is numeric")
            #     if prevclass == 0
            #         # we shall add previous word

            #     elseif prevclass == 2
            #         # push!(res, "sdf" => [1, 2])
            #         res[prev] = res[prev] + 1
            #         prev = c
            #         prevclass = 1
            #     end

        elseif isletter(c) || isnumeric(c) || (c == ''' && length(prev) > 0)
            print(c * " is ascii")
            if prevclass == 0
                prev = prev * c
            elseif prevclass == 1
                throw("not implemented")
            end

        end
    end
    # adding remaining word to dictionary
    if length(prev) > 0
        tmp = get!(res, prev, 0)
        push!(res, prev => tmp + 1)
        # res[prev] = res[prev] + 1
    end

    return res
end

max_stars_repo_name: stepinski/julia-init

max_issues_repo_name: stepinski/julia-init

max_forks_repo_name: stepinski/julia-init


text:
using FixedEffects, Random, Statistics
Random.seed!(1234)
N = 10000000
K = 100
id1 = rand(1:div(N, K), N)
id2 = rand(1:K, N)
fes = [FixedEffect(id1), FixedEffect(id2)]
x = rand(N)

# simple problem
X = deepcopy(x)
@time solve_residuals!(X, fes; tol = 1e-6)
#  0.807994 seconds (848 allocations: 311.411 MiB, 10.14% gc time)
X = [x x x x x x x x x x]
@time solve_residuals!(X, fes)
#  6.610115 seconds (8.83 k allocations: 311.572 MiB, 1.54% gc time)
X = [x x x x x x x x x x]
@time solve_residuals!(X,fes, method = :lsmr_threads)
#  5.120222 seconds (9.19 k allocations: 1.210 GiB, 8.59% gc time)

# More complicated problem
N = 800000 # number of observations
M = 400000 # number of workers
O = 50000 # number of firms
Random.seed!(1234)
pid = rand(1:M, N)
fid = [rand(max(1, div(x, 8)-10):min(O, div(x, 8)+10)) for x in pid]
x = rand(N)
fes = [FixedEffect(pid), FixedEffect(fid)]
@time solve_residuals!([x x x x], fes; double_precision = false)
# 12.690011 seconds (191.26 k allocations: 67.887 MiB, 0.03% gc time)
@time solve_residuals!([x x x x], fes, method = :lsmr_threads, maxiter = 300)
# 10.660428 seconds (239.89 k allocations: 189.271 MiB, 0.80% gc time)

# cluster (381 iterations)
@time solve_residuals!([x x x x], fes; double_precision = false)
#  11.760708 seconds (292.08 k allocations: 51.847 MiB, 0.95% gc time)
@time solve_residuals!([x x x x], fes, method = :lsmr_gpu, double_precision = false)
#  2.342490 seconds (1.15 M allocations: 90.584 MiB, 5.05% gc time)


text:
function _right_sweep_SVD!(ψ::CuMPS, Dcut::Int=typemax(Int))
    Σ = V = CUDA.ones(eltype(ψ), 1, 1)

    for i ∈ eachindex(ψ)
        A = ψ[i]
        C = Diagonal(Σ) * V'

        # attach
        @cutensor M[x, σ, y] := C[x, α] * A[α, σ, y]
        @cast   M̃[(x, σ), y] |= M[x, σ, y]

        # decompose
        U, Σ, V = _truncate_svd(CUDA.svd(M̃), Dcut)

        # create new
        d = physical_dim(ψ, i)
        @cast A[x, σ, y] |= U[(x, σ), y] (σ:d)
        ψ[i] = A
    end
end

function _left_sweep_SVD!(ψ::CuMPS, Dcut::Int=typemax(Int))
    Σ = U = CUDA.ones(eltype(ψ), 1, 1)

    for i ∈ length(ψ):-1:1
        B = ψ[i]
        C = U * Diagonal(Σ)

        # attach
        @cutensor M[x, σ, y]   := B[x, σ, α] * C[α, y]
        @cast   M̃[x, (σ, y)] |= M[x, σ, y]

        # decompose
        U, Σ, V = _truncate_svd(CUDA.svd(M̃), Dcut)

        # create new
        d = physical_dim(ψ, i)
        VV = conj.(transpose(V)) #hack - @cast does fail with allowscalar and Adjoint type
        @cast B[x, σ, y] |= VV[x, (σ, y)] (σ:d)
        ψ[i] = B
    end
end

function _left_sweep_var!!(ϕ::CuMPS, env::Vector{<:CuMatrix}, ψ::CuMPS, Dcut::Int)
    S = eltype(ϕ)

    # overwrite the overlap
    env[end] = CUDA.ones(S, 1, 1)

    for i ∈ length(ψ):-1:1
        L = env[i]
        R = env[i+1]

        # optimize site
        M = ψ[i]
        @cutensor M̃[x, σ, y] := L[x, β] * M[β, σ, α] * R[α, y] order = (α, β)

        # right canonize it
        @cast MM[x, (σ, y)] |= M̃[x, σ, y]
        Q = _truncate_qr(CUDA.qr(MM'), Dcut)'

        Q = CuMatrix(Q) #hack @cast does not work with Adjoint
        d = size(M, 2)
        @cast B[x, σ, y] |= Q[x, (σ, y)] (σ:d)

        # update ϕ and right environment
        ϕ[i] = B
        A = ψ[i]

        @cutensor RR[x, y] := A[x, σ, α] * R[α, β] * conj(B[y, σ, β]) order = (β, α, σ)
        env[i] = RR
    end
end

function _right_sweep_var!!(ϕ::CuMPS, env::Vector{<:CuMatrix}, ψ::CuMPS, Dcut::Int)
    S = eltype(ϕ)

    # overwrite the overlap
    env[1] = CUDA.ones(S, 1, 1)

    for i ∈ eachindex(ψ)
        L = env[i]
        R = env[i+1]

        # optimize site
        M = ψ[i]
        @cutensor M̃[x, σ, y] := L[x, β] * M[β, σ, α] * R[α, y] order = (α, β)

        # left canonize it
        @cast B[(x, σ), y] |= M̃[x, σ, y]
        Q = _truncate_qr(CUDA.qr(B), Dcut)

        d = size(ϕ[i], 2)
        @cast A[x, σ, y] |= Q[(x, σ), y] (σ:d)

        # update ϕ and left environment
        ϕ[i] = A
        B = ψ[i]

        @cutensor LL[x, y] := conj(A[β, σ, x]) * L[β, α] * B[α, σ, y] order = (α, β, σ)
        env[i+1] = LL
    end
    return real(tr(env[end]))
end

function _qr_fix!(Q::CuMatrix, R::CuMatrix)
    d = diag(R)
    ph = d./abs.(d)
    transpose(ph) .* Q
end


text:
struct Dummy{T} <: AbstractRegion{T}
end

@testset "AbstractRegion" begin
    @test_throws ErrorException SVector(0.0, 0.0, 0.0) ∈ Dummy{Float64}()
    @test_throws ErrorException boundingbox(Dummy{Float64}())
end


text:
using Test
using ChemistryFeaturization
using SimpleWeightedGraphs

include("../src/layers.jl")
using .layers: AGNConv, AGNPool

@testset "AGNConv" begin
    # create simple line graph, populate it with feature of all ones
    adjmat = Float32.([0 1 0; 1 0 1; 0 1 0])
    dummyfzn = [AtomFeat(:feat, [0]) for i in 1:4]
    ag = AtomGraph(SimpleWeightedGraph{Int32, Float32}(adjmat), ["C", "C", "C"], ones(Float32, 4,3), dummyfzn)
    # create a conv layer, initialize weights with ones
    l = AGNConv(4=>4, initW=ones)

    # test output looks as expected
    output_fea = l(ag).features
    @test output_fea[:,1]==output_fea[:,3]
    @test isapprox(output_fea[:,1].+output_fea[:,3], .-output_fea[:,2])

    # and now for a loop
    adjmat = Float32.([0 1 1; 1 0 1; 1 1 0])
    ag = AtomGraph(SimpleWeightedGraph{Int32, Float32}(adjmat), ["C", "C", "C"], ones(Float32, 4,3), dummyfzn)
    l = AGNConv(4=>4, initW=ones)

    @test all(isapprox.(l(ag).features, zero(Float32)))
end

@testset "pooling" begin
    # keep our little line graph, but give it more features
    adjmat = Float32.([0 1 0; 1 0 1; 0 1 0])
    feat = ones(Float32, 50,3)
    dummyfzn = [AtomFeat(:feat, [0]) for i in 1:50]
    ag = AtomGraph(SimpleWeightedGraph{Int32, Float32}(adjmat), ["C", "C", "C"], feat, dummyfzn)

    # make some pooling layers
    meanpool = AGNPool("mean", 50, 10, 0.1)
    maxpool = AGNPool("max", 50, 10, 0.1)

    # start with the easy stuff
    @test meanpool(ag) == ones(Float32, 10, 1)
    @test maxpool(ag) == ones(Float32, 10, 1)

    # one level up
    feat[:,2] .= 0.0
    add_features!(ag, feat, dummyfzn)
    # they're still the same here because maxpool takes max along features but averages across nodes right now
    @test all(isapprox.(meanpool(ag), 2/3))
    @test all(isapprox.(maxpool(ag), 2/3))

    # and now make them different
    feat[1:5:50,2] .= 4.0
    add_features!(ag, feat, dummyfzn)
    @test all(isapprox.(meanpool(ag), 14/15))
    @test all(maxpool(ag) .== 2.0)

    # make sure it complains when it should
    @test_throws AssertionError AGNPool("mean", 10, 20, 0.5)
end


text:
"""
    Plot(graphics; kw...)

# Description

Decoration for objects that are `GraphicsArea` (or `Canvas`).
It is used to surround the inner `GraphicsArea` object with
additional information such as a title, border, and axis labels.

# Usage

    Plot(graphics; $(keywords(; default = (), add = (:title, :xlabel, :ylabel, :zlabel, :border, :margin, :padding, :compact, :labels))))

    Plot(x, y, canvas; $(keywords()))

# Arguments

$(arguments(
    (; graphics = "the `GraphicsArea` (e.g. a subtype of `Canvas`) that the plot should decorate");
    add = (:x, :y, :canvas)
))

# Methods

- `title!(plot::Plot, title::String)`

- `xlabel!(plot::Plot, xlabel::String)`

- `ylabel!(plot::Plot, xlabel::String)`

- `zlabel!(plot::Plot, zlabel::String)`

- `label!(plot::Plot, where::Symbol, value::String)`

- `label!(plot::Plot, where::Symbol, row::Int, value::String)`

Author(s)

- Christof Stocker (github.com/Evizero)

# See also

[`scatterplot`](@ref), [`lineplot`](@ref),
[`BarplotGraphics`](@ref), [`BrailleCanvas`](@ref),
[`BlockCanvas`](@ref), [`AsciiCanvas`](@ref)
"""
mutable struct Plot{T<:GraphicsArea}
    graphics::T
    title::String
    xlabel::String
    ylabel::String
    zlabel::String
    margin::Int
    padding::Int
    border::Symbol
    compact::Bool
    labels_left::Dict{Int,String}
    colors_left::Dict{Int,JuliaColorType}
    labels_right::Dict{Int,String}
    colors_right::Dict{Int,JuliaColorType}
    decorations::Dict{Symbol,String}
    colors_deco::Dict{Symbol,JuliaColorType}
    labels::Bool
    colormap::Any
    colorbar::Bool
    colorbar_border::Symbol
    colorbar_lim::Tuple{Number,Number}
    autocolor::Int
end

function Plot(
    graphics::T;
    title::AbstractString = KEYWORDS.title,
    xlabel::AbstractString = KEYWORDS.xlabel,
    ylabel::AbstractString = KEYWORDS.ylabel,
    zlabel::AbstractString = KEYWORDS.zlabel,
    border::Symbol = KEYWORDS.border,
    compact::Bool = KEYWORDS.compact,
    margin::Int = KEYWORDS.margin,
    padding::Int = KEYWORDS.padding,
    labels::Bool = KEYWORDS.labels,
    colorbar::Bool = KEYWORDS.colorbar,
    colorbar_border::Symbol = KEYWORDS.colorbar_border,
    colorbar_lim = KEYWORDS.colorbar_lim,
    colormap::Any = nothing,
    ignored...,
) where {T<:GraphicsArea}
    margin >= 0 || throw(ArgumentError("Margin must be greater than or equal to 0"))
    rows = nrows(graphics)
    cols = ncols(graphics)
    labels_left = Dict{Int,String}()
    colors_left = Dict{Int,JuliaColorType}()
    labels_right = Dict{Int,String}()
    colors_right = Dict{Int,JuliaColorType}()
    decorations = Dict{Symbol,String}()
    colors_deco = Dict{Symbol,JuliaColorType}()
    p = Plot{T}(
        graphics,
        title,
        xlabel,
        ylabel,
        zlabel,
        margin,
        padding,
        border,
        compact,
        labels_left,
        colors_left,
        labels_right,
        colors_right,
        decorations,
        colors_deco,
        labels && graphics.visible,
        colormap,
        colorbar,
        colorbar_border,
        colorbar_lim,
        0,
    )
    if compact
        xlabel != "" && label!(p, :b, xlabel)
        ylabel != "" && label!(p, :l, round(Int, nrows(graphics) / 2), ylabel)
    end
    p
end

function Plot(
    X::AbstractVector{<:Number},
    Y::AbstractVector{<:Number},
    ::Type{C} = BrailleCanvas;
    title::AbstractString = KEYWORDS.title,
    xlabel::AbstractString = KEYWORDS.xlabel,
    ylabel::AbstractString = KEYWORDS.ylabel,
    zlabel::AbstractString = KEYWORDS.zlabel,
    xscale::Union{Symbol,Function} = KEYWORDS.xscale,
    yscale::Union{Symbol,Function} = KEYWORDS.yscale,
    width::Union{Int,Nothing} = nothing,
    height::Union{Int,Nothing} = nothing,
    border::Symbol = KEYWORDS.border,
    compact::Bool = KEYWORDS.compact,
    blend::Bool = KEYWORDS.blend,
    xlim = KEYWORDS.xlim,
    ylim = KEYWORDS.ylim,
    margin::Int = KEYWORDS.margin,
    padding::Int = KEYWORDS.padding,
    labels::Bool = KEYWORDS.labels,
    unicode_exponent::Bool = KEYWORDS.unicode_exponent,
    colorbar::Bool = KEYWORDS.colorbar,
    colorbar_border::Symbol = KEYWORDS.colorbar_border,
    colorbar_lim = KEYWORDS.colorbar_lim,
    colormap::Any = nothing,
    grid::Bool = KEYWORDS.grid,
    min_width::Int = 5,
    min_height::Int = 2,
) where {C<:Canvas}
    length(xlim) == length(ylim) == 2 ||
        throw(ArgumentError("xlim and ylim must be tuples or vectors of length 2"))
    length(X) == length(Y) || throw(DimensionMismatch("X and Y must have same length"))

    width === nothing && (width = DEFAULT_WIDTH[])
    height === nothing && (height = DEFAULT_HEIGHT[])

    (visible = width > 0) && (width = max(width, min_width))
    height = max(height, min_height)

    min_x, max_x = extend_limits(X, xlim, xscale)
    min_y, max_y = extend_limits(Y, ylim, yscale)

    p_width = max_x - min_x
    p_height = max_y - min_y

    canvas = C(
        width,
        height,
        blend = blend,
        visible = visible,
        origin_x = min_x,
        origin_y = min_y,
        width = p_width,
        height = p_height,
        xscale = xscale,
        yscale = yscale,
    )
    plot = Plot(
        canvas,
        title = title,
        margin = margin,
        padding = padding,
        border = border,
        compact = compact,
        labels = labels,
        xlabel = xlabel,
        ylabel = ylabel,
        zlabel = zlabel,
        colormap = colormap,
        colorbar = colorbar,
        colorbar_border = colorbar_border,
        colorbar_lim = colorbar_lim,
    )
    base_x = xscale isa Symbol ? get(BASES, xscale, nothing) : nothing
    base_y = yscale isa Symbol ? get(BASES, yscale, nothing) : nothing
    m_x, M_x, m_y, M_y = map(
        v -> compact_repr(roundable(v) ? round(Int, v, RoundNearestTiesUp) : v),
        (min_x, max_x, min_y, max_y),
    )
    if unicode_exponent
        m_x, M_x = map(v -> base_x !== nothing ? superscript(v) : v, (m_x, M_x))
        m_y, M_y = map(v -> base_y !== nothing ? superscript(v) : v, (m_y, M_y))
    end
    base_x_str = base_x === nothing ? "" : base_x * (unicode_exponent ? "" : "^")
    base_y_str = base_y === nothing ? "" : base_y * (unicode_exponent ? "" : "^")
    label!(plot, :l, nrows(canvas), base_y_str * m_y, color = :light_black)
    label!(plot, :l, 1, base_y_str * M_y, color = :light_black)
    label!(plot, :bl, base_x_str * m_x, color = :light_black)
    label!(plot, :br, base_x_str * M_x, color = :light_black)
    if grid
        if min_y < 0 < max_y
            for i in range(
                min_x,
                stop = max_x,
                length = width * x_pixel_per_char(typeof(canvas)),
            )
                points!(plot, i, 0.0, :normal)
            end
        end
        if min_x < 0 < max_x
            for i in range(
                min_y,
                stop = max_y,
                length = height * y_pixel_per_char(typeof(canvas)),
            )
                points!(plot, 0.0, i, :normal)
            end
        end
    end
    plot
end

function next_color!(plot::Plot{<:GraphicsArea})
    cur_color = COLOR_CYCLE[plot.autocolor + 1]
    plot.autocolor = ((plot.autocolor + 1) % length(COLOR_CYCLE))
    cur_color
end

"""
    title(plot) -> String

Returns the current title of the given plot.
Alternatively, the title can be changed with `title!`.
"""
title(plot::Plot) = plot.title

"""
    title!(plot, newtitle)

Sets a new title for the given plot.
Alternatively, the current title can be
queried using `title`.
"""
function title!(plot::Plot, title::AbstractString)
    plot.title = title
    plot
end

"""
    xlabel(plot) -> String

Returns the current label for the x-axis.
Alternatively, the x-label can be changed with `xlabel!`
"""
xlabel(plot::Plot) = plot.xlabel

"""
    xlabel!(plot, newlabel)

Sets a new x-label for the given plot.
Alternatively, the current label can be
queried using `xlabel`
"""
function xlabel!(plot::Plot, xlabel::AbstractString)
    plot.xlabel = xlabel
    plot
end

"""
    ylabel(plot) -> String

Returns the current label for the y-axis.
Alternatively, the y-label can be changed with `ylabel!`
"""
ylabel(plot::Plot) = plot.ylabel

"""
    ylabel!(plot, newlabel)

Sets a new y-label for the given plot.
Alternatively, the current label can be
queried using `ylabel`
"""
function ylabel!(plot::Plot, ylabel::AbstractString)
    plot.ylabel = ylabel
    plot
end

"""
    zlabel(plot) -> String

Returns the current label for the z-axis (colorbar).
Alternatively, the z-label can be changed with `zlabel!`
"""
zlabel(plot::Plot) = plot.zlabel

"""
    zlabel!(plot, newlabel)

Sets a new z-label (colorbar label) for the given plot.
Alternatively, the current label can be
queried using `zlabel`
"""
function zlabel!(plot::Plot, zlabel::AbstractString)
    plot.zlabel = zlabel
    plot
end

"""
    label!(plot, where, value, [color])

    label!(plot, where, row, value, [color])

This method is responsible for the setting all the textual decorations of a plot.

Note that `where` can be any of: `:tl` (top-left),
`:t` (top-center), `:tr` (top-right),
`:bl` (bottom-left), `:b` (bottom-center),
`:br` (bottom-right), `:l` (left), `:r` (right).

If `where` is either `:l`, or `:r`, then `row` can be between 1
and the number of character rows of the plots canvas.
"""
function label!(plot::Plot, loc::Symbol, value::AbstractString, color::UserColorType)
    loc ∉ (:t, :b, :l, :r, :tl, :tr, :bl, :br) &&
        throw(ArgumentError("Unknown location: try one of these :tl :t :tr :bl :b :br"))
    if loc == :l || loc == :r
        for row in 1:nrows(plot.graphics)
            if loc == :l
                if !haskey(plot.labels_left, row) || plot.labels_left[row] == ""
                    plot.labels_left[row] = value
                    plot.colors_left[row] = julia_color(color)
                    break
                end
            elseif loc == :r
                if !haskey(plot.labels_right, row) || plot.labels_right[row] == ""
                    plot.labels_right[row] = value
                    plot.colors_right[row] = julia_color(color)
                    break
                end
            end
        end
    else
        plot.decorations[loc] = value
        plot.colors_deco[loc] = julia_color(color)
    end
    plot
end

function annotate!(plot::Plot, loc::Symbol, value::AbstractString, color::UserColorType)
    Base.depwarn("`annotate!` has been renamed to `label!`", :Plot)
    label!(plot, loc, value, color)
end

label!(plot::Plot, loc::Symbol, value::AbstractString; color::UserColorType = :normal) =
    label!(plot, loc, value, color)

function annotate!(
    plot::Plot,
    loc::Symbol,
    value::AbstractString;
    color::UserColorType = :normal,
)
    Base.depwarn("`annotate!` has been renamed to `label!`", :Plot)
    label!(plot, loc, value, color)
end

function label!(
    plot::Plot,
    loc::Symbol,
    row::Int,
    value::AbstractString,
    color::UserColorType,
)
    if loc == :l
        plot.labels_left[row] = value
        plot.colors_left[row] = julia_color(color)
    elseif loc == :r
        plot.labels_right[row] = value
        plot.colors_right[row] = julia_color(color)
    else
        throw(ArgumentError("Unknown location \"$loc\", try `:l` or `:r` instead"))
    end
    plot
end

function annotate!(
    plot::Plot,
    loc::Symbol,
    row::Int,
    value::AbstractString,
    color::UserColorType,
)
    Base.depwarn("`annotate!` has been renamed to `label!`", :Plot)
    label!(plot, loc, row, value, color)
end

label!(
    plot::Plot,
    loc::Symbol,
    row::Int,
    value::AbstractString;
    color::UserColorType = :normal,
) = label!(plot, loc, row, value, color)

function annotate!(
    plot::Plot,
    loc::Symbol,
    row::Int,
    value::AbstractString;
    color::UserColorType = :normal,
)
    Base.depwarn("`annotate!` has been renamed to `label!`", :Plot)
    label!(plot, loc, row, value, color)
end

"""
    annotate!(plot, x, y, text; kw...)

# Description

Adds text to the plot at the position `(x, y)`.

# Arguments

$(arguments((; text = "a string of text"); default = (), add = (:x, :y, :color)))

# Examples

```julia-repl
julia> plt = lineplot([1, 2, 7], [9, -6, 8], title = "My Lineplot");

julia> annotate!(plt, 5, 5, "My text")
       ⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀My Lineplot
       ┌────────────────────────────────────────┐
    10 │⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
       │⢇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⡠│
       │⠘⡄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⠤⠊⠁⠀│
       │⠀⢣⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀My text⠀⠀⣀⠔⠊⠁⠀⠀⠀⠀│
       │⠀⠈⡆⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣀⠔⠊⠀⠀⠀⠀⠀⠀⠀⠀│
       │⠀⠀⢸⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡠⠔⠉⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
       │⠀⠀⠀⢇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⡠⠒⠉⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
       │⠤⠤⠤⠼⡤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⢤⠤⠶⠥⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤│
       │⠀⠀⠀⠀⢣⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⠤⠊⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
       │⠀⠀⠀⠀⠈⡆⠀⠀⠀⠀⠀⠀⠀⣀⠔⠊⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
       │⠀⠀⠀⠀⠀⢱⠀⠀⠀⠀⡠⠔⠊⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
       │⠀⠀⠀⠀⠀⠀⢇⡠⠔⠉⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
       │⠀⠀⠀⠀⠀⠀⠈⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
       │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
   -10 │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
       └────────────────────────────────────────┘
       ⠀1⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀7
```

# See also

[`Plot`](@ref), [`lineplot`](@ref), [`scatterplot`](@ref),
[`stairs`](@ref), [`BrailleCanvas`](@ref), [`BlockCanvas`](@ref),
[`AsciiCanvas`](@ref), [`DotCanvas`](@ref)
"""
function annotate!(
    plot::Plot{<:Canvas},
    x::Number,
    y::Number,
    text::Union{Char,AbstractString};
    color = :normal,
    kw...,
)
    color = color == :auto ? next_color!(plot) : color
    annotate!(plot.graphics, x, y, text, color; kw...)
    plot
end

function lines!(plot::Plot{<:Canvas}, args...; kw...)
    lines!(plot.graphics, args...; kw...)
    plot
end

function pixel!(plot::Plot{<:Canvas}, args...; kw...)
    pixel!(plot.graphics, args...; kw...)
    plot
end

function points!(plot::Plot{<:Canvas}, args...; kw...)
    points!(plot.graphics, args...; kw...)
    plot
end

function print_title(
    io::IO,
    left_pad::AbstractString,
    title::AbstractString,
    right_pad::AbstractString,
    blank::Char;
    p_width::Int = 0,
    color = :normal,
)
    title == "" && return
    offset = round(Int, p_width / 2 - length(title) / 2, RoundNearestTiesUp)
    pre_pad = repeat(blank, offset > 0 ? offset : 0)
    print(io, left_pad, pre_pad)
    print_color(color, io, title)
    post_pad = repeat(blank, max(0, p_width - length(pre_pad) - length(title)))
    print(io, post_pad, right_pad)
    nothing
end

function print_border(
    io::IO,
    loc::Symbol,
    length::Int,
    left_pad::AbstractString,
    right_pad::AbstractString,
    bmap = BORDERMAP[:solid],
    color::UserColorType = :light_black,
)
    print(io, left_pad)
    print_color(
        color,
        io,
        bmap[Symbol(loc, :l)],
        repeat(bmap[loc], length),
        bmap[Symbol(loc, :r)],
    )
    print(io, right_pad)
    nothing
end

_nocolor_string(str) = replace(string(str), r"\e\[[0-9]+m" => "")

function print_labels(
    io::IO,
    mloc::Symbol,
    p::Plot,
    border_length,
    left_pad::AbstractString,
    right_pad::AbstractString,
    blank::Char,
)
    p.labels || return
    lloc      = Symbol(mloc, :l)
    rloc      = Symbol(mloc, :r)
    left_str  = get(p.decorations, lloc, "")
    left_col  = get(p.colors_deco, lloc, :light_black)
    mid_str   = get(p.decorations, mloc, "")
    mid_col   = get(p.colors_deco, mloc, :light_black)
    right_str = get(p.decorations, rloc, "")
    right_col = get(p.colors_deco, rloc, :light_black)
    if left_str != "" || right_str != "" || mid_str != ""
        left_len  = length(left_str)
        mid_len   = length(mid_str)
        right_len = length(right_str)
        print(io, left_pad)
        print_color(left_col, io, left_str)
        cnt = round(Int, border_length / 2 - mid_len / 2 - left_len, RoundNearestTiesAway)
        pad = cnt > 0 ? repeat(blank, cnt) : ""
        print(io, pad)
        print_color(mid_col, io, mid_str)
        cnt = border_length - right_len - left_len - mid_len + 2 - cnt
        pad = cnt > 0 ? repeat(blank, cnt) : ""
        print(io, pad)
        print_color(right_col, io, right_str)
        print(io, right_pad)
    end
    nothing
end

function Base.show(io::IO, p::Plot)
    c = p.graphics
    🗷 = Char(BLANK)  # blank outside canvas
    🗹 = Char(c isa BrailleCanvas ? BLANK_BRAILLE : 🗷)  # blank inside canvas
    ############################################################
    # 🗷 = 'x'  # debug
    # 🗹 = Char(typeof(c) <: BrailleCanvas ? '⠿' : 'o')  # debug
    ############################################################
    border_length = ncols(c)
    p_width = border_length + 2  # left corner + border + right corner

    bmap = BORDERMAP[p.border === :none && c isa BrailleCanvas ? :bnone : p.border]

    # get length of largest strings to the left and right
    max_len_l = if p.labels && !isempty(p.labels_left)
        maximum([length(_nocolor_string(l)) for l in values(p.labels_left)])
    else
        0
    end
    max_len_r = if p.labels && !isempty(p.labels_right)
        maximum([length(_nocolor_string(l)) for l in values(p.labels_right)])
    else
        0
    end
    if !p.compact && p.labels && p.ylabel != ""
        max_len_l += length(p.ylabel) + 1
    end

    # offset where the plot (incl border) begins
    plot_offset = max_len_l + p.margin + p.padding

    # padding-string from left to border
    plot_padding = repeat(🗷, p.padding)

    if p.colorbar
        min_z, max_z = p.colorbar_lim
        min_z_str = string(isinteger(min_z) ? min_z : float_round_log10(min_z))
        max_z_str = string(isinteger(max_z) ? max_z : float_round_log10(max_z))
        cbar_max_len =
            max(length(min_z_str), length(max_z_str), length(_nocolor_string(p.zlabel)))
        cbar_pad = plot_padding * repeat(🗹, 4) * plot_padding * repeat(🗷, cbar_max_len)
    else
        cbar_pad = ""
    end

    # padding-string between labels and border
    border_left_pad = repeat(🗷, plot_offset)

    # trailing
    border_right_pad = repeat(🗷, max_len_r) * plot_padding * cbar_pad

    # plot the title and the top border
    print_title(
        io,
        border_left_pad,
        p.title,
        border_right_pad * '\n',
        🗹;
        p_width = p_width,
        color = :bold,
    )
    print_labels(
        io,
        :t,
        p,
        border_length - 2,
        border_left_pad * 🗹,
        🗹 * border_right_pad * '\n',
        🗹,
    )
    c.visible &&
        print_border(io, :t, border_length, border_left_pad, border_right_pad * '\n', bmap)

    # compute position of ylabel
    y_lab_row = round(nrows(c) / 2, RoundNearestTiesUp)

    callback = colormap_callback(p.colormap)

    # plot all rows
    for row in 1:nrows(c)
        # print left annotations
        print(io, repeat(🗷, p.margin))
        if p.labels
            # Current labels to left and right of the row and their length
            left_str  = get(p.labels_left, row, "")
            left_col  = get(p.colors_left, row, :light_black)
            right_str = get(p.labels_right, row, "")
            right_col = get(p.colors_right, row, :light_black)
            left_len  = length(_nocolor_string(left_str))
            right_len = length(_nocolor_string(right_str))
            if !get(io, :color, false)
                left_str  = _nocolor_string(left_str)
                right_str = _nocolor_string(right_str)
            end
            if !p.compact && row == y_lab_row
                # print ylabel
                print_color(:normal, io, p.ylabel)
                print(io, repeat(🗷, max_len_l - length(p.ylabel) - left_len))
            else
                # print padding to fill ylabel length
                print(io, repeat(🗷, max_len_l - left_len))
            end
            # print the left annotation
            print_color(left_col, io, left_str)
        end
        if c.visible
            # print left border
            print(io, plot_padding)
            print_color(:light_black, io, bmap[:l])
            # print canvas row
            printrow(io, c, row)
            # print right label and padding
            print_color(:light_black, io, bmap[:r])
        end
        if p.labels
            print(io, plot_padding)
            print_color(right_col, io, right_str)
            print(io, repeat(🗷, max_len_r - right_len))
        end
        # print colorbar
        if p.colorbar
            print(io, plot_padding)
            printcolorbarrow(
                io,
                c,
                row,
                callback,
                p.colorbar_border,
                p.colorbar_lim,
                (min_z_str, max_z_str),
                plot_padding,
                p.zlabel,
                cbar_max_len,
                🗷,
            )
        end
        row < nrows(c) && println(io)
    end

    # draw bottom border and bottom labels
    c.visible &&
        print_border(io, :b, border_length, '\n' * border_left_pad, border_right_pad, bmap)
    if p.labels
        print_labels(
            io,
            :b,
            p,
            border_length - 2,
            '\n' * border_left_pad * 🗹,
            🗹 * border_right_pad,
            🗹,
        )
        p.compact || print_title(
            io,
            '\n' * border_left_pad,
            p.xlabel,
            border_right_pad,
            🗹;
            p_width = p_width,
        )
    end
    nothing
end

"""
    savefig(p, filename; color=false)

Print the given plot `p` to a text file.

By default, it does not write ANSI color codes to the file. To enable this, set the keyword `color=true`.

# Examples
```julia-repl
julia> savefig(lineplot([0, 1]), "foo.txt")

```
"""
function savefig(p::Plot, filename::String; color::Bool = false)
    ext = lowercase(splitext(filename)[2])
    if ext in (".png", ".jpg", ".jpeg", ".tif", ".gif", ".svg")
        @warn "`UnicodePlots.savefig` only support writing to text files"
    end
    open(filename, "w") do io
        print(IOContext(io, :color => color), p)
    end
    nothing
end


text:
module RulesBIFMHelperInTest

using Test
using ReactiveMP
using Random
import ReactiveMP: @test_rules

@testset "rules:BIFMHelper:in" begin

     @testset "Belief Propagation: (m_out::Any, )" begin

        @test_rules [ with_float_conversions = true ] BIFMHelper(:in, Marginalisation) [
            (input = (m_out = PointMass(1.0), ),                 output = PointMass(1.0) ),
            (input = (m_out = NormalMeanVariance(9.5, 3.2), ),   output = NormalMeanVariance(9.5, 3.2) ),
            (input = (m_out = MvNormalWeightedMeanPrecision([1.0, 2.0], [3.0 0; 0 2.0]), ), output = MvNormalWeightedMeanPrecision([1.0, 2.0], [3.0 0; 0 2.0]) )
        ]

    end

end

end


text:
# Julia wrapper for header: /usr/include/scip/cons_quadratic.h
# Automatically generated using Clang.jl wrap_c

function SCIPincludeConshdlrQuadratic(scip)
    ccall((:SCIPincludeConshdlrQuadratic, libscip), SCIP_RETCODE, (Ptr{SCIP_},), scip)
end

function SCIPincludeQuadconsUpgrade(scip, quadconsupgd, priority, active, conshdlrname)
    ccall((:SCIPincludeQuadconsUpgrade, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Cvoid}, Cint, UInt32, Cstring), scip, quadconsupgd, priority, active, conshdlrname)
end

function SCIPcreateConsQuadratic(scip, cons, name, nlinvars, linvars, lincoefs, nquadterms, quadvars1, quadvars2, quadcoeffs, lhs, rhs, initial, separate, enforce, check, propagate, _local, modifiable, dynamic, removable)
    ccall((:SCIPcreateConsQuadratic, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_CONS}}, Cstring, Cint, Ptr{Ptr{SCIP_VAR}}, Ptr{Cdouble}, Cint, Ptr{Ptr{SCIP_VAR}}, Ptr{Ptr{SCIP_VAR}}, Ptr{Cdouble}, Cdouble, Cdouble, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32), scip, cons, name, nlinvars, linvars, lincoefs, nquadterms, quadvars1, quadvars2, quadcoeffs, lhs, rhs, initial, separate, enforce, check, propagate, _local, modifiable, dynamic, removable)
end

function SCIPcreateConsBasicQuadratic(scip, cons, name, nlinvars, linvars, lincoefs, nquadterms, quadvars1, quadvars2, quadcoefs, lhs, rhs)
    ccall((:SCIPcreateConsBasicQuadratic, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_CONS}}, Cstring, Cint, Ptr{Ptr{SCIP_VAR}}, Ptr{Cdouble}, Cint, Ptr{Ptr{SCIP_VAR}}, Ptr{Ptr{SCIP_VAR}}, Ptr{Cdouble}, Cdouble, Cdouble), scip, cons, name, nlinvars, linvars, lincoefs, nquadterms, quadvars1, quadvars2, quadcoefs, lhs, rhs)
end

function SCIPcreateConsQuadratic2(scip, cons, name, nlinvars, linvars, lincoefs, nquadvarterms, quadvarterms, nbilinterms, bilinterms, lhs, rhs, initial, separate, enforce, check, propagate, _local, modifiable, dynamic, removable)
    ccall((:SCIPcreateConsQuadratic2, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_CONS}}, Cstring, Cint, Ptr{Ptr{SCIP_VAR}}, Ptr{Cdouble}, Cint, Ptr{SCIP_QUADVARTERM}, Cint, Ptr{SCIP_BILINTERM}, Cdouble, Cdouble, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32), scip, cons, name, nlinvars, linvars, lincoefs, nquadvarterms, quadvarterms, nbilinterms, bilinterms, lhs, rhs, initial, separate, enforce, check, propagate, _local, modifiable, dynamic, removable)
end

function SCIPcreateConsBasicQuadratic2(scip, cons, name, nlinvars, linvars, lincoefs, nquadvarterms, quadvarterms, nbilinterms, bilinterms, lhs, rhs)
    ccall((:SCIPcreateConsBasicQuadratic2, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_CONS}}, Cstring, Cint, Ptr{Ptr{SCIP_VAR}}, Ptr{Cdouble}, Cint, Ptr{SCIP_QUADVARTERM}, Cint, Ptr{SCIP_BILINTERM}, Cdouble, Cdouble), scip, cons, name, nlinvars, linvars, lincoefs, nquadvarterms, quadvarterms, nbilinterms, bilinterms, lhs, rhs)
end

function SCIPaddConstantQuadratic(scip, cons, constant)
    ccall((:SCIPaddConstantQuadratic, libscip), Cvoid, (Ptr{SCIP_}, Ptr{SCIP_CONS}, Cdouble), scip, cons, constant)
end

function SCIPaddLinearVarQuadratic(scip, cons, var, coef)
    ccall((:SCIPaddLinearVarQuadratic, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}, Ptr{SCIP_VAR}, Cdouble), scip, cons, var, coef)
end

function SCIPaddQuadVarQuadratic(scip, cons, var, lincoef, sqrcoef)
    ccall((:SCIPaddQuadVarQuadratic, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}, Ptr{SCIP_VAR}, Cdouble, Cdouble), scip, cons, var, lincoef, sqrcoef)
end

function SCIPaddQuadVarLinearCoefQuadratic(scip, cons, var, coef)
    ccall((:SCIPaddQuadVarLinearCoefQuadratic, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}, Ptr{SCIP_VAR}, Cdouble), scip, cons, var, coef)
end

function SCIPaddSquareCoefQuadratic(scip, cons, var, coef)
    ccall((:SCIPaddSquareCoefQuadratic, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}, Ptr{SCIP_VAR}, Cdouble), scip, cons, var, coef)
end

function SCIPaddBilinTermQuadratic(scip, cons, var1, var2, coef)
    ccall((:SCIPaddBilinTermQuadratic, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}, Ptr{SCIP_VAR}, Ptr{SCIP_VAR}, Cdouble), scip, cons, var1, var2, coef)
end

function SCIPgetNlRowQuadratic(scip, cons, nlrow)
    ccall((:SCIPgetNlRowQuadratic, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}, Ptr{Ptr{SCIP_NLROW}}), scip, cons, nlrow)
end

function SCIPgetNLinearVarsQuadratic(scip, cons)
    ccall((:SCIPgetNLinearVarsQuadratic, libscip), Cint, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPgetLinearVarsQuadratic(scip, cons)
    ccall((:SCIPgetLinearVarsQuadratic, libscip), Ptr{Ptr{SCIP_VAR}}, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPgetCoefsLinearVarsQuadratic(scip, cons)
    ccall((:SCIPgetCoefsLinearVarsQuadratic, libscip), Ptr{Cdouble}, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPgetNQuadVarTermsQuadratic(scip, cons)
    ccall((:SCIPgetNQuadVarTermsQuadratic, libscip), Cint, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPgetQuadVarTermsQuadratic(scip, cons)
    ccall((:SCIPgetQuadVarTermsQuadratic, libscip), Ptr{SCIP_QUADVARTERM}, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPsortQuadVarTermsQuadratic(scip, cons)
    ccall((:SCIPsortQuadVarTermsQuadratic, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPfindQuadVarTermQuadratic(scip, cons, var, pos)
    ccall((:SCIPfindQuadVarTermQuadratic, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}, Ptr{SCIP_VAR}, Ptr{Cint}), scip, cons, var, pos)
end

function SCIPgetNBilinTermsQuadratic(scip, cons)
    ccall((:SCIPgetNBilinTermsQuadratic, libscip), Cint, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPgetBilinTermsQuadratic(scip, cons)
    ccall((:SCIPgetBilinTermsQuadratic, libscip), Ptr{SCIP_BILINTERM}, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPgetLhsQuadratic(scip, cons)
    ccall((:SCIPgetLhsQuadratic, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPgetRhsQuadratic(scip, cons)
    ccall((:SCIPgetRhsQuadratic, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPgetLinvarMayDecreaseQuadratic(scip, cons)
    ccall((:SCIPgetLinvarMayDecreaseQuadratic, libscip), Cint, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPgetLinvarMayIncreaseQuadratic(scip, cons)
    ccall((:SCIPgetLinvarMayIncreaseQuadratic, libscip), Cint, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPcheckCurvatureQuadratic(scip, cons)
    ccall((:SCIPcheckCurvatureQuadratic, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPisConvexQuadratic(scip, cons)
    ccall((:SCIPisConvexQuadratic, libscip), UInt32, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPisConcaveQuadratic(scip, cons)
    ccall((:SCIPisConcaveQuadratic, libscip), UInt32, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPgetViolationQuadratic(scip, cons, sol, violation)
    ccall((:SCIPgetViolationQuadratic, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}, Ptr{SCIP_SOL}, Ptr{Cdouble}), scip, cons, sol, violation)
end

function SCIPisLinearLocalQuadratic(scip, cons)
    ccall((:SCIPisLinearLocalQuadratic, libscip), UInt32, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPaddToNlpiProblemQuadratic(scip, cons, nlpi, nlpiprob, scipvar2nlpivar, names)
    ccall((:SCIPaddToNlpiProblemQuadratic, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}, Ptr{SCIP_NLPI}, Ptr{SCIP_NLPIPROBLEM}, Ptr{SCIP_HASHMAP}, UInt32), scip, cons, nlpi, nlpiprob, scipvar2nlpivar, names)
end

function SCIPchgLhsQuadratic(scip, cons, lhs)
    ccall((:SCIPchgLhsQuadratic, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}, Cdouble), scip, cons, lhs)
end

function SCIPchgRhsQuadratic(scip, cons, rhs)
    ccall((:SCIPchgRhsQuadratic, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}, Cdouble), scip, cons, rhs)
end

function SCIPgetFeasibilityQuadratic(scip, cons, sol, feasibility)
    ccall((:SCIPgetFeasibilityQuadratic, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}, Ptr{SCIP_SOL}, Ptr{Cdouble}), scip, cons, sol, feasibility)
end

function SCIPgetActivityQuadratic(scip, cons, sol, activity)
    ccall((:SCIPgetActivityQuadratic, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}, Ptr{SCIP_SOL}, Ptr{Cdouble}), scip, cons, sol, activity)
end

function SCIPchgLinearCoefQuadratic(scip, cons, var, coef)
    ccall((:SCIPchgLinearCoefQuadratic, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}, Ptr{SCIP_VAR}, Cdouble), scip, cons, var, coef)
end

function SCIPchgSquareCoefQuadratic(scip, cons, var, coef)
    ccall((:SCIPchgSquareCoefQuadratic, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}, Ptr{SCIP_VAR}, Cdouble), scip, cons, var, coef)
end

function SCIPchgBilinCoefQuadratic(scip, cons, var1, var2, coef)
    ccall((:SCIPchgBilinCoefQuadratic, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}, Ptr{SCIP_VAR}, Ptr{SCIP_VAR}, Cdouble), scip, cons, var1, var2, coef)
end

function SCIPgetNAllBilinearTermsQuadratic(scip)
    ccall((:SCIPgetNAllBilinearTermsQuadratic, libscip), Cint, (Ptr{SCIP_},), scip)
end

function SCIPgetAllBilinearTermsQuadratic(scip, x, y, nbilinterms, nunderests, noverests, maxnonconvexity)
    ccall((:SCIPgetAllBilinearTermsQuadratic, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_VAR}}, Ptr{Ptr{SCIP_VAR}}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble}), scip, x, y, nbilinterms, nunderests, noverests, maxnonconvexity)
end

function SCIPaddBilinearIneqQuadratic(scip, x, y, idx, xcoef, ycoef, constant, success)
    ccall((:SCIPaddBilinearIneqQuadratic, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_VAR}, Ptr{SCIP_VAR}, Cint, Cdouble, Cdouble, Cdouble, Ptr{UInt32}), scip, x, y, idx, xcoef, ycoef, constant, success)
end


text:
#unum-test-division.jl

#some simple tests in Unum{4,5}
x = Unum{4,5}(30.0)
y = Unum{4,5}(1.5)
z = (x / y)
@test is_exact(z)
@test calculate(z) == 20.0

x = Unum{4,6}(30.0)
y = Unum{4,6}(1.5)
z = (x / y)
@test is_exact(z)
@test calculate(z) == 20.0

x = Unum{4,7}(30.0)
y = Unum{4,7}(1.5)
z = (x / y)
@test is_exact(z)
@test calculate(z) == 20.0

x = Unum{4,8}(30.0)
y = Unum{4,8}(1.5)
z = (x / y)
@test is_exact(z)
@test calculate(x / y) == 20.0

x = Unums.make_exact!(Unum{4,5}(1) / Unum{4,5}(3))
@test (x / x) == one(Unum{4,5})

x = Unums.make_exact!(Unum{4,6}(1) / Unum{4,6}(3))
@test (x / x) == one(Unum{4,6})

x = Unum{4,7}(1.0) / Unum{4,7}(3.0)
@test x == Unum{4,7}(0x0000000000000001, UInt64[0x5555555555555555,0x5555555555555555], 0x0001, 0x0002, 0x007f)
Unums.make_exact!(x)
@test (x / x) == one(Unum{4,7})

################################################################################
# test with ubounds
x = Ubound(Unum{4,6}(6), Unum{4,6}(7))
y = Unum{4,6}(3)
@test (x / y) == Ubound(Unum{4,6}(2), Unum{4,6}(7) / Unum{4,6}(3))

################################################################################
# discovered errors.

#as a result of matrix solving.  Jan 2017

x = Unum{3,5}(0x0000000000000001, 0x0000000000000000, 0x0000, 0x0001, 0x0000)
y = Ubound(Unum{3,5}(0x000000000000000F, 0x373FE2B700000000, 0x0001, 0x0003, 0x001F), Unum{3,5}(0x000000000000000F, 0x373FE2D800000000, 0x0001, 0x0003, 0x001E))

@test (x / y) == Ubound(Unum{3,5}(0x0000000000000006, 0xA51D8EDD00000000, 0x0001, 0x0004, 0x001F), Unum{3,5}(0x0000000000000006, 0xA51D8F0C00000000, 0x0001, 0x0004, 0x001F))

x = Ubound(Unum{3,5}(0x0000000000000031, 0x77A419CC00000000, 0x0003, 0x0005, 0x001F), Unum{3,5}(0x0000000000000031, 0x77A4195000000000, 0x0003, 0x0005, 0x001F))
y = Ubound(Unum{3,5}(0x0000000000000018, 0x071A55DD00000000, 0x0003, 0x0004, 0x001F), Unum{3,5}(0x0000000000000018, 0x071A55D200000000, 0x0003, 0x0004, 0x001F))
@test (x / y) ==  Ubound(Unum{3,5}(0x0000000000000018, 0x6D7FFFB800000000, 0x0001, 0x0004, 0x001D), Unum{3,5}(0x0000000000000018, 0x6D80004000000000, 0x0001, 0x0004, 0x001D))

#more parity errors

x = Ubound(Unum{3,5}(0x0000000000000030, 0xDD99D88C00000000, 0x0003, 0x0005, 0x001E), Unum{3,5}(0x0000000000000030, 0xDD99D87800000000, 0x0003, 0x0005, 0x001D))
y = Ubound(Unum{3,5}(0x0000000000000018, 0x1BADF0DB00000000, 0x0001, 0x0004, 0x001F), Unum{3,5}(0x0000000000000018, 0x1BADF0DC00000000, 0x0001, 0x0004, 0x001D))

@test x / y == Ubound(Unum{3,5}(0x000000000000000F, 0xAF00000800000000, 0x0003, 0x0003, 0x001D), Unum{3,5}(0x000000000000000F, 0xAEFFFFF000000000, 0x0003, 0x0003, 0x001F))

max_stars_repo_name: REX-Computing/unumjl

max_issues_repo_name: REX-Computing/unumjl

max_forks_repo_name: REX-Computing/unumjl