text:
# src/qccdevcontrol.jl
# Created by Dirk Oliver Theis, March 19, 2021
# MIT license
# Sub-module QCCDevCtrl

module QCCDDevControl

export QCCDevCtrl

using ..QCCDevDes_Types
using ..QCCDevControl_Types
using ..QCCDev_Utils
using ..QCCDev_Feasible

include("initFunctions.jl")
include("devCtrlAux.jl")
"""
This sub-module provides the type `QCCDevCtrl` and functions for controlling the operation of the
simulated quantum device.

# Exported
- `QCCDevCtrl()`

# Not Exported Interface
- `load()`
- `linear_transport()`, `junction_transport()`,
- `swap()`
- `Rz()`, `Rxy()`, `XX()`, `ZZ()`

# Todo
* Visualization interface
"""

struct QCCDevCtrl__Operation_Not_Allowed_Exception end

####################################################################################################

"""
Function `QCCDevCtrl(::QCCDevDescription ; simulate::Symbol, 𝑜𝑝𝑡𝑖𝑜𝑛𝑠)`

Constructor; initializes an "empty" QCCD as described, with no ions loaded (yet).

# Arguments

- `simulate::Symbol` — one of `:No`, `:PureStates`, `:MixedStates`
- `qnoise_estimate::Bool` — whether estimation of noise takes place

Setting both `simulate=:No` and `qnoise_estimate=false` allows
feasibility check of a schedule.

## Options:
Currently none.  Possible:
- Modify default noise model (in case of `:MixedStates` simulation
- Modify default qnoise parameters
"""
function QCCDevCtrl(qdd             ::QCCDevDescription
                    ;
                    simulate ::Symbol = :No,
                    qnoise_estimate::Bool = false        ) ::QCCDevControl

    @assert simulate        ∈ [:No, :PureStates, :MixedStates]
    @assert qnoise_estimate ∈ [true,false] # 😃

    #-------------------------------------------------------------------#
    # TODO                                                              #
    #                                                                   #
    # Check whether simulation resources are sufficient to accommodate  #
    # the number of qubits (in pure states, mixed states, or tensor     #
    # network (cuQuantum) simulation)                                   #
    #                                                                   #
    #-------------------------------------------------------------------#

    # Initializes devices components
    gateZones = _initGateZone(qdd.gateZone)
    junctions = _initJunctions(qdd.gateZone.gateZones, qdd.auxZone.auxZones, qdd.loadZone.loadZones, qdd.junction.junctions)
    auxZones = _initAuxZones(qdd.auxZone)
    loadingZones = _initLoadingZones(qdd.loadZone)

    # graph = initGraph(qdd)

    # Check errors
    _checkInitErrors(junctions, auxZones, gateZones, loadingZones)

    # Initalizate QCCDevCtrl
    return QCCDevControl(qdd,
                      simulate, qnoise_estimate,
                      gateZones,junctions,auxZones, loadingZones)
end #^ QCCDevCtrl()

####################################################################################################

"""
Function `load()` — loads an ion into the device

# Arguments
* `t::Time_t` — time at which the operation commences.  Must be no earlier than the latest time
  given to previous function calls.

The function returns a named tuple consisting of:
- `new_ion_idx` — the index identifying the new ion
- `t₀` — the time at which the loaded ion will be usable (for transport off the loading zone);
"""
function load(qdc           ::QCCDevControl,
              t             ::Time_t,
              loading_zone  ::Symbol       )  ::@NamedTuple{new_ion_idx::Int,t₀::Time_t}

  # Checks
  isallowed_load(qdc, loading_zone, t)

  # Create new qubit
  local qubit = initQubit(loading_zone)
  qdc.qubits[qubit.id] = deepcopy(qubit)
  qdc.loadingZones[loading_zone].hole = qubit.id

  # Compute and actualize time
  local t₀ = compute_time(qdc, t, OperationTimes[:load])

  return (new_ion_idx=qubit.id, t₀)
end #^ module
# EOF

####################################################################################################

"""
Function `linear_transport()` — moves ions between zones/junctions.

# Arguments
* `t::Time_t` — time at which the operation commences.  Must be no earlier than the latest time
  given to previous function calls.
* `ion_idx`   — index (1-based) of ion to be moved.
* `edge-idx`  — index (1-based) of edge to move along.

The function returns the time at which the operation will be completed.
"""

function linear_transport(qdc           :: QCCDevControl,
                          t             :: Time_t,
                          ion_idx       :: Int,
                          destination_idx      :: Symbol) ::Time_t
  # Checks
  isallowed_linear_transport(qdc, t, ion_idx, destination_idx)

  ion = qdc.qubits[ion_idx]
  origin = giveZone(qdc, ion.position)
  destination = giveZone(qdc, destination_idx)

  # Remove ion from origin, insert it to destination,
  # and check if it has arrived to its destination
  _move_ion(ion, origin, destination)

  # Compute and actualize time
  local t₀ = compute_time(qdc, t, OperationTimes[:linear_transport])

  return t₀
end

####################################################################################################

"""
Function `junction_transport()` — moves around a junction.

# Arguments
* `t::Time_t` — time at which the operation commences.  Must be no earlier than the latest time
  given to previous function calls.
* `ion_idx`   — index (1-based) of ion to be moved.
* `edge_idx` — index (1-based) of edge identifying the edge on which the ion leaves the junction.

The function returns the time at which the operation will be completed.
"""
function junction_transport(qdc           :: QCCDevControl,
                            t             :: Time_t,
                            ion_idx       :: Int,
                            edge_idx      :: Int       ) ::Time_t

end

####################################################################################################

"""
Function `swap()` — physically swaps the positions of two ions

# Arguments
* `t::Time_t` — time at which the operation commences.  Must be no earlier than the latest time
  given to previous function calls.
* `ion𝑖_idx`, 𝑖=1,2, the (1-based) indices of the two ions.  Must be in the same gate zone.

The function returns the time at which the operation will be completed.
"""
function swap(qdc           :: QCCDevControl,
              t             :: Time_t,
              ion1_idx      :: Int,
              ion2_idx      :: Int       ) ::Time_t
  # Checks
  isallowed_swap(qdc, ion1_idx, ion2_idx, t)

  # Swap qubits
  _swap_ions(qdc, ion1_idx, ion2_idx)

  # Compute and actualize time
  local t₀ = compute_time(qdc, t, OperationTimes[:swap])

  return t₀
end

####################################################################################################

"""
Function `split()` — move ion out of gate zone into edge

# Arguments

* `t::Time_t` — time at which the operation commences.  Must be no earlier than the latest time
  given to previous function calls.
* `ion_idx` (1-based) index of the ion.
* `edge_idx` (1-based) index of edge at the end of which the ion will sit after split.

The function returns the time at which the operation will be completed.
"""
function split(qdc           :: QCCDevControl,
               t             :: Time_t,
               ion_idx       :: Int,
               edge_idx      :: Int) ::Time_t

end

####################################################################################################

"""
Function `merge()` — move ion out of gate zone into edge

# Arguments
* `t::Time_t` — time at which the operation commences.  Must be no earlier than the latest time
  given to previous function calls.
* `ion_idx`   — (1-based) index of the ion.
* `edge_idx`  — (1-based) index of edge at the end of which the ion sits before merge.

The function returns the time at which the operation will be completed.
"""
function merge(qdc           :: QCCDevControl,
               t             :: Time_t,
               ion_idx       :: Int,
               edge_idx      :: Int) ::Time_t

end

####################################################################################################

"""
Function `Rz` — single qubit Z-rotation

# Arguments
* `t::Time_t`   — time at which the operation commences.  Must be no earlier than
  the latest time given to previous function calls.
* `ion_idx`     — (1-based) index of the ion.
* `θ`           — rotation angle

The function returns the time at which the operation will be completed.
"""
function Rz(qdc           :: QCCDevControl,
            t             :: Time_t,
            ion_idx       :: Int,
            θ             :: Real      ) ::Time_t

end

####################################################################################################

"""
Function `Rxy` — single qubit XY-plane rotation

# Arguments
* `t::Time_t`   — time at which the operation commences.  Must be no earlier than
  the latest time given to previous function calls.
* `ion_idx`     — (1-based) index of the ion.
* `ϕ`           — rotation axis is cos(phi)cdot sigma_x + sin(phi)cdot sigma_y
* `θ`           — rotation angle

The function returns the time at which the operation will be completed.
"""
function Rxy(qdc           :: QCCDevControl,
             t             :: Time_t,
             ion_idx       :: Int,
             ϕ             :: Real,
             θ             :: Real      ) ::Time_t
    # Attention: Not all values of ϕ may work for the device

end

####################################################################################################

"""
Function `XX` — two qubit XX-rotation

# Arguments
* `t::Time_t`      — time at which the operation commences.  Must be no earlier than
  the latest time given to previous function calls.
* `ion𝑖_idx`, 𝑖=1,2 — the (1-based) indices of the two ions.  Must be in the same gate zone.
* `θ`              — rotation angle

The function returns the time at which the operation will be completed.
"""
function XX(qdc           :: QCCDevControl,
            t             :: Time_t,
            ion1_idx      :: Int,
            ion2_idx      :: Int,
            θ             :: Real      ) ::Time_t
    # Attention: May not work on all devices

end

####################################################################################################

"""
Function `ZZ` — two qubit ZZ-rotation

# Arguments
* `t::Time_t`      — time at which the operation commences.  Must be no earlier than
  the latest time given to previous function calls.
* `ion𝑖_idx`, 𝑖=1,2 — the (1-based) indices of the two ions.  Must be in the same gate zone.
* `θ`              — rotation angle

The function returns the time at which the operation will be completed.
"""
function ZZ(qdc           :: QCCDevControl,
            t             :: Time_t,
            ion1_idx      :: Int,
            ion2_idx      :: Int,
            θ             :: Real      ) ::Time_t
    # Attention: May not work on all devices

end

end #^ module QCCDevCtrl

# EOF

max_stars_repo_name: TartuQC/trapped-ion-qccd-simulator

max_issues_repo_name: TartuQC/trapped-ion-qccd-simulator

max_forks_repo_name: TartuQC/trapped-ion-qccd-simulator


text:
# "smoothness and data costs are set up one by one, individually.
#  grid neighborhood structure is assumed."
#
using GCoptimization

width = 10
height = 5
num_pixels = width * height
num_labels = 7
result = zeros(num_pixels)    # "stores result of optimization"

gco = GCoptimizationGridGraph(width, height, num_labels)

# "first set up data costs individually"
for i = 0:num_pixels-1, l = 0:num_labels-1
    if i < 25
        if l == 0
            setDataCost(gco, i, l, 0)
        else
            setDataCost(gco, i, l, 10)
        end
    else
        if l == 5
            setDataCost(gco, i, l, 0)
        else
            setDataCost(gco, i, l, 10)
        end
    end
end

# "next set up smoothness costs individually"
for l1 = 0:num_labels-1, l2 = 0:num_labels-1
    cost = (l1-l2)^2 <= 4 ? (l1-l2)^2 : 4
    setSmoothCost(gco, l1, l2, cost)
end

println("Before optimization energy is ", compute_energy(gco))
expansion(gco, 2)    # "run expansion for 2 iterations."
println("After optimization energy is ", compute_energy(gco))

for i = 0:num_pixels-1
    result[i+1] = whatLabel(gco, i)
end

showall(result)

ext: jl


text:
# plot residual oil saturation versus contact angle
using Plots
sor_ww=0.15
sor_ow=0.2
sor_mw=0.08
teta_ow=π
teta_ww=0.0
teta_mw=π/2+0.2 # just to shift it a bit
a=(sor_mw-0.5*(sor_ww-sor_ow)*cos(teta_mw)-0.5*(sor_ww+sor_ow))/(cos(teta_mw)^2-1)
b=0.5*(sor_ww-sor_ow)
c=0.5*(sor_ww+sor_ow)-a

teta=linspace(0, π, 100)
sor=a*cos(teta).^2+b*cos(teta)+c

plot(
plot(teta, sor, xlabel="contact angle [rad]", ylabel="Sor"),
plot(cos(teta), sor, xlabel="contact angle [rad]", ylabel="Sor")
)

max_stars_repo_name: simulkade/peteng

max_issues_repo_name: simulkade/peteng

max_forks_repo_name: simulkade/peteng


text:
using DashBootstrapComponents

input_groups = html_div([
    dbc_inputgroup(
        [dbc_inputgrouptext(dbc_radiobutton()), dbc_input()],
        className = "mb-3",
    ),
    dbc_inputgroup([dbc_inputgrouptext(dbc_checkbox()), dbc_input()]),
]);

max_issues_repo_name: sthagen/dash-bootstrap-components

max_forks_repo_name: sthagen/dash-bootstrap-components


text:
module SmithNormalForm

# import Base: show, getindex
using LinearAlgebra
using SparseArrays

export snf, smith

include("bezout.jl")
include("snf.jl")

struct Smith{P,Q<:AbstractMatrix} <: Factorization{P}
    S::Q
    Sinv::Q
    T::Q
    Tinv::Q
    SNF::AbstractVector{P}
    Smith{P,Q}(S::AbstractMatrix{P}, Sinv::AbstractMatrix{P},
               T::AbstractMatrix{P}, Tinv::AbstractMatrix{P},
               A::AbstractVector{P}) where {P,Q} = new(S, Sinv, T, Tinv, A)
end
Smith(S::AbstractMatrix{P}, T::AbstractMatrix{P}, A::AbstractVector{P}) where {P} =
    Smith{P,typeof(A)}(S, similar(S, 0, 0), T, similar(T, 0, 0), A)

function smith(X::AbstractMatrix{P}; inverse=true, debug=false) where {P}
    S, T, A, Sinv, Tinv = snf(X, inverse=inverse, debug=debug)
    return Smith{P, typeof(X)}(S, Sinv, T, Tinv, diag(A))
end

"""Retrive SNF matrix from the factorization"""
function LinearAlgebra.diagm(F::Smith{P,Q}) where {P,Q}
    base = issparse(F.SNF) ? spzeros(P, size(F.S,1), size(F.T,1)) : zeros(P, size(F.S,1), size(F.T,1))
    for i in 1:length(F.SNF)
        base[i,i] = F.SNF[i]
    end
    return base
end

function Base.show(io::IO, F::Smith)
    println(io, "Smith normal form:")
    show(io, diagm(F))
end

end # module


text:
abstract type AbstractChargeCloud end

include("PlatonicSolids/PlatonicSolids.jl")
include("ParticleTypes.jl")

"""
    struct NBodyChargeCloud{T, N, SH} <: AbstractChargeCloud

Struct which defines a charge cloud consisting of multiple point-like charge carriers,
initially distributed around a given center.

## Parametric Types
* `T`: Precision type.
* `N`: Number of shells.
* `SH`: Geometry of the shells.

## Fields
* `locations::Vector{CartesianPoint{T}}`: Positions of the charge carriers that are part of the charge cloud.
* `energies::Vector{T}`: Energies of the respective charge carriers, in the same order as `locations`.
* `shell_structure::SH`: Initial geometry of the charge carriers around the `center` point, relevant for plotting.
"""
struct NBodyChargeCloud{T <: SSDFloat, N, SH} <: AbstractChargeCloud
    locations::Vector{CartesianPoint{T}}
    energies::Vector{T} # in units of eV
    shell_structure::SH
end

"""
    NBodyChargeCloud(center::CartesianPoint{T}, energy::T, particle_type::Type{PT} = Gamma; kwargs...)

Returns an [`NBodyChargeCloud`](@ref) for a given `energy` deposition at a position that defines the `center` of the charge cloud,
given by a center charge surrounded by shells consisting of platonic solids.

## Arguments
* `center::CartesianPoint{T}`: Center position of the [`NBodyChargeCloud`](@ref).
* `energy::RealQuantity`: Deposited energy with units. If no units are given, the value is parsed in units of eV.
* `particle_type`: [`ParticleType`](@ref) of the particle that deposited the energy. Default is `Gamma`.

## Keywords
* `radius::T`: Estimate for the radius of the [`NBodyChargeCloud`](@ref). Default is determined from `particle_type` via `radius_guess`.
* `number_of_shells::Int`: Number of shells around the `center` point. Default is `2`.
* `shell_structure`: Geometry with which the charges are distributed in the shells. Default is `Dodecahedron`.

## Example
```julia
center = CartesianPoint{T}([0,0,0])
energy = 1460u"keV"
NBodyChargeCloud(center, energy, number_of_shells = 3, shell_structure = SolidStateDetectors.Icosahedron)
```

!!! note
    Using values with units for `energy` requires the package [Unitful.jl](https://github.com/PainterQubits/Unitful.jl).
"""
function NBodyChargeCloud(center::CartesianPoint{T}, energy::RealQuantity, particle_type::Type{PT} = Gamma;
        radius::T = radius_guess(T(to_internal_units(energy)), particle_type), number_of_shells::Int = 2, shell_structure = Dodecahedron
    )::NBodyChargeCloud{T, number_of_shells, typeof(shell_structure{T})} where {T, PT <: ParticleType}

    locations::Vector{CartesianPoint{T}} = CartesianPoint{T}[center]
    energies::Vector{T} = T[1]

    n_shell = 1
    while n_shell <= number_of_shells
        shell = shell_structure(center, n_shell * radius).locations
        locations = vcat(locations, shell)
        energies = vcat(energies, [exp(-n_shell^2 / 2) for i in 1:length(shell)])
        n_shell += 1
    end

    return NBodyChargeCloud{T, number_of_shells, typeof(shell_structure{T})}( locations, energies./sum(energies) * T(to_internal_units(energy)), shell_structure{T})
end

function create_regular_sphere(center::CartesianPoint{T}, N::Integer, R::T)::Vector{CartesianPoint{T}} where {T <: SSDFloat}

    # Source: https://www.cmu.edu/biolphys/deserno/pdf/sphere_equi.pdf

    # There might be fluctuations of the order 10 points or 1% (whatever is larger)
    locations = Vector{CartesianPoint{T}}(undef, N + max(10, round(Int,1.01*N)))
    Ncount = 0
    a = 4π / N
    Mθ = round(π/sqrt(a))
    dθ = π/Mθ
    for m in 0:Mθ-1
        θ = π * (m + 0.5)/Mθ
        Mφ = round(2π * sin(θ)*dθ/a)
        for n in 0:Mφ-1
            φ = 2π*n/Mφ
            Ncount += 1
            locations[Ncount] = center + CartesianVector{T}(R*cos(φ)*sin(θ), R*sin(φ)*sin(θ), R*cos(θ))
        end
    end

    return locations[1:Ncount]
end

"""
    NBodyChargeCloud(center::CartesianPoint{T}, energy::T, N::Integer, particle_type::Type{PT} = Gamma; kwargs...)

Returns an [`NBodyChargeCloud`](@ref) for a given `energy` deposition at a position that defines the `center` of the charge cloud,
given by a center charge surrounded by shells of approximately `N` point charges equally distributed on the surface of a sphere.
Find the algorithm to create the shells [here](https://www.cmu.edu/biolphys/deserno/pdf/sphere_equi.pdf).

## Arguments
* `center::CartesianPoint{T}`: Center position of the [`NBodyChargeCloud`](@ref).
* `energy::RealQuantity`: Deposited energy with units. If no units are given, the value is parsed in units of eV.
* `N::Integer`: Approximate number of charges in the [`NBodyChargeCloud`](@ref) (might vary by around 1%).
* `particle_type`: [`ParticleType`](@ref) of the particle that deposited the energy. Default is `Gamma`.

## Keywords
* `radius::T`: Estimate for the radius of the [`NBodyChargeCloud`](@ref). Default is determined from `particle_type` via `radius_guess`.
* `number_of_shells::Int`: Number of shells around the `center` point. Default is `2`.

## Example
```julia
center = CartesianPoint{T}([0,0,0])
energy = 1460u"keV"
NBodyChargeCloud(center, energy, 200, number_of_shells = 3)
```

!!! note
    Using values with units for `energy` requires the package [Unitful.jl](https://github.com/PainterQubits/Unitful.jl).
"""
function NBodyChargeCloud(center::CartesianPoint{T}, energy::RealQuantity, N::Integer, particle_type::Type{PT} = Gamma;
        radius::T = radius_guess(T(to_internal_units(energy)), particle_type), number_of_shells::Int = 2,
    )::NBodyChargeCloud{T, number_of_shells, NTuple{N > 1 ? number_of_shells : 0, Int}} where {T, PT <: ParticleType}

    locations::Vector{CartesianPoint{T}} = CartesianPoint{T}[center]
    energies::Vector{T} = T[1]
    n = Int[]

    if N > 1
        n_shell = 1
        expected_points_in_shells::Int = div(4 * (4^(number_of_shells) - 1), 3)
        while n_shell <= number_of_shells
            shell = create_regular_sphere(center, round(Int, N * (4^n_shell)/expected_points_in_shells), radius * n_shell)
            locations = vcat(locations, shell)
            push!(n, length(shell))
            energies = vcat(energies, [exp(-n_shell^2 / 2) for i in 1:length(shell)])
            n_shell += 1
        end
    end

    return NBodyChargeCloud{T, number_of_shells, NTuple{length(n), Int}}( locations, energies./sum(energies) * T(to_internal_units(energy)), Tuple(n))
end


text:
module Query

import ...Pkg
using ..Types

function requirements(reqs::Dict, fix::Dict, avail::Dict)
    for (p1,f1) in fix
        for p2 in keys(f1.requires)
            haskey(avail, p2) || haskey(fix, p2) || error("unknown package $p2 required by $p1")
        end
        satisfies(p1, f1.version, reqs) ||
            warn("$p1 is fixed at $(f1.version) conflicting with top-level requirement: $(reqs[p1])")
        for (p2,f2) in fix
            satisfies(p1, f1.version, f2.requires) ||
                warn("$p1 is fixed at $(f1.version) conflicting with requirement for $p2: $(f2.requires[p1])")
        end
    end
    reqs = deepcopy(reqs)
    for (p1,f1) in fix
        merge_requires!(reqs, f1.requires)
    end
    for (p,f) in fix
        delete!(reqs, p)
    end
    reqs
end

function dependencies(avail::Dict, fix::Dict = (ByteString=>Fixed)["julia"=>Fixed(VERSION)])
    avail = deepcopy(avail)
    conflicts = (ByteString=>Set{ByteString})[]
    for (fp,fx) in fix
        delete!(avail, fp)
        for (ap,av) in avail, (v,a) in copy(av)
            if satisfies(fp, fx.version, a.requires)
                delete!(a.requires, fp)
            else
                haskey(conflicts, ap) || (conflicts[ap] = Set{ByteString}())
                push!(conflicts[ap], fp)
                delete!(av, v)
            end
        end
    end
    again = true
    while again
        again = false
        deleted_pkgs = ByteString[]
        for (ap,av) in avail
            if isempty(av)
                delete!(avail, ap)
                push!(deleted_pkgs, ap)
                again = true
            end
        end
        for dp in deleted_pkgs
            for (ap,av) in avail, (v,a) in copy(av)
                if haskey(a.requires, dp)
                    haskey(conflicts, ap) || (conflicts[ap] = Set{ByteString}())
                    union!(conflicts[ap], conflicts[dp])
                    delete!(av, v)
                end
            end
        end
    end
    avail, conflicts
end

typealias PackageState Union(Nothing,VersionNumber)

function diff(have::Dict, want::Dict, avail::Dict, fixed::Dict)
    change = Array((ByteString,(PackageState,PackageState)),0)
    remove = Array((ByteString,(PackageState,PackageState)),0)

    for pkg in collect(union(keys(have),keys(want)))
        h, w = haskey(have,pkg), haskey(want,pkg)
        if h && w
            if have[pkg] != want[pkg]
                push!(change, (pkg,(have[pkg], want[pkg])))
            end
        elseif h
            push!(remove, (pkg,(have[pkg],nothing)))
        elseif w
            push!(change, (pkg,(nothing,want[pkg])))
        end
    end
    append!(sort!(change), sort!(remove))
end

function check_requirements(reqs::Requires, deps::Dict{ByteString,Dict{VersionNumber,Available}}, fix::Dict)
    for (p,vs) in reqs
        if !any([(vn in vs) for vn in keys(deps[p])])
            remaining_vs = VersionSet()
            err_msg = "fixed packages introduce conflicting requirements for $p: \n"
            available_list = sort(collect(keys(deps[p])))
            for (p1,f1) in fix
                f1r = f1.requires
                haskey(f1r, p) || continue
                err_msg *= "         $p1 requires versions $(f1r[p])"
                if !any([vn in f1r[p] for vn in available_list])
                    err_msg *= " [none of the available versions can satisfy this requirement]"
                end
                err_msg *= "\n"
                remaining_vs = intersect(remaining_vs, f1r[p])
            end
            if isempty(remaining_vs)
                err_msg *= "       the requirements are unsatisfiable because their intersection is empty"
            else
                err_msg *= "       available versions are $(join(available_list, ", ", " and "))"
            end
            error(err_msg)
        end
    end
end

# Reduce the number of versions by creating equivalence classes, and retaining
# only the highest version for each equivalence class.
# Two versions are equivalent if:
#   1) They appear together as dependecies of another package (i.e. for each
#      dependency relation, they are both required or both not required)
#   2) They have the same dependencies
# Also, if there are explicitly required packages, dicards all versions outside
# the allowed range (checking for impossible ranges while at it).
function prune_versions(reqs::Requires, deps::Dict{ByteString,Dict{VersionNumber,Available}})

    np = length(deps)

    # To each version in each package, we associate a BitVector.
    # It is going to hold a pattern such that all versions with
    # the same pattern are equivalent.
    vmask = (ByteString=>Dict{VersionNumber, BitVector})[]

    # Parse requirements and store allowed versions.
    allowed = (ByteString=>Dict{VersionNumber, Bool})[]
    for (p,vs) in reqs
        @assert !haskey(allowed, p)
        allowed[p] = (VersionNumber=>Bool)[]
        allowedp = allowed[p]
        for vn in keys(deps[p])
            allowedp[vn] = vn in vs
        end
        @assert !isempty(allowedp)
        @assert any(collect(values(allowedp)))
    end

    filtered_deps = (ByteString=>Dict{VersionNumber,Available})[]
    for (p,depsp) in deps
        filtered_deps[p] = (VersionNumber=>Available)[]
        allowedp = get(allowed, p, (VersionNumber=>Bool)[])
        fdepsp = filtered_deps[p]
        for (vn,a) in depsp
            if !isempty(allowedp) && !allowedp[vn]
                continue
            end
            fdepsp[vn] = a
        end
    end

    # For each package, we examine the dependencies of its versions
    # and put together those which are equal.
    # While we're at it, we also collect all dependencies into alldeps
    alldeps = Set{(ByteString,VersionSet)}()
    for (p, fdepsp) in filtered_deps

        # Extract unique dependencies lists (aka classes), thereby
        # assigning an index to each class.
        uniqdepssets = unique(values(fdepsp))

        # Store all dependencies seen so far for later use
        for a in uniqdepssets, r in a.requires
            push!(alldeps, r)
        end

        # If the package has just one version, it's uninteresting
        if length(deps[p]) == 1
            continue
        end

        # Grow the pattern by the number of classes
        luds = length(uniqdepssets)
        @assert !haskey(vmask, p)
        vmask[p] = (VersionNumber=>BitVector)[]
        vmaskp = vmask[p]
        for vn in keys(fdepsp)
            vmaskp[vn] = falses(luds)
        end
        for (vn,a) in fdepsp
            vmind = findfirst(uniqdepssets, a)
            @assert vmind >= 0
            vm = vmaskp[vn]
            vm[vmind] = true
        end
    end

    # Produce dependency patterns.
    for (p,vs) in alldeps
        # packages with just one version, or dependencies
        # which do not distiguish between versions, are not
        # interesting
        if length(deps[p]) == 1 || vs == VersionSet()
            continue
        end

        # Store the dependency info in the patterns
        @assert haskey(vmask, p)
        vmaskp = vmask[p]
        for (vn,vm) in vmaskp
            push!(vm, in(vn, vs))
        end
    end

    # At this point, the vmask patterns are computed. We divide them into
    # classes so that we can keep just one version for each class.
    pruned_vers = (ByteString=>Vector{VersionNumber})[]
    eq_classes = (ByteString=>Dict{VersionNumber,Vector{VersionNumber}})[]
    for (p, vmaskp) in vmask
        vmask0_uniq = unique(values(vmaskp))
        nc = length(vmask0_uniq)
        classes = [ VersionNumber[] for c0 = 1:nc ]
        for (vn,vm) in vmaskp
            c0 = findfirst(vmask0_uniq, vm)
            push!(classes[c0], vn)
        end
        map(sort!, classes)

        # For each nonempty class, we store only the highest version)
        pruned_vers[p] = VersionNumber[]
        prunedp = pruned_vers[p]
        eq_classes[p] = (VersionNumber=>Vector{VersionNumber})[]
        eqclassp = eq_classes[p]
        for cl in classes
            if !isempty(cl)
                vtop = maximum(cl)
                push!(prunedp, vtop)
                @assert !haskey(eqclassp, vtop)
                eqclassp[vtop] = cl
            end
        end
        sort!(prunedp)
    end
    # Put non-allowed versions into eq_classes
    for (p, allowedp) in allowed
        haskey(eq_classes, p) || continue
        eqclassp = eq_classes[p]
        for (vn, a) in allowedp
            a && continue
            eqclassp[vn] = [vn]
        end
    end
    # Put all remaining packages into eq_classes
    for (p, depsp) in deps
        haskey(eq_classes, p) && continue
        eq_classes[p] = (VersionNumber=>Vector{VersionNumber})[]
        eqclassp = eq_classes[p]
        for vn in keys(depsp)
            eqclassp[vn] = [vn]
        end
    end

    # Recompute deps. We could simplify them, but it's not worth it
    new_deps = (ByteString=>Dict{VersionNumber,Available})[]

    for (p,depsp) in deps
        @assert !haskey(new_deps, p)
        if !haskey(pruned_vers, p)
            new_deps[p] = depsp
            continue
        end
        new_deps[p] = (VersionNumber=>Available)[]
        pruned_versp = pruned_vers[p]
        for (vn,a) in depsp
            in(vn, pruned_versp) || continue
            new_deps[p][vn] = a
        end
    end

    #println("pruning stats:")
    #numvers = 0
    #numdeps = 0
    #for (p,d) in deps, (vn,a) in d
    #    numvers += 1
    #    for r in a.requires
    #        numdeps += 1
    #    end
    #end
    #numnewvers = 0
    #numnewdeps = 0
    #for (p,d) in new_deps, (vn,a) in d
    #    numnewvers += 1
    #    for r in a.requires
    #        numnewdeps += 1
    #    end
    #end
    #println("  before: vers=$numvers deps=$numdeps")
    #println("  after: vers=$numnewvers deps=$numnewdeps")
    #println()

    return new_deps, eq_classes
end
prune_versions(deps::Dict{ByteString,Dict{VersionNumber,Available}}) =
    prune_versions((ByteString=>VersionSet)[], deps)

# Build a subgraph incuding only the (direct and indirect) dependencies
# of a given package set
function dependencies_subset(deps::Dict{ByteString,Dict{VersionNumber,Available}}, pkgs::Set{ByteString})

    np = length(deps)

    staged = pkgs
    allpkgs = pkgs
    while !isempty(staged)
        staged_next = Set{ByteString}()
        for p in staged, a in values(deps[p]), rp in keys(a.requires)
            if !(rp in allpkgs)
                push!(staged_next, rp)
            end
        end
        allpkgs = union(allpkgs, staged_next)
        staged = staged_next
    end

    sub_deps = Dict{ByteString,Dict{VersionNumber,Available}}()
    for p in allpkgs
        haskey(sub_deps, p) || (sub_deps[p] = (VersionNumber=>Available)[])
        sub_depsp = sub_deps[p]
        for (vn, a) in deps[p]
            sub_depsp[vn] = a
        end
    end

    return sub_deps
end

function prune_dependencies(reqs::Requires, deps::Dict{ByteString,Dict{VersionNumber,Available}})
    deps = dependencies_subset(deps, Set{ByteString}(keys(reqs)))
    deps, _ = prune_versions(reqs, deps)

    return deps
end

end # module

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia


text:
function retrieve(c::Client, ::Type{AuditLog}, g::AbstractGuild; kwargs...)
    return get_guild_audit_log(c, g.id; kwargs...)
end


text:
using Underscores
using Utils

const Mappings = Dict{String,String}

function parse_file(f)
    @_ split_blocks(f) |> map(Dict(split(item, ':') for item=split(_)), __)
end

function load()
    open(joinpath(@__DIR__, "04_passports.txt"), "r") do f
        parse_file(f)
    end
end

function isvalid(d)
    isempty(setdiff(Set(["byr", "iyr", "eyr", "hgt", "hcl", "ecl", "pid"]), keys(d)))
end

colors = Set(["amb", "blu", "brn", "gry", "grn", "hzl", "oth"])

function fd(s, fromto)
    parsed = tryparse(UInt, s)
    if isnothing(parsed)
        false
    else
        parsed in fromto
    end
end

function isvalid_strict(d)
    try
        fd(d["byr"], 1920:2020) &&
            fd(d["iyr"], 2010:2020) &&
            fd(d["eyr"], 2020:2030) &&
            d["ecl"] in colors &&
            occursin(r"^#[0-9a-f]{6}$", d["hcl"]) &&
            occursin(r"^\d{9}$", d["pid"]) &&
            occursin(r"^(1([5-8]\d|9[0-3])cm|(59|6\d|7[0-6])in)$", d["hgt"])
    catch KeyError
        false
    end
end

function count_valid(isvalid, passports)
    @_ sum(isvalid(_) ? 1 : 0, passports)
end

passports = load()

println("Number of          valid passports: ", count_valid(isvalid, passports))
println("Number of striclty valid passports: ", count_valid(isvalid_strict, passports))

using Test

@testset "Passport Processing" begin
    @testset "example 1" begin
        example_input = """
            ecl:gry pid:860033327 eyr:2020 hcl:#fffffd
            byr:1937 iyr:2017 cid:147 hgt:183cm

            iyr:2013 ecl:amb cid:350 eyr:2023 pid:028048884
            hcl:#cfa07d byr:1929

            hcl:#ae17e1 iyr:2013
            eyr:2024
            ecl:brn pid:760753108 byr:1931
            hgt:179cm

            hcl:#cfa07d eyr:2025 pid:166559648
            iyr:2011 ecl:brn hgt:59in
            """
        expected = [
            Dict(
                "ecl" => "gry",
                "pid" => "860033327",
                "eyr" => "2020",
                "hcl" => "#fffffd",
                "byr" => "1937",
                "iyr" => "2017",
                "cid" => "147",
                "hgt" => "183cm",
            ),
            Dict(
                "iyr" => "2013",
                "ecl" => "amb",
                "cid" => "350",
                "eyr" => "2023",
                "pid" => "028048884",
                "hcl" => "#cfa07d",
                "byr" => "1929",
            ),
            Dict(
                "hcl" => "#ae17e1",
                "iyr" => "2013",
                "eyr" => "2024",
                "ecl" => "brn",
                "pid" => "760753108",
                "byr" => "1931",
                "hgt" => "179cm",
            ),
            Dict(
                "hcl" => "#cfa07d",
                "eyr" => "2025",
                "pid" => "166559648",
                "iyr" => "2011",
                "ecl" => "brn",
                "hgt" => "59in",
            ),
        ]

        example = parse_file(IOBuffer(example_input))

        @test example == expected

        @test isvalid(example[1])
        @test !isvalid(example[2])
        @test isvalid(example[3])
        @test !isvalid(example[4])
    end

    @testset "example 2 valid" begin
        passports = parse_file(IOBuffer("""
            pid:087499704 hgt:74in ecl:grn iyr:2012 eyr:2030 byr:1980
            hcl:#623a2f

            eyr:2029 ecl:blu cid:129 byr:1989
            iyr:2014 pid:896056539 hcl:#a97842 hgt:165cm

            hcl:#888785
            hgt:164cm byr:2001 iyr:2015 cid:88
            pid:545766238 ecl:hzl
            eyr:2022

            iyr:2010 hgt:158cm hcl:#b6652a ecl:blu byr:1944 eyr:2021 pid:093154719
            """))

        @test (@_ all(isvalid_strict(_), passports))
    end

    @testset "example 2 invalid" begin
        passports = parse_file(IOBuffer("""
            eyr:1972 cid:100
            hcl:#18171d ecl:amb hgt:170 pid:186cm iyr:2018 byr:1926

            iyr:2019
            hcl:#602927 eyr:1967 hgt:170cm
            ecl:grn pid:012533040 byr:1946

            hcl:dab227 iyr:2012
            ecl:brn hgt:182cm pid:021572410 eyr:2020 byr:1992 cid:277

            hgt:59cm ecl:zzz
            eyr:2038 hcl:74454a iyr:2023
            pid:3556412378 byr:2007
            """))

        @test !(@_ any(isvalid_strict(_), passports))
    end

    @testset "TDD" begin
        first = passports[1]
        @test length(first) == 6
    end

    @testset "results" begin
        @test count_valid(isvalid, passports) == 222
        @test count_valid(isvalid_strict, passports) == 140
    end
end

max_stars_repo_name: CmdQ/AoC


text:
include("loaders.jl")

export ir_load_brainweb_t1_256


text:
# DON'T USE!!! Just for visualisation
## PARAMS
mutable struct Param
	# don't modify
	range
	value
	unit::String
end
Param() = Param([],0,"")

mutable struct Params
	# variable number, add your own but keep the struct name as Params!
	# contains items of type Param
	#layers
    #Params(layers...) = new(layers)
	temp::Param
	CI::Param
	CQ::Param
	IQ::Param
	VI::Param
	VQtyp::Param
end

## AUX
mutable struct Aux
	# variable number, add your own but keep the struct name as Params!
	# contains items of type Param or Out
	aux1
end

## COMMENTS
mutable struct Comment
	# don't modify
	MyComment
	PicturePath
end

## OUTS
mutable struct Out
	# don't modify
	range
	value
	unit::String
end

mutable struct Outs
	# variable number, add your own but keep the struct name as Outs!
	# contains items of type Out
	DVline::Out
	DVload::Out
	VQmax::Out
	VQmin::Out
	VQmean::Out
	VQstd::Out
end

## WAVEFORMS
mutable struct Signal
	y
	dt
	t
end

mutable struct Signals
	vi::Signal
	vq::Signal
end

## TOP LEVEL
mutable struct MyData
	param::Params
	aux::Aux
	comment::Comment
	out::Outs
	wfm::Signals
end
MyData(x::Signals) = MyData(Params(Param()),)

# OR use a macro and be done with it. basic example

# First layer
# Params
temp    = Param([25,150], 25, "C")
CI    	= Param([1], 1, "uF")
CQ    	= Param([1], 1, "uF")
IQ    	= Param([1,10,100], 1, "mA")
VI      = Param([6,8,10], 8, "V")
VQtyp   = Param([3.3,5], 5, "V")
# Outs
DVline  = Out([], 0, "mV")
DVload  = Out([], 0, "mV")
VQmax   = Out([], 0, "V")
VQmin   = Out([], 0, "m")
VQmean  = Out([], 0, "V")
VQstd   = Out([], 0, "mV")
# Waveforms
vi      = Signal(ch1, dt, t)
vq      = Signal(ch2, dt, t)
# Second layer
param   = Params(temp, CI, CQ, IQ, VI, VQtyp)
aux     = Aux("")
comment = Comment("", "")
out     = Outs(DVline, DVload, VQmax, VQmin, VQmean, VQstd)
wfm     = Signals(vi, vq)
# Third layer
data    = MyData(param, aux, comment, out, wfm)


text:
using Distributions, StatsBase, Parameters, LinearAlgebra, DataStructures, PyPlot, Random

import Base: isless

include("parameters.jl")
include("rand_engine_model.jl")
include("event_controllers.jl")
include("views.jl")


text:
#=
P11 (*) Modified run-length encoding.
Modify the result of problem P10 in such a way that if an element has no duplicates it is simply copied into the result list. Only elements with duplicates are transferred as (N E) lists.

Example:
* (encode-modified '(a a a a b c c a a d e e e e))
((4 A) B (2 C) (2 A) D (4 E))
=#

# From P09
function _pack(coll, current_group, finished_groups)
    if isempty(coll)
        vcat(finished_groups, Vector{}[current_group])
    else
        if coll[1] == current_group[end]
            _pack(coll[2:end], vcat(current_group, coll[1]), finished_groups)
        else
            _pack(coll[2:end], [coll[1]], vcat(finished_groups, Vector{}[current_group]))
        end
    end
end

# From P09
function pack(coll)
    first, rest = coll[1], coll[2:end]
    _pack(rest, [first], [])
end

function encode_modified(coll)
    packed = pack(coll)
    map(group -> length(group) > 1 ? [length(group), group[1]] : group[1], packed)
end

encoded = encode_modified([:a, :a, :a, :a, :b, :c, :c, :a, :a, :d, :e, :e, :e, :e])
@assert encoded[1] == [4, :a]
@assert encoded[2] == :b
@assert encoded[3] == [2, :c]
@assert encoded[4] == [2, :a]
@assert encoded[5] == :d
@assert encoded[6] == [4, :e]

println("Tests passed: JL-11.jl")


text:
using PyPlot, transD_GP, Random, LinearAlgebra, Statistics

## model
Random.seed!(3)
nlayers, dz = 50, 5.
ρbg        = 1000
zfixed     = [-1e5,   (0:dz:(nlayers-1)*dz )...]
rho        = [1e12,   10 .^(log10(ρbg) .+ randn(nlayers))...]
ianom      = 5
rho[ianom] = 10.
δ          = 1e-4  # fd step size in log10

##  geometry and modeling parameters
rRx = 13.
zRx = -42.0
zTx = -40.0
freqlow = 1e-3
include("electronics_halt.jl")
calcjacobian = true

Fhm = transD_GP.AEM_VMD_HMD.HFieldDHT(
                      lowpassfcs = lowpassfcs,
                      nmax   = length(zfixed),
                      times  = HM_times,
                      ramp   = HM_ramp,
                      zTx    = zTx,
                      rRx    = rRx,
                      rTx    = rTx,
                      zRx    = zRx,
                      freqlow = freqlow,
                      calcjacobian = calcjacobian,
                      )
transD_GP.AEM_VMD_HMD.getfieldTD!(Fhm, zfixed, rho)

## plot HM
f, ax = plt.subplots(1,3, sharey=true, figsize=(8,6))
delz, delt = diff(zfixed[2:end]), diff(log10.(Fhm.times))
tlabel = [(log10.(Fhm.times) .- delt[1]/2)..., log10.(Fhm.times)[end]+delt[end]/2]
zlab = [zfixed[2:end]..., zfixed[end]+delz[end]]
ax[1].step(log10.(rho[2:end]), zfixed[2:end])
im2 = ax[2].pcolormesh(tlabel, zlab, Fhm.dBzdt_J[2:end,:], shading="auto")
cb2 = colorbar(im2, ax=ax[2])
im3 = ax[3].pcolormesh(tlabel, zlab, log10.(abs.(Fhm.dBzdt_J[2:end,:])))
cb2 = colorbar(im3, ax=ax[3])
ax[1].set_xlabel("log10 Ω-m")
ax[1].invert_xaxis()
ax[2].invert_yaxis()
ax[2].set_xlabel("log10 s")
ax[3].set_xlabel("log10 s")
ax[1].set_ylabel("Depth m")
ax[2].set_title("∂f/∂(log10σ)")
ax[3].set_title("∂f/∂(log10σ)")

## approximate partials with central differences and compare
J = copy(Fhm.dBzdt_J)

function getJ_TD!(Fin, zin, ρ, Δ)
    # sensitivities computed in log10σ
    # same as in AEM_VMD_HMD code
    Jout = similar(Fin.dBzdt_J)
    modelfd = copy(ρ)
    jflag = Fin.calcjacobian
    Fin.calcjacobian = false
    for i = 2:length(zin)
        modelfd[:] = ρ
        modelfd[i] = 10 .^-(-log10(ρ[i]) .+ Δ/2)
        transD_GP.AEM_VMD_HMD.getfieldTD!(Fin, zin, modelfd)
        HTD_z1 = copy(Fin.dBzdt)
        modelfd[i] = 10 .^-(-log10(ρ[i]) .- Δ/2)
        transD_GP.AEM_VMD_HMD.getfieldTD!(Fin, zin, modelfd)
        Jout[i,:] = (HTD_z1 - Fin.dBzdt) / Δ
    end
    Fin.calcjacobian = jflag
    Jout
end
J_fd = getJ_TD!(Fhm, zfixed, rho, δ)

## comparison plots
f, ax = plt.subplots(1,3, sharey=true, figsize=(10,6))
ax[1].step(log10.(rho[2:end]), zfixed[2:end])
im2 = ax[2].pcolormesh(tlabel, zlab, J[2:end,:])
cb2 = colorbar(im2, ax=ax[2])
im3 = ax[3].pcolormesh(tlabel, zlab, J_fd[2:end,:])
cb2 = colorbar(im3, ax=ax[3])
ax[1].set_xlabel("log10 Ω-m")
ax[1].invert_xaxis()
ax[2].invert_yaxis()
ax[2].set_xlabel("log10 s")
ax[3].set_xlabel("log10 s")
ax[1].set_ylabel("Depth m")
ax[2].set_title("∂f/∂(log10σ) analytic")
ax[3].set_title("∂f/∂(log10σ) FD")

## more detailed
figure()
nlayer = 40
semilogx(HM_times, J_fd[nlayer,:], "--k")
semilogx(HM_times, J[nlayer,:], "-")

## Benchmark timings
using BenchmarkTools
@info "analytic timings"
@btime transD_GP.AEM_VMD_HMD.getfieldTD!($Fhm, $zfixed, $rho)
@info "Cent Diff timing"
@btime J_fd = getJ_TD!($Fhm, $zfixed, $rho, $δ);


text:
# SPDX-License-Identifier: X11
# 2020-12-08
# Day 4, Part 2

const REQ_FIELDS = [
  "byr",
  "iyr",
  "eyr",
  "hgt",
  "hcl",
  "ecl",
  "pid"
]

const REG = r"^(cid|byr|iyr|eyr|hgt|hcl|ecl|pid):(.*)$"

const FIELD_REG = Dict(
  "yr"  => r"^\d{4}$",
  "hgt" => r"^(\d*)(cm|in)$",
  "hcl" => r"^#[0-9a-f]{6}$",
  "ecl" => r"^(amb|blu|brn|gry|grn|hzl|oth)$",
  "pid" => r"^\d{9}$"
)

validatebyr(v::AbstractString)::Bool = let x = parse(Int, v)
  occursin(FIELD_REG["yr"], v) && x >= 1920 && x <= 2002
end
validateiyr(v::AbstractString)::Bool = let x = parse(Int, v)
  occursin(FIELD_REG["yr"], v) && x >= 2010 && x <= 2020
end
validateeyr(v::AbstractString)::Bool = let x = parse(Int, v)
  occursin(FIELD_REG["yr"], v) && x >= 2020 && x <= 2030
end
function validatehgt(v::AbstractString)::Bool
  m = match(FIELD_REG["hgt"], v)
  m === nothing && return false
  x = parse(Int, m.captures[1])
  return m.captures[2] == "cm" ? x >= 150 && x <= 193 : x >= 59 && x <= 76
end
validatehcl(v::AbstractString)::Bool = occursin(FIELD_REG["hcl"], v)
validateecl(v::AbstractString)::Bool = occursin(FIELD_REG["ecl"], v)
validatepid(v::AbstractString)::Bool = occursin(FIELD_REG["pid"], v)

const FIELD_VALF = Dict(
  "byr" => validatebyr,
  "iyr" => validateiyr,
  "eyr" => validateeyr,
  "hgt" => validatehgt,
  "hcl" => validatehcl,
  "ecl" => validateecl,
  "pid" => validatepid
)

function main()
  input = split.(split(readchomp(stdin), "\n\n"; keepempty=false))

  (println ∘ length ∘ filter)(validate, input)
end

function validate(passport::Vector{SubString{String}})::Bool
  fields = (Dict ∘ zip)(REQ_FIELDS::Array{String},
                        (falses ∘ length)(REQ_FIELDS::Array{String}))

  for f ∈ passport
    m = match(REG::Regex, f)
    m.captures[1] == "cid" && continue

    !validatefield(m.captures[1], m.captures[2]) && return false
    fields[m.captures[1]] = true
  end

  return (all ∘ values)(fields)
end

validatefield(field::AbstractString, val::AbstractString)::Bool =
  FIELD_VALF[field](val)

main()

max_stars_repo_name: jthvai/competitive-programming

max_issues_repo_name: jthvai/competitive-programming

max_forks_repo_name: jthvai/competitive-programming


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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  # ramified

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

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

  # dyadic

  k = norm(G, l)

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

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

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

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

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

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

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

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

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

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

# TODO: caching
# TODO: better documentation

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

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

See [Kir16, Definition 8.3.1].
"""
function genus(L::HermLat, p)
  sym = _genus_symbol(L, p)
  G = genus(HermLat, nf(base_ring(L)), p, sym)
  # Just for debugging
  if isdyadic(G) && isramified(G)
    GG = _genus_symbol_kirschmer(L, p)
    for i in 1:length(G)
      @assert GG[i][4] == G.ni[i]
    end
    #
  end
  return G
end

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

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

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

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

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

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

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

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

  t = length(G1)

  p = prime(G1)

  # We now check the Jordan type

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

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

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

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

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

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

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

  # Dyadic ramified case

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

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

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

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

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

    d = dL1prod * dL2prod

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

  return true
end

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

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

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

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

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

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

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

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

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

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

  return true
end

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

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

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

# This could be made more efficient

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

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

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

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

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

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

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

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

#  Bad:= IsHermitian(L1) select BadPrimes(L1) join { p: p in MySupport( Discriminant(BaseRing(L1)) ) | Minimum(p) eq 2 }
#                          else BadPrimes(L1 : Even);
# return forall{ p: p in Bad | IsLocallyIsometric(L1, L2, p) };

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

base_field(G::GenusHerm) = G.E

rank(G::GenusHerm) = G.rank

function representative(G::GenusHerm)
  P = _non_norm_primes(G.LGS)
  E = base_field(G)
  V = hermitian_space(E, _hermitian_form_with_invariants(base_field(G), rank(G), P, G.signatures))
  M = maximal_integral_lattice(V)
  for g in G.LGS
    p = prime(g)
    L = representative(g)
    M = find_lattice(M, L, p)
  end
  return M
end

#  V = m.HermitianFormWithInvariants(E, self._rank, Pm, self._signatures)
#  if rational:
#      return V.ChangeRing(Em).sage()
#  M = V.MaximalIntegralHermitianLattice()
#  #M = 27*M
#  #M = M.MaximalIntegralLattice()
#  for sym in self._local_symbols:
#      p = sym._prime
#      g = sym.gram_matrix()
#      g = m(g).ChangeRing(E)
#      L = m.HermitianLattice(g)
#      pm = self._ideal_to_magma(Km, p)
#      if Km.sage() == QQ:
#          pm = gcd([ZZ(a) for a in p.gens()])
#      M = m.FindLattice(M, L, pm)
#  return M
#

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

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

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

  K = number_field(order(p))

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

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

  if !is_ramified
    # I add the 0 to make the compiler happy
    return [ genus(HermLat, E,p, Tuple{Int, Int, Int, Int}[(b..., 1, 0) for b in g]) for g in scales_rks]
  end

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

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

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

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

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

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

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

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

  local_symbols = []

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

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

  return res
end

################################################################################
#
#  First attempt, which mirrors Markus' Magma code
#
################################################################################

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

prime(G::LocalGenusSymbol) = G.P

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

data(G::LocalGenusSymbol) = G.data

base_field(G::LocalGenusSymbol) = G.E

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

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

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

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

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

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

#@doc Markdown.doc"""
#    change_uniformizer(G::LocalGenusSymbol, a::NfOrdElem) -> LocalGenus
#
#Returns an equivalent? genus symbol with uniformizer $a$.
#"""
function change_uniformizer(G::LocalGenusSymbol{QuadLat}, unif::NfOrdElem)
  if unif == uniformizer(G)
    return G
  end
  P = prime(G)
  @assert isodd(minimum(P))
  @assert valuation(unif, P) == 1
  _, mF = ResidueField(order(P), P)
  mFF = extend(mF, nf(order(P)))
  b,_ = issquare(mFF(unif//uniformizer(G)))
  if b
    return LocalGenusSymbol{QuadLat}(P, G.data, unif, G.iseven, G.E, nothing, nothing)
  else
    e = G.data[1]
    return LocalGenusSymbol{QuadLat}(P, (e[1], e[2], isodd(e[1] * e[2]) ? -e[3] : e[3]), unif, G.iseven, G.E, nothing, nothing)
  end
end

function Base.:(==)(G1::LocalGenusSymbol{QuadLat}, G2::LocalGenusSymbol{QuadLat})
  if uniformizer(G1) != uniformizer(G2)
    error("Uniformizers of the genus symbols must be equal")
  end
  return data(G1) == data(G2)
end

function uniformizer(G::LocalGenusSymbol{HermLat})
  E = base_field(G)
  K = base_field(E)
  if isramified(G)
    lP = prime_decomposition(maximal_order(E), prime(G))
    @assert length(lP) == 1 && lP[1][2] == 2
    Q = lP[1][1]
    pi = uniformizer(Q)
    A = automorphisms(E)
    uni = A[1](elem_in_nf(pi)) * A[2](elem_in_nf(pi))
    @assert iszero(coeff(uni, 1))
    @assert islocal_norm(E, coeff(uni , 0), prime(G))
    return coeff(uni, 0)
  else
    return uniformizer(prime(G))
  end
end

@doc Markdown.doc"""
    det(G::LocalGenusSymbol) -> NumFieldElem

Return the norm class of the determinant of $G$.
"""
function det(G::LocalGenusSymbol{HermLat})
  z = G.data
  d = prod(b[3] for b in z)
  v = sum(b[1] * b[2] for b in z)
  if isramified(maximal_order(G.E), G.P)
    v = div(v, 2)
  end
  if d == 1
    u = base_field(G.E)(1)
  else
    @assert isramified(G)
    u = _non_norm_rep(G)
  end
  return u * uniformizer(G)^v
end

function rank(G::LocalGenusSymbol{HermLat})
  z = G.data
  return sum(b[2] for b in z)::Int
end


text:
# Exercise 4-1
using ThinkJulia

turtle = Turtle()
@svg begin
    forward(turtle, 100)
    turn(turtle, -90)
    forward(turtle, 100)
    turn(turtle, -90)
    forward(turtle, 100)
    turn(turtle, -90)
    forward(turtle, 100)
end

# Exercise 4-2
function square(t)
    @svg begin
        for i in 1:4
            forward(turtle, 100)
            turn(turtle, -90)
        end
    end
end

turtle = Turtle()
square(turtle)

# Exercise 4-3
turtle = Turtle()
@svg begin
    square(turtle)
end

# Exercise 4-4
function square(t, len)
    for i in 1:4
        forward(turtle, len)
        turn(turtle, -90)
    end
end

turtle = Turtle()
@svg begin
    square(turtle, 200)
end

# Exercise 4-5
function polygon(t, n, len)
    angle = 360 / n
    for i in 1:n
        forward(t, len)
        turn(t, -angle)
    end
end

turtle = Turtle()
@svg begin
    polygon(turtle, 100, 8)
end

# Exercise 4-6
function circle(t, r)
    polygon(t, r, 30)
end

turtle = Turtle()
@svg begin
    circle(turtle, 10)
end

# Exercise 4-7
function arc(t, r, angle)
    circ = 2 * pi * r
    frac = angle / 360
    arclength = circ * frac
    n = 30
    len = arclength / n
    for i in 1:n
        forward(turtle, len)
        turn(turtle, angle/n)
    end
end

turtle = Turtle()
@svg begin
    arc(turtle, 100, 235)
end


text:
@enum H5PhPType begin
  VELOCITY = 1
  DENSITY = 2
  QUALITYFACTOR = 3
  EPSILON = 4
  DELTA = 5
  THETA = 6
  PHI = 7
end

function rotate_Nx2_array(arr, ϕ)
  return arr * [cos(ϕ) -sin(ϕ); sin(ϕ) cos(ϕ)]
end

function rotate_2xN_array(arr, ϕ)
  return [cos(ϕ) -sin(ϕ); sin(ϕ) cos(ϕ)] * arr
end

function h5geo2judiTraceHeaderName(h5geoHdrName::String)
  h5geo = pyimport("h5geopy._h5geo")
  _, h5geoShortHeaderNames = h5geo.getTraceHeaderNames()

  ind = findall(x->x==h5geoHdrName, h5geoShortHeaderNames)
  if length(ind) < 1
    return
  end

  judiHeaderNames = judiTraceHeaders()
  if ind[1] > length(judiHeaderNames)
    return
  end

  return judiHeaderNames[ind[1]]
end

function judi2h5geoTraceHeaderName(judiHdrName::String)
  h5geo = pyimport("h5geopy._h5geo")
  judiHeaderNames = judiTraceHeaders()

  ind = findall(x->x == judiHdrName, judiHeaderNames)
  if length(ind) < 1
    return
  end

  _, h5geoShortHeaderNames = h5geo.getTraceHeaderNames()
  if ind[1] > length(h5geoShortHeaderNames)
    return
  end

  return h5geoShortHeaderNames[ind[1]]
end

function judiTraceHeaders()
  return [
    "TraceNumWithinLine"
    "TraceNumWithinFile"
    "FieldRecord"
    "TraceNumber"
    "EnergySourcePoint"
    "CDP"
    "CDPTrace"
    "TraceIDCode"
    "NSummedTraces"
    "NStackedTraces"
    "DataUse"
    "Offset"
    "RecGroupElevation"
    "SourceSurfaceElevation"
    "SourceDepth"
    "RecDatumElevation"
    "SourceDatumElevation"
    "SourceWaterDepth"
    "GroupWaterDepth"
    "ElevationScalar"
    "RecSourceScalar"
    "SourceX"
    "SourceY"
    "GroupX"
    "GroupY"
    "CoordUnits"
    "WeatheringVelocity"
    "SubWeatheringVelocity"
    "UpholeTimeSource"
    "UpholeTimeGroup"
    "StaticCorrectionSource"
    "StaticCorrectionGroup"
    "TotalStaticApplied"
    "LagTimeA"
    "LagTimeB"
    "DelayRecordingTime"
    "MuteTimeStart"
    "MuteTimeEnd"
    "ns"
    "dt"
    "GainType"
    "InstrumentGainConstant"
    "InstrumntInitialGain"
    "Correlated"
    "SweepFrequencyStart"
    "SweepFrequencyEnd"
    "SweepLength"
    "SweepType"
    "SweepTraceTaperLengthStart"
    "SweepTraceTaperLengthEnd"
    "TaperType"
    "AliasFilterFrequency"
    "AliasFilterSlope"
    "NotchFilterFrequency"
    "NotchFilterSlope"
    "LowCutFrequency"
    "HighCutFrequency"
    "LowCutSlope"
    "HighCutSlope"
    "Year"
    "DayOfYear"
    "HourOfDay"
    "MinuteOfHour"
    "SecondOfMinute"
    "TimeCode"
    "TraceWeightingFactor"
    "GeophoneGroupNumberRoll"
    "GeophoneGroupNumberTraceStart"
    "GeophoneGroupNumberTraceEnd"
    "GapSize"
    "OverTravel"
    "CDPX"
    "CDPY"
    "Inline3D"
    "Crossline3D"
    "ShotPoint"
    "ShotPointScalar"
    "TraceValueMeasurmentUnit"
    "TransductionConstnatMantissa"
    "TransductionConstantPower"
    "TransductionUnit"
    "TraceIdentifier"
    "ScalarTraceHeader"
    "SourceType"
    "SourceEnergyDirectionMantissa"
    "SourceEnergyDirectionExponent"
    "SourceMeasurmentMantissa"
    "SourceMeasurementExponent"
    "SourceMeasurmentUnit"
    "Unassigned1"
    "Unassigned2"
    ]
end

# function h5geoTraceHeaders()
#   return ["SEQWL", "SEQWR", "FFID", "TRCFLD", "SP", "CDP", "TRCNUM", "TRCID", "NVST", "NHST", "DU", "DSREG", "RGE", "SES", "SDBS", "DERG", "DES", "WDS", "WGD", "SAED", "SAC", "SRCX", "SRCY", "GRPX", "GRPY", "UNITS", "WVEL", "SVEL", "UTSRC", "UTGRP", "SECSCOR", "GRPSCOR", "TSA", "LAGTA", "LAGTB", "DELRECT", "MTSTART", "MTEND", "NSMP", "SI", "GTFI", "IG", "IGC", "CORREL", "SFSTART", "SFEND", "SLEN", "STYP", "SSTRLS", "SSTLE", "TTYP", "AFF", "AFS", "NFF", "NFS", "LOCF", "HOCF", "LOCS", "HICS", "YEAR", "DAY", "HOUR", "MINUTE", "SCE", "TMBS", "TWF", "GGNSW", "GGN1ST", "GGNLST", "GAPSZ", "OAWT", "CDP_X", "CDP_Y", "INLINE", "XLINE", "SPN", "SPS", "TVMU"]
# end


text:
## NotANumber

struct NotANumber <: Number end

"""
   Infinity()

represents infinite cardinality. Note that `Infinity <: Integer` to support
being treated as an index.
"""
struct Infinity <: Integer end

const ∞ = Infinity()

show(io::IO, y::Infinity) = print(io, "∞")

promote_rule(::Type{Infinity}, ::Type{II}) where II<:Integer = Union{Infinity,II}

convert(::Type{Float64}, ::Infinity) = Inf64
convert(::Type{Float32}, ::Infinity) = Inf32
convert(::Type{Float16}, ::Infinity) = Inf16
convert(::Type{AF}, ::Infinity) where AF<:AbstractFloat = convert(AF, Inf)

sign(y::Infinity) = 1
angle(x::Infinity) = 0

==(x::Infinity, y::Infinity) = true

isinf(::Infinity) = true
isfinite(::Infinity) = false

==(x::Infinity, y::Number) = isinf(y) && angle(y) == angle(x)
==(y::Number, x::Infinity) = x == y

isless(x::Infinity, y::Infinity) = false
isless(x::Real, y::Infinity) = isfinite(x) || sign(y) == -1
isless(x::AbstractFloat, y::Infinity) = isless(x, convert(typeof(x), y))
isless(x::Infinity, y::AbstractFloat) = false
isless(x::Infinity, y::Real) = false

+(::Infinity, ::Infinity) = ∞
+(::Number, y::Infinity) = ∞
+(::Infinity, ::Number) = ∞
-(::Infinity, ::Number) = ∞
-(x::Number, ::Infinity) = x + (-∞)

+(::Integer, y::Infinity) = ∞
+(::Infinity, ::Integer) = ∞
-(::Infinity, ::Integer) = ∞
-(x::Integer, ::Infinity) = x + (-∞)
+(::Complex, y::Infinity) = ∞
+(::Infinity, ::Complex) = ∞
-(::Infinity, ::Complex) = ∞
-(x::Complex, ::Infinity) = x + (-∞)

# ⊻ is xor
*(::Infinity) = ∞
*(::Infinity, ::Infinity) = ∞

for OP in (:fld,:cld,:div)
  @eval $OP(::Infinity, ::Number) = ∞
end

min(::Infinity, ::Infinity) = ∞
max(::Infinity, ::Infinity) = ∞
min(x::Real, ::Infinity) = x
max(::Real, ::Infinity) = ∞
min(::Infinity, x::Real) = x
max(::Infinity, ::Real) = ∞

≤(::Infinity, ::Infinity) = true
<(::Infinity, ::Infinity) = false
≥(::Infinity, ::Infinity) = true
>(::Infinity, ::Infinity) = false

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

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

# angle is π*a where a is (false==0) and (true==1)
struct OrientedInfinity{T<:Real} <: Number
    angle::T
end

OrientedInfinity{T}() where T = OrientedInfinity(zero(T))
OrientedInfinity() = OrientedInfinity{Bool}()
OrientedInfinity{T}(::Infinity) where T<:Real = OrientedInfinity{T}()
OrientedInfinity(::Infinity) = OrientedInfinity()
-(::Infinity) = OrientedInfinity(true)
+(::Infinity) = ∞

isinf(::OrientedInfinity) = true
isfinite(::OrientedInfinity) = false

promote_rule(::Type{Infinity}, ::Type{OrientedInfinity{T}}) where T = OrientedInfinity{T}
convert(::Type{OrientedInfinity{T}}, ::Infinity) where T = OrientedInfinity{T}()
convert(::Type{OrientedInfinity}, ::Infinity) = OrientedInfinity()

sign(y::OrientedInfinity{<:Integer}) = mod(y.angle,2) == 0 ? 1 : -1
angle(x::OrientedInfinity) = π*x.angle
mod(::OrientedInfinity{<:Integer}, ::Integer) = NotANumber()

function show(io::IO, y::OrientedInfinity{B}) where B<:Integer
    if sign(y) == 1
        print(io, "+∞")
    else
        print(io, "-∞")
    end
end

show(io::IO, x::OrientedInfinity) = print(io, "$(exp(im*π*x.angle))∞")

==(x::OrientedInfinity, y::Infinity) = x.angle == 0
==(y::Infinity, x::OrientedInfinity) = x.angle == 0
==(x::OrientedInfinity, y::OrientedInfinity) = x.angle == y.angle

==(x::OrientedInfinity, y::Number) = isinf(y) && angle(y) == angle(x)
==(y::Number, x::OrientedInfinity) = x == y

isless(x::OrientedInfinity{Bool}, y::OrientedInfinity{Bool}) = x.angle && !y.angle
isless(x::Number, y::OrientedInfinity{Bool}) = !y.angle && x ≠ ∞
isless(x::OrientedInfinity{Bool}, y::Number) = x.angle && y ≠ -∞

-(y::OrientedInfinity{B}) where B<:Integer = sign(y) == 1 ? OrientedInfinity(one(B)) : OrientedInfinity(zero(B))

function +(x::OrientedInfinity, y::OrientedInfinity)
    x == y || throw(ArgumentError("Angles must be the same to add ∞"))
    promote_type(typeof(x),typeof(y))(x.angle)
end

+(x::OrientedInfinity, y::Infinity) = x+OrientedInfinity(y)
+(x::Infinity, y::OrientedInfinity) = OrientedInfinity(x)+y
+(::Number, y::OrientedInfinity) = y
+(y::OrientedInfinity, ::Number) = y
-(y::OrientedInfinity, ::Number) = y
-(::Number, y::OrientedInfinity) = -y

# ⊻ is xor
*(a::OrientedInfinity{Bool}, b::OrientedInfinity{Bool}) = OrientedInfinity(a.angle ⊻ b.angle)
*(a::OrientedInfinity, b::OrientedInfinity) = OrientedInfinity(a.angle + b.angle)
*(a::Infinity, b::OrientedInfinity) = OrientedInfinity(a)*b
*(a::OrientedInfinity, b::Infinity) = a*OrientedInfinity(b)

*(a::Real, y::OrientedInfinity) = a > 0 ? y : (-y)
*(y::OrientedInfinity, a::Real) = a*y

*(a::Number, y::OrientedInfinity) = OrientedInfinity(y.angle+angle(a)/π)
*(y::OrientedInfinity, a::Number) = a*y

*(a::Number,y::Infinity) = a*OrientedInfinity(y)
*(y::Infinity, a::Number) = OrientedInfinity(y)*a

*(a::Integer,y::Infinity) = a*OrientedInfinity(y)
*(a::Complex,y::Infinity) = a*OrientedInfinity(y)
*(y::Infinity, a::Complex) = OrientedInfinity(y)*a
*(y::Infinity, a::Integer) = OrientedInfinity(y)*a

for OP in (:fld,:cld,:div)
  @eval $OP(y::OrientedInfinity, a::Number) = y*(1/sign(a))
end

min(x::OrientedInfinity{B}, y::OrientedInfinity{B}) where B<:Integer = sign(x) == -1 ? x : y
max(x::OrientedInfinity{B}, ::OrientedInfinity{B}) where B<:Integer = sign(x) == 1 ? x : y
min(x::Real, y::OrientedInfinity{B}) where B<:Integer = sign(y) == 1 ? x : y
min(x::OrientedInfinity{B}, y::Real) where B<:Integer = min(y,x)
max(x::Real, y::OrientedInfinity{B}) where B<:Integer = sign(y) == 1 ? y : x
max(x::OrientedInfinity{B}, y::Real) where B<:Integer = max(y,x)

for OP in (:<,:≤)
    @eval begin
        $OP(x::Real, y::OrientedInfinity{B}) where B<:Integer = sign(y) ==  1
        $OP(y::OrientedInfinity{B}, x::Real) where B<:Integer = sign(y) == -1
    end
end

for OP in (:>, :≥)
    @eval begin
        $OP(x::Real, y::OrientedInfinity{B}) where B<:Integer = sign(y) == -1
        $OP(y::OrientedInfinity{B}, x::Real) where B<:Integer = sign(y) == 1
    end
end


text:
"""
appropriate_dims(n1, n2, N)
# Description
Create an array of size n1, with the value N at the location n2, and ones elswhere

# Arguments
- `n1`: int | size of the array
- `n2`: int | location of modification
- `N` : int | value at the modification index n2

# Return
A tuple with all 1's except at location n2, where it is N
"""
function appropriate_dims(n1, n2, N)
    y = ones(Int, n1)
    y[n2] = N
    return Tuple(y)
end

max_stars_repo_name: upiterbarg/QuasiGeostrophy

max_issues_repo_name: sandreza/QuasiGeostrophy

max_forks_repo_name: upiterbarg/QuasiGeostrophy


text:
mutable struct RandomPositionsXXZWithDegeneracyLifted <: Model
    posdatadescriptor
    interaction
    δs
    basis
end

# small constructor - useful for loading
RandomPositionsXXZWithDegeneracyLifted(posdatadescriptor, interaction) = RandomPositionsXXZWithDegeneracyLifted(posdatadescriptor, interaction, nothing, nothing)

function Base.getproperty(xxz::RandomPositionsXXZWithDegeneracyLifted, p::Symbol)
    if hasfield(typeof(xxz), p)
        return getfield(xxz, p)
    else
        return getproperty(getfield(xxz, :posdatadescriptor), p)
    end
end

model_name(::RandomPositionsXXZWithDegeneracyLifted) = "Jxxz_with_degeneracy_lifted"
model_fileprefix(xxz::RandomPositionsXXZWithDegeneracyLifted) = @sprintf("%s_%id_alpha_%.1f_N_%02i_lifted", xxz.geometry, xxz.dimension, xxz.interaction.α, xxz.system_size)

struct PreparedJXXZWithDegeneracyLifted{JT,SB}
    J::JT
    δs::Vector{Float64}
    basis::SB
end

# shot, field, rho
parameter_dims(model::PreparedJXXZWithDegeneracyLifted) = (size(model.J)[3], length(model.δs), size(model.J)[4])

function initialize_model!(model::RandomPositionsXXZWithDegeneracyLifted, array_initializer)
    posdata = load(model.posdatadescriptor)
    shots = posdata.shots
    N = posdata.system_size
    ρs = posdata.ρs
    J = fill!(array_initializer(Float64,N,N,shots, length(ρs)), 0)
    # see Jxxz_x.jl
    _compute_interactions!(J, posdata, model.interaction, nothing)

    return PreparedJXXZWithDegeneracyLifted(J, model.δs, model.basis)
end

_flat_to_indices(index, nshots) = fldmod1(index, nshots)

function split_workload(model::PreparedJXXZWithDegeneracyLifted, parts)
    nshots = size(model.J, 3)
    nρs = size(model.J, 4)
    return _workload_chunks(nshots*nρs, parts)
end

function do_parameters(diag_callback, model::PreparedJXXZWithDegeneracyLifted, parameter_chunk)
    N = size(model.J, 1)
    nshots = size(model.J, 3)
    for index in parameter_chunk
        logmsg("Starting index $(index) of $(parameter_chunk)")
        ρindex, shot = _flat_to_indices(index, nshots)

        J = @view model.J[:,:, shot, ρindex]
        H_int = symmetrize_operator(xxzmodel(J, -0.73), model.basis)

        for (k, δ) in enumerate(model.δs)
            logmsg("Doing #rho=$ρindex #shot=$shot #delta=$k")
            parameterI = CartesianIndex(shot, k, ρindex)
            z_fields = dropdims(sum(J; dims=1); dims=1)
            diag_callback(parameterI, H_int + δ * symmetrize_operator(z_field(z_fields), model.basis))
        end
    end
end


text:
module Threadicare

if VERSION < v"1.3"
    @error "Threadicare requires Julia v1.3 or newer."
end

using KernelAbstractions,
    BenchmarkTools,
    LinearAlgebra

import Base:
    +, -, *, /,
    size, length, eltype,
    iterate, similar, show,
    getindex, lastindex, setindex!,
    push!

include("test.jl")

export test_args

end # module

max_stars_repo_name: sandreza/Threadicare