text:
struct Cost{T}
    #TODO: types for methods
    val
    grad
    hess
    sp::Vector{Vector{Int}}
    val_cache::Vector{T}
    grad_cache::Vector{T}
    hess_cache::Vector{T}
    idx::Vector{Int}
end

function Cost(f::Function, nx::Int, nu::Int, nw::Int, idx::Vector{Int}; eval_hess=false)
    #TODO: option to load/save methods
    @variables x[1:nx], u[1:nu], w[1:nw]
    val = f(x, u, w)
    grad = Symbolics.gradient(val, [x; u])
    val_func = eval(Symbolics.build_function([val], x, u, w)[2])
    grad_func = eval(Symbolics.build_function(grad, x, u, w)[2])
    if eval_hess
        hess = Symbolics.sparsehessian(val, [x; u])
        hess_func = eval(Symbolics.build_function(hess.nzval, x, u, w)[2])
        sparsity = [findnz(hess)[1:2]...]
        nh = length(hess.nzval)
    else
        hess_func = Expr(:null)
        sparsity = [Int[]]
        nh = 0
    end

    return Cost(val_func, grad_func, hess_func, sparsity,
        zeros(1), zeros(nx + nu), zeros(nh), idx)
end

Objective{T} = Vector{Cost{T}} where T

function eval_obj(obj::Objective, x, u, w)
    J = 0.0
    for cost in obj
        for t in cost.idx
            cost.val(cost.val_cache, x[t], u[t], w[t])
            J += cost.val_cache[1]
        end
    end
    return J
end

function eval_obj_grad!(grad, idx, obj::Objective, x, u, w)
    for cost in obj
        for t in cost.idx
            cost.grad(cost.grad_cache, x[t], u[t], w[t])
            @views grad[idx[t]] .+= cost.grad_cache
            # fill!(cost.grad_cache, 0.0) # TODO: confirm this is necessary
        end
    end
end

#TODO: test
function eval_obj_hess!(hess, idx, obj::Objective, x, u, w, σ)
    i = 1
    for cost in obj
        for t in cost.idx
            cost.hess(cost.hess_cache, x[t], u[t], w[t])
            cost.hess_cache .*= σ
            @views hess[idx[i]] .+= cost.hess_cache
            i += 1
            fill!(cost.hess_cache, 0.0) # TODO: confirm this is necessary
        end
    end
end

function sparsity_hessian(obj::Objective, nx::Vector{Int}, nu::Vector{Int})
    row = Int[]
    col = Int[]
    for cost in obj
        if !isempty(cost.sp[1])
            for t in cost.idx
                shift = (t > 1 ? (sum(nx[1:t-1]) + sum(nu[1:t-1])) : 0)
                push!(row, (cost.sp[1] .+ shift)...)
                push!(col, (cost.sp[2] .+ shift)...)
            end
        end
    end
    return collect(zip(row, col))
end

function hessian_indices(obj::Objective, key::Vector{Tuple{Int,Int}}, nx::Vector{Int}, nu::Vector{Int})
    idx = Vector{Int}[]
    for cost in obj
        if !isempty(cost.sp[1])
            for t in cost.idx
                row = Int[]
                col = Int[]
                shift = (t > 1 ? (sum(nx[1:t-1]) + sum(nu[1:t-1])) : 0)
                push!(row, (cost.sp[1] .+ shift)...)
                push!(col, (cost.sp[2] .+ shift)...)
                rc = collect(zip(row, col))
                push!(idx, [findfirst(x -> x == i, key) for i in rc])
            end
        end
    end
    return idx
end


text:
using Test
using LyapunovExponents
using LyapunovExponents: PhaseRelaxer,
    get_integrator, get_tangent_integrator, van_der_pol
using LyapunovExponents.TestTools: @test_isapprox_elemwise

@time @testset "Integrator time continuity" begin
    prob = van_der_pol().prob
    relaxer = PhaseRelaxer(get_integrator(prob.phase_prob),
                           prob.t_tran)
    dt = 1.0
    step!(relaxer.integrator, dt, true)
    @test relaxer.integrator.t > prob.phase_prob.tspan[1]

    integrator = get_tangent_integrator(prob, relaxer)
    @test integrator.t == relaxer.integrator.t

    step!(relaxer.integrator, dt, true)
    step!(integrator, dt, true)
    @test integrator.t == relaxer.integrator.t

    @test_isapprox_elemwise(integrator.u[:, 1],
                            relaxer.integrator.u,
                            rtol = 0.05)
end

ext: jl


text:
@testset "SDE Adaptive Distribution Tests" begin

using StochasticDiffEq, StatsBase, Distributions, HypothesisTests
using Random, DiffEqProblemLibrary

prob = prob_sde_linear
Random.seed!(200)
N = 100
M= 5
ps = Vector{Float64}(undef,M)
T = prob.tspan[2]

for j = 1:M
  Wends = Vector{Float64}(undef,N)
  for i = 1:N
    sol =solve(prob,SRI(),dt=1/2^(4),abstol=1e-2,reltol=0,adaptivealg=:RSwM1)
    Wends[i] = sol.W.W[end]
  end
  kssol = ApproximateOneSampleKSTest(Wends/sqrt(T), Normal())
  ps[j] = pvalue(kssol) #Should be not significant (most of the time)
end

@test sum(ps .> 0.05) > length(ps)/2 ### Make sure more passes than fails

for j = 1:M
  Wends = Vector{Float64}(undef,N)
  for i = 1:N
    sol =solve(prob,SRI(),dt=1/2^(4),abstol=1e-2,reltol=0,adaptivealg=:RSwM2)
    Wends[i] = sol.W.W[end]
  end
  kssol = ApproximateOneSampleKSTest(Wends/sqrt(T), Normal())
  ps[j] = pvalue(kssol) #Should be not significant (most of the time)
end

@test sum(ps .> 0.05) > length(ps)/2 ### Make sure more passes than fails

for j = 1:M
  Wends = Vector{Float64}(undef,N)
  for i = 1:N
    sol =solve(prob,SRI(),dt=1/2^(4),abstol=1e-2,reltol=0,adaptivealg=:RSwM3)
    Wends[i] = sol.W.W[end]
  end
  kssol = ApproximateOneSampleKSTest(Wends/sqrt(T), Normal())
  ps[j] = pvalue(kssol) #Should be not significant (most of the time)
end

@test sum(ps .> 0.05) > length(ps)/2 ### Make sure more passes than fails

end


text:
module libm

import Base.Math.@horner

macro horner_fma(x, p...)
    ex = esc(p[end])
    for i = length(p)-1:-1:1
        ex = :(fma(t,$ex,$(esc(p[i]))))
    end
    Expr(:block, :(t = $(esc(x))), ex)
end

macro horner_muladd(x, p...)
    ex = esc(p[end])
    for i = length(p)-1:-1:1
        ex = :(muladd(t,$ex,$(esc(p[i]))))
    end
    Expr(:block, :(t = $(esc(x))), ex)
end

half(x::Float64) = 0.5
half(x::Float32) = 0.5f0

immutable Double{T}
    hi::Float64
    lo::Float64
end

typealias Double64 Double{Float64}

trunc32(x::Float64) = reinterpret(Float64,reinterpret(Uint64,x) & 0xffff_ffff_0000_0000)

# include("frexp.jl")
# include("exponent.jl")
# include("scalbn.jl")

# include("trunc.jl")
# include("round.jl")

include("log.jl")
# include("simd.jl")
# include("log_simd.jl")
# include("log_fma.jl")
include("log_tang.jl")
include("log2_tang.jl")
include("log_yeppp.jl")
# include("log_table2.jl")
include("expm1.jl")

include("pow.jl")
include("pow1p.jl")

end # module


text:
struct LinearInterpolant{T,N}
    y::Array{T,N}
    Δx::Float64
end

function (f::LinearInterpolant)(a::Float64)
    Δx = f.Δx
    s = a/(Δx);
    n₁ = floor(Int64,s);
    n₂ = ceil(Int64,s);
    Δn = (s - n₁)
    (1-Δn)*f.y[n₁+1] + Δn*f.y[n₂+1]
end

function qnms(;l=0,m=0,n=0,s=0,amax=0.99, ϵ = 0.01)
    as, ωs, Almss, Cllss, P, ϵ = GetModes(l,m,n,s;amax = amax,ϵ = ϵ)
    ωfunc = LinearInterpolant(ωs,ϵ)
    Afunc = LinearInterpolant(Almss,ϵ)
    Cllss = LinearInterpolant(Cllss,ϵ)
    ωfunc, Afunc, Cllss
end

function qnm(;l=0,m=0,n=0,s=0,a=0.0, ϵ = 0.01)
     _, ωs, Almss, Cllss, P, _ = GetModes(l,m,n,s;amax = a,ϵ = ϵ)
    #ωs[end], Almss[end], Cllss[:,end], P
    Alms, Cll = ComputeAₗₘ(P.s,P.m,P.a*P.ω, P.Alm,P.lmax)
    ωs[end], Alms, Cll, P
end


text:
################################################################################
#
#  Abstract Jordan decomposition
#
################################################################################

# This holds invariants of a local Jordan decomposition
#
# L = L_1 \perp ... \perp L_r
#
# In the non-dyadic case we store
# - ranks
# - scales
# - determinant (classes)
# of the L_i
#
# In the dyadic case we store
# - norm generators of L_i
# - (valuation of ) weights
# - determinant (classes)
# - Witt invariants

mutable struct JorDec{S, T, U}
  K::S
  p::T
  is_dyadic::Bool
  ranks::Vector{Int}
  scales::Vector{Int}

  # dyadic things
  normgens::Vector{U}
  weights::Vector{Int}
  dets::Vector{U}
  witt::Vector{Int}

  JorDec{S, T, U}() where {S, T, U} = new{S, T, U}()
end

function JorDec(p, scales::Vector{Int}, ranks::Vector{Int}, dets::Vector{nf_elem})
  K = nf(order(p))
  _weight = Vector{Int}()
  _normgen = Vector{elem_type(K)}()
  _f = Int[]
  witt = Int[]
  z = JorDec{typeof(K), typeof(p), elem_type(K)}()
  z.is_dyadic = is_dyadic(p)
  z.K = K
  z.p = p
  z.ranks = ranks
  z.scales = scales
  z.normgens = _normgen
  z.weights = _weight
  z.witt = witt
  z.dets = dets
  return z
end

function Base.show(io::IO, J::JorDec)
  p = J.p
  if !is_dyadic(p)
    for i in 1:length(J)
      print(io, "(", J.scales[i], ", ", J.ranks[i], ", ", J.dets[i], ")")
    end
  else
    for i in 1:length(J)
      print(io, "(", J.scales[i], ", ", J.ranks[i], ", ", J.normgens[i], ", ", J.weights[i], ", ", J.dets[i], ", ", J.witt[i], ")")
    end
  end
end

function Base.show(io::IO, ::MIME"text/plain", J::JorDec)
  p = J.p
  if !get(io, :compact, false)
    print(IOContext(io, :compact => true), "Abstract Jordan decomposition at ", p)
    if !is_dyadic(p)
      print(io, "\n(scale, rank, determinant class)")
    else
      print(io, "\n(scale, rank, norm generator, weight, det, Witt)")
    end
  end
  print(io, "\n")
  if !is_dyadic(p)
    for i in 1:length(J)
      print(io, "(", J.ranks[i], ", ", J.scales[i], ", ", J.dets[i], ")")
    end
  else
    for i in 1:length(J)
      print(io, "(", J.scales[i], ", ", J.ranks[i], ", ", J.normgens[i], ", ", J.weights[i], ", ", J.dets[i], ", ", J.witt[i], ")")
    end
  end
end

length(J::JorDec) = length(J.ranks)

function JorDec(p, sc::Vector{Int}, rks::Vector{Int}, normgens::Vector{nf_elem}, weights::Vector{Int}, dets::Vector{nf_elem}, witts::Vector{Int})
  K = nf(order(p))
  z = JorDec{typeof(K), typeof(p), elem_type(K)}()
  z.is_dyadic = is_dyadic(p)
  z.K = K
  z.p = p
  z.ranks = rks
  z.scales = sc
  z.normgens = normgens
  z.weights = weights
  z.witt = witts
  z.dets = dets
  return z
end

# Construction of abstract Jordan decomposition from a concrete Jordan decomposition
function JorDec(J, G, E, p)
  O = order(p)
  K = nf(O)
  t = length(G)
  _, _h = ResidueField(O, p)
  h = extend(_h, K)
  if !is_dyadic(p)
    dets = elem_type(K)[det(G[i]) for i in 1:t]
    _weight = Vector{Int}()
    _normgen = Vector{elem_type(K)}()
    _f = Int[]
    witt = Int[]
    #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], is_square(Gs[i])[1] ? 1 : -1) for i in 1:length(J)]
    #return LocalGenusSymbol{QuadLat}(p, x, unif, false, base_field(L), nothing, nothing)
  else
    sL = Int[ minimum(Union{PosInf, Int}[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
    _sym = Vector{Tuple{Int, Int, Int}}(undef, t)
    e = ramification_index(p)
    aL = elem_type(K)[]
    uL = Int[]
    wL = Int[]
    witt = Int[]
    dets = elem_type(K)[]
    unif = elem_in_nf(uniformizer(p))
    for i in 1:t
      _sym[i] = (nrows(J[i]), sL[i], 0)
      push!(dets, det(G[i]))
      push!(witt, witt_invariant(quadratic_space(K, G[i]), p))
      #GG = diagonal_matrix(eltype(G)[ j < i ? unif^(2*(sL[i] - sL[j])) * G[j] : G[j] for j in 1:t])
      D = diagonal(G[i])
      m, pos = findmin(Union{PosInf, Int}[iszero(d) ? inf : 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(Union{PosInf, Int}[uL[i] + quadratic_defect(d//aL[i], p) for d in D])))
    end
    _normgen = aL
    _weight = wL
  end

  z = JorDec{typeof(K), typeof(p), elem_type(K)}()
  z.is_dyadic = is_dyadic(p)
  z.K = K
  z.p = p
  z.ranks = Int[nrows(J[i]) for i in 1:length(J)]
  z.scales = E
  z.normgens = _normgen
  z.weights = _weight
  z.witt = witt
  z.dets = dets
  return z
end

function JorDec(L::QuadLat, p)
  O = base_ring(L)
  K = nf(O)
  J, G, E = jordan_decomposition(L, p)
  return JorDec(J, G, E, p)
end

@doc Markdown.doc"""
    gram_matrix(J::JorDec, i::Int) -> MatElem

Given an abstract Jordan decomposition $J$ and $i$, return the Gram matrix of
of the $i$-th block of $J$.
"""
function gram_matrix(J::JorDec, i::Int)
  @req 1 <= i <= length(J) "Index ($i) must be in range 1:$(length(J))"
  p = J.p
  if J.is_dyadic
    r, s, w = J.ranks[i], J.scales[i], J.weights[i]
    d, a, wi = J.dets[i], J.normgens[i], J.witt[i]
    # Not that these are the invariants of G(L_i)
    # We now write G(L_i) = pi^s * G(L_i') with L_i' = L_i^(pi^-s) and L_i' is
    # unimodular. We first have to translate the invariants of L_i to L_i'
    # and then use the classification of unimodular lattices.
    pi = elem_in_nf(uniformizer(J.p))
    _w = w - s
    _d = pi^(-s * r) * d
    @assert valuation(_d, p) == 0
    _a = a//pi^s
    # Adjust the Witt invariant
    ha = _witt_hasse(wi, r, d, p)
    # this is the hasse invariant
    # this is the hasse invariant of the scaled lattice
    hanew = ha * hilbert_symbol((-1)^(div(r*(r-1), 2)) * d^(r - 1), inv(pi)^s, p)
    winew = _witt_hasse(hanew, r, _d, p)
    # this is the new witt invariant
    z = _quadratic_unimodular_lattice_dyadic(p, r, _w, _d, _a, winew)
    @assert begin L = quadratic_lattice(nf(order(p)), gram = z); witt_invariant(L, p) == winew end

    @assert valuation(det(z), p) == 0
    zz = pi^s * z
    if wi < 2
      @assert begin L = quadratic_lattice(nf(order(p)), gram = zz); witt_invariant(L, p) == wi end
    end
    return zz
  else
    r = J.ranks[i]
    s = J.scales[i]
    z = identity_matrix(J.K, r)
    pi = elem_in_nf(uniformizer(p))
    q = J.dets[i]//pi^(r * s)
    @assert valuation(q, p) == 0
    if is_local_square(pi^(r * s) * J.dets[i], p)
      z[r, r] = one(parent(pi))
    else
      z[r, r] = _non_square(J.K, p)
    end
    return pi^s * z
  end
end

@doc Markdown.doc"""
    gram_matrix(J::JorDec) -> MatElem

Given an abstract Jordan decomposition, return the Gram matrix of a lattice
admitting this Jordan decompositon.
"""
function gram_matrix(J::JorDec)
  K = J.K
  if length(J) > 0
    D = diagonal_matrix(dense_matrix_type(K)[gram_matrix(J, i) for i in 1:length(J)])
  else
    D = zero_matrix(K, 0, 0)
  end
  return D
end

@doc Markdown.doc"""
    lattice(J::JorDec) -> MatElem

Given an abstract Jordan decomposition, return a lattice admitting this Jordan
decompositon.
"""
function lattice(J::JorDec)
  return quadratic_lattice(J.K, gram = gram_matrix(J))
end

@doc Markdown.doc"""
    genus(J::JorDec) -> LocalGenusQuad

Given an abstract Jordan decomposition, return the local genus to which the
corresponding matrix belongs.
"""
function genus(J::JorDec)
  r = length(J.ranks)
  sca = J.scales
  p = J.p
  pi = elem_in_nf(uniformizer(p))
  if !is_dyadic(p)
    detclass = Int[]
    for i in 1:length(J.ranks)
      push!(detclass, is_local_square(J.dets[i]//pi^(J.ranks[i] * J.scales[i]), p) ? 1 : -1)
    end
    z = genus(QuadLat, p, pi, J.ranks, J.scales, detclass)
    z.dets = J.dets
    z.jordec = J
    return z
  end
  K = nf(order(p))
  normgens = J.normgens
  new_normgens = Vector{eltype(normgens)}(undef, r)
  new_weights = Vector{Int}(undef, r)
  e = ramification_index(p)
  for i in 1:r
    # find the minimal norm
    # L^(s_i) = \sum p^(max(0, s_i - s_j)) L_j
    new_normgen = zero(K)
    new_norm = inf # valuation(normgens[i], p)
    for j in 1:r
      poss_normgen = pi^max(0, 2*(sca[i] - sca[j])) * normgens[j]
      poss_n = max(0, 2*(sca[i] - sca[j])) + valuation(normgens[j], p)
      if poss_n < new_norm
        new_norm = poss_n
        new_normgen = poss_normgen
      end
    end

    # Now determine w(L^(s_i)) using this norm generator
    # We first determine the diagonal of a gram matrix of L^(s_i)
    D = elem_type(K)[]
    for j in 1:r
      # TODO: We could determine just the diagonal, that would be easier.
      # Also, there is some redundancy here, the z is acutally independent of i
      z = gram_matrix(J, j)
      for l in 1:nrows(z)
        push!(D, pi^(max(0, 2 * (sca[i] - sca[j]))) * z[l, l])
      end
    end

    # Lemma 3.3.9

    #@show new_normgen
    #@show new_norm
    #@show D

    # Now the scale of L^(s_i) is si
    weight = minimum(Union{PosInf, Int}[new_norm + quadratic_defect(divexact(d, new_normgen), p) for d in D])
    #@show weight
    weight = min(weight, e + sca[i])
    #@show weight
    new_weights[i] = weight
    new_normgens[i] = new_normgen
  end
  g =  genus(QuadLat, p, J.ranks, J.scales, new_weights, J.dets, new_normgens, J.witt)
  g.jordec = J
  return g
end

function orthogonal_sum(J1::JorDec{S, T, U}, J2::JorDec{S, T, U}) where {S, T, U}
  @req J1.p === J2.p "Jordan decompositions must be over same prime"
  if !(J1.is_dyadic)
    i1 = 1
    i2 = 1
    _sca = Int[]
    _rk = Int[]
    _dets = U[]
    while i1 <= length(J1) && i2 <= length(J2)
      if J1.scales[i1] < J2.scales[i2]
        push!(_sca, J1.scales[i1])
        push!(_rk, J1.ranks[i1])
        push!(_dets, J1.dets[i1])
        i1 += 1
      elseif J2.scales[i2] < J1.scales[i1]
        push!(_sca, J2.scales[i2])
        push!(_rk, J2.ranks[i2])
        push!(_dets, J2.dets[i2])
        i2 += 1
      else
        push!(_sca, J1.scales[i1])
        push!(_rk, J1.ranks[i1] + J2.ranks[i2])
        push!(_dets, J1.dets[i1] * J2.dets[i2])
        i1 += 1
        i2 += 1
      end
    end

    if i1 <= length(J1)
      append!(_sca, J1.scales[i1:length(J1)])
      append!(_rk, J1.ranks[i1:length(J1)])
      append!(_dets, J1.dets[i1:length(J1)])
    end

    if i2 <= length(J2)
      append!(_sca, J2.scales[i2:length(J2)])
      append!(_rk, J2.ranks[i2:length(J2)])
      append!(_dets, J2.dets[i2:length(J2)])
    end

    return JorDec(J1.p, _sca, _rk, _dets)
  else
    # Lazy
    return JorDec(orthogonal_sum(lattice(J1), lattice(J2))[1], J1.p)
  end
end

################################################################################
#
#  Construct unimodular lattices with given invariants
#
################################################################################

# p = prime ideal
# r = rank
# w = weight (valuation)
# d = determinant (class)
# alpha = norm generator
# wi = Witt invariants
#
# This is [Kirschmer2016, Proposition 3.3.11]
function _quadratic_unimodular_lattice_dyadic(p, r, w, d, alpha, wi)
  m = r
  K = nf(order(p))
  e = ramification_index(p)
  pi = elem_in_nf(uniformizer(p))
  if m == 1
    return matrix(K, 1, 1, elem_type(K)[d])
  end

  D = kummer_generator_of_local_unramified_quadratic_extension(p)
  rho = divexact(1 - D, 4)
  @assert valuation(rho, p) == 0
  @assert quadratic_defect(D, p) == valuation(4, p)

  if isodd(m)
    r = div(m - 1, 2)
    @assert (w == e) || (w < e && isodd(w))
    mats = dense_matrix_type(K)[]
    if wi == 1
      push!(mats, _Amatrix(K, pi^w, 0))
      push!(mats, matrix(K, 1, 1, [(-1)^r * d]))
      for i in 1:(r - 1)
        push!(mats, _Amatrix(K, 0, 0))
      end
    else
      @assert wi == -1
      push!(mats, _Amatrix(K, pi^w, 4 * rho * pi^(-w)))
      push!(mats, matrix(K, 1, 1, [inv(D) * (-1)^r * d]))
      for i in 1:(r - 1)
        push!(mats, _Amatrix(K, 0, 0))
      end
    end
    return diagonal_matrix(mats)
  else
    r = div(m, 2)
    @assert valuation(d, p) == 0
    mmod4 = mod(m, 4)
    if mmod4 == 0 || mmod4 == 1
      gamma = _find_special_class(d, p) - 1
      @assert is_local_square(d//(1 + gamma), p)
      @assert valuation(d, p) == valuation(1 + gamma, p)
    else
      gamma = _find_special_class(-d, p) - 1
      @assert is_local_square(-d//(1 + gamma), p)
      @assert valuation(-d, p) == valuation(1 + gamma, p)
    end
    @assert quadratic_defect(1 + gamma, p) == (iszero(gamma) ? inf : valuation(gamma, p))
    if iseven(valuation(alpha, p) + w)
      @assert w == e
      mats = dense_matrix_type(K)[]
      @assert iszero(gamma) || valuation(gamma, p) >= w + valuation(alpha, p)
      push!(mats, _Amatrix(K, alpha, -gamma * inv(alpha)))
      @assert valuation(det(mats[end]), p) == 0
      for i in 1:(r - 1)
        push!(mats, _Amatrix(K, 0, 0))
      end
      return diagonal_matrix(mats)
    elseif r == 1
      @assert isodd(valuation(alpha, p) + w)
      @assert w == e || (e > w && w == valuation(gamma, p) - valuation(alpha, p))
      return _Amatrix(K, alpha, -gamma * inv(alpha))
    elseif r >= 2 && isodd(valuation(alpha, p) + w)
      mats = dense_matrix_type(K)[]
      # Kirschmer Remark 3.3.12
      if wi == hilbert_symbol(alpha, 1 + gamma, p)
        push!(mats, _Amatrix(K, alpha, -gamma * inv(alpha)))
        push!(mats, _Amatrix(K, pi^w, 0))
        for i in 1:(r - 2)
          push!(mats, _Amatrix(K, 0, 0))
        end
      else
        push!(mats, _Amatrix(K, alpha, -(gamma - 4 * rho) * inv(alpha)))
        push!(mats, _Amatrix(K, pi^w, 4 * rho * pi^(-w)))
        for i in 1:(r - 2)
          push!(mats, _Amatrix(K, 0, 0))
        end
      end
      z = diagonal_matrix(mats)
      @assert valuation(det(z), p) == 0
      return z
    else
      throw(error("This should not happen"))
    end
  end
end

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

# This holds invariants of a local Genus symbol
#
# L = L_1 \perp ... \perp L_r
#
# In the non-dyadic case we store
# - ranks
# - scales
# - determinant (classes)
# of the L_i = L^(s_i)
#
# In the dyadic case we store
# - norm generators of L^(s_i)
# - (valuation of ) weights of L^(s_i)
# - determinant (classes) of L^(s_i)
# - Witt invariants of L_i

mutable struct LocalGenusQuad{S, T, U}
  K::S
  p::T
  is_dyadic::Bool
  witt_inv::Int
  hass_inv::Int
  det::U
  rank::Int

  uniformizer::U

  ranks::Vector{Int}
  scales::Vector{Int}
  detclasses::Vector{Int}

  # dyadic things
  normgens::Vector{U}
  weights::Vector{Int}
  f::Vector{Int}
  dets::Vector{U}
  witt::Vector{Int}
  norms::Vector{Int}

  # Sometimes we know a jordan decomposition
  jordec::JorDec{S, T, U}

  function LocalGenusQuad{S, T, U}() where {S, T, U}
    z = new{S, T, U}()
    z.rank = 0
    z.witt_inv = 0
    z.hass_inv = 0
    return z
  end
end

local_genus_quad_type(K) = LocalGenusQuad{typeof(K), ideal_type(order_type(K)), elem_type(K)}

function in(L::QuadLat, G::LocalGenusQuad)
  return genus(L, prime(G)) == G
end

function local_quadratic_genus_type(K)
  return LocalGenusQuad{typeof(K),
                        ideal_type(order_type(K)),
                        elem_type(K)}
end

# Access
prime(G::LocalGenusQuad) = G.p

length(G::LocalGenusQuad) = length(G.ranks)

scales(G::LocalGenusQuad) = G.scales

ranks(G::LocalGenusQuad) = G.ranks

dets(G::LocalGenusQuad) = G.detclasses

weights(G::LocalGenusQuad) = G.weights

scale(G::LocalGenusQuad, i::Int) = scales(G)[i]

rank(G::LocalGenusQuad, i::Int) = ranks(G)[i]

function witt_invariant(G::LocalGenusQuad)
  if G.witt_inv != 0
    return G.witt_inv
  end

  p = prime(G)

  @assert is_dyadic(p)

  w, d, n = G.witt[1], G.dets[1], G.ranks[1]

  for i in 2:length(G)
    d, w, n = _witt_of_orthgonal_sum(d, w, n, G.dets[i], G.witt[i], G.ranks[i], p)
  end

  G.witt_inv = w

  return w
end

function rank(G::LocalGenusQuad)
  if G.rank != 0
    return G.rank
  end

  rk = sum(G.ranks)

  G.rank = rk
  return rk
end

function det(G::LocalGenusQuad)
  if isdefined(G, :det)
    return G.det
  end

  if isdefined(G, :dets)
    d = prod(G.dets, init = one(nf(order(G.p))))
    G.det = d
  else
    pi = uniformizer(G)
    d = reduce(*, nf_elem[pi^(G.ranks[i] * G.scales[i]) * G.detclasses[i] for i in 1:length(G)], init = one(nf(order(G.p))))
    G.det = d
  end
  return d
end

function det(G::LocalGenusQuad, i::Int)
  #if is_dyadic(G)
    return G.dets[i]
  #else
  #  pi = uniformizer(G)
  #  return pi^(G.ranks[i] * G.scales[i]) * G.detclasses[i]
  #end
end

function hasse_invariant(G::LocalGenusQuad)
  if rank(G) == 0
    return 1
  end
  if is_dyadic(G)
    w = witt_invariant(G)
    return _witt_hasse(w, rank(G), det(G), prime(G))
  else
    p = prime(G)
    l = length(G)
    pi = uniformizer(G)
    h = hilbert_symbol((-1)^divexact(G.ranks[1] * (G.ranks[1] - 1), 2) * det(G, 1)^(G.ranks[1] - 1), pi^G.scales[1], p)
    d = det(G, 1)
    for i in 2:l
      h = h * hilbert_symbol(d, det(G, i), p)
      h = h * hilbert_symbol((-1)^divexact(G.ranks[i] * (G.ranks[i] - 1), 2) * det(G, i)^(G.ranks[i] - 1), pi^G.scales[i], p)
      d = d * det(G, i)
    end
    @assert hasse_invariant(representative(G), p) == h
    return h
  end
end

function uniformizer(G::LocalGenusQuad)
  @req !is_dyadic(G) "Genus symbol must not be dyadic"
  return G.uniformizer
end

is_dyadic(G::LocalGenusQuad) = G.is_dyadic

function norm_generators(G::LocalGenusQuad)
  @req is_dyadic(G) "Genus symbol must be dyadic"
  return G.normgens
end

function norms(G::LocalGenusQuad)
  @req is_dyadic(G) "Genus symbol must be dyadic"
  if !isdefined(G, :norms)
    p = prime(G)
    G.norms = Int[valuation(a, p) for a in norm_generators(G)]
    return G.norms
  else
    return G.norms
  end
end

function jordan_decomposition(g::LocalGenusQuad)
  if isdefined(g,:jordec)
    j = g.jordec
  elseif rank(g) == 0
    j = JorDec(g.p, Int[], Int[], nf_elem[])
    g.jordec = j
    return j
  else
    # We don't know how to determine the Jordan decomposition
    # directly from the fundamental invariants. We just list all
    # possible jordan decompositions.
    # TODO: We could do it in the good case
    possible_jdec = _local_jordan_decompositions(number_field(order(g.p)),
                                                 g.p,
                                                 collect(zip(g.scales, g.ranks)))
    for j in possible_jdec
      if genus(j) == g
        g.jordec = j
        return j
      end
    end
  end
  return j
end

function Base.show(io::IO, G::LocalGenusQuad{S, T, U}) where {S, T, U}
  if !is_dyadic(G)
    for i in 1:length(G)
      print(io, "(", G.scales[i], ", ", G.ranks[i], ", ", G.detclasses[i], ")")
    end
  else
    for i in 1:length(G)
      print(io, "(", G.scales[i], ", ", G.ranks[i], ", ", G.normgens[i], ", ", G.weights[i], ", ", G.dets[i], ", ", G.witt[i], ")")
    end
  end
end

function Base.show(io::IO, ::MIME"text/plain", G::LocalGenusQuad{S, T, U}) where {S, T, U}
  p = prime(G)
  if !get(io, :compact, false)
    print(io, "Local quadratic genus for prime ")
    print(IOContext(io, :compact => true), p)

    if length(G) == 0
      print(io, " of rank zero ")
      return
    end

    if !is_dyadic(p)
      print(io, " with respect to uniformizer ", uniformizer(G))
      print(io, "\n(scale, rank, determinant class)\n")
    else
      print(io, "\n(scale, rank, norm generator, weight, det, Witt)\n")
    end
  end
  if !is_dyadic(G)
    for i in 1:length(G)
      print(io, "(", G.scales[i], ", ", G.ranks[i], ", ", G.detclasses[i], ")")
    end
  else
    for i in 1:length(G)
      print(io, "(", G.scales[i], ", ", G.ranks[i], ", ", G.normgens[i], ", ", G.weights[i], ", ", G.dets[i], ", ", G.witt[i], ")")
    end
  end
end

# Creation of non-dyadic genus symbol

function genus(::Type{QuadLat}, p, pi::nf_elem, ranks::Vector{Int},
                                                scales::Vector{Int},
                                                normclass::Vector{Int})
  @req !is_dyadic(p) "Prime ideal must not be dyadic"
  K = nf(order(p))
  z = LocalGenusQuad{typeof(K), typeof(p), elem_type(K)}()
  z.p = p
  z.uniformizer = pi
  z.is_dyadic = false
  z.ranks = ranks
  z.scales = scales
  z.detclasses = normclass
  return z
end

# Creation of dyadic genus symbols, with the f already computed

function genus(::Type{QuadLat}, p, ranks::Vector{Int}, scales::Vector{Int},
                                   weights::Vector{Int}, normgens::Vector{T},
                                   dets::Vector{T}, witt::Vector{Int},
                                   f::Vector{Int}) where {T}
  @req is_dyadic(p) "Prime ideal must be dyadic"
  K = nf(order(p))
  z = LocalGenusQuad{typeof(K), typeof(p), elem_type(K)}()
  z.is_dyadic = true
  z.p = p
  z.ranks = ranks
  z.scales = scales
  z.weights = weights
  z.normgens = normgens
  z.dets = dets
  z.witt = witt
  z.f = f
  return z
end

# Creation of dyadic genus symbols, without the f already computed
function genus(::Type{QuadLat}, p, ranks::Vector{Int}, scales::Vector{Int},
                                   weights::Vector{Int}, dets::Vector{T},
                                   normgens::Vector{T},
                                   witt::Vector{Int}) where {T}
  @req is_dyadic(p) "Prime ideal must be dyadic"
  K = nf(order(p))
  z = LocalGenusQuad{typeof(K), typeof(p), elem_type(K)}()
  z.is_dyadic = true
  z.p = p
  z.ranks = ranks
  z.scales = scales
  z.weights = weights
  z.normgens = normgens
  z.dets = dets
  z.witt = witt
  t = length(ranks)
  # I should do this only on demand ...
  uL = Int[valuation(normgens[i], p) for i in 1:t]
  wL = weights
  sL = scales
  aL = normgens
  _f = Vector{Int}()
  e = ramification_index(p) # absolute
  for k in 1:(t - 1)
    exp = uL[k] + uL[k + 1]
    push!(_f, (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
  z.f = _f
  return z
end

# creation of rank zero genus symbol
function genus(::Type{QuadLat}, p)
  if is_dyadic(p)
    T = elem_type(nf(order(p)))
    return genus(QuadLat, p, Int[], Int[], Int[], T[], T[], Int[])
  else
    return genus(QuadLat, p, elem_in_nf(uniformizer(p)), Int[], Int[], Int[])
  end
end

################################################################################
#
#  Equality of genus symbols
#
################################################################################

function Base.:(==)(G1::LocalGenusQuad, G2::LocalGenusQuad)
  if G1 === G2
    return true
  end

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

  p = prime(G1)

  # Test if the rational spaces are isometric
  if is_dyadic(G1)
    # Could be sped up for low rank
    w1 = witt_invariant(G1)
    d1 = prod(nf_elem[G1.dets[i] for i in 1:length(G1)])
    n1 = rank(G1)
    #s1 = _witt_hasse(w1, n1, d1, p)

    w2 = witt_invariant(G2)
    d2 = prod(nf_elem[G2.dets[i] for i in 1:length(G2)])
    n2 = rank(G2)
    #s2 = _witt_hasse(w2, n2, d2, p)

    if n1 != n2
      return false
    end
    if w1 != w2
      return false
    end
    if !is_local_square(d1 * d2, p)
      return false
    end
  end

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

  if ranks(G1) != ranks(G2) || scales(G1) != scales(G2)
    return false
  end

  if !is_dyadic(G1.p)
    if G1.uniformizer == G2.uniformizer
      return G1.detclasses == G2.detclasses
    else
      q = divexact(G2.uniformizer, G1.uniformizer)
      fl = is_local_square(q, p)
      if fl
        return G1.detclasses == G2.detclasses
      else
        G2adj = Vector{Int}(undef, length(G2))
        for i in 1:length(G2)
          if isodd(G1.ranks[i] * G1.scales[i])
            G2adj[i] = -G2.detclasses[i]
          else
            G2adj[i] = G2.detclasses[i]
          end
        end
        return G1.detclasses == G2adj
      end
    end
  end

  if ranks(G1) != ranks(G2) || scales(G1) != scales(G2) || weights(G1) != weights(G2)
    return false
  end

  r = length(G1.ranks)

  uL1 = norms(G1)
  uL2 = norms(G2)

  if uL1 != uL2
    return false
  end

  bL = nf_elem[divexact(G1.normgens[i], G2.normgens[i]) for i in 1:r]
  qL = Union{PosInf, Int}[quadratic_defect(bL[i], p) for i in 1:r]

  for k in 1:r
    if qL[k] < G1.weights[k] - uL2[k]
      return false
    end
  end

  e = ramification_index(p)

  aL1 = G1.normgens
  aL2 = G2.normgens

  d1 = nf_elem[prod(nf_elem[G1.dets[i] for i in 1:j]) for j in 1:r]
  d2 = nf_elem[prod(nf_elem[G2.dets[i] for i in 1:j]) for j in 1:r]

  for i in 1:(r - 1)
    detquot = divexact(d1[i], d2[i])
    if valuation(detquot, p) != 0
      return false
    end

    if quadratic_defect(detquot, p) < G1.f[i]
      return false
    end

    if G1.f[i] > 2 * e + uL1[i + 1] - G2.weights[i + 1]
      # We have to test if L1^(i) embeds into L2^(i) \perp aL2[i + 1]
      # So let's compute the Hasse invariant of those fields.
      haG1 = _witt_hasse(G1.witt[1], G1.ranks[1], G1.dets[1], p)
      _d1 = G1.dets[1]
      haG2 = _witt_hasse(G2.witt[1], G2.ranks[1], G2.dets[1], p)
      _d2 = G2.dets[1]
      ra = G1.ranks[1]

      for j in 2:i
        haG1 = haG1 * _witt_hasse(G1.witt[j], G1.ranks[j], G1.dets[j], p) * hilbert_symbol(_d1, G1.dets[j], p)
        _d1 = _d1 * G1.dets[j]
        haG2 = haG2 * _witt_hasse(G2.witt[j], G2.ranks[j], G2.dets[j], p) * hilbert_symbol(_d2, G2.dets[j], p)
        _d2 = _d2 * G2.dets[j]
        ra = ra + G1.ranks[j]
      end

      detsecondspace = _d2 * aL2[i + 1]
      hassesecondspace = haG2 * hilbert_symbol(_d2, aL2[i + 1], p)

      if !_can_locally_embed(r, _d1, haG1, r + 1, detsecondspace, hassesecondspace, p)
        return false
      end
    end

    if G1.f[i] > 2 * e + uL1[i] - G2.weights[i]
      # We have to test if L1^(i) embeds into L2^(i) \perp aL2[i]
      # So let's compute the Hasse invariant of those fields.
      haG1 = _witt_hasse(G1.witt[1], G1.ranks[1], G1.dets[1], p)
      _d1 = G1.dets[1]
      haG2 = _witt_hasse(G2.witt[1], G2.ranks[1], G2.dets[1], p)
      _d2 = G2.dets[1]
      ra = G1.ranks[1]

      for j in 2:i
        haG1 = haG1 * _witt_hasse(G1.witt[j], G1.ranks[j], G1.dets[j], p) * hilbert_symbol(_d1, G1.dets[j], p)
        _d1 = _d1 * G1.dets[j]
        haG2 = haG2 * _witt_hasse(G2.witt[j], G2.ranks[j], G2.dets[j], p) * hilbert_symbol(_d2, G2.dets[j], p)
        _d2 = _d2 * G2.dets[j]
        ra = ra + G1.ranks[j]
      end

      detsecondspace = _d2 * aL2[i]
      hassesecondspace = haG2 * hilbert_symbol(_d2, aL2[i], p)

      if !_can_locally_embed(r, _d1, haG1, r + 1, detsecondspace, hassesecondspace, p)
        return false
      end
    end
  end

  return true
end

function is_locally_isometric(L::QuadLat, M::QuadLat, p::NfOrdIdl)
  fl = genus(L, p) == genus(M, p)
  #@assert fl == is_locally_isometric_kirschmer(L, M, p)
  return fl
end

function genus(L::QuadLat, p)
  if !has_attribute(L, :local_genera)
    local_genera = Dict{ideal_type(base_ring(L)), local_genus_quad_type(base_field(L))}()
    set_attribute!(L, :local_genera, local_genera)
  else
    local_genera = get_attribute(L, :local_genera)
  end

  if haskey(local_genera, p)
    return local_genera[p]::local_genus_quad_type(base_field(L))
  end

  pi, _sym, _weight, _normgen, _f, dets, witt, uL, J, G, E  = _genus_symbol(L, p)
  ranks = Int[d[1] for d in _sym]
  scales = Int[d[2] for d in _sym]

  if is_dyadic(p)
    g = genus(QuadLat, p, ranks, scales, _weight, _normgen, dets, witt,  _f)
  else
    normclass = Int[d[3] for d in _sym]
    g = genus(QuadLat, p, pi, ranks, scales, normclass)
  end
  g.norms = uL
  g.jordec = JorDec(J, G, E, p)
  local_genera[p] = g
  return g::local_genus_quad_type(base_field(L))
end

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

  J, G, E = jordan_decomposition(L, p)

  t = length(G)

  local pi::nf_elem

  if minimum(p) != 2
    _, _h = ResidueField(O, p)
    h = extend(_h, nf(O))
    _sym = Vector{Tuple{Int, Int, Int}}(undef, length(J))
    pi = elem_in_nf(uniformizer(p))
    for i in 1:t
      r = nrows(J[i])
      _sym[i] = (r, E[i], is_local_square(det(G[i])//pi^(r * E[i]), p) ? 1 : - 1)
    end
    _weight = Vector{Int}()
    _normgen = Vector{elem_type(K)}()
    _f = Int[]
    dets = elem_type(K)[det(G[i]) for i in 1:t]
    witt = Int[]
    uL = Int[]
    #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], is_square(Gs[i])[1] ? 1 : -1) for i in 1:length(J)]
    #return LocalGenusSymbol{QuadLat}(p, x, unif, false, base_field(L), nothing, nothing)
  else
    sL = Int[ minimum(Union{PosInf, Int}[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
    _sym = Vector{Tuple{Int, Int, Int}}(undef, t)
    e = ramification_index(p)
    aL = elem_type(K)[]
    uL = Int[]
    wL = Int[]
    witt = Int[]
    dets = elem_type(K)[]
    unif = elem_in_nf(uniformizer(p))
    pi = zero(K)
    for i in 1:t
      _sym[i] = (nrows(J[i]), sL[i], 0)
      push!(dets, det(G[i]))
      push!(witt, witt_invariant(quadratic_space(K, G[i]), p))
      GG = diagonal_matrix(eltype(G)[ j < i ? unif^(2*(sL[i] - sL[j])) * G[j] : G[j] for j in 1:t])
      D = diagonal(GG)
      m, pos = findmin(Union{Int, PosInf}[iszero(d) ? inf : 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(Union{PosInf, Int}[uL[i] + quadratic_defect(d//aL[i], p) for d in D])))
    end
    _f = Int[]
    for k in 1:(t - 1)
      exp = uL[k] + uL[k + 1]
      push!(_f, (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
    _normgen = aL
    _weight = wL
  end
  return pi, _sym, _weight, _normgen, _f, dets, witt, uL, J, G, E
end

################################################################################
#
#  Sum of genus symbols
#
################################################################################

function orthogonal_sum(G1::LocalGenusQuad, G2::LocalGenusQuad)
  @req prime(G1) === prime(G2) "Local genera must have the same prime ideal"
  if !G1.is_dyadic
    p = prime(G1)
    if uniformizer(G1) != uniformizer(G2)
      q = divexact(G2.uniformizer, G1.uniformizer)
      fl = is_local_square(q, p)
      local G2adj::Vector{Int}
      if fl
        G2adj = G2.detclasses
      else
        G2adj = Vector{Int}(undef, length(G2))
        for i in 1:length(G2)
          if isodd(G2.ranks[i] * G2.scales[i])
            G2adj[i] = -G2.detclasses[i]
          else
            G2adj[i] = G2.detclasses[i]
          end
        end
      end
    else
      G2adj = G2.detclasses
    end
    G3 = _direct_sum_easy(G1, G2, G2adj)
  else
    L1 = representative(G1)
    L2 = representative(G2)
    L3, = orthogonal_sum(L1, L2)
    G3 = genus(L3, prime(G1))
  end

  if isdefined(G1, :jordec) && isdefined(G2, :jordec)
    G3.jordec = orthogonal_sum(G1.jordec, G2.jordec)
  end

  return G3
end

function _direct_sum_easy(G1::LocalGenusQuad, G2::LocalGenusQuad, detclassesG2 = G2.detclasses)
  # We do a merge sort
  i1 = 1
  i2 = 1
  _sca = Int[]
  _rk = Int[]
  _detclass = Int[]
  while i1 <= length(G1) && i2 <= length(G2)
    if scale(G1, i1) < scale(G2, i2)
      push!(_sca, scale(G1, i1))
      push!(_rk, rank(G1, i1))
      push!(_detclass, G1.detclasses[i1])
      i1 += 1
    elseif scale(G2, i2) < scale(G1, i1)
      push!(_sca, scale(G2, i2))
      push!(_rk, rank(G2, i2))
      push!(_detclass, G2.detclasses[i2])
      i2 += 1
    else
      push!(_sca, scale(G1, i1))
      push!(_rk, rank(G1, i1) + rank(G2, i2))
      push!(_detclass, G1.detclasses[i1] * G2.detclasses[i2])
      i1 += 1
      i2 += 1
    end
  end

  if i1 <= length(G1)
    append!(_sca, G1.scales[i1:length(G1)])
    append!(_rk, G1.ranks[i1:length(G1)])
    append!(_detclass, G1.detclasses[i1:length(G1)])
  end

  if i2 <= length(G2)
    append!(_sca, G2.scales[i2:length(G2)])
    append!(_rk, G2.ranks[i2:length(G2)])
    append!(_detclass, G2.detclasses[i2:length(G2)])
  end

  return genus(QuadLat, prime(G1), uniformizer(G1), _rk, _sca, _detclass)
end

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

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

function _Amatrix(K, a, b)
  z = zero_matrix(K, 2, 2)
  z[1, 1] = a
  z[2, 2] = b
  z[1, 2] = 1
  z[2, 1] = 1
  return z
end

function _non_square(K, p)
  O = order(p)
  R, mR = ResidueField(O, p)
  u = elem_in_nf(mR\non_square(R))
  @assert !is_local_square(u, p)[1]
  return u
end

function representative(G::LocalGenusQuad)
  K = nf(order(G.p))
  return lattice(quadratic_space(K, gram_matrix(jordan_decomposition(G))))
end

######

mutable struct LocalGenusSymbol{S}
  P
  data
  x
  iseven::Bool
  E
  is_ramified
  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 (scale, rank, det) 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], is_square(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(Union{PosInf, Int}[iszero(d) ? inf : 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

################################################################################
#
#  Enumeration of local Jordan decompositions and genera
#
################################################################################

# Find a_1,...,a_n in p^n such that for y in p^n we have
# y = a_i \mod p^m for some i. (\mod = the "quadratic" equivalence)
function _representatives_for_equivalence_and_witt(p, n, m, c)
  @assert n <= m
  u = elem_in_nf(uniformizer(p))^n
  G, mG = local_multiplicative_group_modulo_squares(p)
  k = ngens(G)
  reps = typeof(u)[ u * mG(g) for g in G if iszero(g[k])]
  finer_reps = Tuple{typeof(u), Int}[(reps[1], hilbert_symbol(reps[1], reps[1] * c, p))]
  for j in 2:length(reps)
    uu = reps[j]
    new = true
    hi = hilbert_symbol(uu, uu * c, p)
    for k in 1:length(finer_reps)
      if _is_isometric_quadratic(uu, finer_reps[k][1], p, m) && finer_reps[k][2] == hi
        new = false
        break
      end
    end
    if new
      push!(finer_reps, (uu, hi))
    end
  end
  return first.(finer_reps)::Vector{typeof(u)}
end

function _is_isometric_quadratic(a, b, p, n)
  ap = valuation(a, p)
  if ap != valuation(b, p)
    return false
  end

  if quadratic_defect(divexact(a, b), p) >= n - ap
    return true
  end

  return false
end

# uni-unimodular, that is, determinant/discriminant equal to 1
function _unimodular_jordan_block(p, m)
  E = nf(order(p))
  e = ramification_index(p)
  @assert is_dyadic(p)
  # weight, normgen, det, witt
  res = Vector{Tuple{Int, nf_elem, nf_elem, Int}}()

  G, mG = local_multiplicative_group_modulo_squares(p)
  pi = elem_in_nf(uniformizer(p))
  k = ngens(G)
  reps_squares = typeof(pi)[ mG(g) for g in G if iszero(g[k])]

  if isodd(m)
    r = div(m - 1, 2)
    # Norm is always 0
    if m == 1
      weights = Int[e]
    else
      weights = Int[ i for i in 0:e if isodd(i)] # I do the sieving because of the Witt invariant later
      if !isodd(e)
        push!(weights, e)
      end
    end

    # The only norm generator is d

    D = kummer_generator_of_local_unramified_quadratic_extension(p)
    invD = inv(D)

    for w in weights
      for d in reps_squares
        push!(res, (w, (-1)^r * d, d, 1))
        if w != e
          # If w != e, then also -1 is possible
          push!(res, (w, invD * (-1)^r * d, d, -1))
        end
      end
    end
  else
    r = div(m, 2)

    norms = collect(0:e)

    _find_special_class_dict = Dict{nf_elem, nf_elem}()

    __find_special_class = (x, p) -> get!(_find_special_class_dict, x, _find_special_class(x, p))

    mmod4 = mod(m, 4)

    for n in norms
      normgens = pi^n .* reps_squares

      if isodd(n) # odd norm
        _r = Int[ j for j in n:e if iseven(j) ]
        if isodd(e) && (length(_r) == 0 || _r[end] != e)
          push!(_r, e)
        end
      else # iseven(s[3])
        _r = Int[ j for j in n:e if isodd(j) ]
        if iseven(e) && (length(_r) == 0 || _r[end] != e)
          push!(_r, e)
        end
      end
      weights = _r

      a = gen(E)

      for w in weights
        for d in reps_squares
          if mmod4 == 0 || mmod4 == 1
            gamma = __find_special_class(d, p) - 1
            #@assert is_local_square(d//(1 + gamma), p)
            #@assert valuation(d, p) == valuation(1 + gamma, p)
            dis = d
          else
            gamma = __find_special_class(-d, p) - 1
            #@assert is_local_square(-d//(1 + gamma), p)
            #@assert valuation(-d, p) == valuation(1 + gamma, p)
            dis = -d
          end
          #@assert quadratic_defect(1 + gamma, p) == (iszero(gamma) ? inf : valuation(gamma, p))
          if !(iszero(gamma) || valuation(gamma, p) >= w + n)
            continue
          end

          if iseven(n + w)
            @assert w == e
          end

          # This is a good candidate for the determinant
          # It is not sufficient to look for
          _normgens = _representatives_for_equivalence_and_witt(p, n, w, -dis)
          for ng in _normgens

            if r == 1 && isodd(w + n)
              if !(w == e || (e > w && !iszero(gamma) && w == valuation(gamma, p) - valuation(ng, p)))
                continue
              end
            end

            if (m == 2 && isodd(n + w))
              # A(\alpha, -\gamma \alpha^-1)
              wi = hilbert_symbol(ng, ng * -dis, p)
              push!(res, (w, ng, d, wi))
            elseif iseven(n + w)
              wi = hilbert_symbol(ng, ng * -dis, p)
              push!(res, (w, ng, d, wi))
            else
              push!(res, (w, ng, d, 1))
              push!(res, (w, ng, d, -1))
            end
          end
        end
      end
    end
  end
  return res
end

@doc Markdown.doc"""
    local_jordan_decomposition(E::NumField, p; rank::Int
                                               det_val::Int
                                               max_scale::Int)

Given a prime ideal $\mathfrak p$, returns all abstract Jordan decompositions
of rank `r` with determinant valuation `det_val` and scales of the blocks
bounded by `max_scale`.
"""
function local_jordan_decompositions(E, p; rank::Int, det_val::Int, max_scale = nothing)
  local _max_scale::Int

  if max_scale === nothing
    _max_scale = det_val
  else
    _max_scale = max_scale
  end

  scales_rks = Vector{Tuple{Int, Int}}[]

  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

  res = JorDec{typeof(E), typeof(p), elem_type(E)}[]

  if !is_dyadic(p)
    ns = _non_square(E, p)
    u = elem_in_nf(uniformizer(p))
    for scalerank in scales_rks
      _local_jordan_decompositions_nondyadic!(res, E, p, scalerank, ns, u)
    end
    return res
  else
    G, mG = local_multiplicative_group_modulo_squares(p)
    pi = elem_in_nf(uniformizer(p))
    k = ngens(G)
    reps_squares = typeof(pi)[ mG(g) for g in G if iszero(g[k])]
    @assert all(valuation(u, p) == 0 for u in reps_squares)

    e = ramification_index(p)

    # collect the possible ranks

    possible_ranks = unique!(reduce(vcat, Vector{Int}[[r[2] for r in s] for s in scales_rks]))

    decs_per_rank = Dict{Int, Vector{Tuple{Int, nf_elem, nf_elem, Int}}}()

    for m in possible_ranks
      decs_per_rank[m] = _unimodular_jordan_block(p, m)
    end

    for sr in scales_rks
      _local_jordan_decompositions_dyadic!(res, E, p, sr, G, mG, pi, k, reps_squares, e, possible_ranks, decs_per_rank)
    end
    return res
    #for srwndw in scales_rks_norms_weights_normgens_dets_witts
    #  push!(res, JorDec(p, Int[s[1] for s in srwndw], Int[s[2] for s in srwndw], nf_elem[pi^s[1] * s[5] for s in srwndw], Int[s[1] + s[4] for s in srwndw], nf_elem[pi^(s[1] * s[2]) * s[6] for s in srwndw], Int[isodd(s[2]) ? s[7] : s[7] * hilbert_symbol(pi^s[1], (-1)^divexact(s[2]*(s[2] - 1), 2) * pi^(s[1] * s[2]) * s[6], p) for s in srwndw]))
    #end
    #return res
  end
end

function _local_jordan_decompositions(E, p, scalerank)
  res = JorDec{typeof(E), typeof(p), elem_type(E)}[]
  if is_dyadic(p)
    _local_jordan_decompositions_dyadic!(res, E, p, scalerank)
  else
    _local_jordan_decompositions_nondyadic!(res, E, p, scalerank)
  end
  return res
end

function _local_jordan_decompositions_nondyadic!(res, E, p, scalerank)
  ns = _non_square(E, p)
  u = elem_in_nf(uniformizer(p))
  return _local_jordan_decompositions_nondyadic!(res, E, p, scalerank, ns, u)
end

function _local_jordan_decompositions_nondyadic!(res, E, p, scalerank, ns, u)
  class1 = elem_type(E)[u^(s[1] * s[2]) for s in scalerank]
  class2 = elem_type(E)[ns * u^(s[1] * s[2]) for s in scalerank]
  l = length(scalerank)
  @assert l <= sizeof(UInt) * 8
  # I need to compute all possiblities to distribute class1/class2
  # among the blocks.
  t = zero(UInt)
  for i in 1:2^l
    v = Vector{elem_type(E)}(undef, l)
    for j in 1:l
      if Bool((t >> (j - 1)) & 1)
        v[j] = class1[j]
      else
        v[j] = class2[j]
      end
    end
    push!(res, JorDec(p, Int[s[1] for s in scalerank], Int[s[2] for s in scalerank], v))
    t += 1
  end
end

function _local_jordan_decompositions_dyadic!(res, E, p, scalerank, G, mG, pi, k, reps_squares, e, possible_ranks, decs_per_rank)
  sr = scalerank

  it = Iterators.product([decs_per_rank[r] for (s, r) in sr]...)
  for local_blocks in it
    # local blocks in form weight, normgen, det, witt
    # JorDec(p, sc::Vector{Int},
    #           rks::Vector{Int},
    #           normgens::Vector{nf_elem},
    #           weights::Vector{Int},
    #           dets::Vector{nf_elem},
    #           witts::Vector{Int})

    l = length(sr)
    J = JorDec(p, Int[s[1] for s in sr],
               Int[s[2] for s in sr],
               nf_elem[pi^sr[i][1] * local_blocks[i][2] for i in 1:l],
               Int[sr[i][1] + local_blocks[i][1] for i in 1:l],
               nf_elem[pi^(sr[i][1] * sr[i][2]) * local_blocks[i][3] for i in 1:l],
               Int[isodd(sr[i][2]) ? local_blocks[i][4] : local_blocks[i][4] * hilbert_symbol(pi^sr[i][1], (-1)^divexact(sr[i][2]*(sr[i][2] - 1), 2) * pi^(sr[i][1] * sr[i][2]) * local_blocks[i][3], p) for i in 1:l])
    push!(res, J)
  end
end

function _local_jordan_decompositions_dyadic!(res, E, p, scalerank)
  G, mG = local_multiplicative_group_modulo_squares(p)
  pi = elem_in_nf(uniformizer(p))
  k = ngens(G)
  reps_squares = typeof(pi)[ mG(g) for g in G if iszero(g[k])]
  @assert all(valuation(u, p) == 0 for u in reps_squares)

  e = ramification_index(p)

  # collect the possible ranks

  possible_ranks = unique!(reduce(vcat, Vector{Int}[[r[2] for r in s] for s in [scalerank]]))

  decs_per_rank = Dict{Int, Vector{Tuple{Int, nf_elem, nf_elem, Int}}}()

  for m in possible_ranks
    decs_per_rank[m] = _unimodular_jordan_block(p, m)
  end

  return _local_jordan_decompositions_dyadic!(res, E, p, scalerank, G, mG, pi, k, reps_squares, e, possible_ranks, decs_per_rank)
end

function local_genera_quadratic(K, p; rank::Int, det_val::Int, max_scale = nothing)
  J = local_jordan_decompositions(K, p, rank = rank, det_val = det_val, max_scale = max_scale)
  res = local_quadratic_genus_type(K)[]
  for j in J
    g = genus(j)
    if !(g in res)
      push!(res, g)
    end
  end
  return res
end

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

mutable struct GenusQuad{S, T, U}
  K::S
  primes::Vector{T}
  LGS::Vector{LocalGenusQuad{S, T, U}}
  rank::Int
  signatures::Dict{InfPlc, Int}
  d::U
  space

  function GenusQuad{S, T, U}(K) where {S, T, U}
    z = new{typeof(K), ideal_type(order_type(K)), elem_type(K)}()
    z.rank = -1
    return z
  end
end

genus_quad_type(K) = GenusQuad{typeof(K), ideal_type(order_type(K)), elem_type(K)}

function GenusQuad(K, d, LGS, signatures)
  z = genus_quad_type(K)(K)
  z.LGS = LGS
  z.signatures = signatures
  z.d = d
  z.rank = rank(LGS[1])
  z.primes = [prime(g) for g in LGS]
  z.K = K
  return z
end

function genus(L::QuadLat{})
  return get_attribute!(L, :genus) do
    bad = bad_primes(L, even = true)
    S = real_places(base_field(L))
    D = diagonal(rational_span(L))
    signatures = Dict{InfPlc, Int}(s => count(d -> is_negative(d, s), D) for s in S)
    G = GenusQuad(base_field(L), prod(D), [genus(L, p) for p in bad], signatures)
    return G::genus_quad_type(base_field(L))
  end
end

function Base.:(==)(G1::GenusQuad, G2::GenusQuad)
  if G1.K != G2.K
    return false
  end

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

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

  for g1 in G1.LGS
    p1 = prime(g1)
    found = false
    for g2 in G2.LGS
      p2 = prime(g2)
      if p1 == p2
        found = true
        if g1 != g2
          return false
        end
        break
      end
    end
    if !found
      return false
    end
  end

  return true
end

primes(G::GenusQuad) = G.primes

function genera_quadratic(K; rank::Int, signatures, det)
  OK = maximal_order(K)

  _max_scale = 2 * det

  primes = support(2 * det)

  local_symbols = Vector{local_quadratic_genus_type(K)}[]

  ms = _max_scale
  ds = det
  for p in primes
    det_val = valuation(ds, p)
    mscale_val = valuation(ms, p)
    push!(local_symbols, local_genera_quadratic(K, p, rank = rank, det_val = det_val, max_scale = mscale_val))
  end

  res = genus_quad_type(K)[]
  it = Iterators.product(local_symbols...)
  for gs in it
    c = collect(gs)::Vector{local_quadratic_genus_type(K)}
    de = _possible_determinants(K, c, signatures)::Vector{nf_elem}
    for d in de
      b = _check_global_quadratic_genus(c, d, signatures)
      if b
        push!(res, GenusQuad(K, d, c, signatures))
      end
    end
  end

  return res
end

function _check_global_quadratic_genus(c, d, signatures)
  # c = vector or tuple of local quadratic genera
  # signatures = dict{pl => n}
  P = [ i for i in 1:length(c) if hasse_invariant(c[i]) == -1]
  if rank(c[1]) == 1 && length(P) == 0
    return false
  end

  if rank(c[1]) == 2
    for i in P
      #@show is_local_square(-d, prime(c[i]))
      if is_local_square(-d, prime(c[i]))
        return false
      end
    end
  end

  s = length([ s for (s, n) in signatures if (n % 4) in [2, 3]])
  return iseven(length(P) + s)
end

function _possible_determinants(K, local_symbols, signatures)
  # This is probably independent of the local symbol?
  C, mC = class_group(K)
  if length(local_symbols) == 0
    OK = maximal_order(K)
  else
    OK = order(prime(local_symbols[1]))
  end
  I = 1 * OK
  for g in local_symbols
    I = I * prime(g)^valuation(det(g), prime(g))
  end
  classofI = mC\(I)
  U, mU = unit_group_fac_elem(OK)
  classes = elem_type(C)[]
  for c in C
    if iszero(classofI + 2 * c)
      push!(classes, c)
    end
  end
  s = ideal_type(OK)[]
  current_classes = elem_type(C)[]
  primes = PrimeIdealsSet(OK, 2, -1)

  if id(C) in classes
    push!(s, 1 * OK)
    push!(current_classes, id(C))
  end

  for P in primes
    if length(classes) == length(current_classes)
      break
    end

    c = mC\(P)
    if c in current_classes
      continue
    end

    if c in classes
      push!(current_classes, c)
      push!(s, P)
    end
  end
  rlp = real_places(K)
  local R::GrpAbFinGen
  R, _exp, _log = infinite_primes_map(OK, rlp, 1 * OK)
  tar = R(Int[isodd(signatures[sigma]) ? 1 : 0 for sigma in rlp])
  gensU = gens(U)
  S, mS = sub(R, elem_type(R)[_log(mU(u)) for u in gensU], false)
  # I need totally positive units / OK^*2
  Q, mQ = quo(U, 2, false)
  f = hom(Q, R, elem_type(R)[_log(mU(mQ\u)) for u in gens(Q)])
  Ker, mKer = kernel(f, false)
  transver = elem_type(K)[ evaluate(mU(mQ\(mKer(k)))) for k in Ker]

  dets = elem_type(K)[]

  for P in s
    J = I * P^2
    fl, u = is_principal(J)
    @assert fl
    # I need to change u such that sign(u, sigma) = (-1)^signatures[sigma]
    v = _log(elem_in_nf(u)) + tar
    fl, y = haspreimage(mS, v)
    if fl
      z = elem_in_nf(u) * prod(elem_type(K)[ evaluate(mU(gensU[i]))^y.coeff[i] for i in 1:length(gensU)])
      @assert z * OK == J
      @assert all(sigma -> sign(z, sigma) == (-1)^signatures[sigma], rlp)
      for t in transver
        possible_det = z * t
        good = true
        for g in local_symbols
          # I need to test if the determinant is locally correct
          if !is_local_square(possible_det * det(g), prime(g))
            good = false
            break
          end
        end
        if good
          push!(dets, possible_det)
        end
      end
    end
  end
  return dets
end

function quadratic_space(G::GenusQuad)
  if isdefined(G, :space)
    return G.space::quadratic_space_type(G.K)
  end

  c = G.LGS
  K = G.K
  P = ideal_type(order_type(K))[ prime(c[i]) for i in 1:length(c) if hasse_invariant(c[i]) == -1]
  d = G.d
  signa = G.signatures
  rk = G.rank
  G.space = quadratic_space(G.K, _quadratic_form_with_invariants(rk, d, P, signa))
  return G.space::quadratic_space_type(G.K)
end

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

function representative(G::GenusQuad)
  K = G.K
  OK = order(primes(G)[1])
  # Let's follow the Lorch paper. This is also how we do it in the Hermitian case.
  V = quadratic_space(G)
  M = maximal_integral_lattice(V)
  for g in G.LGS
    p = prime(g)
    @vprint :Lattice 1 "Finding representative for $g at $(prime(g))...\n"
    L = representative(g)
    M = locally_isometric_sublattice(M, L, p)
    @assert is_locally_isometric(M, L, p)
  end
  return M
end

function locally_isometric_sublattice(M::QuadLat, L::QuadLat, p)
  k, h = ResidueField(order(p), p)
  m = rank(M)
  chain = typeof(L)[ L ]
  local LL::typeof(M)
  ok, LL = is_maximal_integral(L, p)
  E = nf(order(p))
  while !ok
    push!(chain, LL)
    ok, LL = is_maximal_integral(LL, p)
  end
  pop!(chain)
  LL = M
  reverse!(chain)
  for X in chain
    BM = local_basis_matrix(LL, p, type = :submodule)
    pM = _module_scale_ideal(pseudo_matrix(LL), fractional_ideal(order(p), p))
    while true
      local v::Vector{elem_type(k)}
      v = elem_type(k)[ rand(k) for i in 1:m ]
      while all(iszero, v)
        v = elem_type(k)[ rand(k) for i in 1:m ]
      end
      _, KM = kernel(matrix(k, length(v), 1, v), side = :left)
      KM = map_entries(x -> E(h\x)::elem_type(E), KM)
      _new_pmat = _sum_modules(pseudo_matrix(KM * BM), pM)
      LL = lattice(ambient_space(M), _new_pmat)
      if is_locally_isometric(X, LL, p)
        break
      end
    end
  end
  @assert is_locally_isometric(L, LL, p)
  return LL
end

function representatives(G::GenusQuad)
  return genus_representatives(representative(G))
end

################################################################################
#
#  Orthogonal sum of genus symbols
#
################################################################################

function orthogonal_sum(G1::GenusQuad{S, T, U}, G2::GenusQuad{S, T, U}) where {S, T, U}
  @req G1.K === G2.K "Global genus symbols must be defined over the same field"
  K = G1.K
  LGS = local_genus_quad_type(K)[]
  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)::Int
      g1 = G1.LGS[i]
    else
      @assert !is_dyadic(p)
      fl = is_local_square(G1.d, p)
      dcl = fl ? 1 : -1
      g1 = genus(QuadLat, p, [(0, rank(G1), dcl)])
    end

    if p in P2
      i = findfirst(g -> prime(g) == p, G2.LGS)
      g2 = G2.LGS[i]
    else
      @assert !is_dyadic(p)
      fl = is_local_square(G2.d, p)
      dcl = fl ? 1 : -1
      g2 = genus(QuadLat, p, [(0, rank(G2), dcl)])
    end
    g3 = orthogonal_sum(g1, g2)
    push!(LGS, g3)
  end
  sig1 = G1.signatures
  sig2 = G2.signatures
  sig3 = merge(+, sig1, sig2)

  return GenusQuad(K, G1.d * G2.d, LGS, sig3)
end

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

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

@doc Markdown.doc"""
    in(L::QuadLat, G::GenusQuad) -> Bool

Test if the lattice $L$ is contained in the genus $G$.
"""
Base.in(L::QuadLat, G::GenusQuad) = genus(L) == G


text:
# eigen solver types

abstract type AbstractEigenSolver end

struct FullEigen <: AbstractEigenSolver end

function Base.show(io::IO, ::MIME"text/plain", l::FullEigen)
    println(io, "FullEigen")
end

struct ArpackEigen <: AbstractEigenSolver
    n_first::Int
end

function Base.show(io::IO, ::MIME"text/plain", l::ArpackEigen)
    println(io, "ArpackEigen")
end

struct KrylovEigen <: AbstractEigenSolver
    n_first::Int
end

function Base.show(io::IO, ::MIME"text/plain", l::KrylovEigen)
    println(io, "KrylovEigen")
end


text:
using SpineOpt

template = SpineOpt.template()

template["object_classes"] = [x for x in template["object_classes"] if x[1] != "output"]

for url in ARGS
	SpineOpt.import_data(url, template, "Load SpineOpt template")
end

max_stars_repo_name: Energy-Reform/Mahdavi-test-systems

max_issues_repo_name: Energy-Reform/Mahdavi-test-systems

max_forks_repo_name: Energy-Reform/Mahdavi-test-systems


text:
using Mixers,
      Test,
      MacroTools,
      Parameters,
      Unitful

# @premix macro

q = quote @premix struct Premixdrinks{M,B}
       milkshake::M
       beer::B
    end
end

Mixers.defmacro(q, true)

@premix struct Premixdrinks{M,B}
   milkshake::M
   beer::B
end

@Premixdrinks mutable struct Drinks{J}
   juice::J
end

@test fieldnames(Drinks) == (:milkshake, :beer, :juice)
d = Drinks(1.9, 13, 7)
@test d.milkshake == 1.9

# Redeclaring doesn't error:
@Premixdrinks mutable struct Drinks{J}
   juice::J
end

# @mix macro
@mix struct Juice{P,S}
   mango::P
   lime::S
end
@mix struct NoJuice end

# Can also use with a parent struct <: Smoothie
abstract type Cocktail end
@Juice struct Daiquiri{S,M} <: Cocktail
    sugar::S
    rum::M
end
@test fieldnames(Daiquiri) == (:sugar, :rum, :mango, :lime)

# Type parameters work as expected
@test_throws MethodError Daiquiri(:none, 1.5, 2, "lots")
dq = Daiquiri(:none, 1.5, 2, :lots)
@test dq.mango == 2
@test dq.lime == :lots
@test dq.sugar == :none
@test dq.rum == 1.5

# Empty things stay empty
@NoJuice struct Empty end
@test fieldnames(Empty) == ()

# Empty things with {} can have type parameters added
abstract type AbstractPunch end
@Juice struct Punch{} <: AbstractPunch end
@test fieldnames(Punch) == (:mango, :lime)

# Inheritance
abstract type AbstractBeverage{G} end
abstract type AbstractGlass end
struct Lowball <: AbstractGlass end

@premix struct Liquid{L}
    liquid::L
end
@Liquid mutable struct Beverage{S,G<:AbstractGlass} <: AbstractBeverage{G}
    salt::S
    glass::G
end
@test fieldnames(Beverage) == (:liquid, :salt, :glass)
@test Beverage <: AbstractBeverage

# macro composition
@mix @with_kw struct Softdrinks{C,L}
    cola::C = 1.5u"L"
    lemonade::L = 2.0u"L"
end

# @with_kw works with no local macro
@Softdrinks struct Fridge{} end

fridge = Fridge()
@test fridge.cola == 1.5u"L"
@test fridge.lemonade == 2.0u"L"

# @with_kw is merged with local macro
@Softdrinks @with_kw struct Esky{B}
    beer::B = 6u"L"
end

esky = Esky()
@test esky.lemonade == 2.0u"L"
@test esky.beer == 6u"L"

# Nested @mix
@Softdrinks @mix struct BagODrinks{J}
    juice::J = 3u"L"
end

@BagODrinks struct Icebox{I}
    ice::I = 5u"kg"
end

@test fieldnames(Icebox) == (:ice, :juice, :cola, :lemonade)

icebox = Icebox()
@test icebox.cola == 1.5u"L"
@test icebox.lemonade == 2.0u"L"
@test icebox.juice == 3u"L"
@test icebox.ice == 5u"kg"

# @pour
@pour hello begin
    "Hello world"
end
function sayhi()
   @hello
end
@test sayhi() == "Hello world"


text:
@with_kw mutable struct CostModMDP{S, A} <: MDP{S,A}
    mdp::MDP{S,A}
    cost_model
end

POMDPs.initialstate(mdp::CostModMDP) = initialstate(mdp.mdp)
POMDPs.actions(mdp::CostModMDP) = actions(mdp.mdp)
POMDPs.actionindex(mdp::CostModMDP, a) = actionindex(mdp.mdp)

POMDPs.isterminal(mdp::CostModMDP, s) = isterminal(mdp.mdp, s)
POMDPs.discount(mdp::CostModMDP) = discount(mdp.mdp)

render(mdp::CostModMDP, s, a = nothing; kwargs...) = render(mdp.mdp, s, a; kwargs...)
render(mdp::CostModMDP; kwargs...) = render(mdp.mdp; kwargs...)

function POMDPs.gen(mdp::CostModMDP, s, a, rng::AbstractRNG = Random.GLOBAL_RNG; info=Dict())
    sp, r = gen(mdp.mdp, s, a, rng, info=info)
    cnew = mdp.cost_model((s=s, a=a, sp=sp, r=r, info=info))
    info["cost"] = cnew
    return (sp=sp, r=r)
end

@with_kw mutable struct CostModPOMDP{S, A, O} <: POMDP{S,A,O}
    pomdp::POMDP{S,A,O}
    cost_model
end

POMDPs.initialstate(mdp::CostModPOMDP) = initialstate(mdp.pomdp)
POMDPs.initialobs(mdp::CostModPOMDP, s) = initialobs(mdp.pomdp, s)

POMDPs.actions(mdp::CostModPOMDP) = actions(mdp.pomdp)
POMDPs.actionindex(mdp::CostModPOMDP, a) = actionindex(mdp.pomdp)

POMDPs.isterminal(mdp::CostModPOMDP, s) = isterminal(mdp.pomdp, s)
POMDPs.discount(mdp::CostModPOMDP) = discount(mdp.pomdp)

POMDPs.observation(mdp::CostModPOMDP, s) = observation(mdp.pomdp, s)

render(mdp::CostModPOMDP, s, a = nothing; kwargs...) = render(mdp.pomdp, s, a; kwargs...)
render(mdp::CostModPOMDP; kwargs...) = render(mdp.pomdp; kwargs...)

function POMDPs.gen(mdp::CostModPOMDP, s, a, rng::AbstractRNG = Random.GLOBAL_RNG; info=Dict())
    sp, o, r = gen(mdp.pomdp, s, a, rng, info=info)
    cnew = mdp.cost_model((s=s, a=a, sp=sp, r=r, info=info))
    info["cost"] = cnew
    return (sp=sp, o=o, r=r)
end

CostMod(mdp::T, cost_model) where T <: MDP = CostModMDP(mdp, cost_model)
CostMod(pomdp::T, cost_model) where T <: POMDP = CostModPOMDP(pomdp, cost_model)


text:
"""
    KLocalOperator

A representation of a local operator, acting on certain `sites`, each with
`hilb_dims` dimension. The local operator is written in matrix form in this
basis as `mat`. For every row of `mat` there are several non-zero values contained
in `mel`, and to each of those, `to_change` contains the sites that must change the basis
new value, contained in `new_value`
"""
struct KLocalOperator{SV,M,Vel,Vti,Vtc,Vtv,OC} <: AbsLinearOperator
    # list of sites on which this operator acts
    sites::SV

    # local hilbert space dimensions for each site
    hilb_dims::SV

    # Operator (as a matrix) represented
    mat::M

    # Matrix elements, per row
    mel::Vector{Vel}

    # List of indices that should be changed
    to_change::Vector{Vtc}

    # List of values of that should be changed
    new_values::Vector{Vtv}

    # Connections
    op_conns::Vector{OC}

    # List of state indices corresponding to mel
    # This is redundant, and equivalent to specifying to_change and new_values
    new_indices::Vector{Vti}
end

"""
    KLocalOperatorRow(sites::Vector, hilb_dims::Vector, operator)

Creates a KLocalOperator where connections are stored by row for the operator
    `operator` acting on `sites` each with `hilb_dims` local dimension.
"""
KLocalOperatorRow(sites::AbstractVector, hilb_dims::AbstractVector, operator) =
    KLocalOperatorRow(eltype(operator), sites, hilb_dims, operator)

function KLocalOperatorRow(T::Type{<:Number}, sites::AbstractVector, hilb_dims::AbstractVector, operator)
    # TODO Generalize to arbistrary spaces and not only uniform
    st  = NAryState(real(T), first(hilb_dims), length(sites))
    st1 = NAryState(real(T), first(hilb_dims), length(sites))

    mel         = Vector{Vector{T}}()
    new_indices = Vector{Vector{Int}}()
    to_change   = Vector{Vector{Vector{eltype(Int)}}}()
    new_values  = Vector{Vector{Vector{eltype(st)}}}()
    op_conns    = Vector{OpConnection{Vector{T},
                        Vector{eltype(Int)}, Vector{eltype(st)}}}()
    #SCT = StateChanges{eltype(to_change_els), eltype(new_values_els)}

    for (r, row) = enumerate(eachrow(operator))
        # create the containers for this row
        mel_els        = T[]
        mel_indices    = Int[]
        to_change_els  = Vector{Vector{Int}}()
        new_values_els = Vector{Vector{T}}()
        conns_els      = eltype(op_conns)()

        if abs(row[r]) > 10e-6
            push!(mel_els, row[r])
            push!(mel_indices, r)
            push!(to_change_els, Int[])
            push!(new_values_els, eltype(st)[])
            push!(conns_els, (row[r], Int[], eltype(st)[]))
        end

        set_index!(st, r)
        for (c, val) = enumerate(row)
            r == c && continue
            abs(val) < 10e-6 && continue

            set_index!(st1, c)
            cngd = Int[]
            nwvls = eltype(st)[]
            for (i, (v, vn)) = enumerate(zip(config(st), config(st1)))
                if v != vn
                    push!(cngd, sites[i])
                    push!(nwvls, vn)
                end
            end
            push!(mel_els, val)
            push!(mel_indices, c)
            push!(to_change_els, cngd)
            push!(new_values_els, nwvls)
            push!(conns_els, (val, cngd, nwvls))
        end

        push!(mel, mel_els)
        push!(new_indices, mel_indices)
        push!(to_change, to_change_els)
        push!(new_values, new_values_els)
        push!(op_conns, conns_els)
    end
    KLocalOperator(sites, hilb_dims, convert(Matrix{T}, operator |> collect),
                        mel, to_change, new_values, op_conns, new_indices)
end

## Accessors
"""
    sites(op::KLocalOperator)

Returns the vector of `Int` labelling the lattice sites on which this operator
acts, in no particular order.
"""
sites(op::KLocalOperator) = op.sites

operators(op::KLocalOperator) = (op,)

conn_type(top::Type{KLocalOperator{SV,M,Vel,Vti,Vtc,Vtv,OC}}) where {SV, M, Vel, Vti, Vtc, Vtv, OC} =
    OpConnection{Vel, eltype(Vtc), eltype(Vtv)}
conn_type(op::KLocalOperator{SV,M,Vel,Vti,Vtc,Vtv,OC}) where {SV, M, Vel, Vti, Vtc, Vtv, OC} =
    OpConnection{Vel, eltype(Vtc), eltype(Vtv)}

# Copy
function duplicate(op::KLocalOperator)
    KLocalOperator(deepcopy(op.sites), deepcopy(op.hilb_dims), deepcopy(op.mat),
                    deepcopy(op.mel), deepcopy(op.to_change),
                    deepcopy(op.new_values), deepcopy(op.op_conns),
                    deepcopy(op.new_indices))
end

##
function row_valdiff!(conn::OpConnection, op::KLocalOperator, v::State)
    # Find row index
    r = local_index(v, sites(op))

    # Return values
#mel = op.mel[r]
#    tc  = op.to_change[r]
#    nv  = op.new_values[r]

#    append!(conn, (mel, tc, nv))
    append!(conn, op.op_conns[r])
end

function row_valdiff_index!(conn::OpConnectionIndex, op::KLocalOperator, v::State)
    # Find row index
    r = local_index(v, sites(op))

    # Return values
    mel = op.mel[r]
    ids = op.new_indices[r]

    append!(conn, (mel, ids))
end

# sum
function sum_samesite!(op_l::KLocalOperator, op_r::KLocalOperator)
    @assert op_l.sites == op_r.sites
    @assert length(op_l.mel) == length(op_r.mel)

    op_l.mat .+= op_r.mat

    n = length(op_l.mel)
    for (i, (mel_vals_l, mel_vals_r)) = enumerate(zip(op_l.mel, op_r.mel))

        for (mel, ids, tch, nv) = zip(op_r.mel[i], op_r.new_indices[i],
                                      op_r.to_change[i], op_r.new_values[i])
            found = 0
            # if mel_val_l is empty, then for sure there is nothing in here
            if !isempty(mel_vals_l)
                ids_tch = findall(isequal(tch), op_l.to_change[i])
                for id=ids_tch
                    if op_l.new_values[i][id] == nv
                        found = id
                        break
                    end
                end
            end

            if found == 0
                push!(op_l.mel[i], mel)
                push!(op_l.new_indices[i], ids)
                push!(op_l.to_change[i], tch)
                push!(op_l.new_values[i], tch)
            else
                op_l.mel[i][found] += mel
            end
        end
    end

    for (i, (connsᵢ_l , connsᵢ_r)) = enumerate(zip(op_l.op_conns, op_r.op_conns))
        add!(connsᵢ_l, connsᵢ_r)
    end
    return op_l
end

sum_samesite(op_l::KLocalOperator, op_r::KLocalOperator) = sum_samesite!(duplicate(op_l), op_r)

Base.transpose(op::KLocalOperator) = KLocalOperatorRow(deepcopy(op.sites),
                                                       deepcopy(op.hilb_dims),
                                                       transpose(op.mat)|>collect)

function Base.conj!(op::KLocalOperator)
    conj!(op.mat)
    conj!(op.op_conns)
    for el=op.mel
        conj!.(el)
    end
    return op
end

Base.conj(op::KLocalOperator) = conj!(duplicate(op))

Base.adjoint(op::KLocalOperator) = conj(tranpose(op))

Base.eltype(::T) where {T<:KLocalOperator} = eltype(T)
Base.eltype(T::Type{KLocalOperator{SV,M,Vel,Vti,Vtc,Vtv,OC}}) where {SV,M,Vel,Vti,Vtc,Vtv,OC} =
    eltype(eltype(Vel))


text:
# Generic case for most constraints, excluding those on single variables.
function _info(::Optimizer, key::MOI.ConstraintIndex{MOI.VariableIndex, <:Any})
    throw(MOI.InvalidIndex(key))
    return
end

function _info(model::Optimizer, key::MOI.ConstraintIndex)
    if haskey(model.constraint_info, key)
        return model.constraint_info[key]
    end
    throw(MOI.InvalidIndex(key))
    return
end

function MOI.is_valid(
    model::Optimizer,
    c::MOI.ConstraintIndex{F, S},
) where {F <: MOI.AbstractFunction, S <: MOI.AbstractSet}
    info = get(model.constraint_info, c, nothing)
    return info !== nothing && typeof(info.set) == S
end

function MOI.add_constraint(
    model::Optimizer,
    f::F,
    s::S,
) where {F <: MOI.AbstractFunction, S <: MOI.AbstractSet}
    index = MOI.ConstraintIndex{F, S}(length(model.constraint_info) + 1)
    constr = _build_constraint(model, f, s)
    cpo_java_add(model.inner, constr)
    model.constraint_info[index] = ConstraintInfo(index, constr, f, s)
    return index
end

function MOI.get(
    model::Optimizer,
    ::MOI.ConstraintFunction,
    c::MOI.ConstraintIndex{F, S},
) where {F <: MOI.AbstractFunction, S <: MOI.AbstractSet}
    return _info(model, c).f
end

function MOI.get(
    model::Optimizer,
    ::MOI.ConstraintSet,
    c::MOI.ConstraintIndex{F, S},
) where {F <: MOI.AbstractFunction, S <: MOI.AbstractSet}
    return _info(model, c).set
end


text:
cd(@__DIR__); include("setups/rgc100.jl")
gr(dpi = 200)

plt = gplot(W, X; width = 1); scatter_gplot!(X; marker = 𝚽[:, 1142], ms = 4)
        plot!(xlims = [-180, 220], xlabel = "x(μm)", ylabel = "y(μm)",
        clims = (-0.12, 0.12), frame = :box, size = (550, 500), right_margin = 5mm)
savefig(plt, "../figs/RGC100_twotypes_semi_oscillation.png")

plt = gplot(W, X; width = 1); scatter_gplot!(X; marker = 𝚽[:, 1143], ms = 4,
        plotOrder = :s2l); plot!(xlims = [-180, 220], xlabel = "x(μm)",
        ylabel = "y(μm)", clims = (-0.3, 0.3), frame = :box, size = (550, 500),
        right_margin = 5mm)
savefig(plt, "../figs/RGC100_twotypes_localized.png")


text:
using BinDeps

@BinDeps.setup

QT_ROOT = get(ENV, "QT_ROOT", "")

loadqml = library_dependency("loadqml", aliases=["libloadqml"])

used_homebrew = false
if QT_ROOT == ""
  if is_apple()
    std_hb_root = "/usr/local/opt/qt5"
    if(!isdir(std_hb_root))
      if Pkg.installed("Homebrew") === nothing
        error("Homebrew package not installed, please run Pkg.add(\"Homebrew\")")
      end
      using Homebrew
      if !Homebrew.installed("qt5")
        Homebrew.add("qt5")
      end
      if !Homebrew.installed("cmake")
        Homebrew.add("cmake")
      end
      used_homebrew = true
      QT_ROOT = joinpath(Homebrew.prefix(), "opt", "qt5")
    else
      QT_ROOT = std_hb_root
    end
  end

  if is_linux()
    try
      run(pipeline(`cmake --version`, stdout=DevNull, stderr=DevNull))
      try
        run(pipeline(`qmake-qt5 --version`, stdout=DevNull, stderr=DevNull))
      catch
        run(pipeline(`qmake --version`, stdout=DevNull, stderr=DevNull))
      end
      run(pipeline(`qmlscene $(joinpath(dirname(@__FILE__), "imports.qml"))`, stdout=DevNull, stderr=DevNull))
    catch
      println("Installing Qt and cmake...")

      function printrun(cmd)
        println("Running install command, if this fails please run manually:\n$cmd")
        run(cmd)
      end

      if BinDeps.can_use(AptGet)
        printrun(`sudo apt-get install cmake cmake-data qtdeclarative5-dev qtdeclarative5-qtquick2-plugin qtdeclarative5-dialogs-plugin qtdeclarative5-controls-plugin qtdeclarative5-quicklayouts-plugin qtdeclarative5-window-plugin qmlscene qt5-default`)
      elseif BinDeps.can_use(Pacman)
        printrun(`sudo pacman -S --needed qt5-quickcontrols2`)
      elseif BinDeps.can_use(Yum)
        printrun(`sudo yum install cmake qt5-qtbase-devel qt5-qtquickcontrols qt5-qtquickcontrols2-devel`)
      end
    end
  end
end

cmake_prefix = QT_ROOT

prefix=joinpath(BinDeps.depsdir(loadqml),"usr")
loadqml_srcdir = joinpath(BinDeps.depsdir(loadqml),"src","loadqml")
loadqml_builddir = joinpath(BinDeps.depsdir(loadqml),"builds","loadqml")
lib_prefix = @static is_windows() ? "" : "lib"
lib_suffix = @static is_windows() ? "dll" : (@static is_apple() ? "dylib" : "so")

makeopts = ["--", "-j", "$(Sys.CPU_CORES+2)"]

println("Using Qt from $cmake_prefix")

qml_steps = @build_steps begin
	`cmake -DCMAKE_INSTALL_PREFIX="$prefix" -DCMAKE_BUILD_TYPE="Debug" -DCMAKE_PREFIX_PATH="$cmake_prefix" $loadqml_srcdir`
	`cmake --build . --config Debug --target install $makeopts`
end

# If built, always run cmake, in case the code changed
if isdir(loadqml_builddir)
  BinDeps.run(@build_steps begin
    ChangeDirectory(loadqml_builddir)
    qml_steps
  end)
end

provides(BuildProcess,
  (@build_steps begin
    CreateDirectory(loadqml_builddir)
    @build_steps begin
      ChangeDirectory(loadqml_builddir)
      FileRule(joinpath(prefix,"lib", "$(lib_prefix)loadqml.$lib_suffix"),qml_steps)
    end
  end),loadqml)

@BinDeps.install Dict([(:loadqml, :_l_loadqml)])

if used_homebrew
  # Not sure why this is needed on homebrew Qt, but here goes:
  envfile_path = joinpath(dirname(@__FILE__), "env.jl")
  plugins_path = joinpath(QT_ROOT, "plugins")
  qml_path = joinpath(QT_ROOT, "qml")
  open(envfile_path, "w") do envfile
    println(envfile, "ENV[\"QT_PLUGIN_PATH\"] = \"$plugins_path\"")
    println(envfile, "ENV[\"QML2_IMPORT_PATH\"] = \"$qml_path\"")
  end
end


text:
using Test
using CartesianMesh
using PolynomialBasis
using ImplicitDomainQuadrature
# using Revise
using CutCell

function displacement(alpha, x::V) where {V<:AbstractVector}
    u1 = alpha * x[2] * sin(pi * x[1])
    u2 = alpha * (x[1]^3 + cos(pi * x[2]))
    return [u1, u2]
end

function displacement(alpha, x::M) where {M<:AbstractMatrix}
    npts = size(x)[2]
    return hcat([displacement(alpha, x[:, i]) for i = 1:npts]...)
end

function body_force(lambda, mu, alpha, x)
    b1 = alpha * (lambda + 2mu) * pi^2 * x[2] * sin(pi * x[1])
    b2 =
        -alpha * (6mu * x[1] + (lambda + mu) * pi * cos(pi * x[1])) +
        alpha * (lambda + 2mu) * pi^2 * cos(pi * x[2])
    return [b1, b2]
end

function displacement_nodeids(mesh)
    bn = CutCell.bottom_boundary_node_ids(mesh)
    tn = CutCell.top_boundary_node_ids(mesh)
    ln = CutCell.left_boundary_node_ids(mesh)
    return unique!(vcat(bn, tn, ln))
end

function is_traction_boundary(x)
    if x[1] ≈ 1.0
        return true
    else
        return false
    end
end

function right_face_traction_force(lambda, mu, alpha, x)
    s11 =
        (lambda + 2mu) * alpha * pi * x[2] * cos(pi * x[1]) -
        lambda * alpha * pi * sin(pi * x[2])
    s12 = alpha * mu * (3x[1]^2 + sin(pi * x[1]))

    return [s11, s12]
end

function add_cell_error_squared!(
    err,
    interpolater,
    exactsolution,
    cellmap,
    quad,
    detjac,
)
    for (p, w) in quad
        numsol = interpolater(p)
        exsol = exactsolution(cellmap(p))
        err .+= (numsol - exsol) .^ 2 * detjac * w
    end
end

function mesh_L2_error(
    nodalsolutions,
    exactsolution,
    nodalconnectivity,
    cellmaps,
    basis,
    errorquad,
)

    ndofs = size(nodalsolutions)[1]
    detjac = CutCell.determinant_jacobian(cellmaps[1])
    err = zeros(ndofs)
    interpolater = InterpolatingPolynomial(ndofs, basis)
    for (cellid, cellmap) in enumerate(cellmaps)
        nodeids = nodalconnectivity[:, cellid]
        elementsolutions = nodalsolutions[:, nodeids]
        update!(interpolater, elementsolutions)
        add_cell_error_squared!(
            err,
            interpolater,
            exactsolution,
            cellmap,
            errorquad,
            detjac,
        )
    end
    return sqrt.(err)
end

function apply_displacement_boundary_condition!(
    matrix,
    rhs,
    displacement_function,
    femesh,
)

    boundarynodeids = displacement_nodeids(femesh)
    nodalcoordinates = CutCell.nodal_coordinates(femesh)
    boundarynodecoordinates = nodalcoordinates[:, boundarynodeids]
    boundarydisplacement = displacement_function(boundarynodecoordinates)
    CutCell.apply_dirichlet_bc!(
        matrix,
        rhs,
        boundarynodeids,
        boundarydisplacement,
    )
end

function linear_system(
    basis,
    quad,
    facequads,
    stiffness,
    femesh,
    bodyforcefunc,
    tractionfunc,
    istractionboundary,
)

    cellmaps = CutCell.cell_maps(femesh)
    nodalconnectivity = CutCell.nodal_connectivity(femesh)
    cellconnectivity = CutCell.cell_connectivity(femesh)

    sysmatrix = CutCell.SystemMatrix()
    sysrhs = CutCell.SystemRHS()

    cellmatrix = CutCell.bilinear_form(basis, quad, stiffness, cellmaps[1])
    CutCell.assemble_bilinear_form!(sysmatrix, cellmatrix, nodalconnectivity, 2)
    CutCell.assemble_body_force_linear_form!(
        sysrhs,
        bodyforcefunc,
        basis,
        quad,
        cellmaps,
        nodalconnectivity,
    )
    CutCell.assemble_traction_force_linear_form!(
        sysrhs,
        tractionfunc,
        basis,
        facequads,
        cellmaps,
        nodalconnectivity,
        cellconnectivity,
        istractionboundary,
    )

    K = CutCell.make_sparse(sysmatrix, femesh)
    R = CutCell.rhs(sysrhs, femesh)

    return K, R
end

function error_for_num_elements(numelements, polyorder, numqp)
    x0 = [0.0, 0.0]
    widths = [1.0, 1.0]
    nelements = [numelements, numelements]
    mesh = UniformMesh(x0, widths, nelements)
    stiffness = plane_strain_voigt_hooke_matrix(lambda, mu)

    basis = TensorProductBasis(2, polyorder)
    quad = tensor_product_quadrature(2, numqp)
    facequads = CutCell.face_quadratures(numqp)
    errorquad = tensor_product_quadrature(2, numqp + 2)

    femesh = CutCell.Mesh(mesh, basis)

    matrix, rhs = linear_system(
        basis,
        quad,
        facequads,
        stiffness,
        femesh,
        x -> body_force(lambda, mu, alpha, x),
        x -> right_face_traction_force(lambda, mu, alpha, x),
        is_traction_boundary,
    )
    apply_displacement_boundary_condition!(
        matrix,
        rhs,
        x -> displacement(alpha, x),
        femesh,
    )

    sol = matrix \ rhs
    nodalsolutions = reshape(sol, 2, :)

    nodalconnectivity = CutCell.nodal_connectivity(femesh)
    cellmaps = CutCell.cell_maps(femesh)

    err = mesh_L2_error(
        nodalsolutions,
        x -> displacement(alpha, x),
        nodalconnectivity,
        cellmaps,
        basis,
        errorquad,
    )
    return err
end

function required_quadrature_order(polyorder)
    ceil(Int, 0.5 * (2polyorder + 1))
end

function convergence(nelements, polyorder)
    numqp = required_quadrature_order(polyorder)
    err = zeros(length(nelements), 2)
    for (idx, nelement) in enumerate(nelements)
        err[idx, :] = error_for_num_elements(nelement, polyorder, numqp)
    end
    dx = 1.0 ./ nelements
    u1err = err[:,1]
    u2err = err[:,2]

    u1rate = sum(diff(log.(u1err)) ./ diff(log.(dx)))/(length(nelements)-1)
    u2rate = sum(diff(log.(u2err)) ./ diff(log.(dx)))/(length(nelements)-1)

    return round(u1rate,digits=3),round(u2rate,digits=3)
end

const lambda = 1.0
const mu = 2.0
const alpha = 0.01

numelements = [2,4,8,16,32,64]
u1rate, u2rate = convergence(numelements,1)
# println("Convergence of linear elements : ", u1rate, "    ", u2rate)
@test isapprox(u1rate,2.0,atol=0.05)
@test isapprox(u2rate,2.0,atol=0.05)

u1rate, u2rate = convergence(numelements,2)
# println("Convergence of quadratic elements : ", u1rate, "    ", u2rate)
@test isapprox(u1rate,3.0,atol=0.05)
@test isapprox(u2rate,3.0,atol=0.05)

# u1rate, u2rate = convergence(numelements,3)
# # println("Convergence of cubic elements : ", u1rate, "    ", u2rate)
# @test isapprox(u1rate,4.0,atol=0.05)
# @test isapprox(u2rate,4.0,atol=0.05)
#
#
# u1rate, u2rate = convergence(numelements,4)
# # println("Convergence of quartic elements : ", u1rate, "    ", u2rate)
# @test isapprox(u1rate,5.0,atol=0.1)
# @test isapprox(u2rate,5.0,atol=0.1)


text:
using RobustDataPipelines
using Test

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


text:
module LaTeXWriterTests

using Test
import Documenter.Writers.LaTeXWriter

@testset "file ordering" begin
    # Single page returns a single page
    @test LaTeXWriter.files(["a.md"]) == [("", "a.md", 1)]

    # Multiple pages at the top-level
    @test LaTeXWriter.files(["a.md", "b.md"]) ==
        [("", "a.md", 1), ("", "b.md", 1)]

    # Single header
    @test LaTeXWriter.files(["A" => "a.md"]) == [("A", "a.md", 1)]

    # Single page and a header
    @test LaTeXWriter.files(["a.md", "B" => "b.md"]) ==
        [("", "a.md", 1), ("B", "b.md", 1)]

    # Single page and a vector
    @test LaTeXWriter.files(["a.md", "B" => ["b.md", "c.md"]]) ==
        [("", "a.md", 1), ("B", "", 1), ("", "b.md", 2), ("", "c.md", 2)]

    # Multiple layers of nesting
    @test LaTeXWriter.files(["a.md", "B" => ["b.md", "C" => ["c.md"]]]) == [
        ("", "a.md", 1),
        ("B", "", 1),
        ("", "b.md", 2),
        ("C", "", 2),
        ("", "c.md", 3),
    ]
end

end


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    AbstractPath

A path on a spatial domain of type `D`.
"""
abstract type AbstractPath{D<:AbstractDomain} end

"""
    Base.iterate(path, state=1)

Iterate `path` from a given `state`.
"""
Base.iterate(::AbstractPath, state=1) = @error "not implemented"

"""
    Base.length(path)

Return the length of a `path`.
"""
Base.length(p::AbstractPath) = npoints(p.domain)

"""
    Base.eltype(::Type{P})

Return element type of path type `P`.
"""
Base.eltype(::Type{P}) where {P<:AbstractPath} = Int

#------------------
# IMPLEMENTATIONS
#------------------
include("paths/linear.jl")
include("paths/random.jl")
include("paths/source.jl")
include("paths/shifted.jl")


text:
pushfirst!(LOAD_PATH, joinpath(@__DIR__, "../.."))

using OceanBoundaryLayerParameterEstimation
using ParameterEstimocean
using ParameterEstimocean.EnsembleKalmanInversions: NaNResampler, FullEnsembleDistribution
using Oceananigans.TurbulenceClosures.CATKEVerticalDiffusivities: CATKEVerticalDiffusivity
using LinearAlgebra, CairoMakie

examples_path = joinpath(pathof(ParameterEstimocean), "../../examples")
include(joinpath(examples_path, "intro_to_inverse_problems.jl"))

###
### Pick a parameter set defined in `./utils/parameters.jl`
###

parameter_set = CATKEParametersRiDependent
closure = closure_with_parameters(CATKEVerticalDiffusivity(Float64;), parameter_set.settings)

###
### Pick the secret `true_parameters`
###

true_parameters = (Cᵟu = 0.5, CᴷRiʷ = 1.0, Cᵂu★ = 2.0, CᵂwΔ = 1.0, Cᴷeʳ = 5.0, Cᵟc = 0.5, Cᴰ = 2.0, Cᴷc⁻ = 0.5, Cᴷe⁻ = 0.2, Cᴷcʳ = 3.0, Cᴸᵇ = 1.0, CᴷRiᶜ = 1.0, Cᴷuʳ = 4.0, Cᴷu⁻ = 0.8, Cᵟe = 0.5)
true_closure = closure_with_parameters(closure, true_parameters)

###
### Generate and load synthetic observations
###

kwargs = (tracers = (:b, :e), Δt = 10.0, stop_time = 1day)

transformation = (b = Transformation(normalization=ZScore()),
                  u = Transformation(normalization=ZScore()),
                  v = Transformation(normalization=ZScore()),
                  e = Transformation(normalization=RescaledZScore(1e-1)))

# Note: if an output file of the same name already exist, `generate_synthetic_observations` will return the existing path and skip re-generating the data.
observ_path1 = generate_synthetic_observations("perfect_model_observation1"; Qᵘ = 3e-5, Qᵇ = 7e-9, f₀ = 1e-4, closure = true_closure, kwargs...)
observ_path2 = generate_synthetic_observations("perfect_model_observation2"; Qᵘ = -2e-5, Qᵇ = 3e-9, f₀ = 0, closure = true_closure, kwargs...)
observation1 = SyntheticObservations(observ_path1, field_names = (:b, :e, :u, :v), transformation)
observation2 = SyntheticObservations(observ_path2, field_names = (:b, :e, :u, :v), transformation)
observations = [observation1, observation2]

###
### Build Ensemble Simulation
###

ensemble_simulation, closure = build_ensemble_simulation(observations; Nensemble = 100, architecture = GPU())

###
### Build Inverse Problem
###

build_prior(name) = ScaledLogitNormal(bounds=bounds(name, parameter_set))
free_parameters = FreeParameters(named_tuple_map(names(parameter_set), build_prior))

# Pack everything into Inverse Problem `calibration`
track_times = Int.(floor.(range(1, stop = length(observations[1].times), length = 3)))
calibration = InverseProblem(observations, ensemble_simulation, free_parameters, output_map = ConcatenatedOutputMap(track_times));

###
### Run EKI
###

noise_covariance = 0.01
resampler = NaNResampler(; abort_fraction=0.5, distribution=FullEnsembleDistribution())

eki = EnsembleKalmanInversion(calibration; noise_covariance = noise_covariance,
                                           resampler = resampler)
params = iterate!(eki; iterations = 5)

###
### Summary Plots
###

using CairoMakie
using LinearAlgebra

directory = "calibrate_catke_to_perfect_model/"
isdir(directory) || mkpath(directory)

visualize!(calibration, true_parameters;
    field_names = (:u, :v, :b, :e),
    directory,
    filename = "perfect_model_visual_true_params.png"
)

visualize!(calibration, params;
    field_names = [:u, :v, :b, :e],
    directory,
    filename = "perfect_model_visual_calibrated.png"
)
@show params

plot_parameter_convergence!(eki, directory, true_parameters, n_columns=3)
plot_pairwise_ensembles!(eki, directory, true_parameters)
plot_error_convergence!(eki, directory, true_parameters)


text:
if Base.libllvm_version >= v"7.0"
    include(joinpath("gcn_intrinsics", "math.jl"))
end
include(joinpath("gcn_intrinsics", "indexing.jl"))
include(joinpath("gcn_intrinsics", "assertion.jl"))
include(joinpath("gcn_intrinsics", "synchronization.jl"))
include(joinpath("gcn_intrinsics", "extras.jl"))


text:
#**************************************************************************************
# GeoEfficiency.jl package
# ========================
#
# represent a fast and flexible tool to calculate in batch or individually the geometrical efficiency
# for a set of common radiation detectors shapes (cylindrical, Bore-hole, Well-type) as seen form a sources.
# the source can be a point, a disc or even a cylinder.
#
#**************************************************************************************

module GeoEfficiency
print_with_color(:black,"""\n
    \t *************************************************
    \t **            -=) GeoEfficiency (=-            **
    \t **  Accurate Geometrical Efficiency Calculator **
    \t **                                             **
    \t *************************************************

  @author: Mohamed Krar
  @Profile: https://www.researchgate.net/profile/Mohamed_Krar3
  @repositry: https://github.com/DrKrar/GeoEfficiency.jl/
  @documentation: http://geoefficiencyjl.readthedocs.org
  @version: v"0.8.3"
  Created on Fri Aug 14 20:12:01 2015

  loading the Package:-
  ====================
  Pkg.add("GeoEfficiency")
  using GeoEfficiency

  Quick usage:-
  ============
  geoEff()	: calculate once and return the value.
  calc()	: calculate once and display full information on the console.
  calcN()	: calculate untill you quit.
  batch()	: calculate using data in the "GeoEfficiency" folder.
""")

export 	input,
		getfloat,
		getDetectors,

		Point,
		source,
		RadiationDetector,
		CylDetector,
		BoreDetector,
		WellDetector,
		detectorFactory,

		geoEff,

		calc,
		calcN,
		batch

include("Physical_model.jl")
include("Input_Interface.jl")
include("calculations.jl")
include("Output_Interface.jl")

end #module


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

mutable struct IoK8sApiCoreV1SecretProjection <: SwaggerModel
    items::Any # spec type: Union{ Nothing, Vector{IoK8sApiCoreV1KeyToPath} } # spec name: items
    name::Any # spec type: Union{ Nothing, String } # spec name: name
    optional::Any # spec type: Union{ Nothing, Bool } # spec name: optional

    function IoK8sApiCoreV1SecretProjection(;items=nothing, name=nothing, optional=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1SecretProjection, Symbol("items"), items)
        setfield!(o, Symbol("items"), items)
        validate_property(IoK8sApiCoreV1SecretProjection, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        validate_property(IoK8sApiCoreV1SecretProjection, Symbol("optional"), optional)
        setfield!(o, Symbol("optional"), optional)
        o
    end
end # type IoK8sApiCoreV1SecretProjection

const _property_map_IoK8sApiCoreV1SecretProjection = Dict{Symbol,Symbol}(Symbol("items")=>Symbol("items"), Symbol("name")=>Symbol("name"), Symbol("optional")=>Symbol("optional"))
const _property_types_IoK8sApiCoreV1SecretProjection = Dict{Symbol,String}(Symbol("items")=>"Vector{IoK8sApiCoreV1KeyToPath}", Symbol("name")=>"String", Symbol("optional")=>"Bool")
Base.propertynames(::Type{ IoK8sApiCoreV1SecretProjection }) = collect(keys(_property_map_IoK8sApiCoreV1SecretProjection))
Swagger.property_type(::Type{ IoK8sApiCoreV1SecretProjection }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1SecretProjection[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1SecretProjection }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1SecretProjection[property_name]

function check_required(o::IoK8sApiCoreV1SecretProjection)
    true
end

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


text:
function interpolate(obj::MultivariateInterpolator, arg0::Vector{Vector{jdouble}}, arg1::Vector{jdouble})
    return jcall(obj, "interpolate", MultivariateFunction, (Vector{Vector{jdouble}}, Vector{jdouble}), arg0, arg1)
end


text:
module IcyAntics

export ContinuumIceModel, ViscoElasticRheology

include("ContinuumIceModels.jl")

using .ContinuumIceModels

end # module


text:
__precompile__(true)

module Discreet

using StatsBase
using Random

export
    entropy,
    estimate_entropy,
    estimate_joint_entropy,
    mutual_information_contingency,
    mutual_information

### Source files
include("entropy.jl")
include("mutual_information.jl")

end # module