text:
module RingLists
using Random
import Random.shuffle
export RingList, insertafter!, insertbefore!, previous, next, shuffle

import Base: ==, length, getindex, keys, haskey, insert!, eltype
import Base: Vector, show, hash, reverse, first, delete!, Set, collect, copy

struct RingList{T}
    data::Dict{T,T}
    revdata::Dict{T,T}
end

"""
`RingList{T}()` creates a new, empty `RingList` holding elements of type `T`.
Example: `RingList{Int}()`.

`RingList(list)` creates a new `RingList` from the elements in the
one-dimensional array `list`. Example: `RingList([1,2,3])`.

`RingList(x...)` creates a new `RingList` from the arguments.
Example: `RingList(1,2,3)`.
"""
function RingList(T::Type = Any)
    return RingList(Dict{T,T}(), Dict{T,T}())
end

function _rev(d::Dict{T,T}) where {T}
    rd = Dict{T,T}()
    for k in keys(d)
        v = d[k]
        rd[v] = k
    end
    return rd
end

function RingList(vals::Vector{T}) where {T}
    d = Dict{T,T}()
    n = length(vals)
    if n != length(unique(vals))  # check there are no repeats
        error("List of values may not have a repeat")
    end
    if n == 0
        return RingList(d, d)
    end
    if n == 1
        a = vals[1]
        d[a] = a
        return RingList(d, d)
    end
    for i = 1:n-1
        a = vals[i]
        b = vals[i+1]
        d[a] = b
    end
    a = vals[end]
    b = vals[1]
    d[a] = b
    return RingList(d, _rev(d))
end

function RingList(x...)
    return RingList(collect(x))
end

function RingList{T}() where {T}
    d = Dict{T,T}()
    return RingList(d, d)
end

copy(a::RingList) = RingList(collect(a))

==(a::RingList, b::RingList) = a.data == b.data
length(a::RingList) = length(a.data)
keys(a::RingList) = keys(a.data)

"""
`next(a::RingList, x)` returns the element that follows `x` in `a`,
or throws an error if no such element exists. Abbreviation: `a[x]`.
    See also `previous`.
"""
next(a::RingList, x) = a.data[x]
getindex(a::RingList, x) = next(a, x)

"""
`previous(a::RingList,x)` returns the element `y` so that `a[y]==x`, or throws an
error if no such element exists. Abbreviation: `a(x)`. See also `next`.
"""
previous(a::RingList, x) = a.revdata[x]
(a::RingList)(x) = previous(a, x)

haskey(a::RingList, x) = haskey(a.data, x)
eltype(a::RingList{T}) where {T} = T

"""
`first(a::RingList, true_first::Bool = true)` returns
an element of `a`. If `true_first` is `true`, then try to
return the smallest element held in `a`.
"""
function first(a::RingList, true_first::Bool = true)
    if true_first
        try
            x = minimum(keys(a.data))
            return x
        catch
        end
    end
    return first(a.data)[1]
end

"""
`insert!(a::RingList,x)` adds the element `x` into `a` after the
first element of `a` (which might not be the smallest element in `a`).
Ergo, no guarantee where it will go. See also
`insertafter!`
"""
function insert!(a::RingList{T}, x::T) where {T}
    if length(a) == 0
        a.data[x] = x
        a.revdata[x] = x
        return nothing
    end
    if haskey(a, x)
        error("$x already in this RingList")
    end
    y = first(a)  # get the other elements
    insertafter!(a, x, y)
end

"""
`insertafter!(a::RingList,x,y)` inserts `x` into `a` after `y`.
"""
function insertafter!(a::RingList, x, y)
    if haskey(a, x)
        error("Element $x already in this RingList")
    end
    if !haskey(a, y)
        error("Element $y not in this RingList, cannot insert after")
    end
    # who is currently after y?
    z = a[y]
    # we have y --> z
    # and change to y --> x --> z

    a.data[y] = x
    a.data[x] = z

    a.revdata[x] = y
    a.revdata[z] = x

    nothing
end

"""
`insertbefore!(a::RingList,x,y)` inserts `x` into `a` before `y`.
"""
function insertbefore!(a::RingList, x, y)
    if haskey(a,x)
        error("Element $x already in this RingList")
    end
    if !haskey(a, y)
        error("Element $y not in this RingList, cannot insert before")
    end
    # who is currently before y?
    z = a(y)
    # we have z-->y
    # we want to change to z-->x-->y
    a.data[z] = x
    a.data[x] = y

    a.revdata[y] = x
    a.revdata[x] = z

    nothing
end

"""
`delete!(a::RingList,x)` removes `x` from `a`.
"""
function delete!(a, x)
    if !haskey(a, x)
        error("$x not in the RingList")
    end
    if length(a) == 1
        delete!(a.data, x)
        return
    end

    # prev --> x --> next
    next = a[x]
    prev = previous(a, x)
    delete!(a.data, x)
    delete!(a.revdata, x)
    a.data[prev] = next
    a.revdata[next] = prev
    return
end

"""
`Vector(a::RingList)` converts `a` into an ordinary list (`Vector`).
Also: `collect(a)`.
"""
Vector(a::RingList) = [x for x in a]
collect(a::RingList) = Vector(a)

"""
`reverse(a::RingList)` returns a new `RingList` containing the
same elements as `a` but in reverse order.
"""
reverse(a::RingList{T}) where {T} = RingList{T}(copy(a.revdata), copy(a.data))

"""
`Set(a::RingList)` returns the elements of `a` as a `Set`.
"""
Set(a::RingList{T}) where {T} = Set{T}(keys(a))

hash(a::RingList) = hash(a.data)
hash(a::RingList, h::UInt64) = hash(a.data, h)

"""
`shuffle(a::RingList)` returns a new `RingList` with the same
elements, but in a randomized order.
"""
function shuffle(a::RingList)
    n = length(a)
    p = randperm(n)
    d = collect(a)
    return RingList(d[p])
end

include("iter.jl")
include("show.jl")

end # module


text:
################################################################################
#
#  Type from field
#
################################################################################

quadratic_space_type(K::S) where {S <: Field} =
    QuadSpace{S, dense_matrix_type(elem_type(S))}

################################################################################
#
#  Constructors
#
################################################################################

@doc Markdown.doc"""
    quadratic_space(K::NumField, n::Int) -> QuadSpace

Create the quadratic space over `K` with dimension `n` and Gram matrix
equal to the identity matrix.
"""
function quadratic_space(K::Field, n::Int)
  @req n >= 0 "Dimension ($n) must be positive"
  G = identity_matrix(K, n)
  return QuadSpace(K, G)
end

@doc Markdown.doc"""
    quadratic_space(K::NumField, G::Int) -> QuadSpace

Create the quadratic space over `K` with Gram matrix `G`.
The matrix `G` must be square and symmetric.
"""
function quadratic_space(K::Field, G::MatElem; check::Bool = true)
  if check
    @req issquare(G) "Gram matrix must be square ($(nrows(G)) x $(ncols(G))"
    @req issymmetric(G) "Gram matrix must be symmetric"
  end
  return QuadSpace(K, G)
end

################################################################################
#
#  Predicates
#
################################################################################

isquadratic(V::QuadSpace) = true

ishermitian(V::QuadSpace) = false

_base_algebra(V::QuadSpace) = V.K

################################################################################
#
#  Properties
#
################################################################################

involution(V::QuadSpace) = identity

fixed_field(V::QuadSpace) = base_ring(V)

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

function Base.show(io::IO, V::QuadSpace)
  print(io, "Quadratic space over\n")
  println(io, base_ring(V))
  println(io, "with Gram matrix")
  print(io, gram_matrix(V))
end

################################################################################
#
#  Inner product
#
################################################################################

# TODO: Make this non-allocating using an additonal temporary vector
function _inner_product(V, v, w)
  mv = matrix(base_ring(V), 1, nrows(V), v)
  mw = matrix(base_ring(V), ncols(V), 1, w)
  return (mv * V * mw)[1, 1]
end

inner_product(V::QuadSpace, v::Vector, w::Vector) = _inner_product(gram_matrix(V), v, w)

################################################################################
#
#  Diagonalization
#
################################################################################

function diagonal(V::QuadSpace)
  D, _ = _gram_schmidt(gram_matrix(V), involution(V))
  return diagonal(D)
end

################################################################################
#
#  Hasse and Witt invariant
#
################################################################################

# Auxiliary function which works with a diagonal
function _hasse_invariant(D::Vector, p)
  h = 1
  n = length(D)
  for i in 1:n
    for j in (i + 1):n
      h = h * hilbert_symbol(D[i], D[j], p)
    end
  end
  return h
end

@doc Markdown.doc"""
    hasse_invariant(V::QuadSpace, p::Union{InfPlc, NfOrdIdl}) -> Int

Returns the Hasse invariant of the quadratic space `V` at `p`. This is equal
to the product of local Hilbert symbols $(a_i, a_j)_p$, $i < j$, where $V$ is
isometric to $\langle a_1,\dotsc,a_n\rangle$.
"""
function hasse_invariant(V::QuadSpace, p)
  return _hasse_invariant(diagonal(V), p)
end

# This can be refactored to operate on the diagonal of a gram schmidt basis and
# the gram matrix.
# (Probably only on the diagonal of a gram schmidt basis)
function witt_invariant(L::QuadSpace, p::NfOrdIdl)
  h = hasse_invariant(L, p)
  F = gram_matrix(L)
  dett = det(F)
  K = base_ring(L)
  ncolsFmod8 = mod(ncols(F), 8)
  if ncolsFmod8 == 3 || ncolsFmod8 == 4
    c = -dett
  elseif ncolsFmod8 == 5 || ncolsFmod8 == 6
    c = K(-1)
  elseif ncolsFmod8 == 7 || ncolsFmod8 == 0
    c = dett
  else
    c = K(1)
  end
  return h * hilbert_symbol(K(-1), c, p)
end

# di = determinant
# wi = witt invariant
# ni = rank
# Lam p. 117
function _witt_of_orthgonal_sum(d1, w1, n1, d2, w2, n2, p)
  _n1 = mod(n1, 4)
  if _n1 == 0 || _n1 == 1
    disc1 = d1
  else
    disc1 = -d1
  end

  _n2 = mod(n2, 4)
  if _n2 == 0 || _n2 == 1
    disc2 = d2
  else
    disc2 = -d2
  end

  if n1 % 2 == n2 % 2
    w3 = w1 * w2 * hilbert_symbol(disc1, disc2, p)
  elseif n1 % 2 == 1
    w3 = w1 * w2 * hilbert_symbol(-disc1, disc2, p)
  else
    @assert n2 % 2 == 1
    w3 = w1 * w2 * hilbert_symbol(disc1, -disc2, p)
  end
  return d1 * d2, w3, n1 + n2
end

# n = rank, d = det
function _witt_hasse(s, n, d, p)
  nmod8 = mod(n, 8)
  K = parent(d)
  if nmod8 == 3 || nmod8 == 4
    c = -d
  elseif nmod8 == 5 || nmod8 == 6
    c = K(-1)
  elseif nmod8 == 7 || nmod8 == 0
    c = d
  else
    c = K(1)
  end
  return s * hilbert_symbol(K(-1), c, p)
end

function witt_invariant(L::QuadSpace, p::InfPlc)
  if iscomplex(p)
    return 1
  end

  h = hasse_invariant(L, p)
  F = gram_matrix(L)
  dett = det(F)
  K = base_ring(L)
  ncolsFmod8 = mod(ncols(F), 8)
  if ncolsFmod8 == 3 || ncolsFmod8 == 4
    c = -dett
  elseif ncolsFmod8 == 5 || ncolsFmod8 == 6
    c = K(-1)
  elseif ncolsFmod8 == 7 || ncolsFmod8 == 0
    c = dett
  else
    c = K(1)
  end
  @assert !iszero(c)
  if isnegative(c, p)
    return -h
  else
    return h
  end
end

@doc Markdown.doc"""
    witt_invariant(V::QuadSpace, p::Union{InfPlc, NfOrdIdl}) -> Int

Returns the Witt invariant of the quadratic space `V` at `p`.

See [Definition 3.2.1, Kir16].
"""
witt_invariant(V::QuadSpace, p)

################################################################################
#
#  Local equivalence
#
################################################################################

function isequivalent(L::QuadSpace, M::QuadSpace, p::NfOrdIdl)
  GL = gram_matrix(L)
  GM = gram_matrix(M)
  if GL == GM
    return true
  end

  return rank(GL) == rank(GM) &&
         islocal_square(det(GL) * det(GM), p) &&
         hasse_invariant(L, p) == hasse_invariant(M, p)
end

function isequivalent(L::QuadSpace, M::QuadSpace, p::InfPlc)
  if rank(L) != rank(M)
    return false
  end

  if iscomplex(p)
    return true
  end

  DL = diagonal(L)
  DM = diagonal(M)
  return count(x -> isnegative(x, p), DL) == count(x -> isnegative(x, p), DM)
end

################################################################################
#
#  Quadratic form with given invariants
#
################################################################################

function _quadratic_form_invariants(M::fmpq_mat; minimal = true)
  G, _ = _gram_schmidt(M, identity)
  D = diagonal(G)
  sup = fmpz[]
  for i in 1:length(D)
    for (p, e) in factor(numerator(D[i]))
      if isodd(e)
        push!(sup, p)
      end
    end
    for (p, e) in factor(denominator(D[i]))
      if isodd(e)
        push!(sup, p)
      end
    end
  end
  push!(sup, fmpz(2))
  sup = unique!(sup)
  F = Dict{fmpz, Int}()
  for p in sup
    e = _hasse_invariant(D, p)
    if e == -1 | !minimal
      F[p] = e
    end
  end
  I = [ (inf, count(x -> x < 0, D)) ]
  nu = numerator(prod(D))
  de = denominator(prod(D))
  return squarefree_part(de * nu), F, I
end

function _quadratic_form_invariants(M; minimal = true)
  return _quadratic_form_invariants(M, maximal_order(base_ring(M)), minimal = minimal)
end

function _quadratic_form_invariants(M, O; minimal = true)
  G, _ = _gram_schmidt(M, identity)
  D = diagonal(G)
  K = base_ring(M)
  sup = Dict{ideal_type(O), Bool}()
  for i in 1:length(D)
    f = factor(D[i] * O)
    for (P, e) in f
      if isodd(e)
        sup[P] = true
      end
    end
  end
  for (P, e) in prime_decomposition(O, 2)
    sup[P] = true
  end
  F = Dict{ideal_type(O), Int}()
  for P in keys(sup)
    e = _hasse_invariant(D, P)
    if e == -1 || !minimal
      F[P] = e
    end
  end
  I = [ (P, count(x -> isnegative(x, P), D)) for P in real_places(K) ];
  return reduce(*, D, init = one(K)), F, I
end

@doc Markdown.doc"""
    invariants(M::QuadSpace)
          -> FieldElem, Dict{NfOrdIdl, Int}, Vector{Tuple{InfPlc, Int}}

Returns a triple `(d, H, I)` of invariants of `M`, which determine the
equivalence class completely. The element `d` is the determinant of a Gram
matrix, `H` contains the non-trivial Hasse invariants and `I` contains for
each real place the negative index of inertia.

Note that `d` is determined only modulo squares.
"""
invariants(V::QuadSpace) = _quadratic_form_invariants(gram_matrix(V))

################################################################################
#
#  Global equivalence
#
################################################################################

@doc Markdown.doc"""
    isequivalent(M::QuadSpace, L::QuadSpace) -> Bool

Tests if `M` and `L` are equivalent.
"""
function isequivalent(M::QuadSpace, L::QuadSpace)
  if gram_matrix(M) == gram_matrix(L)
    return true
  end
  d1, H1, I1 = invariants(M)
  d2, H2, I2 = invariants(L)
  return I1 == I2 && H1 == H2 && issquare(d1 * d2)[1]
end

################################################################################
#
#  Quadratic form with given invariants
#
################################################################################

# The following is over Q
function _quadratic_form_with_invariants(dim::Int, det::fmpz,
                                         finite::Vector{fmpz}, negative::Int)
#{Computes a quadratic form of dimension Dim and determinant Det that has Hasse invariants -1 at the primes in Finite.
 #The number of negative entries of the real signature is given in Negative}
  @assert dim >= 1
  @assert !iszero(det)
  @assert negative in 0:dim

  sign(det) != (-1)^(negative % 2) && throw(error("Real place information does not match the sign of the determinant"))

  if dim == 1
    !isempty(finite) && throw(error("Impossible Hasse invariants"))
    return matrix(FlintQQ, 1, 1, fmpz[det])
  end

  finite = unique(finite)
  @assert all(isprime(p) for p in finite)

  if dim == 2
    ok = all(!islocal_square(-det, p) for p in finite)
    if !ok
      q = [p for p in finite if islocal_square(-det, p)][1]
      throw(error("A binary form with determinant $det must have Hasse invariant +1 at the prime $q"))
    end
  end

  # product formula check

  !iseven((negative % 4 >= 2 ? 1 : 0) + length(finite)) && throw(error("The number of places (finite or infinite) with Hasse invariant -1 must be even"))

  # reduce the number of bad primes
  det = squarefree_part(det)

  dim0 = dim
  det0 = det
  finite0 = copy(finite)
  negative0 = negative

#  // Pad with ones
  k = max(0, dim - max(3, negative))
  D = ones(Int, k)
  dim = dim - k

#  // Pad with minus ones
  if dim >= 4
    @assert dim == negative
    k = dim - 3
    d = (-1)^k
    f = (k % 4 >= 2) ? Set(fmpz[2]) : Set(fmpz[])
    PP = append!(fmpz[p for (p, e) in factor(2 * det)], finite)
    PP = unique!(PP)
    finite = fmpz[ p for p in PP if hilbert_symbol(d, -det, p) * (p in f ? -1 : 1) * (p in finite ? -1 : 1) == -1]
    finite = unique!(finite)
    D = append!(D, Int[-1 for i in 1:k])
    det = isodd(k) ? -det : det
    dim = 3
    negative = 3
  end

  # ternary case
  if dim == 3
#    // The primes at which the form is anisotropic
    PP = append!(fmpz[p for (p, e) in factor(2 * det)], finite)
    PP = unique!(PP)
    PP = filter!(p -> hilbert_symbol(-1, -det, p) != (p in finite ? -1 : 1), PP)
#    // Find some a such that for all p in PP: -a*Det is not a local square
#    // TODO: Find some smaller a?! The approach below is very lame.
    a = prod(det % p == 0 ? one(FlintZZ) : p for p in PP)
    if negative == 3
      a = -a
      negative = 2
    end

    PP = append!(fmpz[p for (p, e) in factor(2 * det * a)], finite)
    PP = unique!(PP)
    finite = fmpz[ p for p in PP if hilbert_symbol(a, -det, p) * (p in finite ? -1 : 1) == -1]
    det = squarefree_part(det * a)
    dim = 2
    push!(D, a)
  end

#  // The binary case
  a = _find_quaternion_algebra(fmpq(-det), finite, negative == 2 ? PosInf[inf] : PosInf[])
  Drat = map(FlintQQ, D)
  Drat = append!(Drat, fmpq[a, squarefree_part(FlintZZ(det * a))])

  M = diagonal_matrix(Drat)

  d, f, n = _quadratic_form_invariants(M)

  @assert dim0 == length(Drat)
  @assert d == det0
  @assert issetequal(collect(keys(f)), finite0)
  @assert n[1][2] == negative0
  return M
end

function _quadratic_form_with_invariants(dim::Int, det::fmpq,
                                         finite::Vector{fmpz}, negative::Int)
  _det = numerator(det) * denominator(det)
  return _quadratic_form_with_invariants(dim, _det, finite, negative)
end

#{Computes a quadratic form of dimension Dim and determinant Det that has Hasse invariants -1 at the primes in Finite.
# The number of negative entries of the i-th real signature is given in Negative[i]}
function _quadratic_form_with_invariants(dim::Int, det::nf_elem, finite::Vector, negative::Dict{InfPlc, Int})
  @assert dim >= 1
  @assert !iszero(det)
  K = parent(det)
  inf_plcs = real_places(K)
  @assert length(inf_plcs) == length(negative)
  # All real places must be present
  @assert all(Bool[0 <= c <= dim for (_, c) in negative])
  # Impossible negative entry at plc
  @assert all(sign(det, p) == (-1)^(negative[p]) for p in inf_plcs)
  # Information at the real place plc does not match the sign of the determinant

  if dim == 1
    @assert isempty(finite) # Impossible Hasse invariants
    return matrix(K, 1, 1, nf_elem[det])
  end

  local OK::order_type(K)

  if !isempty(finite)
    OK = order(finite[1])
    @assert ismaximal(OK)
  else
    OK = maximal_order(K)
  end

  finite = unique(finite)

  # Finite places check

  if dim == 2
    ok = all(!islocal_square(-det, p) for p in finite)
    if !ok
      q = eltype(finite)[p for p in finite if islocal_square(-det, p)][1]
      throw(error("A binary form with determinant $det must have Hasse invariant +1 at the prime $q"))
    end
  end

  @assert iseven(length([ p for (p, n) in negative if n % 4 >= 2]) + length(finite))
 #   "The number of places (finite or infinite) with Hasse invariant -1 must be even";

 # // OK, a space with these invariants must exist.
 # // For final testing, we store the invariants.

  dim0 = dim
  det0 = det
  finite0 = copy(finite)
  finite = copy(finite)
  negative = copy(negative)
  negative0 = copy(negative)

  # det = _reduce_modulo_squares(det)

  k = max(0, dim - max(3, maximum(values(negative))))
  D = elem_type(K)[one(K) for i in 1:k]
  dim = dim - k
  local D2::Vector{nf_elem}
  local D::Vector{nf_elem}

  if dim >= 4
#    // Pad with minus ones
    k = min(dim - 3, minimum(values(negative)))
    D2 = elem_type(K)[-one(K) for i in 1:k]
    dim = dim - k
    for (p, n) in negative
      negative[p] = n - k
    end
#    // Pad with other entries
    while dim >= 4
      V = InfPlc[]
      _signs = Int[]
      for (p, n) in negative
        if n == 0
          push!(V, p)
          push!(_signs, +1)
        elseif n == dim
          push!(V, p)
          push!(_signs, -1)
        end
      end

      x = _weak_approximation(V, _signs)::nf_elem
      s = signs(x)
      #@assert all(Bool[sign(x, V[i]) == _signs[i] for i in 1:length(V)])
      let negative = negative, dim = dim
        k = minimum(vcat(Int[dim - 3], Int[s[p] == 1 ? (dim - c) : c for (p, c) in negative]))
      end
      D2 = append!(D2, elem_type(K)[x for i in 1:k])
      dim = dim - k
      for (p, n) in negative
        if s[p] == -1
          negative[p] = negative[p] - k
        end
      end
    end

    local _d::nf_elem
    local _f::Dict{NfAbsOrdIdl{AnticNumberField,nf_elem},Int64}
    _d, _f = _quadratic_form_invariants(diagonal_matrix(D2))

    PP = append!(support(K(2), OK), finite)
    PP = unique!(PP)
    local _finite::Vector{ideal_type(OK)}
    let finite = finite
      _finite = ideal_type(OK)[ p for p in PP if hilbert_symbol(_d, -det, p) * (haskey(_f, p) ? -1 : 1) * (p in finite ? -1 : 1) == -1]
    end
    finite = _finite

    D = append!(D, D2)

    det::nf_elem = det * _d
#    # TODO: reduce det modulo squares
  end

#  // The ternary case
  if dim == 3
    PP = append!(support(K(2), OK), finite)
    append!(PP, support(det, OK))
    PP = unique!(PP)
    PP = ideal_type(OK)[p for p in PP if hilbert_symbol(K(-1), -det, p) != (p in finite ? -1 : 1)]
#    // The primes at which the form is anisotropic

#    // Find some a such that for all p in PP: -a*Det is not a local square
#    // TODO: Find some smaller a?! The approach below is very lame.
#    // We simply make sure that a*Det has valuation 1 at each prime in PP....

    if length(PP) == 0
      a = one(K)
    else
      a = approximate(Int[(1 + valuation(det, p)) % 2 for p in PP], PP)
    end
#    // Fix the signs of a if necessary.
    s = signs(a)
    idx = InfPlc[ p for (p, n) in negative if n in [0, 3]]
    S = Int[ negative[p] == 0 ? s[p] : -s[p] for p in idx]
    if length(PP) > 0
      b = _weak_approximation_coprime(idx, S, prod(PP))
      @assert iscoprime(b * OK, prod(PP))
    else
      b = _weak_approximation_coprime(idx, S, 1 * OK)
    end
    a = a * b

#    // Adjust invariants for the last time:
    s = signs(a)
    for p in InfPlc[p for (p,c) in negative if s[p] < 0]
      negative[p] = negative[p] - 1
    end
    PP = support(K(2))
    append!(PP, support(det, OK))
    append!(PP, support(a, OK))
    append!(PP, finite)
    PP = unique!(PP)
    finite = ideal_type(OK)[p for p in PP if hilbert_symbol(a, -det, p) * (p in finite ? -1 : 1) == -1]
    det = det * a
    # TODO: reduce det
    push!(D, a)
  end

#  // The binary case
  a = _find_quaternion_algebra(-det::nf_elem, finite::Vector{NfOrdIdl}, InfPlc[p for (p, n) in negative if n == 2])
  push!(D, a)
  push!(D, det * a)
  M = diagonal_matrix(D)

  d, f, n = _quadratic_form_invariants(M, OK)
  @assert dim0 == length(D)
  @assert issquare(d * det0)[1]
  @assert issetequal(collect(keys(f)), finite0)
  @assert issetequal(n, collect((p, n) for (p, n) in negative0))

  return M
end

################################################################################
#
#  Isotropic
#
################################################################################

isisotropic(V::QuadSpace, p::InfPlc) = _isisotropic(V, p)

function _isisotropic(D::Array, p)
  n = length(D)
  if n == 0
    return false
  end
  K = parent(D[1])
  d = reduce(*, D, init = one(K))
  if d == 0
    return true
  elseif n <= 1
    return false
  elseif n == 2
    return islocal_square(-d, p)
  elseif n == 3
    return _hasse_invariant(D, p) == hilbert_symbol(K(-1), -d, p)
  elseif n == 4
    return !islocal_square(d, p) || (_hasse_invariant(D, p) == hilbert_symbol(K(-1), K(-1), p))
  else
    return true
  end
end

function isisotropic(V::QuadSpace, p)
  @assert base_ring(V) == nf(order(p))
  d = det(V)
  n = rank(V)
  K = base_ring(V)
  if d == 0
    return true
  elseif n <= 1
    return false
  elseif n == 2
    return islocal_square(-d, p)
  elseif n == 3
    return hasse_invariant(V, p) == hilbert_symbol(K(-1), K(-1), p)
  elseif n == 4
    return !islocal_square(d, p) || (hasse_invariant(V, p) == hilbert_symbol(K(-1), K(-1), p))
  else
    return true
  end
end

################################################################################
#
#  Embeddings
#
################################################################################

# This is O'Meara, 63:21
#
# n, a, ha = dimension, determinant (class) and Hasse symbol of first space
# Similar for m, a, hb
# p is the prime idela
function _can_locally_embed(n::Int, da, ha::Int, m::Int, db, hb::Int, p)
  de = m - n
  if de == 0
    return islocal_square(da * db, p) && ha == hb
  elseif de == 1
    return ha * hilbert_symbol(da * db, da, p) == hb
  elseif de == 2 && islocal_square(-da * db, p)
    # Test if U \perp H \cong V
    # H has Hasse invariant 1
    return ha * hilbert_symbol(da, -1, p) == hb
  else
    return true
  end
end

@doc Markdown.doc"""
    islocally_represented_by(U::QuadSpace, V::QuadSpace, p)

Return whether $U$ is represented by $V$ locally at $\mathfrak p$.
"""
function islocally_represented_by(U::QuadSpace, V::QuadSpace, p)
  n, da, ha = rank(U), det(U), hasse_invariant(U, p)
  m, db, hb = rank(V), det(V), hasse_invariant(V, p)
  return _can_locally_embed(n, da, ha, m, db, hb, p)
end

# We are using O'Meara p. 160 ff.
#
# We need to test if U is represented by V at all places of K.
#
# If the place is complex, there is only the trivial restriction (ranks)
# If the place is real, we have to check if the signatures fit.
#
# If p does not divide 2 * dU * dV and h(U) = 1 = h(V) (Hasse-invariant),
# then U is represented by V locally at p. This follows from the local
# characterization. But the Hasse invariant is zero outside the support
# of the diagonal. Thus we get only finitely many conditions.

@doc Markdown.doc"""
    isrepresented_by(U::QuadSpace, V::QuadSpace)

Return whether $U$ is represented by $V$, that is, whether $U$ embeds into $V$.
"""
function isrepresented_by(U::QuadSpace, V::QuadSpace)
  v = rank(V) - rank(U)
  if v < 0
    return false
  end

  if v == 0
    return isequivalent(U, V)
  end

  K = base_ring(U)

  rlp = real_places(K)

  dU = diagonal(U)
  dV = diagonal(V)

  rkU = rank(U)
  rkV = rank(V)

  negU = Int[ count(x -> isnegative(x, P), dU) for P in rlp ]
  signU = Tuple{Int, Int}[ (i, rkU - i) for i in negU]

  negV = Int[ count(x -> isnegative(x, P), dV) for P in rlp ]
  signV = Tuple{Int, Int}[ (i, rkV - i) for i in negV]

  for i in 1:length(rlp)
    if signU[i][1] > signV[i][1] || signU[i][2] > signV[i][2]
      return false
    end
  end

  OK = maximal_order(K)

  ds = elem_type(OK)[]

  for d in dU
    push!(ds, OK(numerator(d)))
    push!(ds, OK(denominator(d)))
  end

  for d in dV
    push!(ds, OK(numerator(d)))
    push!(ds, OK(denominator(d)))
  end

  push!(ds, OK(2))

  lp = coprime_base(ds)

  # lp is a list of coprime integral ideals, such that elements in ds factor
  # over lp. But these ideals are not necessarily prime.

  for a in lp
    for p in support(a)
      if !islocally_represented_by(U, V, p)
        return false
      end
    end
  end

  return true
end

################################################################################
#
#  Isometry computation
#
################################################################################

function _solve_conic(a::Integer, b::Integer, c::Integer)
  _solve_conic(fmpq(a), fmpq(b), fmpq(c))
end

function _solve_conic(a, b, c, u, v)

  K = parent(a)
	@assert !iszero(a)
	@assert !iszero(b)
	@assert !iszero(c)

	fl, z = ispower(-b//a, 2)
	if fl
    x, y, z = z, K(1), K(0)
    @goto finish
  end

	fl, z = ispower(-c//a, 2)
	if fl
    x, y, z = z, K(0), K(1)
    @goto finish
  end

  Kx, x = PolynomialRing(K, "x", cached = false)
  d = -b//a
  den = denominator(d)
  L, y = number_field(x^2 - d * den^2)
  _Ls, _LstoL = simplify(L)
  fl, _n = isnorm(L, -c//a)
  if L isa AnticNumberField
    n = evaluate(_n)
  else
    n = _n
  end
  if fl
    x, y, z = coeff(n, 0), coeff(n, 1) * den, K(1)
    @goto finish
  end

  return false, a, a, a, u, u, u

  @label finish

  @assert x^2 * a + y^2 * b + z^2 * c == 0

  # Cremona, Conic paper
  # x = Q1(U, V) = ax0U^2 + 2by0UV − bx0V^2
  # y = Q2(U, V) = −ay0U^2 + 2ax0UV + by0V^2
  # z = Q3(U, V) = az0U^2 + bz0V^2

  q1 = a * x * u^2 + 2 * b * y * u * v - b * x * v^2
  q2 = -a * y * u^2 + 2*a*x*u*v + b*y*v^2
  q3 = a*z*u^2 + b*z*v^2

  @assert a * q1^2 + b * q2^2 + c * q3^2 == 0

  return true, x, y, z, q1, q2, q3
end

function _isisometric_with_isometry(a1, a2, b1, b2)
  # I assume that they are isometric ...
  #
  # I want to find an isometry from (a1, a2) to (b1, b2)
  # Let us call the matrix (a b; c d)
  # Then a^2 a_1 + b^2 a_2 = z1^2 * b1 and
  #

  K = parent(a1)
  Kuv, (u, v) = PolynomialRing(K, ["u", "v"], cached = false)

  fl, _aa, _bb, _z1, a, b, z1 = _solve_conic(a1, a2, -b1, u, v)
  @show _aa, _bb, _z1
  @assert fl

  # a^2 a_1 + b^2 a_2 = z2^2 b2 and
  fl, _cc, _dd, _z2, c, d, z2 = _solve_conic(a1, a2, -b2, u, v)
  @show _cc, _dd, _z2
  @assert fl

  @show _aa * _cc * a1 + _bb * _dd * a2

  @show a
  @show b
  @show c
  @show d

  # a * c * a1 + b * d * a2 = 0

  @show z1, z2

  s =  a * c * a1 + b * d * a2
  if s == 0
    return _aa, _bb, _cc, _dd, _z1, _z2
  end
  _a, _b, _c = coeff(s, u^4), coeff(s, u^2 * v^2), coeff(s, v^4)
  @show _a, _b, _c
  @show s
  if 4 * _a * _c == _b^2
    @assert 4*_c*s == (_b * u^2 + 2 * _c * v^2)^2
    # u^2//v^2 = -b/c
    fl, z = ispower(-(2 * _c)//_b, 2)
    # (u/v)^2 == -2c/b
    @assert fl
    v = one(K)
    u = z
    @assert b * u^2 + 2 * c * v^2 == 0
    @assert s(u, v) == 0
  end

  # This should be a parabola?
end

function _solve_conic_affine(A, B, a)
  # Solve Au^2 + B*w^2 = a
  # Gives one solutation

  # a = u^2 + B/A v^2 = (u - sqrt(B/A)v)(u + sqrt(B/A)) = N(u + v sqrt(B/A))

  K = parent(A)

  Kz, z = PolynomialRing(K, "z", cached = false)
  D = -B//A
  de = denominator(D)
  L, _ = number_field(z^2 - de^2 * D)
  fl, _n = isnorm(L, a//(A) * de^2)

  if !fl
    return false, zero(K), zero(K)
  end

  if L isa AnticNumberField
    n = evaluate(_n)
  else
    n = _n
  end

  @assert norm(n) == a//(A) * de^2

  u1, w1 = coeff(n, 0)//de, coeff(n, 1)

  @assert u1^2 * A + w1^2 * B == a

  return true, u1, w1
end

function _solve_conic_affine(A, B, a, t)
  # Solve Au^2 + B*w^2 = a
  # Gives one solutation and a parametrization
  # This assumes that a solution exists!

  # a = u^2 + B/A v^2 = (u - sqrt(B/A)v)(u + sqrt(B/A)) = N(u + v sqrt(B/A))

  K = parent(A)

  Kz, z = PolynomialRing(K, "z", cached = false)
  D = -B//A
  de = denominator(D)
  L, _ = number_field(z^2 - de^2 * D)
  fl, _n = isnorm(L, a//(A) * de^2)

  @assert fl

  if L isa AnticNumberField
    n = evaluate(_n)
  else
    n = _n
  end

  @assert norm(n) == a//(A) * de^2

  u1, w1 = coeff(n, 0)//de, coeff(n, 1)

  @assert u1^2 * A + w1^2 * B == a
  u = (-A * u1 + B * t^2 * u1 - 2 * B * t * w1)//(A + B * t^2)
  w = (-2 * A * t * u1   + A * w1 - B * t^2 * w1)//(A + B * t^2)

  @assert u^2 * A + w^2 * B == a

  return true, u1, w1, u, w
end

# Return true, T such that T * [A 0; 0 B] T^t = [a 0; 0 b] or false, 0 if no such T exists.
function _isisometric_with_isometry_dan(A, B, a, b)
  K = parent(A)

  Kkt, (k, t) = PolynomialRing(K, ["k", "t"], cached = false)

  fl, u1, w1, u, w = _solve_conic_affine(A, B, a, t)
  if !fl
    return false, zero_matrix(K, 0, 0)
  end

  fl, s3, v3, s, v = _solve_conic_affine(B, A, b, k)
  if !fl
    return false, zero_matrix(K, 0, 0)
  end

  lin = ((2 * (-2 * A^2 * B * s3 * t * u1 + A^3 * u1 * v3 - A^2 * B * t^2 * u1 * v3 + A^2 * B * s3 * w1 - A * B^2 * s3 * t^2 * w1 + 2 * A^2 * B * t * v3 * w1))) * k - (-2 * A^2 * B * s3 * u1 +  2 * A * B^2 * s3 * t^2 * u1 - 4 * A^2 * B * t * u1 * v3 - 4 * A * B^2 * s3 * t * w1 + 2 * A^2 * B * v3 * w1 - 2 * A * B^2 * t^2 * v3 * w1)
  sq = 4 * A * B * (A + B * t^2)^2 * (B * s3^2 + A * v3^2) * (A * u1^2 + B * w1^2)

  junk = 4 * (-2 * A^2 * B * s3 * t * u1 + A^3 * u1 * v3 - A^2 * B * t^2 * u1 * v3 + A^2 * B * s3 * w1 -  A * B^2 * s3 * t^2 * w1 + 2 * A^2 * B * t * v3 * w1) * (B + A * k^2) * (A + B * t^2)

  t0 = K(1)
  @assert !iszero(A + B * t0^2)

  middle = A * u * v + B * s * w

  @assert lin^2 - sq == junk * middle

  _sq = sq(0, t0)

  fl, rt = ispower(_sq, 2)

  if !fl
    return false, zero_matrix(K, 0, 0)
  end

  k0 = (rt + (-2 * A^2 * B * s3 * u1 +  2 * A * B^2 * s3 * t^2 * u1 - 4 * A^2 * B * t * u1 * v3 - 4 * A * B^2 * s3 * t * w1 + 2 * A^2 * B * v3 * w1 - 2 * A * B^2 * t^2 * v3 * w1))//((2 * (-2 * A^2 * B * s3 * t * u1 + A^3 * u1 * v3 - A^2 * B * t^2 * u1 * v3 + A^2 * B * s3 * w1 - A * B^2 * s3 * t^2 * w1 + 2 * A^2 * B * t * v3 * w1)))

  kk = numerator(k0)(0, t0)//denominator(k0)(0, t0)

  @assert !iszero(junk(kk, t0))
  @assert !iszero(B + A * kk^2)

  uu = numerator(u)(kk, t0)//denominator(u)(kk, t0)
  ww = numerator(w)(kk, t0)//denominator(w)(kk, t0)
  ss = numerator(s)(kk, t0)//denominator(s)(kk, t0)
  vv = numerator(v)(kk, t0)//denominator(v)(kk, t0)

  T = matrix(K, 2, 2, [uu, ww, vv, ss])
  D1 = diagonal_matrix([A, B])
  D2 = diagonal_matrix([a, b])
  @assert T * D1 * transpose(T) == D2

  return true, T
end

@doc Markdown.doc"""
    isisometric_with_isometry(V::QuadSpace, W::QuadSpace)

Returns wether $V$ and $W$ are isometric together with an isometry in case it
exists. The isometry is given as an invertible matrix $T$ such that
$T G_W T^t = G_V$, where $G_V$, $G_W$ are the Gram matrices.
"""
function isequivalent_with_isometry(V::QuadSpace, W::QuadSpace)
  if !isequivalent(V, W)
    return false, zero_matrix(base_ring(V), 0, 0)
  end

  @req max(rank(V), rank(W)) <= 2 "Rank must be <= 2"

  K = base_ring(V)

  GV = gram_matrix(V)
  GW = gram_matrix(W)

  DV, MV = _gram_schmidt(gram_matrix(V), involution(V))
  DW, MW = _gram_schmidt(gram_matrix(W), involution(W))

  A, B = DV[1, 1], DV[2, 2]
  a, b = DW[1, 1], DW[2, 2]

  @assert MV * GV * transpose(MV) == diagonal_matrix([A, B])
  @assert MW * GW * transpose(MW) == diagonal_matrix([a, b])

  fl, T = _isisometric_with_isometry_dan(A, B, a, b)
  @assert fl

  @assert T * DV * transpose(T) == DW

  # T * DV * T^t == DW
  # T * MV * GV * (T * MV)^t == MW * GW * MW^t
  # GV = MV^-1 * T^-1 * MW * GW * (MV^-1 * T^-1 * MW)^t

  T = inv(MV) * inv(T) * MW
  @assert T * GW * transpose(T) == GV
  return true,  T
end

################################################################################
#
#  Isotropic vector
#
################################################################################

_to_gf2(x) = x == 1 ? 0 : 1
#ToGF2:= func< a,b,p | HilbertSymbol(a,b,p) eq 1 select 0 else 1 >;
#SignGF2:= func< x, p | Evaluate(x, p) lt 0 select 1 else 0 >;
#MyFact:= func< R, d | Type(R) eq RngInt select FactorizationOfQuotient(Rationals() ! d) else Factorization(R*d) >;

# F must be symmetric

function _isisotropic_with_vector(F::MatrixElem)
  K = base_ring(F)
  _D, T = _gram_schmidt(F, identity, false)
  D = diagonal(_D)
  i = findfirst(==(zero(K)), D)
  if i isa Int
    return true, elem_type(K)[T[i, j] for j in 1:ncols(T)]
  end

  if length(D) <= 1
    return false, elem_type(K)[]
  end

  for i in 1:length(D)
    for j in (i + 1):length(D)
      if D[i] == -D[j]
        return true, elem_type(K)[T[i, k] + T[j, k] for k in 1:ncols(T)]
      end
    end
  end

  fl, y = issquare_with_root(-D[1]//D[2])
  if fl
    return true, elem_type(K)[T[1, k] + y * T[2, k] for k in 1:ncols(T)]
  elseif length(D) == 2
    return false, elem_type(K)[]
  end

  if length(D) == 3
    fl, a, b = _solve_conic_affine(D[1], D[2], -D[3])
    if fl
      v = elem_type(K)[a, b, one(K)]
      vma = matrix(K, 1, length(v), v) * T
      @assert vma * F * transpose(vma) == 0
      return true, elem_type(K)[vma[1, i] for i in 1:ncols(vma)]
    else
      return false, elem_type(K)[]
    end
  elseif length(D) == 4
    fl, v = _isisotropic_with_vector(diagonal_matrix(D[3], D[4]))
    if fl
      return true, elem_type(K)[v[1] * T[3, k] + v[2] * T[4, k] for k in 1:ncols(T)]
    end

    for v in real_places(K)
      if !_isisotropic(D, v)
        return false, elem_type(K)[]
      end
    end

    R = maximal_order(K)
    P = ideal_type(R)[]
    for d in append!(elem_type(K)[K(2)], D)
      for (p, _) in factor(d * R)
        if p in P
          continue
        end
        if !_isisotropic(D, p)
          return false, elem_type(K)[]
        else
          push!(P, p)
        end
      end
    end

    # At this point we know that the space is isotropic.
    # But we need to determine the vector.

    # Find x != 0 such that <D[1], D[2]> and <-D[3], -D[4]> both represent x.

    rlp = real_places(K)

    _target = append!(Int[_to_gf2(hilbert_symbol(D[1], D[2], p)) for p in P], Int[_to_gf2(hilbert_symbol(-D[3], -D[4], p)) for p in P])

    I = eltype(rlp)[]

    for p in rlp
      s = sign(D[1], p)
      if s == sign(D[2], p)
        push!(I, p)
        push!(_target, _to_gf2(s))
      else
        s = sign(-D[3], p)
        if s == sign(-D[4], p)
          push!(I, p)
          push!(_target, _to_gf2(s))
        end
      end
    end

    V = abelian_group(Int[2 for i in 1:length(_target)])
    target = V(_target)
    # Find x such that target equals the vector
    # [ _to_gf2_(hilbert_symbol(-D[1] * D[2], x, p)) for p in ] vcat
    # [ _to_gf2(hilbert_symbol(-D[3]*D[4], x, p)) for p in P ] vcat
    # [ _to_gf2(sign(x, p) for p in I ]
    if iszero(target)
      x = one(K)
    else
      found = false
      S, mS = sub(V, elem_type(V)[], false)
      basis = elem_type(K)[]
      signs = elem_type(V)[]
      L, mL = sunit_group_fac_elem(P)
      Q, mQ = quo(L, 2, false)
      for q in gens(Q)
        x = evaluate(mL(mQ\q))
        _v = append!(Int[_to_gf2(hilbert_symbol(-D[1] * D[2], x, p)) for p in P], Int[_to_gf2(hilbert_symbol(-D[3] * D[4], x, p)) for p in P])
        _v = append!(_v, Int[_to_gf2(sign(x, p)) for p in I])
        s = V(_v)
        fl, _ = haspreimage(mS, s)
        if !fl
          push!(signs, s)
          push!(basis, x)
          S, mS = sub(V, signs, false)
          if haspreimage(mS, target)[1]
            found = true
            break
          end
        end
      end

      # Still not found

      Cl, mCl = class_group(R)
      A = abelian_group(fill(0, length(P)))
      hh = hom(A, Cl, [mCl\(p) for p in P])
      SS, mSS = image(hh, false)
      Q, mQ = quo(Cl, [mSS(SS[i]) for i in 1:ngens(SS)])

      p  = 2
      while !found
        p = next_prime(p)
        lp = prime_decomposition(R, p)
        for (q, _) in lp
          if q in P
            continue
          end
          o = order(mQ(mCl\(q)))
          c = -(hh\(o * (mCl\(q))))
          fl, _x = isprincipal(q * prod(P[i]^Int(c.coeff[i]) for i in 1:length(P)))
          x = elem_in_nf(_x)
          _v = append!(Int[_to_gf2(hilbert_symbol(-D[1] * D[2], x, p)) for p in P], Int[_to_gf2(hilbert_symbol(-D[3] * D[4], x, p)) for p in P])
          _v = append!(_v, Int[_to_gf2(sign(x, p)) for p in I])
          s = V(_v)
          if haspreimage(mS, s + target)[1]
            push!(basis, x)
            push!(signs, s)
            found = true
            break
          end
        end
      end

      FF = GF(2, cached = false)
      fl, expo = can_solve_with_solution(matrix(FF, length(signs), length(_target), [ s.coeff[1, i] for s in signs, i in 1:length(_target)]), matrix(FF, 1, length(_target), _target), side = :left)
      @assert fl

      x = evaluate(FacElem(basis, map(fmpz, [lift(expo[1, i]) for i in 1:length(basis)])))
    end
    ok, v = _isisotropic_with_vector(diagonal_matrix([D[1], D[2], -x]))
    @assert ok
    ok, w = _isisotropic_with_vector(diagonal_matrix([D[3], D[4],  x]))
    @assert ok
    v = inv(v[3]) .* v
    w = inv(w[3]) .* w
    v = matrix(K, 1, 4, [v[1], v[2], w[1], w[2]]) * T
    @assert v * F * v' == 0
    return true, elem_type(K)[v[1, i] for i in 1:4]
  else
    # Dimension >= 5, we need to only take care of the real places
    ok = all(v -> _isisotropic(D, v), real_places(K))
    if !ok
      return false, elem_type(K)[]
    end

    # We need D[3..5] to yield both signs at every real place
    rlp = real_places(K)
    found = false
    for i in 1:length(D), j in (i + 1):length(D)
      if all(p -> sign(D[i], p) != sign(D[j], p), rlp)
        TT = identity_matrix(K, nrows(F))
        found = true
        if i != 3
          swap_cols!(TT, 3, i)
        end
        if j != 4
          swap_cols!(TT, 4, j)
        end
        T = TT * T
        _D = (T * F * transpose(T))
        @assert isdiagonal(_D)
        D = diagonal(_D)
        break
      end
    end
    if !found
      fix = Int[]
      signs = Int[]
      for i in 1:length(rlp)
        s = sign(D[4], rlp[i])
        if s != sign(D[5], rlp[i])
          continue
        end
        if s == sign(D[4], rlp[i])
          _a = _real_weak_approximation(rlp[i], rlp[fix])
          a = inv(_a)
          j = findfirst(Bool[sign(D[j], rlp[i]) != s for j in 1:length(D)])::Int
          r = 0
          while true
            r += 1
            t = D[4] + a^(2*r)*D[j]
            if sign(t, rlp[i]) != s && all(sign(t, rlp[fix[k]]) == signs[k] for k in 1:length(fix))
              break
            end
          end
          b = -a^r * D[j]//D[4]
          v = [T[4, k] for k in 1:ncols(T)]
          for k in 1:ncols(T)
            T[4, k] = T[4, k] + a^r * T[j, k]
          end
          for k in 1:ncols(T)
            T[j, k] = T[j, k] + b * v[k]
          end
        end
        push!(fix, i)
        push!(signs, -s)
        _D = T * F * transpose(T)
        @assert isdiagonal(_D)
        D = diagonal(_D)
      end
    end

    ok, v = _isisotropic_with_vector(diagonal_matrix(D[3:5]))
    if ok
      res = Vector{elem_type(K)}(undef, ncols(T))
      for k in 1:ncols(T)
        res[k] = v[1] * T[3, k] + v[2] * T[4, k] + v[3] * T[5, k]
      end
      return true, res
    end

    R = maximal_order(K)
    P = ideal_type(R)[]
    X = Tuple{elem_type(K), elem_type(K)}[]
    M = ideal_type(R)[]
    for p in Set([ p for d in append!(nf_elem[K(2)], D) for p in support(d, R)])
      if _isisotropic(D[3:5], p)
        continue
      end

      if _isisotropic([D[3], D[4], D[5], D[1]], p)
        x = one(K)
        y = zero(K)
      elseif _isisotropic([D[3], D[4], D[5], D[2]], p)
        x = zero(K)
        y = one(K)
      else
        # now D[1] and D[2] represent necessarily the same class
        # leaving this class is enough
        V1 = valuation(D[1], p)
        V2 = valuation(D[2], p)
        V = max(V1, V2)
        pi = elem_in_nf(uniformizer(p))
        k, h = ResidueField(order(p), p)
        hext = extend(h, K)
        y = pi^(div(V - V2, 2))
        yy = pi^(div(V - V1, 2))
        cnt = 1
        while true
          cnt += 1
          @assert cnt <= 1000
          x = (hext\rand(k)) * yy
          if _isisotropic(elem_type(K)[D[3], D[4], D[5], x^2 * D[1] + y^2 * D[2]], p)
            break
          end
        end
      end
      push!(X, (x, y))
      push!(P, p)
      V = valuation(x^2 * D[1] + y^2 * D[2], p) + 1
      if isdyadic(p)
        V = V + 2 * ramification_index(p)
      end
      @show p, V
      push!(M, p^V)
    end
    @assert length(P) != 0

    xx = elem_in_nf(crt([R(x[1]) for x in X], M))
    yy = elem_in_nf(crt([R(x[2]) for x in X], M))
    t = xx^2 * D[1] + yy^2 * D[2]
    ok, w = _isisotropic_with_vector(diagonal_matrix(elem_type(K)[D[3], D[4], D[5], t]))
    @assert ok
    @assert w[1]^2 * D[3] + w[2]^2 * D[4] + w[3]^2 * D[5] + w[4]^2 * t == 0
    w = inv(w[4]) .* w
    v = matrix(K, 1, ncols(T), append!(elem_type(K)[xx, yy, w[1], w[2], w[3]], [zero(K) for i in 1:(nrows(T) - 5)])) * T
    v = lcm(fmpz[denominator(v[1, i]) for i in 1:ncols(v)]) * v
    @assert v * F * transpose(v) == 0
    return true, elem_type(K)[v[1, i] for i in 1:ncols(v)]
  end
end

function _quadratic_form_decomposition(F::MatrixElem)
  # Decompose F into an anisotropic kernel, an hyperbolic space and its radical
  @req issymmetric(F) "Matrix must be symmetric"
  K = base_ring(F)
  r, Rad = left_kernel(F)
  @assert nrows(Rad) == r
  RadComp = _find_direct_sum(Rad)
  newF = RadComp * F * RadComp'
  H = similar(F, 0, ncols(F))
  CurBas = RadComp

  while true
    fl, HH = _find_hyperbolic_subspace(newF)
    if fl
      @assert iszero(sub(HH, 1:1, 1:ncols(HH)) * newF  * sub(HH, 1:1, 1:ncols(HH))')
      H = vcat(H, HH * CurBas)
      CurBas = _orthogonal_complement(newF, HH) * CurBas
      newF = CurBas * F * CurBas'
    else
      break
    end
  end

  @assert iseven(nrows(H))
  if nrows(H) > 0
    D = diagonal_matrix([matrix(K, 2, 2, [1, 0, 0, -1]) for i in 1:div(nrows(H), 2)])
    @assert isequivalent(quadratic_space(K, H * F * H'), quadratic_space(K, D))
  end

  @assert iszero(Rad * F * Rad')

  #H * F * H', CurBas * F * CurBas', Rad * F * Rad'
  return CurBas, H, Rad
end

function _find_hyperbolic_subspace(F)
  K = base_ring(F)
  iso, v = _isisotropic_with_vector(F)
  if !iso
    return false, F
  end
  vv = matrix(K, 1, length(v), v)

  # Find basis vector which has non-trivial product with v
  o = F * vv'
  i = findfirst(j -> !iszero(o[j, 1]), 1:nrows(o))
  @assert i isa Int
  H = vcat(vv, zero_matrix(base_ring(F), 1, ncols(F)))
  H[2, i] = inv(o[i, 1])
  GG = H * F * H'

  if !iszero(GG[2, 2])
    al = -GG[2, 2]//2
    for i in 1:ncols(H)
      H[2, i] = al * H[1, i] + H[2, i]
    end
    GG = H * F * H'
  end

  @assert iszero(GG[1, 1])
  @assert iszero(GG[2, 2])
  @assert isone(GG[1, 2])
  @assert isone(GG[2, 1])

  return true, H
end

function _find_direct_sum(B::MatrixElem)
  # I am very lazy
  C = B
  r = nrows(B)
  K = base_ring(B)
  piv = Int[]
  for i in 1:ncols(B)
    z = zero_matrix(K, 1, ncols(B))
    z[1, i] = one(K)
    CC = vcat(C, z)
    _r = rank(CC)
    if _r == r
      continue
    else
      r = _r
      push!(piv, i)
      C = CC
    end
  end
  l = length(piv)
  @assert nrows(B) + l == ncols(B)
  Z = zero_matrix(K, l, ncols(B))
  for i in 1:length(piv)
    Z[i, piv[i]] = one(K)
  end
  return Z
end

function _orthogonal_complement(F::MatrixElem, B::MatrixElem)
  r, R = left_kernel(F * transpose(B))
  @assert nrows(R) == r
  if !iszero(det(F))
    @assert nrows(R) + nrows(B) == nrows(F)
  end
  return R
end

function _maximal_isotropic_subspace(F::MatElem)
  _, H, R = _quadratic_form_decomposition(F)
  return vcat(sub(H, collect(1:2:nrows(H)), collect(1:ncols(H))), R)
end

function _isisometric_with_isometry(F, G)
  K = base_ring(F)

  if nrows(F) != nrows(G)
    return false, F
  end

  A1, H1, R1 = _quadratic_form_decomposition(F)
  A2, H2, R2 = _quadratic_form_decomposition(G)

  if nrows(H1) != nrows(H2) || nrows(R1) != nrows(R2)
    return false, F
  end

  _d1, _H1, _I1 = _quadratic_form_invariants(A1 * F * A1')
  _d2, _H2, _I2 = _quadratic_form_invariants(A2 * G * A2')

  if !(_I1 == _I2 && _H1 == _H2 && issquare(_d1 * _d2)[1])
    return false, F
  end

  @assert nrows(H1) == nrows(H2) && nrows(R1) == nrows(R2)

  X = zero_matrix(K, 0, ncols(F))
  Y = zero_matrix(K, 0, ncols(F))

  while nrows(A1) > 0
    GA1 = A1 * F * A1'
    GA2 = A2 * G * A2'
    ok, v = _isisotropic_with_vector(diagonal_matrix(GA1, -GA2))
    @assert ok
    n = div(length(v), 2)
    _x = matrix(K, 1, n, v[1:n])
    _y = matrix(K, 1, n, v[(n+1):2*n])
    @assert _x * GA1 * _x' == _y * GA2 *  _y'
    x = matrix(K, 1, n, v[1:n]) * A1
    y = matrix(K, 1, n, v[(n+1):2*n]) * A2
    @assert x * F * x' == y * G * y'
    X = vcat(X, x)
    Y = vcat(Y, y)
    @assert X * F * X' == Y * G * Y'
    _A1 = _orthogonal_complement(GA1, matrix(K, 1, n, v[1:n]))
    A1 = _A1 * A1
    _A2 = _orthogonal_complement(GA2, matrix(K, 1, n, v[(n + 1):2*n]))
    A2 = _A2 * A2
  end

  @assert X * F * X' == Y * G * Y'
  @assert H1 * F * H1' == H2  *  G * H2'

  M = inv(vcat(X, H1, R1)) * vcat(Y, H2, R2)

  @assert M * G * M' == F

  return true, M
end

function isisometric_with_isometry(V::QuadSpace, W::QuadSpace)
  GV = gram_matrix(V)
  GW = gram_matrix(W)
  return _isisometric_with_isometry(GV, GW)
end

function _real_weak_approximation(s, I)
  K = NumberField(s)
  a = gen(K)
  while true
    x = simplest_inside(real(evaluate(a, s, 10)))
    a = 2 * (a - x)
    if all(t -> t == s || abs(evaluate(a, t)) >= 2, I)
      break
    end
  end
  @assert abs(evaluate(a, s)) < 1//2
  return a
end

################################################################################
#
#  Is isotropic
#
################################################################################

function _isisotropic_with_vector_finite(M)
  n = ncols(M)
  k = base_ring(M)
  @assert k isa Field && characteristic(k) != 2
  if n == 0
    ;
  elseif n == 1
    if iszero(M[1, 1])
      return true, elem_type(k)[one(k)]
    end
  else
    if n <= 3
      G, T = _gram_schmidt(M, identity, false) # might be non-degenerate
    else
      G, T = _gram_schmidt(sub(M, 1:3, 1:3), identity, false) # might be non-degenerate
      B = zero_matrix(k, 3, n)
      B[1, 1] = 1
      B[2, 2] = 1
      B[3, 3] = 1
      T = T * B
    end
    for i in 1:ncols(G)
      if iszero(G[i, i])
        return true, elem_type(k)[T[i, j] for j in 1:ncols(T)]
      end
    end

    if n == 2
      ok, s = issquare_with_square_root(-divexact(G[1, 1], G[2, 2]))
      if ok
        return true, elem_type(k)[T[1, i] + s*T[2, i] for i in 1:ncols(T)]
      end
    else
      while true
        x = rand(k)
        y = rand(k)
        ok, z = issquare( -x^2 * G[1, 1] - y^2 * divexact(G[2, 2], G[3, 3]))
        if (ok && (!iszero(x) || !iszero(y)))
          return true,  elem_type(k)[x*T[1, i] + y*T[2, i] + z * T[3, i] for i in 1:ncols(T)]
        end
      end
    end
  end
  return false, elem_type(k)[]
end

ext: jl


text:
module AbstractNFFTs

using Graphics: @mustimplement
using LinearAlgebra
using Printf

# interface
export AbstractFTPlan, AbstractRealFTPlan, AbstractComplexFTPlan,
       AbstractNFFTPlan, AbstractNFCTPlan, AbstractNFSTPlan, AbstractNNFFTPlan,
       plan_nfft, plan_nfct, plan_nfst, mul!, size_in, size_out, nodes!

# optional
export deconvolve!, deconvolve_transpose!, convolve!, convolve_transpose!

# derived
export nfft, nfft_adjoint, ndft, ndft_adjoint,
       nfct, nfct_transpose, ndct, ndct_transpose,
       nfst, nfst_transpose

# misc
export TimingStats, accuracyParams, reltolToParams, paramsToReltol,
       PrecomputeFlags, LINEAR, FULL, TENSOR, POLYNOMIAL

include("misc.jl")
include("interface.jl")
include("derived.jl")

end


text:
"""
    CacheImage <: Augmentor.ImageOperation

Description
--------------

Write the current state of the image into the working memory.
Optionally a user has the option to specify a preallocated
`buffer` to write the image into. Note that if a `buffer` is
provided, then it has to be of the correct size and eltype.

Even without a preallocated `buffer` it can be beneficial in some
situations to cache the image. An example for such a scenario is
when chaining a number of affine transformations after an elastic
distortion, because performing that lazily requires nested
interpolation.

Usage
--------------

    CacheImage()

    CacheImage(buffer)

Arguments
--------------

- **`buffer`** : Optional. A preallocated `AbstractArray` of the
    appropriate size and eltype.

See also
--------------

[`augment`](@ref)

Examples
--------------

```julia
using Augmentor

# make pipeline that forces caching after elastic distortion
pl = ElasticDistortion(3,3) |> CacheImage() |> Rotate(-10:10) |> ShearX(-5:5)

# cache output of elastic distortion into the allocated
# 20x20 Matrix{Float64}. Note that for this case this assumes that
# the input image is also a 20x20 Matrix{Float64}
pl = ElasticDistortion(3,3) |> CacheImage(zeros(20,20)) |> Rotate(-10:10)

# convenience syntax with the same effect as above.
pl = ElasticDistortion(3,3) |> zeros(20,20) |> Rotate(-10:10)
```
"""
struct CacheImage <: ImageOperation end

applyeager(op::CacheImage, img::AbstractArray, param) = maybe_copy(img)

function showconstruction(io::IO, op::CacheImage)
    print(io, typeof(op).name.name, "()")
end

function Base.show(io::IO, op::CacheImage)
    if get(io, :compact, false)
        print(io, "Cache into temporary buffer")
    else
        print(io, "Augmentor.")
        showconstruction(io, op)
    end
end

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

"""
    CacheImageInto <: Augmentor.ImageOperation

see [`CacheImage`](@ref)
"""
struct CacheImageInto{T<:Union{AbstractArray,Tuple}} <: ImageOperation
    buffer::T
end
CacheImage(buffer::AbstractArray) = CacheImageInto(buffer)
CacheImage(buffers::AbstractArray...) = CacheImageInto(buffers)
CacheImage(buffers::NTuple{N,AbstractArray}) where {N} = CacheImageInto(buffers)

@inline supports_lazy(::Type{<:CacheImageInto}) = true

applyeager(op::CacheImageInto, img::AbstractArray, param) = applylazy(op, img)
applyeager(op::CacheImageInto, img::Tuple) = applylazy(op, img)

function applylazy(op::CacheImageInto, img::Tuple)
    throw(ArgumentError("Operation $(op) not compatiable with given image(s) ($(summary(img))). This can happen if the amount of images does not match the amount of buffers in the operation"))
end

function applylazy(op::CacheImageInto{<:AbstractArray}, img::AbstractArray, param)
    copy!(match_idx(op.buffer, axes(img)), img)
end

function applylazy(op::CacheImageInto{<:Tuple}, imgs::Tuple)
    map(op.buffer, imgs) do buffer, img
        copy!(match_idx(buffer, axes(img)), img)
    end
end

function _showconstruction(io::IO, array::AbstractArray)
    print(io, "Array{")
    _showcolor(io, eltype(array))
    print(io, "}(")
    print(io, join(map(i->string(length(i)), axes(array)), ", "))
    print(io, ")")
end

function showconstruction(io::IO, op::CacheImageInto{<:AbstractArray})
    print(io, "CacheImage(") # shows exported API
    _showconstruction(io, op.buffer)
    print(io, ")")
end

function showconstruction(io::IO, op::CacheImageInto{<:Tuple})
    print(io, "CacheImage(")
    for (i, buffer) in enumerate(op.buffer)
        _showconstruction(io, buffer)
        i < length(op.buffer) && print(io, ", ")
    end
    print(io, ")")
end

function Base.show(io::IO, op::CacheImageInto{<:AbstractArray})
    if get(io, :compact, false)
        print(io, "Cache into preallocated ")
        print(io, summary(op.buffer))
    else
        print(io, typeof(op).name, "(")
        Base.showarg(io, op.buffer, false)
        print(io, ")")
    end
end

function Base.show(io::IO, op::CacheImageInto{<:Tuple})
    if get(io, :compact, false)
        print(io, "Cache into preallocated ")
        print(io, "(", join(map(summary, op.buffer), ", "), ")")
    else
        print(io, typeof(op).name, "((")
        for (i, buffer) in enumerate(op.buffer)
            Base.showarg(io, buffer, false)
            i < length(op.buffer) && print(io, ", ")
        end
        print(io, "))")
    end
end


text:
using Dates
LOCK = 1
EMPT = 0
BOARD_HEIGHT = 9
BOARD_WIDTH = 6
FILE_NAME = "julia-logfile.txt"
bk1 = [
    LOCK EMPT EMPT
    LOCK LOCK LOCK
    LOCK EMPT EMPT
    LOCK EMPT EMPT
]
bk2 = [
    LOCK EMPT
    LOCK LOCK
    LOCK EMPT
    LOCK EMPT
    LOCK EMPT
]
bk3 = [
    LOCK LOCK
    LOCK EMPT
    LOCK EMPT
    LOCK EMPT
]
bk4 = [
    LOCK LOCK
    LOCK EMPT
    LOCK EMPT
    LOCK EMPT
    LOCK EMPT
]
bk5 = [
    LOCK LOCK LOCK
    LOCK EMPT EMPT
    LOCK EMPT EMPT
    LOCK EMPT EMPT
]
bk6 = [
    LOCK EMPT
    LOCK LOCK
    EMPT LOCK
    EMPT LOCK
]
bk7 = [
    LOCK EMPT LOCK
    LOCK LOCK LOCK
    LOCK EMPT EMPT
]
bk8 = [
    LOCK EMPT EMPT
    LOCK LOCK LOCK
    EMPT LOCK EMPT
    EMPT LOCK EMPT
]
bk9 = [
    LOCK EMPT
    LOCK LOCK
    LOCK LOCK
]

force_check = 0

function markDateBoard()
    date = Dates.now()
    nboard = zeros(Int8, BOARD_HEIGHT, BOARD_WIDTH)
    wdpos = [[8,4] ,[8,5],[8,6],[9,3] ,[9,4],[9,5] ,[9,6]]
    if Dates.month(date) > BOARD_WIDTH
        nboard[2,Dates.month(date) - BOARD_WIDTH] = LOCK
    else
        nboard[1,Dates.month(date)] = LOCK
    end
    if Dates.month(date) == 1
        global force_check = 1
    end
    r, c = divrem(Dates.day(date), BOARD_WIDTH)
    nboard[r + 3,c] = LOCK
    i, j = wdpos[Dates.dayofweek(date)]
    nboard[i,j] = LOCK
    return nboard
end

function pos(bk,x,y)
    r,c = size(bk)
    if x > 0
        bk = vcat(zeros(Int8,x,c),bk)
    end
    bk = vcat(bk,zeros(Int8,BOARD_HEIGHT - r - x,c))
    r,c = size(bk)
    if y > 0
        bk = hcat(zeros(Int8,r,y),bk)
    end
    bk = hcat(bk,zeros(Int8,r,BOARD_WIDTH - c - y))
    return bk
end

function rotate(bk, times = 0)
    if times < 0
        bk = transpose(bk)
        times += 1
    end
    return rotl90(bk, times)
end

function validBoard(currBoard)
    if maximum(currBoard) == LOCK
        fullfill = filter(x -> x == LOCK,currBoard)
        if length(fullfill) > 21
            fillrows = fdiv(length(fullfill), 8)
            return all(currBoard[1:fillrows,:] .== LOCK)
        else
            global force_check
            return currBoard[1,force_check] == LOCK
        end
    end
    return false
end

function placeBlock(currBoard, blk)
    for i in 0:BOARD_HEIGHT
        for j in 0:BOARD_WIDTH
            row,column = size(blk)
            if row + i < 10 && column + j < 7
                posBlk = pos(blk,i,j)
                curr = (+).(currBoard,posBlk)
                if validBoard(curr)
                    run(`clear`)
                    display(curr)
                    return true,curr,posBlk
                end
            end
        end
    end
    return false,currBoard,blk
end

function enumFill(currBoard,leftBlk,placedBlk)
    curr = currBoard
    lefts = leftBlk
    placed = placedBlk
    if isempty(lefts)
        return true,curr,lefts,placed
    end
    for i in 1:ndims(lefts)
        for j in -5:5
            chk1,curr1,pos1 = placeBlock(curr,rotate(lefts[i],j))
            if chk1
                lefts1 = copy(lefts)
                deleteat!(lefts1,i)
                placed1 = cat(placed,pos1;dims=3)
                chk2,curr2,lefts2,placed2 = enumFill(curr1,lefts1,placed1)
                if chk2
                    return true,curr2,lefts2,placed2
                end
            end
        end
    end
    return false,curr,lefts,placed
end
function main()
    blocks = [bk3,bk2,bk4,bk6,bk1,bk5,bk7,bk8,bk9]
    board = markDateBoard()
    placedBlocks = []
    chk,curr,lefts,placed = enumFill(board,blocks,placedBlocks)
    if chk
        run(`clear`)
        display(placed)
    end
end

main()
# display(rotate(bk1, -1))
# test = placeBlock(board,rotate(bk1, -1))
# display(test)

max_stars_repo_path: kata-julia/A Puzzle A Day.jl

max_stars_repo_name: crimsonskyrem/codewars

max_issues_repo_path: kata-julia/A Puzzle A Day.jl

max_issues_repo_name: crimsonskyrem/codewars

max_forks_repo_path: kata-julia/A Puzzle A Day.jl

max_forks_repo_name: crimsonskyrem/codewars


text:
module TestMeta
    using Test
    using CBSOData3

    include("mock_getodata.jl")

    meta_80137ned = CBSOData3.get_meta("80137ned")
    @test "DataProperties" in keys(meta_80137ned)
    @test Set(["TableInfos", "CategoryGroups", "BedrijfstakkenSBI93", "Persoonskenmerken", "DataProperties", "Perioden"]) == Set(keys(meta_80137ned))

end


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

mutable struct ConnectionMonitorTcpConfiguration <: SwaggerModel
    port::Any # spec type: Union{ Nothing, Int32 } # spec name: port
    disableTraceRoute::Any # spec type: Union{ Nothing, Bool } # spec name: disableTraceRoute

    function ConnectionMonitorTcpConfiguration(;port=nothing, disableTraceRoute=nothing)
        o = new()
        validate_property(ConnectionMonitorTcpConfiguration, Symbol("port"), port)
        setfield!(o, Symbol("port"), port)
        validate_property(ConnectionMonitorTcpConfiguration, Symbol("disableTraceRoute"), disableTraceRoute)
        setfield!(o, Symbol("disableTraceRoute"), disableTraceRoute)
        o
    end
end # type ConnectionMonitorTcpConfiguration

const _property_map_ConnectionMonitorTcpConfiguration = Dict{Symbol,Symbol}(Symbol("port")=>Symbol("port"), Symbol("disableTraceRoute")=>Symbol("disableTraceRoute"))
const _property_types_ConnectionMonitorTcpConfiguration = Dict{Symbol,String}(Symbol("port")=>"Int32", Symbol("disableTraceRoute")=>"Bool")
Base.propertynames(::Type{ ConnectionMonitorTcpConfiguration }) = collect(keys(_property_map_ConnectionMonitorTcpConfiguration))
Swagger.property_type(::Type{ ConnectionMonitorTcpConfiguration }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ConnectionMonitorTcpConfiguration[name]))}
Swagger.field_name(::Type{ ConnectionMonitorTcpConfiguration }, property_name::Symbol) =  _property_map_ConnectionMonitorTcpConfiguration[property_name]

function check_required(o::ConnectionMonitorTcpConfiguration)
    true
end

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


text:
@testset "Protting" begin
  @testset "DistPlot1D" begin
    # TODO
  end

  @testset "DistPlot2D" begin
    # TODO
  end

  @testset "CornerPlot" begin
    # TODO
  end
end


text:
using Test
using WaveProp
using WaveProp.Geometry
using WaveProp.Integration
using WaveProp.Mesh

@testset "Circle" begin
    geo = Circle()
    bnd = boundary(geo)
    @test length(bnd) == 1
    Ω    = Domain([geo])
    mesh =  Mesh.meshgen(Ω;h=0.1)
    @test ambient_dimension(mesh) == 2
end

max_stars_repo_name: WaveProp/WaveProp

max_issues_repo_name: WaveProp/WaveProp

max_forks_repo_name: WaveProp/WaveProp


text:
function circ(h::Int64,H::UnitRange{Int64})::Int64
    (h==minimum(H)) ? idx = maximum(H) : idx = h-1
    return idx
end

function RMP0!(m::JuMP.Model,
              ms::ms_struct,
              mp::mp_struct,
              hc::hc_struct)::JuMP.Model

    # */ ------------ variables ------------ /* #
    @variable(m, f )
    @variable(m, xi[ms.P,ms.Ii] >= .0 )
    @variable(m, xo[1:hc.nx,ms.Io])
    # */ ----------------------------------- /* #

    # */ ----------- obj function ---------- /* #
    @objective(m, Min, f )
    # */ ----------------------------------- /* #

    # */ ----------- constraints ----------- /* #
    @constraint(m, c1, f >= sum(mp.ci[p,ii]*xi[p,ii] for p in ms.P for ii in ms.Ii) + sum(mp.cf[p]*xo[p,io] for p in ms.P for io in ms.Io))
    @constraint(m, c2[p=ms.P,io=ms.Io], xo[p,io] == mp.xh[io,p] + sum(xi[p,ii] for ii in mp.map[p,io]) )
    @constraint(m, c3[p=ms.P,io=ms.Io], xo[p,io] >= .0 )
    @constraint(m, c4[p=ms.P,io=ms.Io], xo[p,io] <= mp.xm[p] )
    @constraint(m, c5[       io=ms.Io], xo[hc.nx0+1,io] == hc.h[io,1] )
    @constraint(m, c6[       io=ms.Io], xo[hc.nx0+2,io] == hc.h[io,2] )
    # */ ----------------------------------- /* #

    return m
end

function LM!(m::JuMP.Model,
            ms::ms_struct,
            mp::mp_struct,
            hc::hc_struct)::JuMP.Model

    RMP0!(m,ms,mp,hc)
    # */ ------------ variables ------------ /* #
    @variable(m, β[ms.Io] >= 0., container=Array)
    @variable(m, ν[ms.Io] >= 0., container=Array)
    @variable(m, ϕ[ms.Io] >= 0., container=Array)
    # */ ----------- obj function ---------- /* #
    @objective(m, Min, objective_function(m) + sum(β[io] for io in ms.Io))
    # */ --------- constraints ---------- /* #
    @constraint(m, cl01[io=ms.Io], ν[io] + ϕ[io]*hc.c[io,1] <= β[io] )

    return m
end

function UM!(m::JuMP.Model,
            ms::ms_struct,
            mp::mp_struct,
            hc::hc_struct)::JuMP.Model

    RMP0!(m,ms,mp,hc)
    # */ ------------ variables ------------ /* #
    @variable(m, β[ms.Io] >= 0. )
    # */ ----------- obj function ---------- /* #
    @objective(m, Min, objective_function(m) + sum(β[io] for io in ms.Io))
    # */ --------- constraints ---------- /* #
    @constraint(m, cu01[io=ms.Io], 0. <= β[io] )
    @constraint(m, cu02[io=ms.Io,ix=1:hc.nx], 0. <= m[:xo][ix,io] )
    @constraint(m, cu03[io=ms.Io], 0. == 1. )

    return m
end

# */ ---------------------------------------------------------------------- /* #

function SP!(m::JuMP.Model,
            ps::ps_struct,
            pp::pp_struct,
            hc::hc_struct)::JuMP.Model

    # */ ------------ variables ------------ /* #
    @variable(m, ϕ[1:hc.nc] )
    @variable(m, c0 )
    @variable(m, yG[ps.G,ps.S,ps.H] >= 0  )
    @variable(m, yI[ps.B,ps.S,ps.H] >= 0  )
    @variable(m, yO[ps.B,ps.S,ps.H] >= 0  )
    @variable(m, yL[ps.B,ps.S,ps.H] >= 0  )
    @variable(m, yR[ps.R,ps.S,ps.H] >= 0  )
    @variable(m, yS[     ps.S,ps.H] >= 0  )
    @variable(m, x[1:hc.nx] )
    # */ ----------- obj function ---------- /* #
    @objective(m, Min, c0 )
    # */ --------- constraints ---------- /* #
    @constraint(m, c01, c0 == sum(sum(sum((pp.cvg[g]+pp.cfg[g]/pp.ηg[g])*yG[g,s,h] for g in ps.G) + pp.cs*yS[s,h] for h in ps.H) for s in ps.S))
    @constraint(m, c02, ϕ[1] == sum(sum(sum(pp.eg[g]/pp.ηg[g]*yG[g,s,h] for g in ps.G) for h in ps.H) for s in ps.S))
    @constraint(m, c03[b=ps.B,s=ps.S,h=ps.H], yL[b,s,h] - yL[b,s,circ(h,ps.H)] == pp.ηb[b]*yI[b,s,h] - yO[b,s,h])
    @constraint(m, c04[g=ps.G,s=ps.S,h=ps.H], yG[g,s,h] - yG[g,s,circ(h,ps.H)] <= pp.rg[g]*x[g])
    @constraint(m, c05[g=ps.G,s=ps.S,h=ps.H], yG[g,s,circ(h,ps.H)] - yG[g,s,h] <= pp.rg[g]*x[g])
    @constraint(m, c06[g=ps.G,s=ps.S,h=ps.H], yG[g,s,h] <= x[g])
    @constraint(m, c07[b=ps.B,s=ps.S,h=ps.H], yI[b,s,h] <= pp.pb[b]*x[ps.b0+b])
    @constraint(m, c08[b=ps.B,s=ps.S,h=ps.H], yO[b,s,h] <= pp.pb[b]*x[ps.b0+b])
    @constraint(m, c09[b=ps.B,s=ps.S,h=ps.H], yL[b,s,h] <= x[ps.b0+b])
    @constraint(m, c10[r=ps.R,s=ps.S,h=ps.H], yR[r,s,h] <= pp.pr[r,s,h]*x[ps.r0+r])
    @constraint(m, c11[s=ps.S,h=ps.H], sum(yG[g,s,h] for g in ps.G) + sum(yO[b,s,h]-yI[b,s,h] for b in ps.B) + sum(yR[r,s,h] for r in ps.R) + yS[s,h] == -x[hc.nx0+1]*pp.pd[s,h] )
    @constraint(m, c12, ϕ[1] <= x[hc.nx0+2]*1.0e6 )

    return m
end

function create_SPo(ps::ps_struct,
                    pp::pp_struct,
                    hc::hc_struct)::SPs

    gurobi_env = @suppress Gurobi.Env()
    m = JuMP.Model(JuMP.with_optimizer(Gurobi.Optimizer,gurobi_env,OutputFlag=0,Method=2))
    SP!(m,ps,pp,hc)
    x,λ = zeros(hc.nx),zeros(hc.nx)
    c,ϕ = zeros(hc.nc),zeros(hc.nc)
    θ   = 0.
    return SPs(m,x,c,θ,λ,ϕ)

end

function update_and_solve_SPo!(spo::SPs,
                                 x::Array{Float64,1},
                                 c::Array{Float64,1})::SPs

    spo.x .= x
    spo.c .= c
    fix.(spo.m[:x],spo.x;force=true)
    @objective(spo.m, Min, spo.m[:c0] + sum(spo.c[ic]*spo.m[:ϕ][ic] for ic in 1:length(spo.c)))
    optimize!(spo.m)
    spo.θ = objective_value(spo.m)
    spo.λ .= dual.(FixRef.(spo.m[:x]))
    spo.ϕ .= value.(spo.m[:ϕ])

    return spo

end

# */ ---------------------------------------------------------------------- /* #

function create_RMPo(ms::ms_struct,
                     mp::mp_struct,
                     ps::ps_struct,
                     pp::pp_struct,
                     hc::hc_struct,
                     al::al_params)::RMPs

    gurobi_env = @suppress Gurobi.Env()
    rmL = JuMP.Model(JuMP.with_optimizer(Gurobi.Optimizer,gurobi_env,OutputFlag=0))
    rmU = JuMP.Model(JuMP.with_optimizer(Gurobi.Optimizer,gurobi_env,OutputFlag=0))
    rmB = JuMP.Model(JuMP.with_optimizer(Gurobi.Optimizer,gurobi_env,OutputFlag=0))
    LM!(rmL,ms,mp,hc)
    UM!(rmU,ms,mp,hc)
    UM!(rmB,ms,mp,hc)
    rmLo = MLo(rmL,Array{ConstraintRef,2}(undef,hc.ny,al.J+1))
    rmUo = MUo(rmU,Array{VariableRef  ,2}(undef,hc.ny,al.J+1))
    rmBo = MUo(rmB,Array{VariableRef  ,2}(undef,hc.ny,al.J+1))
    sp = create_SPo(ps,pp,hc)
    xm,cm = vcat(minimum(mp.xh,dims=1)[:],minimum(hc.h,dims=1)[:]),minimum(hc.c,dims=1)[:]
    ls = LSs(zeros(ms.Io[end],ps.G[end]),zeros(ms.Io[end],ps.B[end]),zeros(ms.Io[end],ps.B[end]),zeros(ms.Io[end],ps.R[end]),zeros(ms.Io[end]),zeros(ms.Io[end]))
    return RMPs(rmLo,rmUo,rmBo,sp,xm,cm,hc.nx,hc.nc,zeros(hc.nx,hc.ny),hc.h,hc.c,ms.Io,zeros(hc.ny),zeros(hc.ny),-Inf,Inf,Inf,zeros(0),zeros(0),zeros(0),0,al.ϵ,al.J,ls)

end

function build_cuts_RMP!(rmp::RMPs)::RMPs

    j = rmp.it+1
    for i in rmp.Io
        rmp.mL.c[i,j] = @constraint(rmp.mL.m, rmp.mL.m[:ν][i] + sum(rmp.mL.m[:ϕ][i,ic]*rmp.sp.c[ic] for ic in 1:rmp.nc) >= rmp.sp.θ + sum(rmp.sp.λ[ix]*(rmp.mL.m[:xo][ix,i]-rmp.sp.x[ix]) for ix in 1:rmp.nx) )
        rmp.mU.v[i,j] = @variable(rmp.mU.m,base_name="μ[$(i),$(j)]",lower_bound=0.)
        set_normalized_coefficient(rmp.mU.m[:cu01][i],rmp.mU.v[i,j],rmp.sp.θ + sum(rmp.sp.ϕ[ic]*(rmp.c[i,ic]-rmp.sp.c[ic]) for ic in 1:rmp.nc))
        for ix in 1:rmp.nx set_normalized_coefficient(rmp.mU.m[:cu02][i,ix],rmp.mU.v[i,j],rmp.sp.x[ix]) end
        set_normalized_coefficient(rmp.mU.m[:cu03][i],rmp.mU.v[i,j],1.)
        rmp.mB.v[i,j] = @variable(rmp.mB.m,base_name="μ[$(i),$(j)]",lower_bound=0.)
        set_normalized_coefficient(rmp.mB.m[:cu01][i],rmp.mB.v[i,j],rmp.sp.θ + sum(rmp.sp.ϕ[ic]*(rmp.c[i,ic]-rmp.sp.c[ic]) for ic in 1:rmp.nc))
        for ix in 1:rmp.nx set_normalized_coefficient(rmp.mB.m[:cu02][i,ix],rmp.mB.v[i,j],rmp.sp.x[ix]) end
        set_normalized_coefficient(rmp.mB.m[:cu03][i],rmp.mB.v[i,j],1.)
    end
    return rmp

end

function do_step0_RMP!(rmp::RMPs)::RMPs

    println("")
    println(" PlanModel algorithm started \r")
    println("")
    update_and_solve_SPo!(rmp.sp,rmp.xm,rmp.cm)
    build_cuts_RMP!(rmp)
    return rmp

end

function do_step_RMP!(rmp::RMPs)::RMPs

    rmp.it += 1
    optimize!(rmp.mL.m)
    rmp.x .= value.(rmp.mL.m[:xo])
    rmp.Lb = objective_value(rmp.mL.m)
    for ix in 1:rmp.nx for iy in rmp.Io fix(rmp.mB.m[:xo][ix,iy],rmp.x[ix,iy];force=true) end end
    optimize!(rmp.mB.m)
    rmp.βL .= value.(rmp.mL.m[:β])
    rmp.βU .= value.(rmp.mB.m[:β])
    i0 = findmax(rmp.βU-rmp.βL)[2]
    update_and_solve_SPo!(rmp.sp,rmp.x[:,i0],rmp.c[i0,:])
    build_cuts_RMP!(rmp)
    optimize!(rmp.mU.m)
    rmp.Ub = objective_value(rmp.mU.m)
    rmp.δ = (rmp.Ub-rmp.Lb)/rmp.Ub*100
    println(" j = $(rmp.it), Δ = $(@sprintf("%.2e",rmp.δ))% \r")
    push!(rmp.LB,rmp.Lb)
    push!(rmp.UB,rmp.Ub)
    push!(rmp.Δ,rmp.δ)
    return rmp

end

function solve_RMP!(rmp::RMPs)::RMPs

    do_step0_RMP!(rmp)
    for j in 1:rmp.J
        do_step_RMP!(rmp)
        (rmp.δ <= rmp.ϵ) ? break : nothing
    end
    return rmp

end

function last_step_RMP!(rmp::RMPs,
                         hc::hc_struct,
                         ps::ps_struct)::RMPs

     xo,co =zeros(hc.nx),zeros(hc.nc)
     for io in rmp.Io
         xo .= value.(rmp.mU.m[:xo][:,io])
         co .= rmp.c[io]
         update_and_solve_SPo!(rmp.sp,xo,co)
         for g in ps.G
             rmp.ls.G[io,g] = sum(value.(rmp.sp.m[:yG][g,:,:]))/1.0e6
         end
         for b in ps.B
             rmp.ls.Bi[io,b] = sum(value.(rmp.sp.m[:yI][b,:,:]))/1.0e6
             rmp.ls.Bo[io,b] = sum(value.(rmp.sp.m[:yO][b,:,:]))/1.0e6
         end
         for r in ps.R
             rmp.ls.R[io,r] = sum(value.(rmp.sp.m[:yR][r,:,:]))/1.0e6
         end
         rmp.ls.E[io] = value(rmp.sp.m[:ϕ][1])/1.0e6
         rmp.ls.C[io] = -dual(rmp.sp.m[:c12])
     end
    return rmp

end

# */ ---------------------------------------------------------------------- /* #

function write_output(rmp::RMPs,
                       ms::ms_struct,
                       ps::ps_struct,
                       sd::sd_struct)

    XLSX.openxlsx("$(pwd())/output/results.xlsx", mode="w") do xf

        sheet = xf[1]
        XLSX.rename!(sheet, "Investment")
        sheet["A1"] = "Investment in new capacity"
        sheet["A3"] = "Years"
        sheet["A5", dim=1] = sd.yrs
        sheet["B3", dim=2] = sd.techs
        sheet["B4", dim=2] = repeat(["(MW)"],length(sd.techs))
        xi = zeros(ms.P[end],ms.Ii[end])
        xi .= value.(rmp.mU.m[:xi])
        iO = (sd.bi.*(1:length(sd.bi)))[Bool.(sd.bi)]
        Xi = zeros(ms.Io[end],ms.P[end])
        for ii in ms.Ii Xi[iO[ii],:] .= xi[:,ii] end
        sheet["B5"] = Xi
        sheet[3,ms.P[end]+3] = "Total_cost"
        sheet[4,ms.P[end]+3] = "(10^9 £)"
        sheet[5,ms.P[end]+3] = objective_value(rmp.mU.m)*1.0e-9
        sheet[3,ms.P[end]+4] = "Investment_cost"
        sheet[4,ms.P[end]+4] = "(10^9 £)"
        sheet[5,ms.P[end]+4] = sum(mp.ci[p,ii]*xi[p,ii] for p in ms.P for ii in ms.Ii)*1.0e-9

        XLSX.addsheet!(xf)
        sheet = xf[2]
        XLSX.rename!(sheet, "Operation")
        sheet["A1"] = "Total yearly available capacity"
        sheet["A3"] = "Years"
        sheet["A5", dim=1] = sd.yrs
        sheet["B3", dim=2] = sd.techs
        sheet["B4", dim=2] = repeat(["(MW)"],length(sd.techs))
        Xo = zeros(ms.Io[end],ms.P[end])
        for p in ms.P for io in ms.Io Xo[io,p] = value(rmp.mU.m[:xo][p,io])  end end
        sheet["B5"] = Xo
        sheet[3,ms.P[end]+3] = "FixO&M"
        sheet[4,ms.P[end]+3] = "(10^9 £)"
        sheet[5,ms.P[end]+3] = sum(mp.cf[p]*Xo[io,p] for p in ms.P for io in ms.Io)*1.0e-9
        sheet[3,ms.P[end]+4] = "Operation cost"
        sheet[4,ms.P[end]+4] = "(10^9 £)"
        sheet[5,ms.P[end]+4] = sum(value.(rmp.mU.m[:β]))*1.0e-9

        XLSX.addsheet!(xf)
        sheet = xf[3]
        XLSX.rename!(sheet, "Usage_Conventional")
        sheet["A1"] = "Total yearly usage of conventional technologies"
        sheet["A3"] = "Years"
        sheet["A5", dim=1] = sd.yrs
        sheet["B3", dim=2] = sd.techs[ps.G]
        sheet["B4", dim=2] = repeat(["(TWh)"],ps.G[end])
        sheet["B5"] = RMP.ls.G

        XLSX.addsheet!(xf)
        sheet = xf[4]
        XLSX.rename!(sheet, "Usage_Storage")
        sheet["A1"] = "Total yearly usage of storage technologies"
        sheet["A3"] = "Years"
        sheet["A5", dim=1] = sd.yrs
        sheet["B3", dim=2] = repeat(sd.techs[ps.b0.+ps.B],2) .* vcat(repeat([" (charge)"],ps.B[end]),repeat([" (discharge)"],ps.B[end]))
        sheet["B4", dim=2] = repeat(["(TWh)"],2*ps.B[end])
        sheet[XLSX.encode_column_number(2).*"5"]           = RMP.ls.Bi
        sheet[XLSX.encode_column_number(2+ps.B[end]).*"5"] = RMP.ls.Bo

        XLSX.addsheet!(xf)
        sheet = xf[5]
        XLSX.rename!(sheet, "Usage_Renewable")
        sheet["A1"] = "Total yearly usage of renewable technologies"
        sheet["A3"] = "Years"
        sheet["A5", dim=1] = sd.yrs
        sheet["B3", dim=2] = sd.techs[ps.r0.+ps.R]
        sheet["B4", dim=2] = repeat(["(TWh)"],ps.R[end])
        sheet["B5"] = RMP.ls.R

        XLSX.addsheet!(xf)
        sheet = xf[6]
        XLSX.rename!(sheet, "Emissions")
        sheet["A1"] = "Total yearly CO2 emissions"
        sheet["A3"] = "Years"
        sheet["A5", dim=1] = sd.yrs
        sheet["B3"] = "CO2 emissions"
        sheet["B4"] = "(MtCO2)"
        sheet["B5", dim=1] = RMP.ls.E
        sheet["C3"] = "additional CO2 price"
        sheet["C4"] = "(£/tCO2)"
        sheet["C5", dim=1] = RMP.ls.C

    end

end

# */ ---------------------------------------------------------------------- /* #

max_stars_repo_name: nimazzi/PlanModelDet

max_issues_repo_name: nimazzi/PlanModelDet

max_forks_repo_name: nimazzi/PlanModelDet


text:
if basename(pwd()) == "aoc"
    cd("2016/1")
end

function part1(filename::AbstractString)
    position = CartesianIndex(0, 0)
    dirs = (CartesianIndex(1, 0), CartesianIndex(0, 1), CartesianIndex(-1, 0), CartesianIndex(0, -1))
    i = 1
    for m in eachmatch(r"[LR]|\d+", readchomp(filename))
        if m.match == "L"
            i = isone(i) ? 4 : i - 1
        elseif m.match == "R"
            i = i % 4 + 1
        else
            position += dirs[i] * parse(Int, m.match)
        end
    end
    abs(position[1]) + abs(position[2])
end

@assert part1("example.txt") == 12
part1("input.txt")

function part2(filename::AbstractString)
    position = CartesianIndex(0, 0)
    visited = Set((position,))
    dirs = (CartesianIndex(1, 0), CartesianIndex(0, 1), CartesianIndex(-1, 0), CartesianIndex(0, -1))
    i = 1
    for m in eachmatch(r"[LR]|\d+", readchomp(filename))
        if m.match == "L"
            i = isone(i) ? 4 : i - 1
        elseif m.match == "R"
            i = i % 4 + 1
        else
            for _ in 1:parse(Int, m.match)
                position += dirs[i]
                if position ∈ visited
                    return abs(position[1]) + abs(position[2])
                else
                    push!(visited, position)
                end
            end
        end
    end
    throw(ArgumentError("Didn't loop"))
end

@assert part2("example2.txt") == 4
part2("input.txt")

max_stars_repo_name: stellartux/advent-of-code

max_issues_repo_name: stellartux/advent-of-code

max_forks_repo_name: stellartux/advent-of-code


text:
export eachtraversal
using LinearAlgebra

"""
    hits = eachtraversal(ray, edges)

Compute the traversal of `ray` through a grid that is the produce of `edges`.

# Example

```jldoctest
julia> using VoxelRayTracers

julia> edges = (0:1:10, 4:2:10,);

julia> ray = (position=[0,0], velocity=[1,1]);

julia> for hit in eachtraversal(ray, edges)
           println(hit)
       end
(voxelindex = CartesianIndex(5, 1), entry_time = 4.0, exit_time = 5.0)
(voxelindex = CartesianIndex(6, 1), entry_time = 5.0, exit_time = 6.0)
(voxelindex = CartesianIndex(7, 2), entry_time = 6.0, exit_time = 7.0)
(voxelindex = CartesianIndex(8, 2), entry_time = 7.0, exit_time = 8.0)
(voxelindex = CartesianIndex(9, 3), entry_time = 8.0, exit_time = 9.0)
(voxelindex = CartesianIndex(10, 3), entry_time = 9.0, exit_time = 10.0)
```
"""
function eachtraversal end

firstlast(xs) = (first(xs), last(xs))

function _start_voxelindex(pos, edges)
    map(pos, edges) do pos, walls
        voxelindex = searchsortedlast(walls, pos)
        clamp(voxelindex, firstindex(walls), lastindex(walls) - 1)
    end
end

function interval_entry_exit_time(pos, vel, walls)
    x_left, x_right = firstlast(walls)
    invvel = 1/vel
    if vel == 0
        T = typeof((x_left  - pos) *invvel)
        isinside = x_left <= pos <= x_right
        if isinside
            T(-Inf), T(Inf)
        else
            T(Inf), T(-Inf)
        end
    else
        t_left  = (x_left  - pos) * invvel
        t_right = (x_right - pos) * invvel
        minmax(t_left, t_right)
    end
end

function entry_exit_time(pos, vel, limits)
    ts = map(interval_entry_exit_time, pos, vel, limits)::Tuple
    entry_time = maximum(first, ts)
    exit_time = minimum(last, ts)
    entry_time, exit_time
end

struct WallHit{N,T,V}
    position::NTuple{N,T}
    index::NTuple{N,Int}
    hittimes::NTuple{N,T}
    time::T
    velocity::NTuple{N,T}
    voxel_traversal_times::V
end

abstract type AbstractTraversalTime end

struct TraversalTime_Range{T} <: AbstractTraversalTime
    t::T
end
Base.getindex(o::TraversalTime_Range, i) = o.t
has_finite_travel_time(o::TraversalTime_Range) = isfinite(o.t)

struct TraversalTime_Vector{T,V} <: AbstractTraversalTime
    walls::V
    inv_velocity::T
end
Base.@propagate_inbounds function Base.getindex(o::TraversalTime_Vector, i)
    lo = get(o.walls, i , first(o.walls))
    hi = get(o.walls, i+1, last(o.walls))
    o.inv_velocity * (hi - lo)
end
has_finite_travel_time(o::TraversalTime_Vector) = isfinite(o.inv_velocity)

function inv_abs(x)
    T = typeof(x)
    ifelse(iszero(x), T(Inf), inv(abs(x)))
end
function AbstractTraversalTime(walls, velocity)
    TraversalTime_Vector(walls, inv_abs(velocity))
end
function AbstractTraversalTime(walls::AbstractRange, velocity)
    TraversalTime_Range(step(walls)*inv_abs(velocity))
end

function WallHit(;
    position,
    index,
    hittimes,
    time,
    velocity,
    voxel_traversal_times,
    )
    WallHit(
        position,
        index,
        hittimes,
        time,
        velocity,
        voxel_traversal_times,
    )
end
function WallHit(o::WallHit; position, time, index, hittimes, )
    WallHit(position, index, hittimes, time, o.velocity, o.voxel_traversal_times)
end
function Base.show(io::IO, o::WallHit)
    show(io, typeof(o))
    println(io, "(;")
    for pname in propertynames(o)
        val = getproperty(o, pname)
        println(io, "   $pname = $val,")
    end
    println(io, ")")
end

int_sign(x)::Int = Int(sign(x))
function next_wallhit(o::WallHit)
    @inbounds begin
        dim          = argmin(o.hittimes)
        time_new     = o.hittimes[dim]
        position_new = let Δt = time_new-o.time
            map(o.velocity, o.position) do vel, pos
                vel*Δt + pos
            end
        end
        s = int_sign(o.velocity[dim])
        index_dim    = o.index[dim]+s
        index_new    = Base.setindex(o.index, index_dim, dim)
        hittime_new  = o.hittimes[dim] + o.voxel_traversal_times[dim][index_dim]
        hittimes_new = Base.setindex(o.hittimes, hittime_new, dim)
        WallHit(o,
            position=position_new,
            time=time_new,
            index=index_new,
            hittimes=hittimes_new,
        )
    end
end

function enter(focus, velocity, edges)
    entry_time, exit_time = entry_exit_time(focus, velocity, edges)
    if exit_time >= 0
        time = max(entry_time,zero(entry_time))
        position = let time=time
            map(focus, velocity) do foc, vel
                foc + time * vel
            end
        end
        index = _start_voxelindex(position, edges)
        return (;position, index, time)
    else
        return nothing
    end
end

function grid_entry_time(focus, velocity, edges) # this might not be on a wall
    entry_time, exit_time = entry_exit_time(focus, velocity, edges)
    max.(entry_time, zero(entry_time))
end

struct EachTraversal{N,T,E}
    edges::E
    position::NTuple{N,T}
    velocity::NTuple{N,T}
    function EachTraversal(edges::NTuple{N,Any}, position, velocity) where {N}
        E = typeof(edges)
        x1 = sum(first,edges)
        x2 = sum(position)
        x3 = sum(velocity)
        T = typeof(x1 + x2 / x3)
        pos = _makeNTuple(NTuple{N,T}, position)
        vel = _makeNTuple(NTuple{N,T}, velocity)
        new{N,T,E}(edges, pos, vel)
    end
end

function eachtraversal(ray, edges)
    EachTraversal(edges, ray.position, ray.velocity)
end

function _makeNTuple(::Type{NTuple{1,T}}, xs)::Tuple{T} where {T}
    x1, = xs
    (T(x1),)
end
function _makeNTuple(::Type{NTuple{2,T}}, xs)::NTuple{2,T} where {T}
    x1,x2 = xs
    (T(x1),T(x2))
end
function _makeNTuple(::Type{NTuple{3,T}}, xs)::NTuple{3,T} where {T}
    x1,x2,x3 = xs
    (T(x1),T(x2), T(x3))
end
function _makeNTuple(::Type{NTuple{4,T}}, xs)::NTuple{4,T} where {T}
    x1,x2,x3,x4 = xs
    (T(x1),T(x2), T(x3), T(x4))
end
function _makeNTuple(::Type{NTuple{N,T}}, xs)::NTuple{N,T} where {N,T}
    NTuple{N,T}(xs)
end

const VoxelTraversal{N,T} = NamedTuple{(:voxelindex, :entry_time, :exit_time), Tuple{CartesianIndex{N}, T, T}}

Base.eltype(::Type{<:EachTraversal{N,T}}) where {N,T} = VoxelTraversal{N,T}
Base.IteratorSize(::Type{<:EachTraversal}) = Base.SizeUnknown()

function is_inside_walls(x, walls)
    lo, hi = firstlast(walls)
    lo <= x <= hi
end
function is_inside_edges(pos, edges)
    all(is_inside_walls.(pos,edges))
end

function Base.iterate(o::EachTraversal)
    entry_time = grid_entry_time(o.position, o.velocity, o.edges)
    isfinite(entry_time) || return nothing
    position = let t = entry_time
        map(o.position, o.velocity) do pos, vel
            pos + t * vel
        end
    end
    is_inside_edges(position, o.edges) || return nothing
    index = _start_voxelindex(position, o.edges)
    exit_state = _first_wallhit(position, o.velocity, o.edges, index, entry_time)
    exit_state === nothing && return nothing
    any(has_finite_travel_time, exit_state.voxel_traversal_times) || return nothing
    item::eltype(o) = (
        voxelindex=CartesianIndex(index),
        entry_time=entry_time,
        exit_time=exit_state.time,
    )
    item, exit_state
end

function _first_wallhit(position, velocity, edges, index, time)
    limits = map(edges, index) do r, i
        @inbounds (r[i], r[i+1])
    end
    hittimes = let time = time
        map(position, velocity, limits) do pos, vel, lims
            interval_hit_times(pos,vel,lims) + time
        end
    end
    voxel_traversal_times = map(AbstractTraversalTime, edges,velocity)
    h = WallHit(;
        position,
        index,
        hittimes,
        time,
        velocity,
        voxel_traversal_times,
    )
    next_wallhit(h)
end

function interval_hit_times(pos::T, vel::T, lims) where {T}
    lo,hi = firstlast(lims)
    if iszero(vel)
        T(Inf)
    else
        max((hi - pos) / vel, (lo - pos) / vel)
    end
end

function is_inbound_voxelindex(o::EachTraversal, index::Tuple)
    all(
        map(o.edges,index) do r,i
            firstindex(r) <= i < lastindex(r)
        end
    )
end

function Base.iterate(o::EachTraversal, entry_state::WallHit)
    is_inbound_voxelindex(o, entry_state.index) || return nothing
    exit_state = next_wallhit(entry_state)
    item::eltype(o) = (
        voxelindex=CartesianIndex(entry_state.index),
        entry_time=entry_state.time,
        exit_time=exit_state.time,
    )
    item, exit_state
end

################################################################################
##### Adapt
################################################################################
import Adapt
Adapt.@adapt_structure EachTraversal
Adapt.@adapt_structure WallHit
Adapt.@adapt_structure TraversalTime_Vector
Adapt.@adapt_structure TraversalTime_Range


text:
using Random

include("../src/sorting.jl")

function testsort(sortfunc, testpairs)
    for (arr, ans) in testpairs
        temparr = copy(arr)
        sortfunc(temparr)
        @test temparr == ans
    end
end

function testrandom(sortfunc, seednum)
        Random.seed!(seednum)
        testdata = randn(100)

        d1 = copy(testdata)
        sort!(d1)
        d2 = copy(testdata)
        sortfunc(d2)

        @test d1 == d2
end

@testset "Sort" begin

    @testset "Basic cases" begin
        testpairs = [
                ([], []),
                ([3, 1, 2], [1, 2, 3]),
                ([4, 2, 3, 1], [1, 2, 3, 4]),
                ([4, 3, 3, 1, -1], [-1, 1, 3, 3, 4]),
                ([1.28, -2.591, -2.2266], [-2.591, -2.2266, 1.28]),
                ([1], [1]),
                (['b', 'c', 'a'], ['a', 'b', 'c']),
                (['E', 'A', 'S', 'Y', 'Q', 'U',
                  'E', 'S', 'T', 'I', 'O', 'N'],
                ['A', 'E', 'E', 'I', 'N', 'O',
                 'Q', 'S', 'S', 'T', 'U', 'Y'])
            ]

        testsort(insertsort!, testpairs)
        testsort(mergesort!, testpairs)
        testsort(quicksort!, testpairs)
    end

    @testset "Compare with Base.sort" begin
        for seed = [20190128, 20200828, 20211228]
            testrandom(insertsort!, seed)
            testrandom(mergesort!, seed)
            testrandom(quicksort!, seed)
        end
    end

end

max_stars_repo_name: Jetafull/AlgorithmsInJulia

max_issues_repo_name: Jetafull/algorithms-in-julia

max_forks_repo_name: Jetafull/algorithms-in-julia


text:
export iauZr
"""
**  Initialize an r-matrix to the null matrix.
**
**  This function is part of the International Astronomical Union's
**  SOFA (Standards Of Fundamental Astronomy) software collection.
**
**  Status:  vector/matrix support function.
**
**  Returned:
**     r        double[3][3]    r-matrix
**
**  This revision:  2013 June 18
**
**  SOFA release 2018-01-30
**
**  Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# void iauZr(double r[3][3])

function iauZr()
   r = zeros(Float64, 3, 3)

   ccall((:iauZr, libsofa_c), Cvoid,
           (Ptr{Cdouble},),
           r)

    # Transpose since C call return row-major operation
   return SMatrix{3,3}(r')
end


text:
struct CommunitiesCrime <: Name end

# mandatory methods
format(::Type{CommunitiesCrime}) = TabularData
task(::Type{CommunitiesCrime}) = Regression
source(::Type{CommunitiesCrime}) = "https://archive.ics.uci.edu/ml/datasets/communities+and+crime"

function downloadlink(::Type{CommunitiesCrime})
    return "https://archive.ics.uci.edu/ml/machine-learning-databases/communities/communities.data"
end

function preprocess(N::Type{CommunitiesCrime})
    return path -> preprocess_csv(N, path, :train; istarget = 1)
end

nattributes(::Type{CommunitiesCrime}) = (128, )
ninstances(::Type{CommunitiesCrime}) = (1994, 0, 0)

# optional methods
function checksum(::Type{CommunitiesCrime})
    return "09e0b5c07eae24c1efab19b2edee05e160e7f5743b6f31e31eec3d73624da2ea"
end

name(::Type{CommunitiesCrime}) = "Communities and Crime Data Set"
author(::Type{CommunitiesCrime}) = ["Michael Redmond"]

max_stars_repo_name: VaclavMacha/DatasetProvider

max_issues_repo_name: VaclavMacha/DatasetProvider

max_forks_repo_name: VaclavMacha/DatasetProvider


text:
using RecursiveArrayTools: ArrayPartition
using ComponentArrays: ComponentVector
using ..RingPolymers: RingPolymers

get_positions(u::ComponentVector) = u.r
get_velocities(u::ComponentVector) = u.v
get_flat_positions(u::ComponentVector) = @view get_positions(u)[:]
get_flat_velocities(u::ComponentVector) = @view get_velocities(u)[:]

# These must be kept since the symplectic integrators (VelocityVerlet) in DiffEq use ArrayPartitions.
get_positions(u::ArrayPartition) = u.x[2]
get_velocities(u::ArrayPartition) = u.x[1]

function RingPolymers.transform_from_normal_modes!(p::RingPolymers.RingPolymerParameters, u)
    RingPolymers.transform_from_normal_modes!(p, DynamicsUtils.get_velocities(u))
    RingPolymers.transform_from_normal_modes!(p, DynamicsUtils.get_positions(u))
    return nothing
end

function RingPolymers.transform_to_normal_modes!(p::RingPolymers.RingPolymerParameters, u)
    RingPolymers.transform_to_normal_modes!(p, DynamicsUtils.get_velocities(u))
    RingPolymers.transform_to_normal_modes!(p, DynamicsUtils.get_positions(u))
    return nothing
end

function get_quantum_subsystem end

get_mapping_positions(u::ComponentVector) = u.qmap
get_mapping_momenta(u::ComponentVector) = u.pmap
get_mapping_positions(u::ComponentVector, i) = @view get_mapping_positions(u)[:,i]
get_mapping_momenta(u::ComponentVector, i) = @view get_mapping_momenta(u)[:,i]


text:
function codebook_trans3To2(aq_codebook::AbstractArray{<:AbstractFloat, 3})
    (M,K,D) = size(aq_codebook)
    aq_codebook_2d = zeros(M*K,D)
    for i in 1:M
        for j in 1:K
            aq_codebook_2d[(i-1)*K + j,:] = aq_codebook[i,j,:]
        end
    end
    return aq_codebook_2d
end

function codebook_trans2To3(aq_codebook::AbstractArray{<:AbstractFloat, 2},M::Integer)
    (MK,D) = size(aq_codebook)
    K = Int(MK / M)

    aq_codebook_3d = zeros(M, K, D)
    for i in 1:M
        for j in 1:K
            aq_codebook_3d[i,j,:] = aq_codebook[(i-1)*K + j,:]
        end
    end
    return aq_codebook_3d
end

function aq_codesTo_codesB(aq_codes::AbstractArray{<:Integer,2}, K::Integer)
    # aq_codes = (n, M)
    # self.aq_codes_2 = aq_codes + np.arange(M) * K
    rangeM = (1:M) .- 1
    Arr1 = convert(Array{Int64,1}, rangeM.*K)
    codesB = broadcast(+,aq_codes,Arr1[[CartesianIndex()],:])
    return codesB
end


text:
module BlochSim

using LinearAlgebra: Diagonal, diagm

include("helpers.jl")
include("Spin.jl")
include("sequences.jl")

export GAMMA, GAMBAR
export AbstractSpin, Spin, SpinMC
export freeprecess, freeprecess!
export excitation, excitation!
export spoil, spoil!
export combine, applydynamics!

export mese!, mese, spgr!

end


text:
# Michelle Medina
# August 5, 2020

using PyPlot

# x = collect(range(0, stop=1, length = 1025))
# level = 4
# for l in 0:level, j in 0:2^l
# 	plot(x, basis.(l, j, x), "m")
# 	plot(x, derivBasis.(l, j, x))
# 	xlabel("x")
# 	ylabel("phi'(x)")
# 	show()
# end

# plot(x, derivBasis.(0, 0, x))
# plot(x, derivBasis.(0, 1, x))
# plot(x, derivBasis.(1, 0, x))
# plot(x, derivBasis.(1, 1, x))
# plot(x, derivBasis.(1, 2, x))
# show()

# xVal = collect(range(0.0, stop=1.0, length = 129))
# y = calcNodal(x -> sin(4*pi*x), xVal)
# cosx = calcNodal(x -> 4*pi*cos(4*pi*x), xVal)
# modal = Nodal_2_H(y)
#
# dy = zeros(length(xVal))
# for index in CartesianIndices(xVal)
# 	dy[index] = derivEvaluate(modal, xVal[index])
# end

# plot(xVal,y.values)
# plot(xVal,dy, "g .")
# plot(xVal, cosx.values, "r")
#

error = []
Num = []

level = 13
for l in 0:level
	N = 2^l + 1

	xOld = collect(range(0.0, stop=1.0, length = N))
	U = calcNodal(x -> exp(x), xOld)
	dU = calcNodal(x -> exp(x), xOld)
	modal_U = Nodal_2_H(U)

	# Compute the interpolation error
	le = 10 * (2^l+1)
	h = 1/le
	xNew = collect(range(0.0 + h/2, stop=1.0 - h/2, length = le))
	dUNew = calcNodal(x -> exp(x), xNew)
	dU_from_modal_U = calcNodal(x -> derivEvaluate(modal_U, x), xNew)

	error_up_to = sqrt(sum((dUNew.values - dU_from_modal_U.values).^2)/length(dUNew.values))

	append!(error, error_up_to)
	append!(Num, N)

	# @show l, N, error_up_to
end

loglog(Num, error)
loglog(Num, 1 ./Num)
loglog(Num, 1 ./(Num.^2))
show()