text:
export genus, representative, rank, det, uniformizer, det_representative,
       gram_matrix, genera_hermitian, local_genera_hermitian, rank,
       orthogonal_sum, isinert, scales, ranks, dets, issplit, isramified,
       isdyadic, norms, primes, signatures

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

# 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
  issplit::Bool                       # p split in E
  non_norm_rep::FieldElem             # 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

local_genus_herm_type(E) = LocalGenusHerm{typeof(E), ideal_type(order_type(base_field(E)))}

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

function Base.show(io::IO, ::MIME"text/plain", G::LocalGenusHerm)
  compact = get(io, :compact, false)
  if !compact
    if isdyadic(G) && isramified(G)
       print(io, "Local genus symbol (scale, rank, det, norm) at ")
    else
      print(io, "Local genus symbol (scale, rank, det) at ")
    end
    if length(G) > 0
      print(IOContext(io, :compact => true), prime(G), ":")
      print(io, "\n")
    else
      print(IOContext(io, :compact => true), prime(G), " of rank zero")
    end
  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 local genus symbol `g` for hermitian lattices over $E/K$ at a prime $\mathfrak
p$ of $\mathcal O_K$, return the $\mathfrak P$-valuation of the scale of the `i`th
Jordan block of `g`, where $\mathfrak P$ is a prime ideal of $\mathcal O_E$ lying
over $\mathfrak p$.
"""
scale(G::LocalGenusHerm, i::Int) = G.data[i][1]

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

Given a local genus symbol `g` for hermitian lattices over $E/K$ at a prime $\mathfrak
p$ of $\mathcal O_K$, return the $\mathfrak P$-valuation of the scales of the Jordan
blocks of `g`, where $\mathfrak P$ is a prime ideal of $\mathcal O_E$ 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) -> Int

Given a local genus symbol `g` for hermitian lattices, 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 local genus symbol `g` for hermitian lattices over $E/K$ at a prime ideal
$\mathfrak p$ of $\mathcal O_K$, return the rank of any hermitian lattice whose
$\mathfrak p$-adic completion has local genus symbol `g`.
"""
function rank(G::LocalGenusHerm)
  return reduce(+, (rank(G, i) for i in 1:length(G)), init = Int(0))
end

@doc Markdown.doc"""
    ranks(g::LocalGenusHerm) -> Vector{Int}

Given a local genus symbol `g` for hermitian lattices, return the ranks of the
Jordan blocks 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 local genus symbol `g` for hermitian lattices over $E/K$, return the determinant
of the `i`th Jordan block of `g`.

The returned value is $1$ or $-1$ depending on whether the determinant is a local norm in `K`.
"""
det(G::LocalGenusHerm, i::Int) = G.data[i][3]

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

Given a local genus symbol `g` for hermitian lattices over $E/K$ at a prime ideal
$\mathfrak p$ of $\mathcal O_K$, return the determinant of a hermitian lattice
whose $\mathfrak p$-adic completion has local genus symbol `g`.

The returned value is $1$ or $-1$ depending on whether the determinant is a local
norm in `K`.
"""
function det(G::LocalGenusHerm)
  return reduce(*, (det(G, i) for i in 1:length(G)), init = Int(1))
end

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

Given a local genus symbol `g` for hermitian lattices over $E/K$, return the determinants
of the Jordan blocks of `g`.

The returned values are $1$ or $-1$ depending on whether the respective determinants are
are local norms in `K`.
"""
dets(G::LocalGenusHerm) = map(i -> det(G, i), 1:length(G))::Vector{Int}

@doc Markdown.doc"""
    discriminant(g::LocalGenusHerm, i::Int) -> Int

Given a local genus symbol `g` for hermitian lattices over $E/K$, return the discriminant
of the `i`th Jordan block of `g`.

The returned value is $1$ or $-1$ depending on whether the discriminant is a local norm in `K`.
"""
function discriminant(G::LocalGenusHerm, i::Int)
  d = det(G, i)
  r = rank(G, i) % 4
  if !isramified(G) || r == 0 || r == 1
    return d
  end
  E = base_field(G)
  fl = islocal_norm(E, base_field(E)(-1), prime(G))
  if fl
    return d
  else
    return -d
  end
end

@doc Markdown.doc"""
    discriminant(g::LocalGenusHerm) -> Int

Given a local genus symbol `g` for hermitian lattices over $E/K$ at a prime ideal
$\mathfrak p$ of $\mathcal O_K$, return the discriminant of a hermitian lattice
whose $\mathfrak p$-adic completion has local genus symbol `g`.

The returned value is $1$ or $-1$ depending on whether the discriminant is a local
norm in `K`.
"""
function discriminant(G::LocalGenusHerm)
  d = det(G)
  r = rank(G) % 4
  if !isramified(G) || r == 0 || r == 1
    return d
  end
  E = base_field(G)
  fl = islocal_norm(E, base_field(E)(-1), prime(G))
  if fl
    return d
  else
    return -d
  end
end

# this only works if it is ramified and dyadic
@doc Markdown.doc"""
    norm(g::LocalGenusHerm, i::Int) -> Int

Given a ramified dyadic local genus symbol `g` for hermitian lattices over $E/K$ at a
prime ideal $\mathfrak p$ of $\mathcal O_K$, return the $\mathfrak p$-valuation of
the norm of the `i`th Jordan block of `g`.
"""
norm(G::LocalGenusHerm, i::Int) = begin @assert isdyadic(G) && isramified(G); G.norm_val[i] end

# this only works if it is ramified and dyadic
@doc Markdown.doc"""
    norms(g::LocalGenusHerm) -> Vector{Int}

Given a ramified dyadic local genus symbol `g` for hermitian lattices over $E/K$ at a
prime ideal $\mathfrak p$ of $\mathcal O_K$, return the $\mathfrak p$-valuations of the
norms of the Jordan blocks of `g`.
"""
norms(G::LocalGenusHerm) = begin @assert isdyadic(G) && isramified(G); G.norm_val end

@doc Markdown.doc"""
    isramified(g::LocalLenusHerm) -> Bool

Given a local genus symbol `g` for hermitian lattices over $E/K$ at a prime ideal
$\mathfrak p$ of $\mathcal O_K$, return whether $\mathfrak p$ is ramified in
$\mathcal O_E$.
"""
isramified(g::LocalGenusHerm) = g.isramified

@doc Markdown.doc"""
    issplit(g::LocalGenusHerm) -> Bool

Given a local genus symbol `g` for hermitian lattices over $E/K$ at a prime ideal
$\mathfrak p$ of $\mathcal O_K$, return whether $\mathfrak p$ is split in
$\mathcal O_E$.
"""
issplit(g::LocalGenusHerm) = g.issplit

@doc Markdown.doc"""
    isinert(g::LocalGenusHerm) -> Bool

Given a local genus symbol `g` for hermitian lattices over $E/K$ at a prime ideal
$\mathfrak p$ of $\mathcal O_K$, return whether $\mathfrak p$ is inert in
$\mathcal O_E$.
"""
isinert(g::LocalGenusHerm) = !g.isramified && !g.issplit

@doc Markdown.doc"""
    isdyadic(g::LocalGenusHerm) -> Bool

Given a local genus symbol `g` for hermitian lattices over $E/K$ at a prime ideal
$\mathfrak p$ of $\mathcal O_K$, return whether $\mathfrak p$ is dyadic.
"""
isdyadic(G::LocalGenusHerm) = G.isdyadic

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

Given a local genus symbol `g` for hermitian lattices, return the number of Jordan blocks
of `g`.
"""
length(G::LocalGenusHerm) = length(G.data)

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

Given a local genus symbol `g` for hermitian lattices over $E/K$, return `E`.
"""
base_field(G::LocalGenusHerm) = G.E

@doc Markdown.doc"""
    prime(g::LocalGenusHerm) -> NfOrdIdl

Given a local genus symbol `g` for hermitian lattices over $E/K$ at a prime ideal
$\mathfrak p$ of $\mathcal O_K$, return $\mathfrak p$.
"""
prime(G::LocalGenusHerm) = G.p

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

# If G is defined over E/K at a prime p of O_K, this returns an unit in K which is
# not a local norm at p.
function _non_norm_rep(G::LocalGenusHerm)
  if isdefined(G, :non_norm_rep)
    return G.non_norm_rep::elem_type(base_field(base_field(G)))
  end

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

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

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

Given a local genus symbol `g` for hermitian lattices over $E/K$ at a prime ideal
$\mathfrak p$ of $\mathcal O_K$, return a generator for the largest ideal of $\mathcal O_E$
containing $\mathfrak p$ and invariant under the action of the non-trivial involution
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 K(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
#
################################################################################

# If p does not ramify in E, then the determinant of G is pi^v where pi is a
# uniformizer of p. Otherwise pi is a local norm and the class of the
# determinant of G is the same as u*pi^(v//2) where u is 1 if d == 1, otherwise
# it is a non local norm. Indeed the valuation is with respect to p = P^2
# but the scale is with respect to P and thus the determinant of a G is
# represented P^(scale*rank) = p^(scale*rank/2) times u.
@doc Markdown.doc"""
    det_representative(g::LocalGenusHerm) -> NumFieldElem

Given a local genus symbol `g` for hermitian lattices over $E/K$, return a
representative of the norm class of the determinant of `g` in $K^{\times}$.
"""
function det_representative(G::LocalGenusHerm)
  d = det(G)
  v = sum(scale(G, i) * rank(G, i) for i in 1:length(G); init = 0)
  if !isramified(G)
    return uniformizer(G)^v
  end

  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

@doc Markdown.doc"""
    det_representative(g::LocalGenusHerm, i::Int) -> NumFieldElem

Given a local genus symbol `g` for hermitian lattices over $E/K$, return a
representative of the norm class of the determinant of the `i`th Jordan block of
`g` in $K^{\times}$.
"""
function det_representative(G::LocalGenusHerm, i::Int)
  d = det(G, i)
  v = scale(G, i) * rank(G, i)
  if !isramified(G)
    return uniformizer(G)^v
  end

  v = div(v, 2)

  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) -> MatElem

Given a local genus symbol `g` for hermitian lattices over $E/K$ at a prime ideal
$\mathfrak p$ of $\mathcal O_K$, return a Gram matrix `M` of `g`, with coefficients
in `E`.`M` is such that any hermitian lattice over $E/K$ with Gram matrix `M` satisfies
that the local genus symbol of its completion at $\mathfrak p$ is `g`.
"""
function gram_matrix(G::LocalGenusHerm)
  if rank(G) == 0
    return zero_matrix(base_field(G), 0, 0)
  end
  return diagonal_matrix(dense_matrix_type(base_field(G))[gram_matrix(G, i) for i in 1:length(G)])
end

@doc Markdown.doc"""
    gram_matrix(g::LocalGenusHerm, i::Int) -> MatElem

Given a local genus symbol `g` for hermitian lattices over $E/K$ at a prime ideal
$\mathfrak p$ of $\mathcal O_K$, return a Gram matrix `M` of the `i`th Jordan block
of `g`, with coefficients in `E`. `M` is such that any hermitian lattice over $E/K$
with Gram matrix `M` satisfies that the local genus symbol of its completion at
$\mathfrak p$ is equal to the `i`th Jordan block of `g`.
"""
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)))
  conj = involution(E)

  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, elem_type(E)[zero(E), pi^i, conj(pi)^i, zero(E)])

  if !isdyadic(G)
    # non-dyadic
    if iseven(i)
      # According to Kir16, there the last exponent should be i/2 * (1 - m)
      lastexp = div(i, 2) * (1 - m)
      drep = det_representative(G, l)
      return diagonal_matrix(push!(elem_type(E)[E(p)^div(i, 2) for j in 1:(m - 1)], drep * E(p)^(lastexp)))
    else
      return diagonal_matrix(dense_matrix_type(E)[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)
      # hyperbolic
      @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 # not hyperbolic
      @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 symbol `g` for hermitian lattices over $E/K$ at a prime ideal
$\mathfrak p$ of $\mathcal O_K$, return a hermitian lattice over $E/K$ whose
completion at $\mathfrak p$ admits `g` as local genus symbol.
"""
function representative(G::LocalGenusHerm)
  E = base_field(G)
  L = lattice(hermitian_space(E, gram_matrix(G)))
  S = ideal_type(base_ring(base_ring(L)))
  GType = local_genus_herm_type(E)
  symbols = get_attribute!(L, :local_genus) do
    Dict{S, GType}()
  end::Dict{S, GType}

  get!(symbols, prime(G), G)

  return L
end

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

@doc Markdown.doc"""
    genus(HermLat, E::NumField, p::NfOrdIdl, data::Vector; type::Symbol = :det,
		                                           check::Bool = true)
                                                                 -> LocalGenusHerm

Construct the local genus symbol `g` for hermitian lattices over the algebra `E`,
with base field $K$, at the prime ideal `p` of $\mathcal O_K$. Its invariants are
specified by `data`.

- If the prime ideal is good (not ramified-and-dyadic), the elements of `data` must
  be `(s, r, d)::Tuple{Int, Int, Int}` where `s` refers to the scale $\mathfrak P$-valuation,
  `r` the rank and `d` the determinant/discriminant class of a Jordan block of `g`, where
  $\mathfrak P$ is a prime ideal of $\mathcal O_E$ lying over `p`.

  In the unramified case, `d` is determined by `s` and `r` and can be omitted.
  Hence also `(s, r)::Tuple{Int, Int}` is allowed.

- If the prime ideal is bad (ramified and dyadic), the elements of `data` must
  be `(s, r, d, n)::Tuple{Int, Int, Int, Int}`, where in addition `n`
  refers to the norm `p`-valuation.

Additional comments:
- `d` must be in $\{[1, -1\}$;
- If `type == :disc`, the parameter `d` is interpreted as the discriminant;
- Sanity checks can be disabled by setting `check == false`.
"""
genus(::Type{HermLat}, E, p, data; type, check)

# rank zero genus
function genus(::Type{HermLat}, E::S, p::T) where {S, T}
  if isramified(maximal_order(E), p) && isdyadic(p)
    return genus(HermLat, E, p, Vector{Tuple{Int, Int, Int, Int}}())
  else
    return genus(HermLat, E, p, Vector{Tuple{Int, Int, Int}}())
  end
end

# Some comments
#
# We distinguish between bad and good case
# First the good case

# This is the internal function which requires the decomposition behavior of
# the prime to be already determined and does not do any internal checks.
function _genus(::Type{HermLat}, E::S, p::T, data::Vector{Tuple{Int, Int, Int}}, is_dyadic, is_ramified, is_split) where {S, T}
  z = LocalGenusHerm{S, T}()
  z.E = E
  z.p = p
  @hassert :Lattice 1 !(is_dyadic && is_ramified)
  z.isdyadic = is_dyadic
  z.isramified = is_ramified
  z.issplit = is_split
  z.data = data
  return z
end

# This is one of the two user facing functions for the good case, namely the
# unramified case. Here the determinant/discriminant class need not be supplied.
function genus(::Type{HermLat}, E, p, data::Vector{Tuple{Int, Int}}; check::Bool = true)
  @req !isramified(maximal_order(E), p) "For dyadic primes the norm valuation has to be specified"
  if check
    @req all(data[i][2] >= 0 for i in 1:length(data)) "Ranks must be positive"
    @req all(data[i][1] < data[i + 1][1] for i in 1:length(data)-1) "Scales must be strictly increasing"
  end

  is_dyadic = isdyadic(p)

  lp = prime_decomposition(maximal_order(E), p)
  if length(lp) == 2
    is_split = true
    is_ramified = false
  else
    is_split = false
    if lp[1][2] == 1
      is_ramified = false
    else
      is_ramified = true
    end
  end

  cdata = Vector{Tuple{Int, Int, Int}}(undef, length(data))

  l = length(data)

  if !is_split
    # inert case
    for i in 1:l
      # The determinant class is the class of E(p)^(r * s) and E(p)
      # has class -1.
      if isodd(data[i][1]) && isodd(data[i][2])
        cdata[i] = (data[i][1], data[i][2], -1)
      else
        cdata[i] = (data[i][1], data[i][2], 1)
      end
    end
  else
    for i in 1:l
      # split case
      cdata[i] = (data[i][1], data[i][2], 1)
    end
  end
  return _genus(HermLat, E, p, cdata, is_dyadic, is_ramified, is_split)
end

# This is the second user facing functions for the general good case. Here the
# determinant/discriminant class must be supplied. We also do some sanity
# checks in the unramified case concerning the det/disc.
function genus(::Type{HermLat}, E::S, p::T, data::Vector{Tuple{Int, Int, Int}}; type = :det, check::Bool = true) where {S <: NumField, T}
  @req type === :det || type === :disc "type :$type must be :disc or :det"

  if check
    @req all(data[i][2] >= 0 for i in 1:length(data)) "Ranks must be positive"
    @req all(data[i][1] < data[i + 1][1] for i in 1:length(data)-1) "Scales must be strictly increasing"
    @req all(abs(data[i][3]) == 1 for i in 1:length(data)) "Norm classes must be +/-1"
  end

  # Determine the prime decomposition
  lp = prime_decomposition(maximal_order(E), p)
  if length(lp) == 2
    is_split = true
    is_ramified = false
  else
    is_split = false
    if lp[1][2] == 1
      is_ramified = false
    else
      is_ramified = true
    end
  end

  is_dyadic = isdyadic(p)

  @req !(is_dyadic && is_ramified) "For dyadic primes the norm valuation has to be specified"

  l = length(data)

  cdata = copy(data)

  if type === :disc
    # We need to swap the sign depending on whether some rank is 2, 3 mod 4
    # and -1 is not a local norm. The latter can only happen in the ramified case
    if any(data[i][2] % 4 in [2, 3] for i in 1:l) && is_ramified
      fl = islocal_norm(E, base_field(E)(-1), p)
      if !fl
        for i in 1:l
          r = data[i][2] % 4
          if r == 2 || r == 3
            cdata[i] = (cdata[i][1], cdata[i][2], -data[i][3])
          end
        end
      end
    end
  end

  # Now check that the determinant class fits with the scale and rank.
  # There is a restriction only in the unramified case
  if check
    if !is_ramified
      if !is_split
        # inert case
        for i in 1:l
          # The determinant class is the class of E(p)^(r * s) and E(p)
          # has class -1.
          if isodd(cdata[i][1]) && isodd(cdata[i][2])
            @req cdata[i][3] == -1 "$(type === :disc ? "Discriminant" : "Determinant") class does not fit scale and rank"
          else
            @req cdata[i][3] == 1 "$(type === :disc ? "Discriminant" : "Determinant") class does not fit scale and rank"
          end
        end
      else
        for i in 1:l
          # In the split case the class must always be 1.
          @req cdata[i][3] == 1 "$(type === :disc ? "Discriminant" : "Determinant") classes must be 1"
        end
      end
    else
      # Non-dyadic ramified
      # If the scale is odd, then the rank m must be even and
      # the lattice is H(i)^(m/2), hence has determinant class
      # [-1]^(m/2)
      #
      fl = islocal_norm(E, base_field(E)(-1), p)
      for i in 1:l
        if isodd(cdata[i][1])
          @req iseven(cdata[i][2]) "Rank must be even for blocks of odd scale"
          m2 = div(cdata[i][2], 2)
          if fl
            @req cdata[i][3] == 1 "Determinant must be 1 for blocks of odd scale"
          else
            @req cdata[i][3] == (iseven(m2) ? 1 : -1) "Determinant mismatch in block $i"
          end
        end
      end
    end
  end
  # Now call the internal function
  return _genus(HermLat, E, p, cdata, is_dyadic, is_ramified, is_split)
end

# Now comes the bad case.
#
# First the internal function, which has as additonal argument the vector of norm valuations.
function _genus(::Type{HermLat}, E::S, p::T, data::Vector{Tuple{Int, Int, Int}}, norms::Vector{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)
  z.issplit = false
  # We test the cheap thing
  @req z.isdyadic && z.isramified "Prime must be dyadic and ramified"
  z.norm_val = norms
  z.data = data
  z.ni = _get_ni_from_genus(z)
  return z
end

# The user facing function in the bad case.
function genus(::Type{HermLat}, E::S, p::T, data::Vector{Tuple{Int, Int, Int, Int}}; type = :det, check::Bool = true) where {S <: NumField, T}
  is_dyadic = isdyadic(p)
  is_ramified = isramified(maximal_order(E), p)
  @req is_dyadic && is_ramified "Prime must be dyadic and ramified"
  @req type === :det || type === :disc "type :$type must be :disc or :det"

  cdata = Tuple{Int, Int, Int}[Base.front(v) for v in data]
  norm_val = Int[v[end] for v in data]

  if check
    @req all(data[i][1] < data[i + 1][1] for i in 1:length(data)-1) "Scales must be strictly increasing"
    @req all(abs(data[i][3]) == 1 for i in 1:length(data)) "Norm classes must be +/-1"
  end

  l = length(cdata)

  if check
    for i in 1:l
      # If the rank is odd, then n(L) * O_E = s(L), so n = 2 * s,
      # since n is the valuation in K and the extension is ramified.
      v = cdata[i]
      if isodd(v[2])
        @req 2 * norm_val[i] == v[1] """Not a valid local genus in block $(i):
                                        Scale ($(v[1])) must be twice the norm ($(norm_val[i]))"""
      end
      # TODO: We should also check using e, the valuation of the different
    end
  end

  if type === :disc && any(data[i][2] % 4 in [2, 3] for i in 1:l)
    # We need to swap the sign depending on whether some rank is 2, 3 mod 4
    # and -1 is not a local norm.
    fl = islocal_norm(E, base_field(E)(-1), p)
    if !fl
      for i in 1:l
        r = cdata[i][2] % 4
        if r == 0 || r == 1
          continue
        else
          cdata[i] = (data[i][1], data[i][2], (-1) * data[i][3])
        end
      end
    end
  end
  # Now call the internal function, no checks
  return _genus(HermLat, E, p, cdata, norm_val)
end

################################################################################
#
#  Local genus symbol of a lattice
#
################################################################################

# TODO: better documentation

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

Return the local genus symbol `g` for hermitian lattices over $E/K$ of the completion
of the hermitian lattice `L` at the prime ideal `p` of $\mathcal O_K$.
"""
function genus(L::HermLat, q)
  if typeof(q) === ideal_type(base_ring(base_ring(L)))
    # yippii, correct ideal type
    return _genus_correct_ideal_type(L, q)
  else
    if q isa fmpq || q isa Int
      # we allow this in case base_ring(base_ring(L)) == ZZ
      @req base_ring(base_ring(L)) isa FlintIntegerRing "Smaller field must be QQ"
      qq = ideal(base_ring(base_ring(L)), q)::ideal_type(base_ring(base_ring(L)))
      return _genus_correct_ideal_type(L, qq)
    end
  end
end

function _genus_correct_ideal_type(L, q)
  @assert typeof(q) === ideal_type(base_ring(base_ring(L)))
  S = ideal_type(base_ring(base_ring(L)))
  GType = local_genus_herm_type(base_field(L))
  symbols = get_attribute!(L, :local_genus) do
    return Dict{S, GType}()
  end::Dict{S, GType}

  return get!(symbols, q) do
    _genus(L, q)
  end::GType
end

function _genus(L::HermLat, p)
  bad, sym = _genus_symbol(L, p)
  if bad
    G = genus(HermLat, base_field(L), p, sym)
  else
    G = genus(HermLat, base_field(L), p, [Base.front(v) for v in sym])
  end
  # Just for debugging
  @hassert :Lattice 1 begin
    if isdyadic(G) && isramified(G)
      GG = _genus_symbol_kirschmer(L, p)
      all(let G = G; i -> GG[i][4] == G.ni[i]; end, 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

Return whether `g` and the local genus symbol of the completion of the hermitian
lattice `L` at `prime(g)` agree. Note that `L` being in `g` requires both `L` and
`g` to be defined over the same extension $E/K$.
"""
Base.in(L::HermLat, G::LocalGenusHerm) = base_field(L) === base_field(G) && genus(L, prime(G)) == G

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

@doc Markdown.doc"""
    ==(g1::LocalGenusHerm, g2::LocalGenusHerm) -> Bool

Given two local genus symbols `g1` and `g2`, return whether their respective Jordan
decompositions are of the same Jordan type. Note that equality requires `g1` and `g2`
to be defined over the same extension $E/K$ and at the same prime ideal $\mathfrak p$
of $\mathcal O_K$.
"""
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 !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 orthogonal 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), [i for i in 1:t if iseven(scale(G1,i))])
  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, 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)
  local sym::Vector{Tuple{Int, Int, Int, Int}}
  bad = true
  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)]
    bad = false
  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(dense_matrix_type(E)[pi^(max(0, S[i] - S[j])) * G[j] for j in 1:length(B)])
      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
    bad = true
  end
  return bad, sym
end

################################################################################
#
#  Sum of local genera
#
################################################################################

@doc Markdown.doc"""
    orthogonal_sum(g1::LocalGenusHerm, g2::LocalGenusHerm) -> LocalGenusHerm

Given two local genus symbols `g1` and `g2` for hermitian lattices over $E/K$
at the same prime ideal $\mathfrak p$ of $\mathcal O_K$, return their orthogonal
sum. It corresponds to the local genus symbol of the $\mathfrak p$-adic completion
of the orthogonal sum of respective representatives of `g1` and `g2`.
"""
function orthogonal_sum(G1::LocalGenusHerm, G2::LocalGenusHerm)
  @req prime(G1) == prime(G2) "Local genera must have the same prime ideal"
  if !G1.isdyadic || !G2.isramified
    return _direct_sum_easy(G1, G2)
  else
    L1 = representative(G1)
    L2 = representative(G2)
    L3, = orthogonal_sum(L1, L2)
    return genus(L3, prime(G1))
  end
end

function _direct_sum_easy(G1::LocalGenusHerm, G2::LocalGenusHerm)
  # We do a merge sort
  i1 = 1
  i2 = 1
  z = Tuple{Int, Int, Int}[]
  while i1 <= length(G1) && i2 <= length(G2)
    if scale(G1, i1) < scale(G2, i2)
      push!(z, G1.data[i1])
      i1 += 1
    elseif scale(G2, i2) < scale(G1, i1)
      push!(z, G2.data[i2])
      i2 += 1
    else
      @assert scale(G1, i1) == scale(G2, i2)
      push!(z, (scale(G1, i1), rank(G1, i1) + rank(G2, i2), det(G1, i1) * det(G2, i2)))
      i1 += 1
      i2 += 1
    end
  end

  if i1 <= length(G1)
    append!(z, G1.data[i1:length(G1)])
  end

  if i2 <= length(G2)
    append!(z, G2.data[i2:length(G2)])
  end

  return genus(HermLat, base_field(G1), prime(G1), z)
end

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

mutable struct GenusHerm{S, T, U, V}
  E::S
  primes::Vector{T}
  LGS::Vector{U}
  rank::Int
  signatures::V

  function GenusHerm(E::S, r, LGS::Vector{U}, signatures::V) where {S, U, V}
    K = base_field(E)
    primes = Vector{ideal_type(order_type(K))}(undef, length(LGS))

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

genus_herm_type(E) = GenusHerm{typeof(E), ideal_type(order_type(base_field(E))), local_genus_herm_type(E), Dict{InfPlc, Int}}

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

function Base.show(io::IO, ::MIME"text/plain", G::GenusHerm)
  print(io, "Global genus symbol over ")
  print(io, G.E)
  print(io, "\n", "with local genus symbols",)
  for g in G.LGS
    print(io, "\n")
    print(IOContext(io, :compact => true), prime(g), " => ", g)
  end
  print(io, "\n", "and signatures")
  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

function Base.show(io::IO, 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

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

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

Given a global genus symbol `G` for hermitian lattices over $E/K$, return `E`.
"""
base_field(G::GenusHerm) = G.E

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

Given a global genus symbol `G` for hermitian lattices over $E/K$, return
the list of prime ideals of $\mathcal O_K$ at which `G` has a local genus symbol.
"""
primes(G::GenusHerm) = G.primes

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

Given a global genus symbol `G` for hermitian lattices over $E/K$, return
the signatures at the infinite places of `K`. For each real place, it is
given by the negative index of inertia of the Gram matrix of the rational span of
a hermitian lattice whose global genus symbol is `G`.

The output is given as a dictionary with keys the infinite places of `K` and value
the corresponding signatures.
"""
signatures(G::GenusHerm) = G.signatures

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

Return the rank of any hermitian lattice with global genus symbol `G`.
"""
rank(G::GenusHerm) = G.rank

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

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

Given two global genus symbols `G1` and `G2` for hermitian lattices, return whether
they share the same local genus symbols. Note that equality requires `G1` and `G2`
to be defined over the same extension $E/K$.
"""
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

################################################################################
#
#  Sum of global genera
#
################################################################################

@doc Markdown.doc"""
    orthogonal_sum(G1::GenusHerm, G2::GenusHerm) -> GenusHerm

Given two global genus symbols `G1` and `G2` for hermitian lattices over $E/K$,
return their orthogonal sum. It corresponds to the global genus symbol of the
orthogonal sum of respective representatives of `G1` and `G2`.
"""
function orthogonal_sum(G1::GenusHerm, G2::GenusHerm)
  @req G1.E === G2.E "Genera must have same base field"
  E = G1.E
  LGS = local_genus_herm_type(G1.E)[]
  prim = eltype(primes(G1))[]
  P1 = Set(primes(G1))
  P2 = Set(primes(G2))
  for p in union(P1, P2)
    if p in P1
      i = findfirst(g -> prime(g) == p, G1.LGS)
      g1 = G1.LGS[i]
    else
      @assert !isramified(maximal_order(E), p)
      g1 = genus(HermLat,E , p, [(0, rank(G1), 1)])
    end

    if p in P2
      i = findfirst(g -> prime(g) == p, G2.LGS)
      g2 = G2.LGS[i]
    else
      @assert !isramified(maximal_order(E), p)
      g2 = genus(HermLat, E, p, [(0, rank(G2), 1)])
    end

    g3 = orthogonal_sum(g1, g2)
    push!(prim, p)
    push!(LGS, g3)
  end
  sig1 = G1.signatures
  sig2 = G2.signatures
  g3 = merge(+, sig1, sig2)
  return genus(LGS, g3)
end

Base.:(+)(G1::GenusHerm, G2::GenusHerm) = orthogonal_sum(G1, G2)

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

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

Return whether `G` and the global genus symbol of the hermitian lattice `L` agree.
"""
Base.in(L::HermLat, G::GenusHerm) = genus(L) == G

# This could be made more efficient

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

@doc Markdown.doc"""
    genus(S::Vector{LocalGenusHerm}, signatures) -> GenusHerm

Given a set of local genus symbols `S`  and a set of signatures `signatures`,
return the corresponding global genus symbol `G`. Note that the local genus symbols
in `S` must have the same rank which is also, therefore, the rank of `G`.

This constructor requires `S` to be non empty, the signatures to be non negative and
that the local invariants respect the product formula for Hilbert symbols.

`signatures` can be a dictionary with keys infinite places and values the
corresponding signatures, or a list of tuples $(::InfPlc, ::Int)$.
"""
genus(S::Vector{LocalGenusHerm}, signatures)

function genus(L::Vector{<:LocalGenusHerm}, 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))
  @req all(g -> rank(g) == r, L) "Local genus symbols must have the same rank"
  E = base_field(first(L))
  @req all(g -> base_field(g) == E, L) "Local genus symbols must be defined over the same extension E/K"
  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

Return the global genus symbol `G` of the hermitian lattice `L`. `G` satisfies:
- its local genus symbols correspond to those of the completions of `L` at the bad
  prime ideals of `L`, i.e. the prime ideals dividing either the scale of `L`, or the
  volume of `L`, or the discriminant of $\mathcal O_E$, and also the dyadic prime
  ideals of $\mathcal O_K$;
- signatures are those of the Gram matrix of the rational span of `L`. They are given
  at the real infinite places of `K` which split into complex places of `E`.
"""
function genus(L::HermLat)
  c = get_attribute(L, :genus)
  S = ideal_type(base_ring(base_ring(L)))
  if c === nothing
    G = _genus(L)
    set_attribute!(L, :genus => G)
    return G
  else
    return c::genus_herm_type(base_field(L))
  end
end

function _genus(L::HermLat)
  bad = bad_primes(L)
  for p in support(discriminant(base_ring(L)))
    if !(p in bad)
      push!(bad, p)
    end
  end

  for p in support(2 * base_ring(base_ring(L)))
    if !(p in bad)
      push!(bad, p)
    end
  end

  SE = infinite_places(base_field(L))
  # Only taking real places of K which split into complex places
  S = unique([r.base_field_place for r in SE if isreal(r.base_field_place) && !isreal(r)])
  D = diagonal(rational_span(L))
  # Only count the places with stay
  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

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

function _check_global_genus(LGS, signatures)
  _non_norm = _non_norm_primes(LGS, ignore_split = true)
  P = length(_non_norm)
  I = length([(s, N) for (s, N) in signatures if isodd(mod(N, 2))])
  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{LocalGenusHerm{S, T}}; ignore_split = false) where {S, T}
  z = T[]
  for g in LGS
    if ignore_split && issplit(g)
      continue
    end
    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)
  i === nothing && throw(error("No local genus symbol at $P"))
  return G.LGS[i]
end

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

function _hermitian_form_with_invariants(E, dim, P, N)
  K = base_field(E)
  R = maximal_order(K)
  @req all(n -> n in 0:dim, values(N)) "Number of negative entries is impossible"
  infinite_pl = [ p for p in real_places(K) if length(infinite_places(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
  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)
  D = elem_type(E)[]
  for i in 1:(dim - 1)
    if length(I) == 0
      push!(D, one(E))
    else
      el = E(_weak_approximation(infinite_pl, [N[p] >= i ? -1 : 1 for p in infinite_pl]))
      push!(D, el)
    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
  v = involution(E)

  @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 length(infinite_places(E, p)) == 1])
  return ncols(M), collect(keys(P)), I
end

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

Given a global genus symbol `G` for hermitian lattices over $E/K$, return a hermitian
lattice over $E/K$ which admits `G` as global genus symbol.
"""
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\n"
  M = maximal_integral_lattice(V)
  lp = G.primes
  for g in G.LGS
    p = prime(g)
    @vprint :Lattice 1 "Finding representative for $g at $(prime(g))...\n"
    L = representative(g)
    @hassert :Lattice 1 genus(L, p) == g
    @vprint :Lattice 1 "Finding sublattice\n"
    M = locally_isometric_sublattice(M, L, p)
  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 genus symbols for hermitian lattices over the algebra `E`, with base
field $K$, at the prime ideal`p` of $\mathcal O_K$. Each of them has rank equal to
`rank`, scale $\mathfrak P$-valuations bounded by `max_scale` and determinant
`p`-valuations equal to `det_val`, where $\mathfrak P$ is a prime ideal of $\mathcal O_E$
lying above `p`.
"""
function local_genera_hermitian(E, p, rank::Int, det_val::Int, max_scale::Int)
  is_ramified = isramified(maximal_order(E), p)
  is_inert = !is_ramified && length(prime_decomposition(maximal_order(E), p)) == 1
  is_split = !is_ramified && !is_inert
  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))
  @req base_field(E) == K "Incompatible arguments: p must be a prime ideal in the base field of E"

  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}[]
      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))
        else
          push!(z, (b[1], b[2], -1))
        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)
  hyperbolic_det = islocal_norm(E, K(-1), p) ? 1 : -1
  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 cartesian_product_iterator(dets, inplace = true)# Iterators.product(dets...)
        g2 = Vector{Tuple{Int, Int, Int}}(undef, length(g))
        for k in 1:n
          # Again the 0 for dummy purposes
          g2[k] = (g[k]..., d[k])
        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
        k = ceil(Int, i//2)
        while 2 * k < i + e
          push!(dn, Int[1, k])
          push!(dn, Int[-1, k])
          k += 1
        end

        if mod(i + e, 2) == 0
          push!(dn, Int[hyperbolic_det^(div(b[2], 2)), k])
        end
        push!(det_norms, dn)
      end
    end

    for dn in cartesian_product_iterator(det_norms, inplace = true)
      g2 = Vector{Tuple{Int, Int, Int, Int}}(undef, length(g))
      for k in 1:n
        g2[k] = (g[k][1], g[k][2], dn[k][1], dn[k][2])
      end
      h = genus(HermLat, E, p, g2)
      if !(h in symbols)
        push!(symbols, h)
      end
    end
  end
  return symbols
end

@doc Markdown.doc"""
    genera_hermitian(E::NumField, rank::Int,
                                  signatures::Dict{InfPlc, Int},
                                  determinant::nf_elem,
                                  max_scale = nothing) -> Vector{GenusHerm}

Return all global genus symbols for hermitian lattices over the algebra`E` with rank
`rank`, signatures given by `signatures`, scale bounded by `max_scale` and determinant
class equal to `determinant`.

If `max_scale == nothing`, it is set to be equal to `determinant`.
"""
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{local_genus_herm_type(E)}[]

  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 = genus_herm_type(E)[]
  it = cartesian_product_iterator(local_symbols, inplace = true)
  for gs in it
    c = copy(gs)
    b = _check_global_genus(c, signatures)
    if b
      push!(res, GenusHerm(E, rank, c, signatures))
    end
  end

  return res
end


text:
using DrWatson
@quickactivate("HyperspectraWithNeXL")
# Load the necessary libraries
using HyperspectraWithNeXL
using NeXLSpectrum
using Gadfly
using Images, FileIO, Unitful
# Load the HyperSpectrum from disk
lt = 0.72*4.0*18.0*3600.0/(1024*1024) # 18.0 hours on 4 detectors
hs = NeXLSpectrum.compress(HyperSpectrum(
    LinearEnergyScale(0.0,10.0),
    Dict{Symbol,Any}(
      :TakeOffAngle => deg2rad(35.0),
      :ProbeCurrent => 1.0,
      :LiveTime => lt,
      :BeamEnergy => 20.0e3,
      :Name => "Mn Nodule"
    ),
    readrplraw(joinpath(datadep"MnNodule","map[15]")),
    fov = [ 4.096u"mm", 4.096u"mm"], offset= [ 0.0u"mm", 0.0u"mm" ]
))
# hs = hs[16:8:1024, 16:8:1024]
# Build the fitting references
refpath = datadep"MnNodule_Standards"
refs = references( [
    reference(n"C", joinpath(refpath, "C std.msa") ),
    reference(n"Ag", joinpath(refpath, "Ag std.msa") ),
    reference(n"Al", joinpath(refpath, "Al std.msa") ),
    reference(n"C", joinpath(refpath, "C std.msa") ),
    reference(n"Ca", joinpath(refpath, "CaF2 std.msa") ),
    reference(n"Ce", joinpath(refpath, "CeO2 std.msa") ),
    reference(n"Cl", joinpath(refpath, "NaCl std.msa") ),
    reference(n"Cr", joinpath(refpath, "Cr std.msa") ),
    reference(n"Cu", joinpath(refpath, "Cu std.msa") ),
    reference(n"Fe", joinpath(refpath, "Fe std.msa") ),
    reference(n"S", joinpath(refpath, "FeS2 std.msa") ),
    reference(n"P", joinpath(refpath, "GaP std.msa") ),
    reference(n"K", joinpath(refpath, "KBr std.msa") ),
    reference(n"Mg", joinpath(refpath, "Mg std.msa") ),
    reference(n"O", joinpath(refpath, "MgO std.msa") ),
    reference(n"Mn", joinpath(refpath, "Mn std.msa") ),
    reference(n"Na", joinpath(refpath, "NaCl std.msa") ),
    reference(n"Ni", joinpath(refpath, "Ni std.msa") ),
    reference(n"Si", joinpath(refpath, "Si std.msa") ),
    reference(n"Ti", joinpath(refpath, "Ti std.msa") ),
    reference(n"Zn", joinpath(refpath, "Zn std.msa") ) ],
    132.0
)
# Fit the reference spectra to the HyperSpectrum unknown
resf = fit_spectrum(hs, refs, mode = :Fast)
# Select the optimal k-ratio (K, L? or M?) for each element
bestks = optimizeks(SimpleKRatioOptimizer(2.0), resf)
# Convert compound k-ratios to pure element equivalents (matrix correction)
pureks = aspure.(bestks)
# Normalize the k-ratios at each point in the HyperSpectrum
bestf = normalizek(pureks)

# Save the resulting images
path = joinpath(plotsdir(),"Table K-ratio Map")
mkpath(path)
for bf in bestf
  FileIO.save(File{format"PNG"}(joinpath(path,"k_norm[$(bf.xrays)].png")),Log3Band.(bf.kratios))
end

li = labeledimages([ "$(bf.xrays)" for bf in bestf],[Log3Band.(bf.kratios)
                for bf in bestf ], ncols=4)

li |> SVG(joinpath(plotsdir(),"Labeled k-ratio maps.svg"), 8inch, 10inch)

# Clean up
hs, resf, bestks, pureks, bestf = nothing, nothing, nothing, nothing, nothing

ext: jl


text:
# Copyright (c) 2020-2021, Bank of Canada
# All rights reserved.

@testset "MV construct" begin
    @test size(MVTSeries(20Q1)) == (0, 0)
    @test size(MVTSeries(20Q1, :a)) == (0, 1)
    @test size(MVTSeries(20Q1, (:a, :b))) == (0, 2)
    @test size(MVTSeries(20Q1, (:a, "b"))) == (0, 2)
    @test size(MVTSeries(20Q1:20Q4, :a)) == (4, 1)
    @test size(MVTSeries(20Q1:20Q4, (:a, :b))) == (4, 2)
    @test size(MVTSeries(20Q1:20Q4, ("a", :b))) == (4, 2)
    @test size(MVTSeries(20Q1:20Q4, ["a", :b])) == (4, 2)
    @test size(MVTSeries(20Q1:20Q4, ["a", :b], undef)) == (4, 2)
    @test size(MVTSeries(20Q1:20Q4, ["a", :b], 5)) == (4, 2)
    @test size(MVTSeries(20Q1:20Q4, ["a", :b], zeros)) == (4, 2)

    let a = similar(zeros(Int, 0, 0), (20Y:22Y, (:a, :b))),
        b = similar(zeros(Int, 0, 0), Complex, (20Y:22Y, (:a, :b))),
        c = fill(5, 20Q1:20Q4, (:a, :b, :c))

        @test a isa MVTSeries
        @test rangeof(a) == 20Y:22Y
        @test tuple(colnames(a)...) == (:a, :b)
        @test eltype(a) == Int
        @test eltype(b) == Complex
        @test eltype(c) == Int
        @test size(c) == (4, 3)
        @test extrema(c) == (5, 5)
        @test typeof(similar(c)) === typeof(c)
    end

    @test_throws ArgumentError MVTSeries(1U, (:a, :b), zeros(10, 3))
    @test_throws ArgumentError MVTSeries(1U, (:a, :b), zeros(10, 1))

    @test_throws ArgumentError MVTSeries(1U:5U, (:a, :b), zeros(10, 2))
    @test_throws ArgumentError MVTSeries(1U:5U, (:a, :b), zeros(4, 2))
    @test (MVTSeries(1U:5U, (:a, :b), zeros(5, 2)); true)

    @test (MVTSeries(20Q1, ("a",), zeros(5,)); true)
    @test (MVTSeries(20Q1, "a", zeros(5,)); true)
    @test (MVTSeries(1U:5U, ("a",), zeros(5,)); true)

    @test (MVTSeries(20Q1, :a, zeros(5,)); true)
    @test (MVTSeries(1U:5U, :a, zeros(5,)); true)

end

@testset "MV Int Ind" begin
    a = MVTSeries(20Q1, (:a, :b), rand(5, 2))
    for i = 1:10
        x = a[i]
        a[i] = i
    end
    @test a[:] == a.values[:]
    @test a[:] == 1:10
    @test a[3:5] == 3:5
    @test a[:, 1] == 1:5
    @test a[1, :] == [1, 6]
    @test a[1:3, :] == [1 6; 2 7; 3 8]
    @test a[:, 1:2] == a.values[:, 1:2]
    @test a[:, :] == a.values
    for i = 1:5
        for j = 1:2
            x = a[i, j]
            a[i, j] = 10i + j
        end
    end
    @test a[:] == [11, 21, 31, 41, 51, 12, 22, 32, 42, 52]

    b = rand(Bool, size(a))
    @test a[b[:, 1], 1] == a.values[b[:, 1], 1]
    @test a[1, b[1, :]] == a.values[1, b[1, :]]
    @test a[b] == a.values[b]

end

@testset "MV dot" begin
    a = MVTSeries(20Q1, [:a, "b"], rand(10, 2))
    @test a.a isa TSeries
    a.a[20Q1] = 22
    @test a[1, 1] == 22
    @test_throws BoundsError a.c
    @test (a.a = TSeries(20Q1, 1:10); a.values[:, 1] == collect(1:10))
end

@testset "MV" begin
    @test_throws ArgumentError MVTSeries(1M10, (:a, :b, :c), rand(10, 2))
    let nms = (:a, :b), dta = rand(20, 2),
        sd = MVTSeries(2000Q1, nms, copy(dta)),
        dta2 = rand(size(dta)...)

        # if one argument is Colon, fall back on single argument indexing
        # getindex
        @test sd[2000Q1,:] == dta[1,:]
        @test all(sd[2000Q1:2000Q4,:].values == dta[1:4,:])
        @test sd[:,:a].values == dta[:,1]
        @test sd[:,(:a,:b)].values == dta[:,1:2]
        @test sd[:,[:a,:b]].values == dta[:,1:2]

        @test firstdate(sd) == 2000Q1
        @test lastdate(sd) == 2000Q1 + 20 - 1
        @test frequencyof(sd) == Quarterly
        @test sd isa AbstractMatrix
        @test size(sd) == size(dta)
        # integer indexing must be identical to dta
        for i in axes(dta, 2)
            @test sd[:, i] == dta[:, i]
        end
        for j in axes(dta, 1)
            @test sd[j, :] == dta[j, :]
        end
        for i in eachindex(dta)
            @test sd[i] == dta[i]
        end
        # set indexing with integers
        for i in eachindex(dta2)
            sd[i] = dta2[i]
        end
        @test sd[:] == dta2[:]
        for i in axes(dta, 1)
            sd[i, :] = dta[i, :]
        end
        @test sd[:] == dta[:]
        for i in axes(dta2, 2)
            sd[:, i] = dta2[:, i]
        end
        @test sd[:] == dta2[:]
        # access by dot notation (for the columns)
        @test propertynames(sd) == nms
        @test sd.a isa TSeries
        @test sd.b isa TSeries
        @test_throws BoundsError sd.c
        @test sd.a.values == dta2[:, 1]
        sd.a[:] = dta[:, 1]
        sd.b[:] = dta[:, 2]
        @test sd[:] == dta[:]
        @test sd[:a] isa TSeries && sd[:a].values == dta[:, 1]
        @test sd["b"] isa TSeries && sd["b"].values == dta[:, 2]
        #
        sd.a = dta[:, 1]
        sd.b = dta[:, 2]
        @test sd[:] == dta[:]
        #
        sd[:] = dta[:]
        sd.a = sd.b
        @test sd[:, 1] == sd[:, 2]
        sd.a = zeros(size(dta, 1))
        @test sum(abs, sd.a.values) == 0
        @test_throws DimensionMismatch sd.a[:] = ones(length(sd.a) + 5)
        @test_throws BoundsError sd.a = ones(length(sd.a) + 5)
        # access to rows by MIT
        sd[:] = dta[:]
        @test sd[2000Q1] isa Vector{Float64}
        for (i, idx) in enumerate(rangeof(sd))
            @test sd[idx] == dta[i, :]
        end
        sd[2000Q1] = zeros(size(dta, 2))
        @test sum(abs, sd[1, :]) == 0
        @test_throws BoundsError sd[1999Q4]
        @test_throws DimensionMismatch sd[2000Q3] = zeros(size(dta, 2) + 1)
        @test_throws ArgumentError sd[2000Y]
        @test_throws BoundsError sd[1999Q4] = [5.0, 1.1]
        @test_throws BoundsError sd[2010Q4] = rand(2)
        #
        @test sd[2001Q1:2002Q4] isa MVTSeries
        @test sd[2001Q1:2002Q4].values == sd[5:12, :]
        sd[2001Q1:2002Q4] = 1:16
        @test sd[5:12, :] == reshape(1.0:16.0, :, 2)
        @test_throws BoundsError sd[1111Q4:2002Q1]
        @test_throws BoundsError sd[2002Q1:2200Q2]
        @test_throws DimensionMismatch sd[2001Q1:2002Q4] = 1:17
        # assign new column
        sd1 = hcat(sd, c = sd.a .+ 3.0)
        @test sd1[nms] == sd
        @test sd1[(:a, :c)].values == sd1[:, [1, 3]]
        # access with 2 args MIT and Symbol
        @test sd[2001Q2, (:a, :b)] isa Vector{eltype(sd)}
        let foo = sd[2001Q2:2002Q1, (:a, :b)]
            @test foo isa MVTSeries
            @test size(foo) == (4, 2)
            @test firstdate(foo) == 2001Q2
        end
        # access with an MIT range and a Symbol returns a TSeries
        let foo = sd[2001Q2:2002Q1, :a]
            @test foo isa TSeries
            @test size(foo) == (4, )
            @test firstdate(foo) == 2001Q2
            @test foo.values == sd[2001Q2:2002Q1, :a].values
        end
        @test_throws BoundsError sd[1999Q1, (:a,)]
        @test_throws BoundsError sd[2001Q1:2001Q2, (:a, :c)]
        @test_throws Exception sd.c = 5
        @test similar(sd) isa typeof(sd)
        @test_throws BoundsError sd[1999Q1:2000Q4] = zeros(8, 2)
        @test_throws BoundsError sd[2004Q1:2013Q4, [:a, :b]]
        # setindex with two
        @test_throws ArgumentError sd[2001Q1, (:b,)] = 3.5
        sd[2001Q1, :b] = 3.5
        @test sd[5, 2] == 3.5
        @test_throws ArgumentError sd[2000Q1:2001Q4, (:b,)] = 3.7
        sd[2000Q1:2001Q4, (:b,)] = fill(3.7, 8)
        @test all(sd[1:8, 2] .== 3.7)
        @test_throws BoundsError sd[1999Q1:2000Q4, (:a, :b)] = 5.7
        @test_throws BoundsError sd[2000Q1:2000Q4, (:a, :c)] = 5.7

        # getindex mixed_freq_error
        @test_throws ArgumentError sd[1U,:a]
        @test_throws ArgumentError sd[1U:5U,:b]

        # setindex mixed_freq_error
        @test_throws ArgumentError sd[1U,:a] = 5
        @test_throws ArgumentError sd[1U:5U,:a] = 5

        # if one argument is Colon, fall back on single argument indexing
        # setindex
        @test (myvar = [2,2]; sd[2000Q1,:] = myvar; sd[2000Q1,:] == myvar)
        @test (myvar = [1 2;3 4]; sd[2000Q1:2000Q2,:] = myvar; sd[2000Q1:2000Q2,:].values == myvar)
        @test (sd[:,:a] = dta[:,1]; sd[:,:a].values == dta[:,1])
        @test (myvar = rand(size(sd)...); sd[:,(:a,:b)] = myvar; sd[:,(:a,:b)].values == myvar)
        @test (myvar = rand(size(sd)...); sd[:,[:a,:b]] = myvar; sd[:,[:a,:b]].values == myvar)

        # with a range of MIT and a single column, we fall back on TSeries assignment
        @test (myvar = [1,2,3,4]; sd[2000Q1:2000Q4,:a] = myvar; sd[2000Q1:2000Q4,:a].values == myvar)

        # setindex from an MVTSeries to an MVTSeries
        let sd2 = MVTSeries(2001Q1, nms, rand(8, length(nms)))
            sd[2001Q1:2001Q4,[:a, :b]] = sd2
            @test (sd[2001Q1:2001Q4,:].values == sd2[2001Q1:2001Q4,:].values)
        end
    end
end

@testset "MVTSeries show" begin
    # test the case when column labels are longer than the numbers
    let io = IOBuffer(), x = MVTSeries(1U, (:verylongandsuperboringnameitellya, :anothersuperlongnamethisisridiculous, :a), rand(20, 3) .* 100)
        show(io, x)
        lines = readlines(seek(io, 0))
        # labels longer than 10 character are abbreviated with '…' at the end
        @test length(split(lines[2], '…')) == 3
    end
    let io = IOBuffer(), x = MVTSeries(1U, (:alpha, :beta), zeros(24, 2))
        show(io, x)
        lines = readlines(seek(io, 0))
        lens = length.(lines)
        # when labels are longer than the numbers, the entire column stretches to fit the label
        @test lens[2] == lens[3]
    end
    nrow = 24
    letters = Symbol.(['a':'z'...])
    for (nlines, fd) in zip([3, 4, 5, 6, 7, 8, 22, 23, 24, 25, 26, 30], Iterators.cycle((2010Q1, 2010M1, 2010Y, 1U)))
        for ncol in [2, 5, 10, 20]
            # display size if nlines × 80
            # data size is nrow × ncol
            # when printing data there are two header lines - summary and column names
            io = IOBuffer()
            x = MVTSeries(fd, tuple(letters[1:ncol]...), rand(nrow, ncol))
            show(IOContext(io, :displaysize => (nlines, 80)), MIME"text/plain"(), x)
            lines = readlines(seek(io, 0))
            @test length(lines) == max(3, min(nrow + 2, nlines - 3))
            @test maximum(length, lines) <= 80
            io = IOBuffer()
            show(IOContext(io, :limit => false), x)
            lines = readlines(seek(io, 0))
            @test length(lines) == nrow + 2
        end
    end
end

@testset "MV bcast" begin
    x = MVTSeries(20Q1, (:a, :b), rand(10,2))

    # we can do dot operations with scalars
    @test (x .+ 2; true)

    # we can do dot operations with multiple MVTSeries of the same dimension
    @test x .+ 2x ./ 2 == 2x

    # we can do dot operations with multiple MVTSeries of different ranges
    @test begin
        q = x .+ x[begin+1:end-1, :] ./ 2
        q == 1.5x[begin+1:end-1, :]
    end

    # we can do dot operations with matrices of identical size
    @test x .+ 3 .* x.values == 4x

    # with matrixes of wrong size we get a DimensionMismatch
    @test_throws DimensionMismatch x .+ x.values[begin+1:end-1, :]
    @test_throws DimensionMismatch x .+ x.values[:, 1:1]

    t = TSeries(20Q1, ones(10))
    @test x .+ t.values  == (x .+ 1)
    @test x .+ t  == (x .+ 1)
    #
    @test x .+ 1 .* 3  == (x .+ 3)
    @test 1 .* 3 .+ x  == (x .+ 3)
    @test x .+ t.values .* 3  == (x .+ 3)
    @test t.values .* 3 .+ x  == (x .+ 3)
    @test x .+ t .* 3  == (x .+ 3)
    @test t .* 3 .+ x  == (x .+ 3)

end


text:
#
# This file is a part of MolecularGraph.jl
# Licensed under the MIT License http://opensource.org/licenses/MIT
#

export
    hasselfloop, hasmultiedge, issimplegraph

@cache function hasselfloop(graph::Union{OrderedGraph,OrderedDiGraph})
    for (u, v) in edgesiter(graph)
        u == v && return true
    end
    return false
end

@cache function hasmultiedge(graph::Union{OrderedGraph,OrderedDiGraph})
    for nbr in neighborsiter(graph)
        length(nbr) == length(Set(values(nbr))) || return true
    end
    return false
end

@cache issimplegraph(graph::Union{OrderedGraph,OrderedDiGraph}
    ) = !hasselfloop(graph) && !hasmultiedge(graph)


text:
using RigidBodyDynamics
using RigidBodyDynamics: Bounds
using MeshCat
using MeshCatMechanisms
using LCPSim
using JuMP, Gurobi
using Base.Test

@testset "box in a box" begin
    urdf = joinpath(@__DIR__, "urdf", "box_in_a_box.urdf")
    mechanism = parse_urdf(Float64, urdf)
    planar_joint = findjoint(mechanism, "floor_to_box")
    planar_joint.position_bounds .= [
        Bounds(-3, 3),
        Bounds(-3, 3),
        Bounds(-2π, 2π)
    ]
    planar_joint.velocity_bounds .= Bounds(-1000, 1000)

    vis = Visualizer()
    if !haskey(ENV, "CI")
        open(vis)
        wait(vis)
    end
    mvis = MechanismVisualizer(mechanism, URDFVisuals(urdf), vis[:robot])
    env = LCPSim.parse_contacts(mechanism, urdf, 1.0, :xz)
    controller = LCPSim.passive_controller()
    Δt = 0.02
    N = 100
    solver = GurobiSolver(Gurobi.Env(), OutputFlag=0)

    for i in 1:50
        @show i
        srand(i)
        x0 = MechanismState{Float64}(mechanism)
        set_velocity!(x0, zeros(num_velocities(x0)))
        set_configuration!(x0, planar_joint, [rand() - 0.5, -1.0, randn()])
        set_velocity!(x0, planar_joint, 2 * randn(3))
        results = LCPSim.simulate(x0, controller, env, Δt, N, solver);
        for r in results
            set_configuration!(mvis, configuration(r.state))
            sleep(0.1 * Δt)
        end
        @test length(results) == N
    end
end


text:
#
# A single execution of this file generates a single pair (a(r),r).
# The pair is then stored in "LS_pureIOU_2000_narrow_saved.jld2", together with
# all previously computed pairs.
#      If running for the first time (no "LS_pureIOU_2000_narrow_saved.jld2" exists)
#      uncomment lines 24-29. Put the #'s back for all subsequent runs.
#
include("LS-method-7s-init.jl");

# Assume the bounds on the wealth generated by the dual method and create the grid.
#
w_bounds=(-1.6292720455478487, 7.18257249498523);
begin
    Lg=2000 #the only place to change the grid size
    grda=[w_bounds[1]+i*(w_bounds[end]-w_bounds[1])/(Lg-1) for i=0:(Lg-1)]
    β=exp(-ρ)
end;

## Uncomment and execute the next block the very *first* time,
## *only* if "LS_pureIOU_2000_narrow_saved.jld2" does not exist in the current path.
## Put the #'s back for all subsequent runs.

#begin
#    a_list=Array{Float64}(undef, 0)
#    r_list=Array{Float64}(undef, 0)
#    LS_pureIOU_2000_narrow_saved=[r_list,a_list]
#    @save "LS_pureIOU_2000_narrow_saved.jld2" LS_pureIOU_2000_narrow_saved
#end;

##
## The next block guesses the next candidate for r.
## It can be bypassed by setting r=<value> explicitly. See line 68.
##

begin
    @load "LS_pureIOU_2000_narrow_saved.jld2" LS_pureIOU_2000_narrow_saved;
    r_list=LS_pureIOU_2000_narrow_saved[1]
    a_list=LS_pureIOU_2000_narrow_saved[2]
    if length(r_list)==0
        r=0.03
        elseif length(r_list)==1
        r=0.02
        else
        r=0.5*(r_list[findall(x->x<0,a_list)[end]]+r_list[findall(x->x>0,a_list)[end]])
    end
    push!(r_list,r)
end;

# the next line could take up to an hour
# arguments 2 and 3are the requested digits of convergence
#
@time Λ1,tbl2=LS_method(r,5,5,nos,hours,wage,β,grda,Lg,CPROB,R);

## Two important checkpoints!
sum(Λ1) # should print 1
[sum(Λ1[i,:]) for i=1:nos]-BSP # should print a list of zeros

begin
    policy=[grda[Int64(tbl2[kk,ii])] for kk=1:nos, ii=1:Lg]
    a=sum(Λ1.*policy) # the aggregate demand a(r)
end;

##
## If running with an explicit r (see line 32)
## be sure that a_list and r_list exist before running the next block.
##

begin
    a_list=push!(a_list,a)
    LS_pureIOU_2000_narrow_saved=[r_list,a_list,Λ1,tbl2]
    @save "LS_pureIOU_2000_narrow_saved.jld2" LS_pureIOU_2000_narrow_saved
end;


text:
#TODO add more annotation
println("This file demos some of the mathematical operations on Associative Array")

A = ReadCSV("A.csv")

A = logical(A)

print(A)

print(sum(A,1))

print(sum(A,2))

Aa = A[:,"a,"]

Ab = A[:,"b,"]

Aab = nocol(Aa) & nocol(Ab)

F =  ( Aab )' * A;

print(F)

Fn = F ./ sum(A,1)

print(Fn)

AtA = sqIn(A)

d = diag(adj(AtA))

AtA.A = AtA.A - sparse(diagm(d))

print(AtA)


text:
## A test might be comprised of several pg file so that each question
## can be displayed on its own page. However, spreading a test out over
## many scripts might be undesirable. This example shows how to use the
## `JuliaWeBWorK.PAGE` function to create connected pg files
## using one basic script.

using JuliaWeBWorK
PAGE = JuliaWeBWorK.PAGE(@__FILE__)

q = numericq(L"What is \({{:a1}} + {{:a2}}\)?", (a,b) -> a+b, (1:4, 2:5))
PAGE("Addition", (q,)) # writes to SCRIPT_BASE_NAME-1.pg

q = numericq(L"What is \({{:a1}} - {{:a2}}\)?", (a,b) -> a-b, (1:4, 2:5))
PAGE("subtraction", (q,)) # writes to SCRIPT_BASE_NAME-2.pg

q = numericq(L"What is \({{:a1}} * {{:a2}}\)?", (a,b) -> a*b, (1:4, 2:5))
PAGE("multiplication", (q,)) # writes to SCRIPT_BASE_NAME-3.pg


text:
using Documenter, DocumenterTools

# The DOCSARGS environment variable can be used to pass additional arguments to make.jl.
# This is useful on CI, if you need to change the behavior of the build slightly but you
# can not change the .travis.yml or make.jl scripts any more (e.g. for a tag build).
if haskey(ENV, "DOCSARGS")
    for arg in split(ENV["DOCSARGS"])
        push!(ARGS, arg)
    end
end

makedocs(
    modules = [Documenter, DocumenterTools],
    format = Documenter.HTML(
        # Use clean URLs, unless built as a "local" build
        prettyurls = !("local" in ARGS),
        canonical = "https://juliadocs.github.io/Documenter.jl/stable/",
        assets = ["assets/favicon.ico"],
        analytics = "UA-136089579-2",
    ),
    clean = false,
    sitename = "Documenter.jl",
    authors = "Michael Hatherly, Morten Piibeleht, and contributors.",
    linkcheck = !("skiplinks" in ARGS),
    pages = [
        "Home" => "index.md",
        "Manual" => Any[
            "Guide" => "man/guide.md",
            "man/examples.md",
            "man/syntax.md",
            "man/doctests.md",
            "man/latex.md",
            hide("man/hosting.md", [
                "man/hosting/walkthrough.md"
            ]),
            "man/other-formats.md",
        ],
        "Library" => Any[
            "Public" => "lib/public.md",
            hide("Internals" => "lib/internals.md", Any[
                "lib/internals/anchors.md",
                "lib/internals/builder.md",
                "lib/internals/cross-references.md",
                "lib/internals/docchecks.md",
                "lib/internals/docmeta.md",
                "lib/internals/docsystem.md",
                "lib/internals/doctests.md",
                "lib/internals/documenter.md",
                "lib/internals/documentertools.md",
                "lib/internals/documents.md",
                "lib/internals/dom.md",
                "lib/internals/expanders.md",
                "lib/internals/markdown2.md",
                "lib/internals/mdflatten.md",
                "lib/internals/selectors.md",
                "lib/internals/textdiff.md",
                "lib/internals/utilities.md",
                "lib/internals/writers.md",
            ])
        ],
        "contributing.md",
    ],
    strict = !("strict=false" in ARGS),
    doctest = ("doctest=only" in ARGS) ? :only : true,
)

deploydocs(
    repo = "github.com/JuliaDocs/Documenter.jl.git",
    target = "build",
)


text:
function write_results(doc::String, it_results::Array{Tuple{Number, Number, String, String}, 1})
    output = "../../results"*doc[21:end-4]*"_results.txt"
    open(output, "w") do o

        println(o, string("makespan time indicator criteria"))
        for it in it_results
            println(o, string(it[1])*" "*string(it[2])*" "*string(it[3])*" "*string(it[4]))
        end

    end
end

function write_time_results(doc::String, exec_time::Number)
    output = "../../results"*doc[21:end-4]*"_time_results.txt"
    open(output, "w") do o
        println(o, string(exec_time))
    end
end

max_stars_repo_licenses: Zlib MIT

max_issues_repo_licenses: Zlib MIT

max_forks_repo_licenses: Zlib MIT


text:
# Create and return VTK grids that are ready to be filled with data (vtu version)
function build_vtk_grids(::Val{:vtu}, coordinates, levels, center_level_0, length_level_0,
                         n_visnodes, verbose, output_directory, is_datafile, filename)
  # Extract number of spatial dimensions
  ndims_ = size(coordinates, 1)

  # Calculate VTK points and cells
  verbose && println("| Preparing VTK cells...")
  if is_datafile
    @timeit "prepare VTK cells (node data)" begin
      vtk_points, vtk_cells = calc_vtk_points_cells(Val(ndims_), coordinates, levels,
                                                    center_level_0, length_level_0, n_visnodes)
    end
  end

  # Prepare VTK points and cells for celldata file
  @timeit "prepare VTK cells (cell data)" begin
    vtk_celldata_points, vtk_celldata_cells = calc_vtk_points_cells(Val(ndims_), coordinates,
                                                                    levels, center_level_0,
                                                                    length_level_0, 1)
  end

  # Determine output file names
  base, _ = splitext(splitdir(filename)[2])
  vtk_filename = joinpath(output_directory, base)
  vtk_celldata_filename = vtk_filename * "_celldata"

  # Open VTK files
  verbose && println("| Building VTK grid...")
  if is_datafile
    @timeit "build VTK grid (node data)" vtk_nodedata = vtk_grid(vtk_filename, vtk_points,
                                                                 vtk_cells)
  else
    vtk_nodedata = nothing
  end
  @timeit "build VTK grid (cell data)" vtk_celldata = vtk_grid(vtk_celldata_filename,
                                                                vtk_celldata_points,
                                                                vtk_celldata_cells)

  return vtk_nodedata, vtk_celldata
end

# Create and return VTK grids that are ready to be filled with data (vti version)
function build_vtk_grids(::Val{:vti}, coordinates, levels, center_level_0, length_level_0,
                         n_visnodes, verbose, output_directory, is_datafile, filename)
    # Extract number of spatial dimensions
    ndims_ = size(coordinates, 1)

    # Prepare VTK points and cells for celldata file
    @timeit "prepare VTK cells" vtk_celldata_points, vtk_celldata_cells = calc_vtk_points_cells(
        Val(ndims_), coordinates, levels, center_level_0, length_level_0, 1)

    # Determine output file names
    base, _ = splitext(splitdir(filename)[2])
    vtk_filename = joinpath(output_directory, base)
    vtk_celldata_filename = vtk_filename * "_celldata"

    # Open VTK files
    verbose && println("| Building VTK grid...")
    if is_datafile
      # Determine level-wise resolution
      max_level = maximum(levels)
      resolution = n_visnodes * 2^max_level

      Nx = Ny = resolution + 1
      dx = dy = length_level_0/resolution
      origin = center_level_0 .- 1/2 * length_level_0
      spacing = [dx, dy]
      @timeit "build VTK grid (node data)" vtk_nodedata = vtk_grid(vtk_filename, Nx, Ny,
                                                          origin=origin,
                                                          spacing=spacing)
    else
      vtk_nodedata = nothing
    end
    @timeit "build VTK grid (cell data)" vtk_celldata = vtk_grid(vtk_celldata_filename,
                                                        vtk_celldata_points,
                                                        vtk_celldata_cells)

  return vtk_nodedata, vtk_celldata
end

# Determine and return filenames for PVD fiels
function pvd_filenames(filenames, pvd, output_directory)
  # Determine pvd filename
  if !isnothing(pvd)
    # Use filename if given on command line
    filename = pvd

    # Strip of directory/extension
    filename, _ = splitext(splitdir(filename)[2])
  else
    filename = get_pvd_filename(filenames)

    # If filename is empty, it means we were not able to determine an
    # appropriate file thus the user has to supply one
    if filename == ""
      error("could not auto-detect PVD filename (input file names have no common prefix): " *
            "please provide a PVD filename name with the keyword argument `pvd=...`")
    end
  end

  # Get full filenames
  pvd_filename = joinpath(output_directory, filename)
  pvd_celldata_filename = pvd_filename * "_celldata"

  return pvd_filename, pvd_celldata_filename
end

# Determine filename for PVD file based on common name
function get_pvd_filename(filenames::AbstractArray)
  filenames = getindex.(splitdir.(filenames), 2)
  bases = getindex.(splitext.(filenames), 1)
  pvd_filename = longest_common_prefix(bases)
  return pvd_filename
end

# Determine longest common prefix
function longest_common_prefix(strings::AbstractArray)
  # Return early if array is empty
  if isempty(strings)
    return ""
  end

  # Count length of common prefix, by ensuring that all strings are long enough
  # and then comparing the next character
  len = 0
  while all(length.(strings) .> len) && all(getindex.(strings, len+1) .== strings[1][len+1])
    len +=1
  end

  return strings[1][1:len]
end

# Convert coordinates and level information to a list of points and VTK cells (2D version)
function calc_vtk_points_cells(::Val{2}, coordinates::AbstractMatrix{Float64},
                               levels::AbstractVector{Int},
                               center_level_0::AbstractVector{Float64},
                               length_level_0::Float64,
                               n_visnodes::Int=1)
  ndim = 2
  # Create point locator
  pl = PointLocator{2}(center_level_0, length_level_0, 1e-12)

  # Create arrays for points and cells
  n_elements = length(levels)
  points = Vector{SVector{2, Float64}}()
  vtk_cells = Vector{MeshCell}(undef, n_elements * n_visnodes^ndim)
  point_ids = Vector{Int}(undef, 2^ndim)

  # Reshape cell array for easy-peasy access
  reshaped = reshape(vtk_cells, n_visnodes, n_visnodes, n_elements)

  # Create VTK cell for each Trixi element
  for element_id in 1:n_elements
    # Extract cell values
    cell_x = coordinates[1, element_id]
    cell_y = coordinates[2, element_id]
    cell_dx = length_level_0 / 2^levels[element_id]

    # Adapt to visualization nodes for easy-to-understand loops
    dx = cell_dx / n_visnodes
    x_lowerleft = cell_x - cell_dx/2 - dx/2
    y_lowerleft = cell_y - cell_dx/2 - dx/2

    # Create cell for each visualization node
    for j = 1:n_visnodes
      for i = 1:n_visnodes
        # Determine x and y
        x = x_lowerleft + i * dx
        y = y_lowerleft + j * dx

        # Get point id for each vertex
        point_ids[1] = insert!(pl, points, (x - dx/2, y - dx/2))
        point_ids[2] = insert!(pl, points, (x + dx/2, y - dx/2))
        point_ids[3] = insert!(pl, points, (x - dx/2, y + dx/2))
        point_ids[4] = insert!(pl, points, (x + dx/2, y + dx/2))

        # Add cell
        reshaped[i, j, element_id] = MeshCell(VTKCellTypes.VTK_PIXEL, copy(point_ids))
      end
    end
  end

  # Convert array-of-points to two-dimensional array
  vtk_points = Matrix{Float64}(undef, ndim, length(points))
  for point_id in 1:length(points)
    vtk_points[1, point_id] = points[point_id][1]
    vtk_points[2, point_id] = points[point_id][2]
  end

  return vtk_points, vtk_cells
end

# Convert coordinates and level information to a list of points and VTK cells (3D version)
function calc_vtk_points_cells(::Val{3}, coordinates::AbstractMatrix{Float64},
                               levels::AbstractVector{Int},
                               center_level_0::AbstractVector{Float64},
                               length_level_0::Float64,
                               n_visnodes::Int=1)
  ndim = 3
  # Create point locator
  pl = PointLocator{3}(center_level_0, length_level_0, 1e-12)

  # Create arrays for points and cells
  n_elements = length(levels)
  points = Vector{SVector{3, Float64}}()
  vtk_cells = Vector{MeshCell}(undef, n_elements * n_visnodes^ndim)
  point_ids = Vector{Int}(undef, 2^ndim)

  # Reshape cell array for easy-peasy access
  reshaped = reshape(vtk_cells, n_visnodes, n_visnodes, n_visnodes, n_elements)

  # Create VTK cell for each Trixi element
  for element_id in 1:n_elements
    # Extract cell values
    cell_x = coordinates[1, element_id]
    cell_y = coordinates[2, element_id]
    cell_z = coordinates[3, element_id]
    cell_dx = length_level_0 / 2^levels[element_id]

    # Adapt to visualization nodes for easy-to-understand loops
    dx = cell_dx / n_visnodes
    x_lowerleft = cell_x - cell_dx/2 - dx/2
    y_lowerleft = cell_y - cell_dx/2 - dx/2
    z_lowerleft = cell_z - cell_dx/2 - dx/2

    # Create cell for each visualization node
    for k = 1:n_visnodes, j = 1:n_visnodes, i = 1:n_visnodes
      # Determine x and y
      x = x_lowerleft + i * dx
      y = y_lowerleft + j * dx
      z = z_lowerleft + k * dx

      # Get point id for each vertex
      point_ids[1] = insert!(pl, points, (x - dx/2, y - dx/2, z - dx/2))
      point_ids[2] = insert!(pl, points, (x + dx/2, y - dx/2, z - dx/2))
      point_ids[3] = insert!(pl, points, (x - dx/2, y + dx/2, z - dx/2))
      point_ids[4] = insert!(pl, points, (x + dx/2, y + dx/2, z - dx/2))
      point_ids[5] = insert!(pl, points, (x - dx/2, y - dx/2, z + dx/2))
      point_ids[6] = insert!(pl, points, (x + dx/2, y - dx/2, z + dx/2))
      point_ids[7] = insert!(pl, points, (x - dx/2, y + dx/2, z + dx/2))
      point_ids[8] = insert!(pl, points, (x + dx/2, y + dx/2, z + dx/2))

      # Add cell
      reshaped[i, j, k, element_id] = MeshCell(VTKCellTypes.VTK_VOXEL, copy(point_ids))
    end
  end

  # Convert array-of-points to two-dimensional array
  vtk_points = Matrix{Float64}(undef, ndim, length(points))
  for point_id in 1:length(points)
    vtk_points[1, point_id] = points[point_id][1]
    vtk_points[2, point_id] = points[point_id][2]
    vtk_points[3, point_id] = points[point_id][3]
  end

  return vtk_points, vtk_cells
end


text:
nmonomers(species::AbstractSpecies) = length(species.mids)

hasmonomer(species::AbstractSpecies, mid::Integer) = mid in species.mids

monomer_fraction(species::AbstractSpecies, mid::Integer) = 0.0

"""
    density!(species)

Update the internal density fields for an `AbstractSpecies`.
"""
density!(species::AbstractSpecies) = nothing

setup!(species::AbstractSpecies, sys::FieldSystem) = nothing

#==============================================================================#

"""
    chain_grid(N_block, Ns)

Calculate the total chain length and block discretization
for a generic polymer chain with variable number of blocks.

The number of contour segments returned for each block needs to be even
for use of Simpson's rule in density/stress calculations.
"""
function chain_grid(N_block::AbstractVector{<:Real}, Ns::Integer)
    N = sum(N_block)
    nblocks = length(N_block)
    ds = N / Ns

    begin_block = Vector{Int}(undef, nblocks+1)
    ds_block = Vector{Float64}(undef, nblocks)

    bgn = 1
    begin_block[1] = bgn

    # Loop over all blocks in the chain
    for iblk = 1:nblocks
        Ns_half = floor(Int, N_block[iblk] / ds / 2.0 + 0.5)

        # Calculate step size for each block
        if Ns_half == 0
            Ns_half = 1
            ds_block[iblk] = N_block[iblk] / 2.0
        else
            ds_block[iblk] = N_block[iblk] / Ns_half / 2.0
        end

        # Calculate location of block start
        bgn = bgn + Ns_half * 2
        begin_block[iblk+1] = bgn
    end

    return N, begin_block, ds_block
end