text:
@doc Markdown.doc"""
    collect_small_blocks(G::MatElem) -> Vector{MatElem}

Given a block diagonal matrix $G$ consisting of blocks of size $1$ or $2$,
returns the blocks.
"""
function collect_small_blocks(G)
  L = _get_small_block_indices(G)
  blocks = typeof(G)[]
  for i in 1:length(L)
    j, nj = L[i]
    l = j + nj - 1
    push!(blocks, G[j:l, j:l])
  end
  return blocks
end

function _unit_part(x, p)
  if iszero(x)
    return x
  end
  v = _val(x, p)
  y = divexact(x, parent(x)(p^v))
  return y
end

function _ispadic_normal_form(G, p)
  if isodd(p)
    return _ispadic_normal_form_odd(G, p)
  else
    return _ispadic_normal_form_dyadic(G, p)
  end
end

function _ispadic_normal_form_odd(G, p)
  _D = G
  B = collect_small_blocks(_D)
  if G != diagonal_matrix(B)
    return false
  end
  curv = valuation(reduce(gcd, _eltseq(B[1])), p)
  blocks = [[divexact(B[1], fmpq(p)^curv)]]
  for i in 2:length(B)
    v = valuation(reduce(gcd, _eltseq(B[i])), p)
    if v == curv
      push!(blocks[end], divexact(B[i], fmpq(p)^curv))
    else
      curv = v
      push!(blocks, [divexact(B[i], fmpq(p)^curv)])
    end
  end

  for b in blocks
    if !(all(x -> nrows(x) == 1, b))
      return false
    end
  end

  o = identity_matrix(QQ, 1)

  F,  = FiniteField(p, 1, cached = false)

  for i in 1:length(blocks)
    if all(==(o), blocks[i])
      continue
    else
      if !(all(==(o), blocks[i][1:end-1]))
        return false
      end
      u = blocks[i][end][1, 1]
      m = FlintZZ(u)
      if issquare(F(m))
        return false
      end
      for i in 1:(m-1)
        if !issquare(F(i))
          return false
        end
      end
    end
  end
  return true
end

function _ispadic_normal_form_dyadic(G, p)
  _D = G
  B = collect_small_blocks(_D)
  if G != diagonal_matrix(B)
    return false
  end
  curv = valuation(reduce(gcd, _eltseq(B[1])), 2)
  blocks = [[divexact(B[1], fmpq(2)^curv)]]
  for i in 2:length(B)
    v = valuation(reduce(gcd, _eltseq(B[i])), 2)
    if v == curv
      push!(blocks[end], divexact(B[i], fmpq(2)^curv))
    else
      curv = v
      push!(blocks, [divexact(B[i], fmpq(2)^curv)])
    end
  end

  U = matrix(FlintQQ,2,2,[0,1,1,0])
  V = matrix(FlintQQ,2,2,[2,1,1,2])
  W1 = matrix(FlintQQ,1,1,[1])
  W3 = matrix(FlintQQ,1,1,[3])
  W5 = matrix(FlintQQ,1,1,[5])
  W7 = matrix(FlintQQ,1,1,[7])

  for B in blocks
    i = 1
    while i <= length(B) && B[i] == U
      i += 1
    end

    if i <= length(B) && B[i] == V
      i += 1
    end

    if i <= length(B) && (B[i] == W1 || B[i] == W3 || B[i] == W5 || B[i] == W7)
      i += 1
    end

    if i <= length(B) && (B[i] == W1 || B[i] == W3 || B[i] == W5 || B[i] == W7)
      i += 1
    end

    if i != length(B) + 1
      return false
    end
  end
  return true
end

@doc Markdown.doc"""
    padic_normal_form(G::MatElem, p::fmpz; prec::Int = -1, partial::Bool = false)
                                              -> MatElem{fmpq}, MatElem{fmpq}

Return the normal `D` and the transformation `T` to the `p`-adic normal form of
the symmetric matrix `G`, such that `d * D = d * B * G * transpose(B)` holds modulo `p^prec`.
If `prec == -1`,

Let `p` be odd and `u` be the smallest non-square modulo `p`.  The normal form
is a block diagonal matrix with blocks `p^k G_k` such that `G_k` is either the
identity matrix or the identity matrix with the last diagonal entry replaced by
`u`.

If `p = 2`, define the `1` by `1` matrices

    W1 = [1;]
    W3 = [3;]
    W5 = [5;]
    W7 = [7;]

and the `2` by `2` matrices

    U = [0 1;
         1 0]
    V = [2 1;
         1 2]

For `p = 2` the partial normal form is a block diagonal matrix with blocks
`2^k G_k` such that $G_k$ is a block diagonal matrix of the form
`[U`, ... , `U`, `V`, `Wa`, `Wb]`
where we allow `V`, `Wa`, `Wb` to be `0 \times 0` matrices.

Further restrictions to the full normal form apply.
We refer to [MirMor2009]_ IV Definition 4.6. for the details.

If `partial` is set, only the partial normal form is returned.
"""
function padic_normal_form(G, p::IntegerUnion; prec::Int = -1, partial::Bool = false)
  return _padic_normal_form(change_base_ring(FlintQQ, G), fmpz(p), prec = prec, partial = partial)
end

# For a definition in the even case, see Definition 4.6 of Miranda, Morrison,
# "Embeddings of Integral Quadratic Forms", 2009.
function _padic_normal_form(G::fmpq_mat, p::fmpz; prec::Int = -1, partial::Bool = false)

  _G = deepcopy(G)
  dd = denominator(G)
  G0 = change_base_ring(FlintZZ, dd * G)
  d = valuation(dd, p)
  n = nrows(G)
  r = rank(G0)

  if r != n
    _, U = hnf_with_transform(G0)
    _ker = U[(r + 1):n, :]
    _nondeg = U[1:r, :]
    ker = change_base_ring(FlintQQ, _ker)
    nondeg = change_base_ring(FlintQQ, _nondeg)
    G = p^d * nondeg * G * transpose(nondeg)
  else
    ker = zero_matrix(FlintQQ, 0, n)
    nondeg = identity_matrix(FlintQQ, n)
    G = p^d * G
  end
  # continue with the non-degenerate part

  if prec == -1
    prec = valuation(det(G)::fmpq, p) + 4
  end

  modu = p^prec
  R = ResidueRing(FlintZZ, modu, cached = false)
  Gmod = map(q -> R(invmod(denominator(q), modu) * numerator(q)), G) # this will probably fail
  D = deepcopy(Gmod)

  n = ncols(Gmod)

  Qp = PadicField(p, prec, cached = false)

  if n == 0
    return (zero_matrix(FlintQQ, n, n), zero_matrix(FlintQQ, n, n))::Tuple{fmpq_mat, fmpq_mat}
  end

  # the transformation matrix is called B

  if p == 2
    D, B = _jordan_2_adic(Gmod)
  else
    D, B = _jordan_odd_adic(Gmod, p)
  end

  @assert B * Gmod * transpose(B) == D

  D, B1 = _normalize(D, p)
  B = B1 * B

  @assert B * Gmod * transpose(B) == D
  # we have reached a normal form for p != 2
  # for p == 2 extra work is necessary
  if p == 2
    D, B1 = _two_adic_normal_forms(D, p, partial = partial)
    B = B1 * B
  end
  @assert B * Gmod * transpose(B) == D

  __nondeg = B * map_entries(x -> R(numerator(x)), nondeg)
  B = vcat(__nondeg, map_entries(x -> R(numerator(x)), ker))
  D = diagonal_matrix([D, zero_matrix(base_ring(D), nrows(ker), nrows(ker))])
  _Gmod = map(q -> R(invmod(denominator(q), modu) * numerator(q)), p^d * _G)
  if p == 2
    @assert B * _Gmod * transpose(B) == diagonal_matrix(collect_small_blocks(D)) == D
  else
    @assert B * _Gmod * transpose(B) == diagonal_matrix(diagonal(D)) == D
  end
  #    if debug:
  @assert _val(det(B), p) == 0     # B is invertible!

  DD = map_entries(x -> FlintQQ(lift(x))//p^d, D)
  BB = map_entries(x -> FlintQQ(lift(x)), B)

  return (DD, BB)::Tuple{fmpq_mat, fmpq_mat}
end

#def _find_min_p(G, cnt, lower_bound=0):
#    r"""
#    Find smallest valuation below and right from ``cnt`` prefering the diagonal.
#
#    INPUT:
#
#    - ``G`` -- a symmetric `n` by `n` matrix in `\QQ_p`
#    - ``cnt`` -- start search from this index
#    - ``lower_bound`` -- an integer (default: ``0``)
#      a lower bound for the valuations used for optimization
#
#    OUTPUT:
#
#    - ``min`` -- minimal valuation
#    - ``min_i`` -- row index of the minimal valuation
#    - ``min_j`` -- column index of the minimal valuation
#
#    EXAMPLES::
#
#        sage: from sage.quadratic_forms.genera.normal_form import _find_min_p
#        sage: G = matrix(Qp(2, show_prec=False), 3, 3, [4,0,1,0,4,2,1,2,1])
#        sage: G
#        [2^2   0   1]
#        [  0 2^2   2]
#        [  1   2   1]
#        sage: _find_min_p(G, 0)
#        (0, 2, 2)
#        sage: G = matrix(Qp(3, show_prec=False), 3, 3, [4,0,1,0,4,2,1,2,1])
#        sage: G
#        [1 + 3     0     1]
#        [    0 1 + 3     2]
#        [    1     2     1]
#        sage: _find_min_p(G, 0)
#        (0, 0, 0)
#    """

function _find_min_p(G, p, cnt, lower_bound = 0)
  n = ncols(G)
  minval = _val(G[cnt, cnt], p)
  min_i = cnt
  min_j = cnt
  # diagonal has precedence
  for i in cnt:n
    v = _val(G[i, i], p)
    if v == lower_bound
      return lower_bound, i, i
    end
    if v < minval
      min_i = i
      min_j = i
      minval = v
    end
  end
  for i in cnt:n
    for j in (i + 1):n
      v = _val(G[i, j], p)
      if v == lower_bound
        return lower_bound, i, j
      end
      if v < minval
        min_i = i
        min_j = j
        minval = v
      end
    end
  end
  return minval, min_i, min_j
end

#def _get_small_block_indices(G):
#    r"""
#    Return the indices of the blocks.
#
#    For internal use in :meth:`collect_small_blocks`.
#
#    INPUT:
#
#    - ``G`` -- a block_diagonal matrix consisting of `1` by `1` and `2` by `2` blocks
#
#    OUTPUT:
#
#    - a list of integers
#
#    EXAMPLES::
#
#        sage: from sage.quadratic_forms.genera.normal_form import _get_small_block_indices
#        sage: W1 = Matrix([1])
#        sage: U = Matrix(ZZ, 2, [0, 1, 1, 0])
#        sage: G = Matrix.block_diagonal([W1, U, U, W1, W1, U, W1, U])
#        sage: _get_small_block_indices(G)
#        [0, 1, 3, 5, 6, 7, 9, 10]

function _get_small_block_indices(G)
  L = Tuple{Int, Int}[]
  n = ncols(G)
  i = 1
  while i < n
    if !iszero(G[i, i + 1])
      push!(L, (i, 2))
      i += 2
    else
      push!(L, (i, 1))
      i += 1
    end
  end

  if i == n
    push!(L, (i, 1))
  end

  return L
end

#    """

#    L = []
#    n = G.ncols()
#    i = 0
#    while i < n-1:
#        L.append(i)
#        if G[i, i+1]!=0:
#            i += 2
#        else:
#            i += 1
#    if i == n-1:
#        L.append(i)
#    return L[:]
#
#def _get_homogeneous_block_indices(G):
#    r"""
#    Return the indices of the homogeneous blocks.
#
#    We call a matrix homogeneous if it is a multiple of an invertible matrix.
#    Sometimes they are also called modular.
#
#    INPUT:
#
#    - ``G`` - a block diagonal matrix over the p-adics
#      with blocks of size at most `2`.
#
#    OUTPUT:
#
#    - a list of integers
#
#    EXAMPLES::
#
#        sage: from sage.quadratic_forms.genera.normal_form import _get_homogeneous_block_indices
#        sage: W1 = Matrix(Zp(2), [1])
#        sage: V = Matrix(Zp(2), 2, [2, 1, 1, 2])
#        sage: G = Matrix.block_diagonal([W1, V, V, 2*W1, 2*W1, 8*V, 8*W1, 16*V])
#        sage: _get_homogeneous_block_indices(G)
#        ([0, 5, 7, 10], [0, 1, 3, 4])
#        sage: G = Matrix.block_diagonal([W1, V, V, 2*W1, 2*W1, 8*V, 8*W1, 16*W1])
#        sage: _get_homogeneous_block_indices(G)
#        ([0, 5, 7, 10], [0, 1, 3, 4])
#    """
function _get_homogeneous_block_indices(G, p)
  L = Int[]
  vals = Int[]
  n = ncols(G)
  i = 1
  val = -5
  while i < n
    if G[i, i + 1] != 0
      m = _val(G[i, i + 1], p)
    else
      m = _val(G[i, i], p)
    end
    if m > val
      push!(L, i)
      val = m
      push!(vals, val)
    end
    if G[i, i + 1] != 0
      i += 1
    end
    i += 1
  end
  if i == n
    m = _val(G[i, i], p)
    if m > val
      push!(L, i)
      val = m
      push!(vals, val)
    end
  end
  return L, vals
end

#def _homogeneous_normal_form(G, w):
#    r"""
#    Return the homogeneous normal form of the homogeneous ``G``.
#
#    INPUT:
#
#    - ``G`` -- a modular symmetric matrix over the `2`-adic integers
#      in partial normal form
#
#    OUTPUT:
#
#    - ``B`` -- an invertible matrix over the basering of ``G``
#      such that ``B*G*B.T`` is in homogeneous normal form
#
#    EXAMPLES::
#
#        sage: from sage.quadratic_forms.genera.normal_form import _homogeneous_normal_form
#        sage: R = Zp(2, type = 'fixed-mod', print_mode='terse', show_prec=False)
#        sage: U = Matrix(R, 2, [0,1,1,0])
#        sage: V = Matrix(R, 2, [2,1,1,2])
#        sage: W1 = Matrix(R, 1, [1])
#        sage: W3 = Matrix(R, 1, [3])
#        sage: W5 = Matrix(R, 1, [5])
#        sage: W7 = Matrix(R, 1, [7])
#        sage: G = Matrix.block_diagonal([V, W1])
#        sage: B = _homogeneous_normal_form(G, 1)[1]
#        sage: B * G * B.T
#        [2 1 0]
#        [1 2 0]
#        [0 0 1]
#        sage: G = Matrix.block_diagonal([V, W1, W3])
#        sage: B = _homogeneous_normal_form(G, 2)[1]
#        sage: B * G * B.T
#        [2 1 0 0]
#        [1 2 0 0]
#        [0 0 1 0]
#        [0 0 0 3]
#        sage: G = Matrix.block_diagonal([U, V, W1, W5])
#        sage: B = _homogeneous_normal_form(G, 2)[1]
#        sage: B * G * B.T
#        [0 1 0 0 0 0]
#        [1 0 0 0 0 0]
#        [0 0 0 1 0 0]
#        [0 0 1 0 0 0]
#        [0 0 0 0 7 0]
#        [0 0 0 0 0 7]
#        sage: G = Matrix.block_diagonal([U, W7, W3])
#        sage: B = _homogeneous_normal_form(G, 2)[1]
#        sage: B * G * B.T
#        [0 1 0 0]
#        [1 0 0 0]
#        [0 0 3 0]
#        [0 0 0 7]
#        sage: G = Matrix.block_diagonal([V, W5, W5])
#        sage: B = _homogeneous_normal_form(G, 2)[1]
#        sage: B * G * B.T
#        [0 1 0 0]
#        [1 0 0 0]
#        [0 0 3 0]
#        [0 0 0 7]
#        sage: G = Matrix.block_diagonal([V, W3, W3])
#        sage: B = _homogeneous_normal_form(G, 2)[1]
#        sage: B * G * B.T
#        [0 1 0 0]
#        [1 0 0 0]
#        [0 0 1 0]
#        [0 0 0 5]
#        sage: G = Matrix.block_diagonal([V, W1, W3])
#        sage: B = _homogeneous_normal_form(G, 2)[1]
#        sage: B * G * B.T
#        [2 1 0 0]
#        [1 2 0 0]
#        [0 0 1 0]
#        [0 0 0 3]
#        sage: G = Matrix.block_diagonal([W3, W3])
#        sage: B = _homogeneous_normal_form(G, 2)[1]
#        sage: B * G * B.T
#        [7 0]
#        [0 7]
#    """
function _homogeneous_normal_form(G, w, p)
  D = deepcopy(G)
  R = base_ring(G)
  n = ncols(G)
  B = identity_matrix(base_ring(G), n)
  if w == 2
    if n > 2 && !iszero(D[end - 2, end - 2])
      v = 2
    else
      v = 0
    end

    if v == 2
      e1 = _unit_part(D[end - 1, end - 1], p)
      e2 = _unit_part(D[end - 0, end - 0], p)
      e = Set([e1, e2])
      E = [Set(R.([1, 3])), Set(R.([1, 7])), Set(R.([5, 7])), Set(R.([3, 5]))]
      if !(e in E)
        B[end - 3:end, :] = _relations(D[(end - 3):end, (end - 3):end], 5, p) * B[(end - 3):end, :]
        D = B * G * transpose(B)
      end
    end

    e1 = _unit_part(D[end - 1, end - 1], p)
    e2 = _unit_part(D[end - 0, end - 0], p)
    e = Set([e1, e2])
    E = [Set(R.([3, 3])), Set(R.([3, 5])), Set(R.([5, 5])), Set(R.([5, 7]))]
    if e in E
      B[end - 1:end, :] = _relations(D[(end - 1):end, end - 1:end], 1, p) * B[end - 1:end, :]
      D = B * G * transpose(B)
    end
    # assert that e1 < e2
    e1 = _unit_part(D[end - 1, end - 1], p)
    e2 = _unit_part(D[end - 0, end - 0], p)
    if lift(e1) > lift(e2)
      swap_rows!(B, n, n - 1)
      swap_rows!(D, n, n - 1)
      swap_cols!(D, n, n - 1)
    end
  end
  return D, B
end

#        # assert that e1 < e2
#        e1 = D[-2,-2].unit_part()
#        e2 = D[-1,-1].unit_part()
#        if ZZ(e1) > ZZ(e2):
#            B.swap_rows(n-1, n-2)
#            D.swap_rows(n-1, n-2)
#            D.swap_columns(n-1, n-2)
#    return D, B
#
#def _jordan_odd_adic(G):
#    r"""
#    Return the Jordan decomposition of a symmetric matrix over an odd prime.
#
#    INPUT:
#
#    - a symmetric matrix over `\ZZ_p` of type ``'fixed-mod'``
#
#    OUTPUT:
#
#    - ``D`` -- a diagonal matrix
#    - ``B`` -- a unimodular matrix such that ``D = B * G * B.T``
#
#    EXAMPLES::
#
#        sage: from sage.quadratic_forms.genera.normal_form import _jordan_odd_adic
#        sage: R = Zp(3, prec=2, print_mode='terse', show_prec=False)
#        sage: A4 = Matrix(R,4,[2, -1, 0, 0, -1, 2, -1, 0, 0, -1, 2, -1, 0, 0, -1, 2])
#        sage: A4
#        [2 8 0 0]
#        [8 2 8 0]
#        [0 8 2 8]
#        [0 0 8 2]
#        sage: D, B = _jordan_odd_adic(A4)
#        sage: D
#        [2 0 0 0]
#        [0 2 0 0]
#        [0 0 1 0]
#        [0 0 0 8]
#        sage: D == B*A4*B.T
#        True
#        sage: B.determinant().valuation() == 0
#        True
#    """
function _jordan_odd_adic(G, p)
  R = base_ring(G)
  D = deepcopy(G)
  n = nrows(G)

  # transformation matrix
  B = identity_matrix(R, n)
  cnt = 1
  minval = 0
  while cnt <= n
    pivot = _find_min_p(D, fmpz(p), cnt, minval)
    piv1 = pivot[2]
    piv2 = pivot[3]
    minval = pivot[1]
    # the smallest valuation is on the diagonal
    if piv1 == piv2
      # move pivot to position [cnt, cnt]
      if piv1 != cnt
        swap_rows!(B, cnt, piv1)
        swap_rows!(D, cnt, piv1)
        swap_cols!(D, cnt, piv1)
      end
      # we are alread orthgonal to the part with i < cnt
      # now make the rest orthgonal too
      for i in (cnt + 1):n
        if !iszero(D[i, cnt])
          c = divexact(D[i, cnt], D[cnt, cnt])
          for j in 1:ncols(B)
            B[i, j] = B[i, j] - c * B[cnt, j]
            D[i, j] = D[i, j] - c * D[cnt, j]
          end
          for j in 1:nrows(D)
            D[j, i] = D[j, i] - c * D[j, cnt]
          end
        end
      end
      cnt += 1
    else
      # the smallest valuation is off the diagonal
      row = pivot[2]
      col = pivot[3]
      for j in 1:ncols(B)
        B[row, j] = B[row, j] + B[col, j]
        D[row, j] = D[row, j] + D[col, j]
      end
      for j in 1:nrows(D)
        D[j, row] = D[j, row] + D[j, col]
      end
      # the smallest valuation is now on the diagonal
    end
  end
  return D, B
end
#                    D[:, i] += - c * D[:, cnt]
#            cnt = cnt + 1
#        else:
#            # the smallest valuation is off the diagonal
#            row = pivot[1]
#            col = pivot[2]
#            B[row, :] += B[col, :]
#            D[row, :] += D[col, :]
#            D[:, row] += D[:, col]
#            # the smallest valuation is now on the diagonal
#    return D, B
#
#def _jordan_2_adic(G):
#    r"""
#    Transform a symmetric matrix over the `2`-adic integers into jordan form.
#
#    Note that if the precision is too low, this method fails.
#    The method is only tested for input over `\ZZ_2` of ``'type=fixed-mod'``.
#
#    INPUT:
#
#    - ``G`` -- symmetric `n` by `n` matrix in `\ZZ_p`
#
#    OUTPUT:
#
#    - ``D`` -- the jordan matrix
#    - ``B`` -- transformation matrix, i.e, ``D = B * G * B^T``
#
#    The matrix ``D`` is a block diagonal matrix consisting
#    of `1` by `1` and `2` by `2` blocks.
#    The `2` by `2` blocks are matrices of the form
#    `[[2a,  b], [b, 2c]] * 2^k`
#    with `b` of valuation `0`.
#
#    EXAMPLES::
#
#        sage: from sage.quadratic_forms.genera.normal_form import _jordan_2_adic
#        sage: R = Zp(2, prec=3, print_mode='terse', show_prec=False)
#        sage: A4 = Matrix(R,4,[2, -1, 0, 0, -1, 2, -1, 0, 0, -1, 2, -1, 0, 0, -1, 2])
#        sage: A4
#        [2 7 0 0]
#        [7 2 7 0]
#        [0 7 2 7]
#        [0 0 7 2]
#        sage: D, B = _jordan_2_adic(A4)
#        sage: D
#        [ 2  7  0  0]
#        [ 7  2  0  0]
#        [ 0  0 12  7]
#        [ 0  0  7  2]
#        sage: D == B*A4*B.T
#        True
#        sage: B.determinant().valuation() == 0
#        True
#    """

function _jordan_2_adic(G)
  D = deepcopy(G)
  n = ncols(G)
  R = base_ring(G)

  # transformation matrix
  B = identity_matrix(base_ring(G), n)

  # indices of the diagonal entrys which are already used
  cnt = 1
  local minval::Int
  while cnt <= n
    pivot = _find_min_p(D, fmpz(2), cnt)
    minval = pivot[1]
    piv1 = pivot[2]
    piv2 = pivot[3]
    # the smallest valuation is on the diagonal
    if piv1 == piv2
      # move pivot to position [cnt, cnt]
      if piv1 != cnt
        swap_rows!(B, cnt, piv1)
        swap_rows!(D, cnt, piv1)
        swap_cols!(D, cnt, piv1)
      end
      # we are already orthgonal to the part with i < cnt
      # now make the rest orthgonal too
      for i in cnt+1:n
        if !iszero(D[i, cnt])
          c = divexact(D[i, cnt], D[cnt, cnt])
          for j in 1:ncols(B)
            B[i, j] = B[i, j] - c * B[cnt, j]
            D[i, j] = D[i, j] - c * D[cnt, j]
          end
          for j in 1:nrows(D)
            D[j, i] = D[j, i] - c * D[j, cnt]
          end
        end
      end
      cnt = cnt + 1
    else
      # the smallest valuation is off the diagonal
      # move this 2 x 2 block to the top left (starting from cnt)
      if piv1 != cnt
        swap_rows!(B, cnt, piv1)
        swap_rows!(D, cnt, piv1)
        swap_cols!(D, cnt, piv1)
      end
      if piv2 != cnt + 1
        swap_rows!(B, cnt + 1, piv2)
        swap_rows!(D, cnt + 1, piv2)
        swap_cols!(D, cnt + 1, piv2)
      end
      # we split off a 2 x 2 block
      # if it is the last 2 x 2 block, there is nothing to do.
      if cnt != n - 1
        content = R(fmpz(2)^minval)
        DD = D[cnt:cnt + 1, cnt:cnt + 1]
        #@show objectid(DD)
        _eqn_mat = Hecke._eltseq(DD)
        eqn_mat = matrix(R, 2, 2, [divexact(e, content) for e in _eqn_mat])
        # calculate the inverse without using division
        eqn_mat_inv = inv(det(eqn_mat)) * matrix(R, 2, 2, [eqn_mat[2, 2], -eqn_mat[1, 2], -eqn_mat[2, 1], eqn_mat[1, 1]])
        @assert isone(eqn_mat * eqn_mat_inv)
        B1 = B[cnt:cnt+1, :]
        #@show objectid(B1)
        DD = D[cnt + 2:nrows(D), cnt:cnt + 1]
        #@show objectid(DD)
        B2 = DD * eqn_mat_inv
        for i in 1:nrows(B2)
          for j in 1:ncols(B2)
            B2[i, j] = divexact(B2[i, j], content)
          end
        end
        Bv = @view B[cnt + 2:nrows(B), :]
        add!(Bv, Bv, -B2 * B1)
        Dv = @view D[cnt:nrows(D),cnt:ncols(D)]
        BB = B[cnt:nrows(B), :]
        mul!(Dv, BB * G, transpose(BB))
      end
      cnt += 2
    end
  end
  return D, B
end

#
#def _min_nonsquare(p):
#    r"""
#    Return the minimal nonsquare modulo the prime `p`.
#
#    INPUT:
#
#    - ``p`` -- a prime number
#
#    OUTPUT:
#
#    - ``a`` -- the minimal nonsquare mod `p`
#
#    EXAMPLES::
#
#        sage: from sage.quadratic_forms.genera.normal_form import _min_nonsquare
#        sage: _min_nonsquare(2)
#        sage: _min_nonsquare(3)
#        2
#        sage: _min_nonsquare(5)
#        2
#        sage: _min_nonsquare(7)
#        3
#    """
#    R = GF(p)
#    for i in R:
#        if not R(i).is_square():
#            return i

function _min_nonsquare(p)
  Rx, x = PolynomialRing(GF(p, cached = false), "x", cached = false)
  for i in 1:p
    if length(factor(x^2 - i)) == 1
      return i
    end
  end
end

function _issquare(d::Nemo.fmpz_mod, p::fmpz)
  @assert valuation(lift(d), p) == 0
  Rx, x = PolynomialRing(parent(d), "x", cached = false)
  rts = roots(x^2 - d, p, valuation(modulus(parent(d)), p))
  return length(rts) > 0
end

function _issquare(d::nmod, p)
  f = ZZ(modulus(parent(d)))
  R = ResidueRing(FlintZZ, f, cached = false)
  g = R(d)
  return _issquare(g, ZZ(p))
end

function _sqrt(d::Nemo.fmpz_mod, p::fmpz)
  @assert valuation(lift(d), p) == 0
  Rx, x = PolynomialRing(parent(d), "x", cached = false)
  rts = roots(x^2 - d, p, valuation(modulus(parent(d)), p))
  r = rts[1][1]
  @assert r^2 == d
  return r
end

#
#def _normalize(G, normal_odd=True):
#    r"""
#    Return the transformation to sums of forms of types `U`, `V` and `W`.
#
#    Part of the algorithm :meth:`p_adic_normal_form`.
#
#    INPUT:
#
#    - ``G`` -- a symmetric matrix over `\ZZ_p` in jordan form --
#      the output of :meth:`p_adic_normal_form` or :meth:`_jordan_2_adic`
#    - ``normal_odd`` -- bool (default: True) if true and `p` is odd,
#      compute a normal form.
#
#    OUTPUT:
#
#    - ``(D, B)`` -- a pair of matrices such that ``D=B*G*B.T``
#      is a sum of forms of types `U`, `V` and `W` for `p=2` or
#      diagonal with diagonal entries equal `1` or `u`
#      where `u` is the smallest non-square modulo the odd prime `p`.
#
#    EXAMPLES::
#
#        sage: from sage.quadratic_forms.genera.normal_form import _normalize
#        sage: R = Zp(3, prec = 5, type = 'fixed-mod', print_mode='series', show_prec=False)
#        sage: G = matrix.diagonal(R, [1,7,3,3*5,3,9,-9,27*13])
#        sage: D, B =_normalize(G)
#        sage: D
#        [    1     0     0     0     0     0     0     0]
#        [    0     1     0     0     0     0     0     0]
#        [    0     0     3     0     0     0     0     0]
#        [    0     0     0     3     0     0     0     0]
#        [    0     0     0     0   2*3     0     0     0]
#        [    0     0     0     0     0   3^2     0     0]
#        [    0     0     0     0     0     0 2*3^2     0]
#        [    0     0     0     0     0     0     0   3^3]
#    """

function _normalize(G, p, normal_odd = true)
  R = base_ring(G)
  D = deepcopy(G)
  n = ncols(G)
  B = identity_matrix(R, n)
  if p != 2
    # squareclasses 1, v
    v = _min_nonsquare(p)
    vv = R(v)
    non_squares = []
    val = 0
    for i in 1:n
      if _val(D[i, i], p) > val
        # a new block starts
        val = _val(D[i, i], p)
        if normal_odd && length(non_squares) != 0
          # move the non-square to the last entry of the previous block
          j = pop!(non_squares)
          swap_rows!(B, j, i - 1)
        end
      end
      d = _unit_part(D[i, i], p)
      if _issquare(d, p)
        D[i, i] = 1
        c = _sqrt(inv(d), p)
        for j in 1:ncols(B)
          B[i, j] = c * B[i, j]
        end
      else
        D[i, i] = vv
        for j in 1:ncols(B)
          B[i, j] = B[i, j] * _sqrt(vv * inv(d), p)
        end
        if normal_odd && length(non_squares) != 0
          # we combine two non-squares to get the 2 x 2 identity matrix
          j = pop!(non_squares)
          trafo = _normalize_odd_twobytwo(D[[i, j], [i, j]], p)
          _BB = trafo * B[[i, j], :]
          B[i, :] = _BB[1, :]
          B[j, :] = _BB[2, :]
          #B[[i, j], :] = trafo * B[[i, j], :]
          D[i, i] = 1
          D[j, j] = 1
        else
          push!(non_squares, i)
        end
      end
    end
    if normal_odd && length(non_squares) != 0
      j = pop!(non_squares)
      swap_rows!(B, j, n)
    end
  else
    # squareclasses 1, 3, 5, 7 modulo 8
    for i in 1:n
      d = _unit_part(D[i, i], p)
      if !iszero(d)
        vv = R(mod(lift(d), 8))
        for j in 1:ncols(B)
          B[i, j] = B[i, j] * _sqrt(vv * inv(d), p)
        end
      end
    end
    D = B * G * transpose(B)
    for i in 1:(n - 1)
      if !iszero(D[i, i + 1])
        # there is a 2 x 2 block here
        block = D[i:i+1, i:i+1]
        trafo = _normalize_twobytwo(block, p)
        Bv = @view B[i:i+1, :]
        T = trafo * Bv
        for i in 1:nrows(T)
          for j in 1:ncols(T)
            Bv[i, j] = T[i, j]
          end
        end
        #mul!(Bv, trafo, Bv)
      end
    end
  end
  D = B * G * transpose(B)
  return D, B
end

#def _normalize_2x2(G):
#    r"""
#    Normalize this indecomposable `2` by `2` block.
#
#    INPUT:
#
#    ``G`` - a `2` by `2` matrix over `\ZZ_p`
#    with ``type = 'fixed-mod'`` of the form::
#
#        [2a  b]
#        [ b 2c] * 2^n
#
#    with `b` of valuation 1.
#
#    OUTPUT:
#
#    A unimodular `2` by `2` matrix ``B`` over `\ZZ_p` with
#    ``B * G * B.transpose()``
#    either::
#
#        [0 1]              [2 1]
#        [1 0] * 2^n  or    [1 2] * 2^n
#
#    EXAMPLES::
#
#        sage: from sage.quadratic_forms.genera.normal_form import _normalize_2x2
#        sage: R = Zp(2, prec = 15, type = 'fixed-mod', print_mode='series', show_prec=False)
#        sage: G = Matrix(R, 2, [-17*2,3,3,23*2])
#        sage: B =_normalize_2x2(G)
#        sage: B * G * B.T
#        [2 1]
#        [1 2]
#
#        sage: G = Matrix(R,2,[-17*4,3,3,23*2])
#        sage: B = _normalize_2x2(G)
#        sage: B*G*B.T
#        [0 1]
#        [1 0]
#
#        sage: G = 2^3 * Matrix(R, 2, [-17*2,3,3,23*2])
#        sage: B = _normalize_2x2(G)
#        sage: B * G * B.T
#        [2^4 2^3]
#        [2^3 2^4]
#    """
#    from sage.rings.all import PolynomialRing
#    from sage.modules.free_module_element import vector

function roots(f::fmpz_mod_poly, p::fmpz, e::Int)
  F = Fac{fmpz}()
  F.unit = one(fmpz)
  F[p] = e
  return roots(f, F)
end
function roots(f::fmpz_mod_poly, fac::Fac{fmpz})
  res = Nemo.fmpz_mod_poly_factor(base_ring(f))
  _fac = Nemo.fmpz_factor()
  for (p, e) in fac
    ccall((:_fmpz_factor_append, libflint), Cvoid, (Ref{Nemo.fmpz_factor}, Ref{fmpz}, UInt), _fac, p, UInt(e))
  end
  ccall((:fmpz_mod_poly_roots_factored, libflint), Cvoid, (Ref{Nemo.fmpz_mod_poly_factor}, Ref{fmpz_mod_poly}, Cint, Ref{Nemo.fmpz_factor}, Ref{Nemo.fmpz_mod_ctx_struct}), res, f, 1, _fac, base_ring(f).ninv)
  _res = Tuple{Nemo.fmpz_mod, Int}[]
  for i in 1:res.num
    g = parent(f)()
    ccall((:fmpz_mod_poly_factor_get_fmpz_mod_poly, libflint), Nothing,
          (Ref{Nemo.fmpz_mod_poly}, Ref{Nemo.fmpz_mod_poly_factor}, Int,
           Ref{Nemo.fmpz_mod_ctx_struct}),
          g, res, i - 1, base_ring(f).ninv)
    e = unsafe_load(res.exp, i)
    push!(_res, (-coeff(g, 0), e))
  end
  return _res
end

function _normalize_twobytwo(G, p)
  # p = fmpz(2)
  R = base_ring(G)
  B = identity_matrix(R, nrows(G))
  P, x = PolynomialRing(R, "x", cached = false)
  odd1 = _val(G[1, 1], p) < _val(G[2, 1], p)
  odd2 = _val(G[2, 2], p) < _val(G[2, 1], p)
  if odd1 || odd2
    error("Not a valid 2 x 2 block.")
  end
  scale = p^(_val(G[1, 2], p))
  D = matrix(R, 2, 2, [divexact(d, R(scale)) for d in G]) # G is symmetric
  # now D is of the form
  # [2a b ]
  # [b  2c]
  # where b has valuation 1.
  G = deepcopy(D)
  # Make sure G[2, 2] has valuation 1.
  if _val(D[2, 2], p) > _val(D[1, 1], p)
    swap_cols!(B, 1, 2)
    swap_cols!(D, 1, 2)
    swap_rows!(D, 1, 2)
  end
  if _val(D[2, 2], p) != 1
    # this works because
    # D[0, 0] has valuation at least 2
    for j in 1:ncols(B)
      B[2, j] = B[2, j] + B[1, j]
    end
    D = B * G * transpose(B)
  end
  @assert _val(D[2, 2], p) == 1

  if mod(lift(det(D)), 8) == 3
#        #  in this case we can transform D to
#        #  2 1
#        #  1 2
#        # Find a point of norm 2
#        # solve: 2 == D[1,1]*x^2 + 2*D[1,0]*x + D[0,0]
    pol = map_coefficients(y -> divexact(y, R(2)), D[2, 2]*x^2 + 2*D[2, 1]*x + D[1, 1] - 2, parent = parent(x))
    sol = roots(pol, p, valuation(modulus(base_ring(G)), p))[1][1]
    @assert pol(sol) == 0
    #pol = divexact(D[2, 2]*x^2 + 2*D[2, 1]*x + D[1, 1] - 2, R(2))
    #sol = _solve_quadratic_dyadic(pol)
    B[1, 2] = sol
    D = B * G * transpose(B)
    # make D[1, 2] = 1
    for j in 1:ncols(B)
      B[2, j] = B[2, j] * inv(D[2, 1])
    end
    D = B * G * transpose(B)

    if D[2, 2] != 2
      v = matrix(P, 1, 2, [x, -2*x + 1])
      pol = map_coefficients(y -> divexact(y, R(2)), (v * change_base_ring(P, D) * transpose(v))[1, 1] - 2, parent = P)
      sol = roots(pol, p, valuation(modulus(R), p))[1][1]
      for j in 1:ncols(B)
        B[2, j] = sol * B[1, j] + (-2 * sol + 1) * B[2, j]
      end
      D = B * G * transpose(B)
    end
    # check the result
    @assert D == matrix(R, 2, 2, [2, 1, 1, 2])
  elseif mod(lift(det(D)), 8) == 7
    # in this case we can transform D to
    #  0 1
    #  1 0
    # Find a point representing 0
    # solve: 0 == D[1,1]*x^2 + 2*D[1,0]*x + D[0,0]
    pol = map_coefficients(y -> divexact(y, R(2)), D[2, 2]*x^2 + 2*D[2,1]*x + D[1,1], parent = P)
    sol = roots(pol, p, valuation(modulus(R), p))[1][1]
    for j in 1:ncols(B)
      B[1, j] = B[1, j] + sol * B[2, j]
    end
    D = B * G * transpose(B)
    # make the second basis vector have 0 square as well.
    for j in 1:ncols(B)
      B[2, j] = B[2, j] - divexact(D[2, 2], 2 * D[1, 2]) * B[1, j]
    end
    D = B * G * transpose(B)
    # rescale to get D[0,1] = 1
    for j in 1:ncols(B)
      B[1, j] = B[1, j] * inv(D[2, 1])
    end
    D = B * G * transpose(B)
    # check the result
    @assert D == matrix(R, 2, 2, [0, 1, 1, 0])
  end
  return B
end

#def _normalize_odd_2x2(G):
#    r"""
#    Normalize this `2` by `2` block.
#
#    INPUT:
#
#    - ``G`` -- a multiple of the `2` by `2` identity_matrix
#      over the `p`-adics for `p` odd.
#
#    OUTPUT:
#
#    - A transformation matrix ``B`` such that
#      ``B * G * B.T`` is the identity matrix
#
#    EXAMPLES::
#
#        sage: from sage.quadratic_forms.genera.normal_form import _normalize_odd_2x2
#        sage: R = Zp(5, type='fixed-mod', print_mode='terse', show_prec=False)
#        sage: G = 2 * Matrix.identity(R, 2)
#        sage: B = _normalize_odd_2x2(G)
#        sage: B*G*B.T
#        [1 0]
#        [0 1]
#    """
function _normalize_odd_twobytwo(G, p)
  @assert G[1, 1] == G[2, 2]
  u = G[1, 1]
  y = zero(base_ring(G))
  while !_issquare(inv(u) - y^2, p)
    y = y + 1
  end
  x = _sqrt(inv(u) - y^2, p)
  B = identity_matrix(base_ring(G), nrows(G))
  B[1, 1] = x
  B[1, 2] = y
  B[2, 1] = y
  B[2, 2] = -x

  D = B * G * transpose(B)
  # check the result
  @assert D == matrix(base_ring(G), 2, 2, [1, 0, 0, 1])
  return B
end

#    B = copy(G.parent().identity_matrix())
#    B[0,0] = x
#    B[0,1] = y
#    B[1,0] = y
#    B[1,1] = -x
#    return B
#
#def _partial_normal_form_of_block(G):
#    r"""
#    Return the partial normal form of the homogeneous block ``G``.
#
#    For internal use in :meth:`_two_adic_normal_forms`.
#
#    INPUT:
#
#    - ``G`` -- a modular symmetric matrix over the `2`-adic integers
#
#    OUTPUT:
#
#    - ``D, B, w`` -- with ``B`` a transformation matrix such that
#      ``B * G * B.T`` is in partial normal form
#      and `w = 0, 1, 2` is the size of the part consisting of forms of type W
#
#    EXAMPLES::
#
#        sage: from sage.quadratic_forms.genera.normal_form import _partial_normal_form_of_block
#        sage: R = Zp(2,prec=4, type = 'fixed-mod',print_mode='terse', show_prec=False)
#        sage: U = Matrix(R, 2, [0,1,1,0])
#        sage: V = Matrix(R, 2, [2,1,1,2])
#        sage: W1 = Matrix(R, 1, [1])
#        sage: W3 = Matrix(R, 1, [3])
#        sage: W5 = Matrix(R, 1, [5])
#        sage: W7 = Matrix(R, 1, [7])
#        sage: G = Matrix.block_diagonal([W1, U, V, W5, V, W3, V, W7])
#        sage: B = _partial_normal_form_of_block(G)[1]
#        sage: B * G * B.T
#        [0 1 0 0 0 0 0 0 0 0 0 0]
#        [1 0 0 0 0 0 0 0 0 0 0 0]
#        [0 0 0 1 0 0 0 0 0 0 0 0]
#        [0 0 1 0 0 0 0 0 0 0 0 0]
#        [0 0 0 0 0 1 0 0 0 0 0 0]
#        [0 0 0 0 1 0 0 0 0 0 0 0]
#        [0 0 0 0 0 0 0 1 0 0 0 0]
#        [0 0 0 0 0 0 1 0 0 0 0 0]
#        [0 0 0 0 0 0 0 0 2 1 0 0]
#        [0 0 0 0 0 0 0 0 1 2 0 0]
#        [0 0 0 0 0 0 0 0 0 0 1 0]
#        [0 0 0 0 0 0 0 0 0 0 0 7]
#        sage: G = Matrix.block_diagonal([W1, U, V, W1, V, W1, V, W7])
#        sage: B = _partial_normal_form_of_block(G)[1]
#        sage: B * G * B.T
#        [0 1 0 0 0 0 0 0 0 0 0 0]
#        [1 0 0 0 0 0 0 0 0 0 0 0]
#        [0 0 0 1 0 0 0 0 0 0 0 0]
#        [0 0 1 0 0 0 0 0 0 0 0 0]
#        [0 0 0 0 0 1 0 0 0 0 0 0]
#        [0 0 0 0 1 0 0 0 0 0 0 0]
#        [0 0 0 0 0 0 0 1 0 0 0 0]
#        [0 0 0 0 0 0 1 0 0 0 0 0]
#        [0 0 0 0 0 0 0 0 0 1 0 0]
#        [0 0 0 0 0 0 0 0 1 0 0 0]
#        [0 0 0 0 0 0 0 0 0 0 3 0]
#        [0 0 0 0 0 0 0 0 0 0 0 7]
#    """
function _partial_normal_form_of_block(G, p)
  D = deepcopy(G)
  n = ncols(G)
  B = identity_matrix(G, n) # transformation matrixj
  blocks = first.(_get_small_block_indices(D))
  # collect the indices of forms of types U, V and W
  U = Int[]
  V = Int[]
  W = Int[]
  for i in blocks
    if i + 1 in blocks || i == n
      push!(W, i)
    else
      if !iszero(D[i, i])
        append!(V, [i, i + 1])
      else
        append!(U, [i, i + 1])
      end
    end
    if length(W) == 3
      # W W W transforms to W U or W V
      #B[W, :] = _relations(D[W, W], 2, p) * B[W, :]
      _BB = _relations(D[W, W], 2, p) * B[W, :]
      for (j, l) in enumerate(W)
        for ll in 1:ncols(B)
          B[l, ll] = _BB[j, ll]
        end
      end
      D = B * G * transpose(B)
      if mod(lift(_unit_part(det(D[W[2:end], W[2:end]]), p)), 8) == 3
        append!(V, W[2:end])
      else
        append!(U, W[2:end])
      end
      W = Int[W[1]]
    end
    if length(V) == 4
      B[V, :] = _relations(D[V, V], 3, p) * B[V, :]
      append!(U, V)
      V = Int[]
      D = B * G * transpose(B)
    end
  end
  # put everything into the right order
  UVW = append!(copy(U), V)
  UVW = append!(UVW, W)
  B = B[UVW, :]
  D = B * G * transpose(B)
  return D, B, length(W)
end
#def _relations(G,n):
#    r"""
#    Return relations of `2`-adic quadratic forms.
#
#    See [MirMor2009]_ IV Prop. 3.2. This function is for internal use only.
#
#    INPUT:
#
#    - ``n`` -- an integer between 1 and 10 -- the number of the relation
#    - ``G`` -- a block diagonal matrix consisting of blocks of types `U, V, W`
#      the left side of the relation. If ``G`` does not match `n` then the
#      results are unpredictable.
#
#    OUTPUT:
#
#    - square matrix ``B`` such that ``B * G * B.T`` is the right side of the
#      relation which consists of blocks of types `U`, `V`, `W` again
function _relations(G, n, p)
  R = base_ring(G)
  if n == 1
    e1 = _unit_part(G[1, 1], p)
    e2 = _unit_part(G[2, 2], p)
    B = matrix(R, 2, 2, [1, 2, 2*e2, -e1])
  elseif n == 2
    e1 = _unit_part(G[1, 1], p)
    e2 = _unit_part(G[2, 2], p)
    e3 = _unit_part(G[3, 3], p)
    B = matrix(R, 3, 3, [1, 1, 1, e2, -e1, 0, e3, 0, -e1])
  elseif n == 3
    B = matrix(R, 4, 4, [1, 1, 1, 0,  1, 1, 0, 1,  1, 0, -1, -1,  0, 1, -1, -1])
  elseif n == 4
    error("Relation 4 is not needed")
  elseif n == 5
    e1 = _unit_part(G[3, 3],  p)
    e2 = _unit_part(G[4, 4],  p)
    if mod(lift(e1), 4) != mod(lift(e2), 4)
      error("W is of the wrong type for relation 5")
    end
    B = matrix(R, 4, 4, [  1,   0,        1,     1,
                           0,   1,        1,     1,
                         -e2, -e2,        0,     3,
                         -e1, -e1, 2*e2 + 3, -2*e1])
  elseif n == 6
    if _val(G[1, 1], p) + 1 != _val(G[2, 2], p)
      error("Wrong scales for relation 6")
    end
    e1 = _unit_part(G[1, 1],  p)
    e2 = _unit_part(G[2, 2],  p)
    B = matrix(R, 2, 2, [1, 1, -2*e2, e1])
  elseif n == 7
    e = _unit_part(G[1, 1], p)
    B = matrix(R, 3, 3, [-3,  e^2,  e^2,  2*e,  1,  0,  2*e,  0,  1])
  elseif n == 8
    e = _unit_part(G[3, 3], p)
    if iszero(G[1, 1])
      B = matrix(R, 3, 3, [e, 0, -1,
                           0, e, -1,
                           2, 2,  1])
    else
      B = matrix(R, 3, 3, [  1,   0,   1,
                             0,   1,   1,
                           2*e, 2*e, - 3])
    end
  elseif n == 9
    e1 = _unit_part(G[1,1], p)
    e2 = _unit_part(G[2,2], p)
    e3 = _unit_part(G[3,3], p)
    B = matrix(R, 3, 3, [1, 0, 1,
                         2*e3, 1,
                        -e1, -2*e2*e3, 2*e1^2*e3 + 4*e1*e3^2, e1*e2])
  elseif n == 10
    e1 = _unit_part(G[1, 1],  p)
    e2 = _unit_part(G[2, 2],  p)
    B = matrix(R, 2, 2, [1, 1, -4*e2, e1])
  else
    error("Relation ($n) must be between 1 and 10")
  end
  D, B1 = _normalize(B * G * transpose(B), p)
  return B1 * B
end

#def _two_adic_normal_forms(G, partial=False):
#    r"""
#    Return the 2-adic normal form of a symmetric matrix.
#
#    INPUT:
#
#    - ``G`` -- block diagonal matrix with blocks of type `U`, `V`, `W`
#    - ``partial`` -- bool (defaul: ``False``)
#
#    OUTPUT:
#
#    - ``D``, ``B`` -- such that ``D = B * G * B.T``
function _two_adic_normal_forms(G, p; partial = false)
  RR = base_ring(G)
  B = identity_matrix(G, nrows(G))
  h, scales = _get_homogeneous_block_indices(G, p)
  push!(h, ncols(B) + 1)
  # UVlist[k] is a list of indices of the block of scale p^k.
  # It contains the indices of the part of types U or V.
  # So it may be empty.
  UVlist = Vector{Int}[Int[], Int[]] # empty lists are appended to avoid special cases.
  Wlist = Vector{Int}[Int[],Int[]]
  # homogeneous normal form for each part
  for k in 0:(scales[end] - scales[1])
    if k + scales[1] in scales
      i = findfirst(isequal(k + scales[1]), scales)
      @assert i !== nothing
      Gk = G[h[i]:h[i+1]-1, h[i]:h[i+1]-1]
      Dk, Bk, wk = _partial_normal_form_of_block(Gk, p)
      B[h[i]:h[i+1]-1, :] = Bk * B[h[i]:h[i+1]-1, :]
      if !partial
        Dk, B1k = _homogeneous_normal_form(Dk, wk, p)
        B[h[i]:h[i+1]-1,:] = B1k * B[h[i]:h[i+1]-1, :]
      end
      push!(UVlist, collect(h[i]:(h[i + 1] - wk - 1)))
      push!(Wlist, collect((h[i + 1] - wk):h[i + 1] - 1))
    else
      push!(UVlist, Int[])
      push!(Wlist, Int[])
    end
  end
  D = B * G * transpose(B)

  if partial
    return D, B
  end

  # use relations descending in k
  # we never leave partial normal form
  # but the homogeneous normal form may be destroyed
  # it is restored at the end.
  for k in length(UVlist):-1:3
    # setup notation
    W = Wlist[k]
    Wm = Wlist[k - 1]
    Wmm = Wlist[k - 2]
    UV = UVlist[k]
    UVm = UVlist[k - 1]
    UVlistk = UVlist[k]
    # UVlistk is either empty or contains at least two elements
    if length(UVlistk) == 0
      V = Int[]
    else
      V = UVlist[k][end-1:end]
    end
    if length(V) != 0 && D[V[1], V[1]] == 0
      V = Int[]
      # condition b)
    end
    if length(Wm) != 0
      if length(V) == 2
        R = append!([Wm[1]], V)
        B[R, :] = _relations(D[R, R], 7, p) * B[R, :]
        V = Int[]
        D = B * G * transpose(B)
      end
      E = [RR(3), RR(7)]
      for w in W
        if _unit_part(D[w, w], p) in E
          R = [Wm[1], w]
          B[R, :] = _relations(D[R, R], 6, p) * B[R, :]
          D = B * G * transpose(B)
        end
      end
    end
    # condition c)
    # We want type a or W = []
    # modify D[w,w] to go from type b to type a
    x = append!([length(V)], [(mod(lift(_unit_part(w, p)),8)) for w in diagonal(D[W,W])])

    if length(x) == 3 && x[2] > x[3]
      x[2], x[3] = x[3], x[2]
    end
    # the first entry of x is either
    # 0 if there is no type V component or
    # 2 if there is a single type V component
    # a = [[0,1], [2,3], [2,5], [0,7], [0,1,1], [2,1,3], [0,7,7], [0,1,7]]
    b = [[0,5], [2,7], [2,1], [0,3], [0,1,5], [2,1,7], [0,3,7], [0,1,3]]
    if x in b
      w = W[end]
      if x == [0,3,7]
        # relation 10 should be applied to 3 to stay in homogeneous normal form
        w = W[1]
      end
      if length(UVm) > 0
        R = push!(UVm[end-1:end], w)
        B[R,:] = _relations(D[R,R], 8, p) * B[R,:]
      elseif length(Wmm) > 0
        R = push!([Wmm[1]], w)
        B[R,:] = _relations(D[R,R],10,p) * B[R,:]
      elseif length(Wm) == 2
        e0 = _unit_part(D[Wm,Wm][1,1], p)
        e1 = _unit_part(D[Wm,Wm][2,2], p)
        if mod(lift(e1-e0),4) == 0
          R = push!(copy(Wm), w)
          B[R,:] = _relations(D[R,R], 9, p) * B[R,:]
        end
      end
    end
    D = B * G * transpose(B)
    # condition a) - stay in homogeneous normal form
    R = append!(copy(UV), W)
    Dk = D[R,R]
    _, Bk = _homogeneous_normal_form(Dk, length(W), p)
    #B[R,:] = Bk * B[R,:]
    _B = Bk * B[R, :]
    for (l, i) in enumerate(R)
      for j in 1:ncols(B)
        B[i, j] = _B[l, j]
      end
    end
    D = B * G * transpose(B)
    # we need to restore the homogeneous normal form of  k-1
    if length(Wm) > 0
      R = append!(copy(UVm), Wm)
      Dkm = D[R,R]
      _, Bkm = _homogeneous_normal_form(Dkm, length(Wm), p)
      #B[R,:] = Bkm * B[R,:]
      BB = Bkm * B[R,:]
      for (i, r) in enumerate(R)
        for j in 1:ncols(B)
          B[r, j] = BB[i, j]
        end
      end

      D = B * G * transpose(B)
    end
  end
  return D, B
end

function Base.setindex!(A::T, B::T, r::Vector{Int}, ::Colon) where {T <: MatElem}
  for (cnt, i) in enumerate(r)
    for j in 1:ncols(A)
      A[i, j] = B[cnt, j]
    end
  end
  return A
end

################################################################################
#
#  Helper function to treat ZZ/p^k as ZZ_p
#
################################################################################

_val(x::Nemo.fmpz_mod, y::IntegerUnion) = iszero(x) ? inf : valuation(lift(x), y)


text:
# Here Engheim  makes use of multiple dispatch so its not
# exactly idiomatic OOP, where the methods are associated with the struct
# Tou can't do "real" OOP in Julia and that is perfectly alright.
# https://docs.julialang.org/en/v1/manual/types/#Composite-Types-1
# https://discourse.julialang.org/t/is-julias-way-of-oop-superior-to-c-python-why-julia-doesnt-use-class-based-oop/52058
# https://arstechnica.com/science/2020/10/the-unreasonable-effectiveness-of-the-julia-programming-language/

abstract type Cipher end

function encrypt(cipher::Cipher, char::Char)
    error("Implement encrypt(::", typeof(cipher), ", char)")
end

function decrypt(cipher::Cipher, char::Char)
    error("Implement decrypt(::", typeof(cipher), ", char)")
end

# once you define the encrypt/decrypt pair for the cipher
# these functions handle strings for you
function encrypt(cipher::Cipher, message::AbstractString)
#    map(ch -> encrypt(cipher, ch), message)
    map(message) do ch
        encrypt(cipher, ch)
    end
end

function decrypt(cipher::Cipher, ciphertext::AbstractString)
#    map(ch -> decrypt(cipher, ch), ciphertext)
    map(ciphertext) do ch
        decrypt(cipher, ch)
    end
end

# Caeser cipher implementation - same as before, except shift is in strucct

struct CaesarCipher <: Cipher
  shift::Int
end

const n = length('A':'Z')

function encrypt(cipher::CaesarCipher, ch::Char)
    ch in 'A':'Z' ? 'A' + mod((ch - 'A') + cipher.shift, n) :
    ch in 'a':'z' ? 'a' + mod((ch - 'a') + cipher.shift, n) :
    #= otherwise =# ch
end

function decrypt(cipher::CaesarCipher, ch::Char)
    ch in 'A':'Z' ? 'A' + mod((ch - 'A') - cipher.shift, n) :
    ch in 'a':'z' ? 'a' + mod((ch - 'a') - cipher.shift, n) :
    #= otherwise =# ch
end

# Substitution cipher

# Engheim did it this way, but then you have to figure out the substitute yourself
# struct SubstitutionCipher <: Cipher
#   substitute::Dict{Char, Char}
#   alphabet::Dict{Char, Char}

#   function SubstitutionCipher(substitute)
#     sub   = Dict(zip('A':'Z', collect(substitute)))
#     alpha = Dict(zip(collect(substitute), 'A':'Z'))
#     new(sub, alpha)
#   end
#end
# function encrypt(cipher::SubstitutionCipher, ch::Char)
#   get(cipher.substitute, ch, ch)
# end

# function decrypt(cipher::SubstitutionCipher, ch::Char)
#     get(cipher.alphabet, ch, ch)
# end

# this alternate way creates and stores the shuffle for you
using Random
struct SubstitutionCipher <: Cipher
    upshuf::Array{Char}
    lowshuf::Array{Char}
    updict::Dict{Char, Char}
    lowdict::Dict{Char, Char}
    rupdict::Dict{Char, Char}
    rlowdict::Dict{Char, Char}

    function SubstitutionCipher()
        # not sure you have to save ths shuffle since you can get them out of the dict
        # still its convenient
        upshuf = shuffle('A':'Z')
        lowshuf = shuffle('a':'z')
        updict = Dict(zip('A':'Z', upshuf))
        lowdict = Dict(zip('a':'z', lowshuf))
        rupdict = Dict(values(updict) .=> keys(updict))
        rlowdict = Dict(values(lowdict) .=> keys(lowdict))
        new(upshuf, lowshuf, updict, lowdict, rupdict, rlowdict)
    end
end

function encrypt(cipher::SubstitutionCipher, ch::Char)
    ch in 'A':'Z' ? get(cipher.updict, ch, ch) :
    ch in 'a':'z' ? get(cipher.lowdict, ch, ch) :
    #= otherwise =# ch
end

function decrypt(cipher::SubstitutionCipher, ch::Char)
    ch in 'A':'Z' ? get(cipher.rupdict, ch, ch) :
    ch in 'a':'z' ? get(cipher.rlowdict, ch, ch) :
    #= otherwise =# ch
end

# example
# julia> cipher = SubstitutionCipher()
# SubstitutionCipher(['U', 'V', 'J', 'L', 'C', 'F', 'D', 'P', 'N', 'X'  …  'Z', 'O', 'M', 'B', 'A', 'I', 'Q', 'Y', 'W', 'T'], ['s', 'm', 'q', 'b', 'u', 'o', 'y', 'g', 'x', 'e'  …  'h', 'j', 'c', 'a', 'w', 'd', 't', 'f', 'k', 'n'], Dict('E' => 'C', 'Z' => 'T', 'X' => 'Y', 'B' => 'V', 'C' => 'J', 'D' => 'L', 'A' => 'U', 'R' => 'O', 'G' => 'D', 'F' => 'F'…), Dict('n' => 'r', 'f' => 'o', 'w' => 't', 'd' => 'b', 'e' => 'u', 'o' => 'p', 'h' => 'g', 'j' => 'e', 'i' => 'x', 'k' => 'i'…), Dict('E' => 'O', 'Z' => 'Q', 'X' => 'J', 'C' => 'E', 'B' => 'T', 'D' => 'G', 'A' => 'U', 'R' => 'K', 'G' => 'N', 'O' => 'R'…), Dict('n' => 'z', 'f' => 'x', 'w' => 'u', 'd' => 'v', 'e' => 'j', 'o' => 'f', 'j' => 'r', 'h' => 'q', 'i' => 'k', 'r' => 'n'…))

# julia> encrypt(cipher, 'I')
# 'N': ASCII/Unicode U+004E (category Lu: Letter, uppercase)

# julia> decrypt(cipher, 'N')
# 'I': ASCII/Unicode U+0049 (category Lu: Letter, uppercase)

# julia> encrypt(cipher, 'i')
# 'x': ASCII/Unicode U+0078 (category Ll: Letter, lowercase)

# julia> decrypt(cipher, 'x')
# 'i': ASCII/Unicode U+0069 (category Ll: Letter, lowercase)

# julia> encrypt(cipher, '!')
# '!': ASCII/Unicode U+0021 (category Po: Punctuation, other)

# julia> encrypt(cipher, "This is the message!")
# "Bgxc xc agu luccsyu!"

# julia> decrypt(cipher, "Bgxc xc agu luccsyu!")
# "This is the message!"

# Cipher agnostic password keeping service
mutable struct Vault
    passwords::Dict{String, String}
    cipher::Cipher
end

function Vault(cipher::Cipher)
    Vault(Dict{String, String}(), cipher)
end

function addlogin!(vault::Vault, login::AbstractString, password::AbstractString)
    vault.passwords[login] = encrypt(vault.cipher, password)
end

function getpassword(vault::Vault, login::AbstractString)
    decrypt(vault.cipher, vault.passwords[login])
end

# example
# julia> vault = Vault(CaesarCipher(23))
# Vault(Dict{String, String}(), CaesarCipher(23))

# julia> addlogin!(vault, "google", "BING")
# "YFKD"

# julia> addlogin!(vault, "amazon", "SECRET")
# "PBZOBQ"

# julia> getpassword(vault, "google")
# "BING"

# julia> getpassword(vault, "amazon")
# "SECRET"

# julia> cipher = SubstitutionCipher()
# SubstitutionCipher(['G', 'M', 'L', 'V', 'A', 'N', 'D', 'S', 'E', 'C'  …  'I', 'K', 'O', 'H', 'B', 'X', 'J', 'T', 'W', 'Z'], ['n', 'a', 'i', 'z', 'o', 'l', 'v', 'g', 'h', 'k'  …  'b', 'u', 'q', 'f', 'e', 'p', 'r', 'w', 'm', 'y'], Dict('E' => 'A', 'Z' => 'Z', 'X' => 'T', 'B' => 'M', 'C' => 'L', 'D' => 'V', 'A' => 'G', 'R' => 'K', 'G' => 'D', 'F' => 'N'…), Dict('n' => 'x', 'f' => 'l', 'w' => 'r', 'd' => 'z', 'e' => 'o', 'o' => 'd', 'h' => 'g', 'j' => 'k', 'i' => 'h', 'k' => 's'…), Dict('E' => 'I', 'Z' => 'Z', 'X' => 'V', 'C' => 'J', 'B' => 'U', 'D' => 'G', 'A' => 'E', 'R' => 'N', 'G' => 'A', 'N' => 'F'…), Dict('n' => 'a', 'f' => 't', 'w' => 'x', 'd' => 'o', 'e' => 'u', 'o' => 'e', 'h' => 'i', 's' => 'k', 'i' => 'c', 'r' => 'w'…))

# julia> vault2 = Vault(cipher)
# Vault(Dict{String, String}(), SubstitutionCipher(['G', 'M', 'L', 'V', 'A', 'N', 'D', 'S', 'E', 'C'  …  'I', 'K', 'O', 'H', 'B', 'X', 'J', 'T', 'W', 'Z'], ['n', 'a', 'i', 'z', 'o', 'l', 'v', 'g', 'h', 'k'  …  'b', 'u', 'q', 'f', 'e', 'p', 'r', 'w', 'm', 'y'], Dict('E' => 'A', 'Z' => 'Z', 'X' => 'T', 'B' => 'M', 'C' => 'L', 'D' => 'V', 'A' => 'G', 'R' => 'K', 'G' => 'D', 'F' => 'N'…), Dict('n' => 'x', 'f' => 'l', 'w' => 'r', 'd' => 'z', 'e' => 'o', 'o' => 'd', 'h' => 'g', 'j' => 'k', 'i' => 'h', 'k' => 's'…), Dict('E' => 'I', 'Z' => 'Z', 'X' => 'V', 'C' => 'J', 'B' => 'U', 'D' => 'G', 'A' => 'E', 'R' => 'N', 'G' => 'A', 'N' => 'F'…), Dict('n' => 'a', 'f' => 't', 'w' => 'x', 'd' => 'o', 'e' => 'u', 'o' => 'e', 'h' => 'i', 's' => 'k', 'i' => 'c', 'r' => 'w'…)))

# julia> addlogin!(vault2, "amazon", "SECRET")
# "OALKAH"

# julia> addlogin!(vault2, "apple", "JONAGOLD")
# "CPRGDPYV"

# julia> getpassword(vault2, "amazon")
# "SECRET"

# julia> getpassword(vault2, "apple")
# "JONAGOLD"

ext: jl

max_stars_repo_name: FourMInfo/JuliaforBeginners

max_issues_repo_name: FourMInfo/JuliaforBeginners

max_forks_repo_name: FourMInfo/JuliaforBeginners


text:
using Documenter, AxisArrays

makedocs(
    modules = [AxisArrays],
    sitename = "AxisArrays",
)

deploydocs(
    deps   = Deps.pip("mkdocs", "python-markdown-math"),
    repo   = "github.com/JuliaArrays/AxisArrays.jl.git"
)


text:
function floyd_warshall(d, n, f, g)
	for k in 1 : n
		for i in 1 : n
			for j in 1 : n
				d[i, j] = f(d[i, j], g(d[i, k], d[k, j]))
			end
		end
	end
	return d
end

function transitive_closure(d, n)
	d .= reshape([i < Inf ? 1 : 0 for i in d], (n, n))
	or(a, b) = (a != 0 || b != 0) ? 1 : 0
	and(a, b) = (a != 0 && b != 0) ? 1 : 0
	return floyd_warshall(d, n, or, and)
end

function create_preference_matrix(b, v, c)
	ind(c)::Int = c - 'A' + 1
	chk_list::Array{Bool, 1} = Array{Bool, 1}(undef, c)
	ret::Array{Int, 2} = fill(0, c, c)
	for i in b
		fill!(chk_list, false)
		for j in i
			chk_list[ind(j)] = true
			for k in 1 : c
				if !chk_list[k]
					ret[ind(j), k] += 1
				end
			end
		end
	end
	return ret
end

function find_strongest_paths(p, c)
	for i in 1 : c - 1
		for j in i + 1 : c
			if p[i, j] < p[j, i]
				p[i, j] = 0
			else
				p[j, i] = 0
			end
		end
	end
	floyd_warshall(p, c, max, min)
	for i in 1 : c
		p[i, i] = 0
	end
	return p
end

find_schulze_ranking(sp, c) = String([i[2] for i in sort([(sum([sp[j, i] > sp[i, j] ? 0 : 1 for i in 1 : c]), Char(j + 64)) for j in 1 : c])])


text:
using SpikingNetworks

parseParameters"""
σ = 10
ρ = 28
β = 8/3"""

parseEquations"""
dx/dt = σ*(y-x)
dy/dt = x*(ρ-z)-y
dz/dt = x*y-β*z"""

net=init"""
x=[1]
y=[1]
z=[1]"""

stateMonitor"[1]"

net.dt=0.01
@time run!(net,20/net.dt)
using MatlabPlot
mplot3(net.xRecord[1,:],net.yRecord[1,:],net.zRecord[1,:])


text:
function cmatchregret(I::Int64, r::Vector{Vector{Float64}}, gs::GameSet)
    ni_stage, na_stage = gs.ni_stage, gs.na_stage
    na = numactions(I, ni_stage, na_stage)
    σ = Vector{Float64}(undef, na)
    for a in 1:na
        denom = sum(max(r[I][b], 0.0) for b in 1:na)
        if denom > 0.0
            σ[a] = max(r[I][a], 0.0) / denom
        else
            σ[a] = 1.0 / na
        end
    end
    return σ
end

function cupdateutility!(player::Int64, I::Int64, h::SVector, σ::Vector{Float64},
        u::Float64, πi::Float64, πo::Float64,
        r::Vector{Vector{Float64}}, s::Vector{Vector{Float64}},
        g::AirDefenseGame, gs::GameSet, depth::Int64,
        α::Number, β::Number, γ::Number, iter::Int, siter::Vector{Int}, discounted::Bool,
        λl::Vector{Vector{Float64}}, λu::Vector{Vector{Float64}})
    A, ni_stage, na_stage = gs.A, gs.ni_stage, gs.na_stage
    na = numactions(I, ni_stage, na_stage)
    m = Vector{Float64}(undef, na)
    nextactions = actions(depth, A)
    # @show na[I]
    for a in 1:na
        ha = setindex(h, nextactions[a], depth)
        if getplayer(depth) == player
            πip = σ[a] * πi
            up = cupdatetree!(ha, player, πip, πo, r, s, g, gs, depth + 1, α, β, γ, iter, siter, discounted, λl, λu)
            if getplayer(depth) == 2 # constrained player
                lagrange = λl[I][a] * -1 + λu[I][a] * 1 # assumes indpendent constraints
            else
                lagrange = 0.0
            end
            m[a] = up - lagrange
            u = u + σ[a] * up
        else
            πop = σ[a] * πo
            up = cupdatetree!(ha, player, πi, πop, r, s, g, gs, depth + 1, α, β, γ, iter, siter, discounted, λl, λu)
            u = u + up # possibly need a lagrange adjustment here? but I don't think so
        end
    end
    return m, u
end

function cupdateregret!(player::Int64, I::Int64, h::SVector, σ::Vector{Float64},
        u::Float64, m::Vector{Float64}, πi::Float64,
        r::Vector{Vector{Float64}}, s::Vector{Vector{Float64}}, gs::GameSet, depth ::Int64,
        α::Number, β::Number, γ::Number, iter::Int, siter::Vector{Int}, discounted::Bool)
    ni_stage, na_stage = gs.ni_stage, gs.na_stage
    if getplayer(depth) == player
        siter[I] += 1
        for a in 1:numactions(I, ni_stage, na_stage)
            rnew = r[I][a] + m[a] - u
            if discounted
                # r[I][a] = rnew >= 0 ? rnew * (iter^α / (iter^α + 1)) : rnew * (iter^β / (iter^β + 1))
                # s[I][a] = (s[I][a] + πi * σ[a]) * (iter / (iter + 1))^γ
                if rnew >= 0
                    r[I][a] = rnew * (siter[I]^α / (siter[I]^α + 1))
                else
                    r[I][a] = rnew * (float(siter[I])^β / (float(siter[I])^β + 1))
                end
                s[I][a] = (s[I][a] + πi * σ[a]) * (siter[I] / (siter[I] + 1))^γ
            else
                r[I][a] = rnew
                s[I][a] = s[I][a] + πi * σ[a]
            end
        end
    end
end

function cupdatetree!(h::SVector, player::Int64, πi::Float64, πo::Float64,
    r::Vector{Vector{Float64}}, s::Vector{Vector{Float64}}, g::AirDefenseGame,
    gs::GameSet, depth::Int64, α::Number, β::Number, γ::Number, iter::Int, siter::Vector{Int}, discounted::Bool,
    λl::Vector{Vector{Float64}}, λu::Vector{Vector{Float64}})
    # πo is short for \vec{π}_{-i}
    A = gs.A
    (terminal(depth)) && (return leafutility(h, player, πi, πo, g, gs))
    if getplayer(depth) == 3
        as = actions(depth, A)
        chance_weights = pweights([chance(depth, a, g, gs) for a in as])
        a = as[sample(chance_weights)]
        ha = setindex(h, a, depth)
        return cupdatetree!(ha, player, πi, πo, r, s, g, gs, depth + 1, α, β, γ, iter, siter, discounted, λl, λu)
    end
    I = infoset(h, depth, g, gs)
    u = 0.0
    σ = cmatchregret(I, r, gs)
    m, u = cupdateutility!(player, I, h, σ, u, πi, πo, r, s, g, gs, depth, α, β, γ, iter, siter, discounted, λl, λu)
    cupdateregret!(player, I, h, σ, u, m, πi, r, s, gs, depth, α, β, γ, iter, siter, discounted)
    return u
end

function strat_ccfr(s::Vector{Vector{Float64}})
    σ = [similar(si) for si in s]
    for i in eachindex(s)
        denom = sum(s[i])
        if denom == 0
            len = length(s[i])
            σ[i] .= fill(1 / len, len)
        else
            σ[i] .= s[i] ./ denom
        end
    end
    return σ
end

function ccfr_full(g::AirDefenseGame, gs::GameSet, rp2lb, rp2ub;
        iterlimit::Int = 100_000,
        timelimit::Number = 600, tol::Float64 = 5e-5,
        α::Number = 1.5, β::Number = 0.5, γ::Number = 2, discounted::Bool = false,
        λmax::Number = 1_000, λscale::Number = 4)
    ni, ni_stage, na_stage, players = gs.ni, gs.ni_stage, gs.na_stage, gs.players
    A, _, na_stage = build_action(g)
    ni, ni1, ni2, ni_stage = build_info(g)
    ns1, ns2, ns1_stage, ns2_stage = build_nseq(g, na_stage, ni_stage)
    _, _, (seqI1, seqI2), (nextseqI1, nextseqI2) = build_Iset(g, A, na_stage, ns1_stage, ns2_stage, ni1, ni2, ns1, ns2)
    seqvals = (ns1, ns2, seqI1, seqI2, nextseqI1, nextseqI2)
    r = [zeros(numactions(i, ni_stage, na_stage)) for i in 1:ni] # regret
    s = [zeros(numactions(i, ni_stage, na_stage)) for i in 1:ni] # cumulative strategies
    λl = [zeros(numactions(i, ni_stage, na_stage)) for i in 1:ni]
    λu = [zeros(numactions(i, ni_stage, na_stage)) for i in 1:ni]
    λmean = zeros(iterlimit ÷ 10)
    siter = zeros(Int64, ni) # infoset visit counts
    u = 0.0
    n_hist = 10 # number of intermediate strategies to record (for best responses)
    Tx = round(Int64, timelimit * 60 / n_hist) # frequency of recording history
    tx = 1 # index for sigma calculations
    iter = 1 # index for number of iter
    lind = 1 # index for lambda tracking
    u1deriv = 0
    u1prev = 0
    converged = false
    u1 = Float64[]
    σ1 = Vector{Vector{Vector{Float64}}}(undef, n_hist)
    println("CCFR: start $timelimit min or $iterlimit iter at $(now())...")
    tstart, t_hist = now(), now()
    # lamcheck = []
    while !converged && tominute(now() - tstart) < timelimit && iter < iterlimit
        ξ = (λmax ^ λscale) / (1 * sqrt(iter)) # λ learning rate β / (G √T)
        Ψ = real_strat_struct_combined(strat_ccfr(s), gs, seqvals) # could be faster if we only update for player 2
        for I in 1:size(rp2lb, 1)
            for a in 1:length(rp2lb[I])
                λl[I][a] = clamp(λl[I][a] + ξ * (-Ψ[I][a] + rp2lb[I][a]), 0.0, λmax)
                λu[I][a] = clamp(λu[I][a] + ξ * (Ψ[I][a] - rp2ub[I][a]), 0.0, λmax)
            end
        end
        # push!(lamcheck, deepcopy(λl))
        # iter % 100 == 0 && @show λl[1], λu[1]
        for player in players
            h0 = SVector(0, 0, 0, 0, 0, 0)
            u = cupdatetree!(h0, player, 1.0, 1.0, r, s, g, gs, 1, α, β, γ, iter, siter, discounted, λl, λu)
            (player == 1) && push!(u1, u)
        end
        if tosecond(now() - t_hist) > Tx # calculate sigma every Tx seconds
            σ1[tx] = strat_ccfr(s)
            tx += 1
            t_hist = now()
        end
        u1deriv = u1deriv * (1 - 1 / iter) + (u1[iter] - u1prev) * (1 / iter)
        converged = abs(u1deriv) < tol
        u1prev = u1[iter]
        # if iter % 10 == 0
        #     lammean = mean(mean.(λu))
        #     lamcount = sum(sum(x .> 0) for x in λu)
        #     λmean[lind] = lammean
        #     print("\rIter: $iter, λ mean: $(round(lammean, digits = 6)), λ count: $lamcount")
        #     lind += 1
        # end
        iter % 10 == 0 && print("\rIter: $iter")
        iter += 1
    end
    println("")
    σ = strat_ccfr(s)
    # return u1, r, s, σ, σ1, converged, iter, λmean
    return u1, r, s, σ, σ1, converged, iter
end

max_stars_repo_name: ajkeith/Cyber-Air-Defense

max_issues_repo_name: ajkeith/Cyber-Air-Defense


text:
import ONNXRunTime
function onnxruntime_infer(f, inputs...)

	reversedims(a::AbstractArray{T,N}) where {T, N} = permutedims(a, N:-1:1)

	mktempdir() do dir
		modelfile = joinpath(dir, "model.onnx")
		save(modelfile, f, size.(inputs)...)

		model = ONNXRunTime.load_inference(modelfile)
		return model(Dict(ONNXRunTime.input_names(model) .=> reversedims.(inputs))) |> values .|> reversedims |> Tuple
	end

end


text:
abstract type TA_Gaugefields_2D{NC} <: TA_Gaugefields{NC,2}
end

include("./TA_gaugefields_2D_serial.jl")

function TA_Gaugefields(NC,NX,NT;mpi=false)
    if mpi
        return TA_Gaugefields_2D_mpi(NC,NX,NT)
        #error("mpi = $mpi is not supoorted")
    else
        return TA_Gaugefields_2D_serial(NC,NX,NT)
    end
end

function clear_U!(U::Array{T,1}) where T <: TA_Gaugefields_2D
    for μ=1:2
        clear_U!(U[μ])
    end
end


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

include("density_with_diff.jl")
include("grad_func.jl")
include("forwarddiff.jl")
include("zygotediff.jl")


text:
# @author b33rNc0d3
# @date 2020-06-10
# Using Turing to implement JAGS model from Coursera Bayesian Statistics: Techniques and Models

using Turing
using Distributions
using StatsPlots
using Random

Random.seed!(50)

@model TuringJagsModel(y) = begin
    μ ~ TDist(1.0)
    σ = 1.0
    N = length(y)

    for i = 1 : N
        y[i] ~ Normal(μ, σ)
    end
end;

y = [1.2, 1.4, -0.5, 0.3, 0.9, 2.3, 1.0, 0.1, 1.3, 1.9]

iterations = 1000
ϵ = 0.05
τ = 10

chain = sample( TuringJagsModel(y), HMC(ϵ, τ), iterations)

StatsPlots.plot(chain)

print(chain)

plot(chain)


text:
module StaticNumbers

# TODO: This module should be renamed StaticNumbers, and the word `static`
# replaced by `static` everywhere.

export Static, static,
       StaticInteger, StaticReal, StaticNumber, StaticOrInt, StaticOrBool,
       @staticnumbers, ofstatictype

const StaticError = ErrorException("Illegal type parameter for Static.")

"""
A `StaticInteger` is an `Integer` whose value is stored in the type, and which
contains no runtime data.
"""
struct StaticInteger{X} <: Integer
    function StaticInteger{X}() where {X}
        X isa Integer && !(X isa Static) && isimmutable(X) || throw(StaticError)
        new{X}()
    end
end
StaticInteger{X}(x::Integer) where {X} = x==X ? StaticInteger{X}() : throw(InexactError(:StaticInteger, StaticInteger{X}, x))
Base.@pure StaticInteger(x::Integer) = StaticInteger{x}()

"""
A `StaticReal` is a `Real` whose value is stored in the type, and which
contains no runtime data.
"""
struct StaticReal{X} <: Real
    function StaticReal{X}() where {X}
        X isa Real && !(X isa Integer) && !(X isa Static) && isimmutable(X) || throw(StaticError)
        new{X}()
    end
end
StaticReal{X}(x::Real) where {X} = x==X ? StaticReal{X}() : throw(InexactError(:StaticReal, StaticReal{X}, x))
Base.@pure StaticReal(x::Real) = StaticReal{x}()

"""
A `StaticNumber` is a `Number` whose value is stored in the type, and which
contains no runtime data.
"""
struct StaticNumber{X} <: Number
    function StaticNumber{X}() where {X}
        X isa Number && !(X isa Real) && !(X isa Static) && isimmutable(X) || throw(StaticError)
        new{X}()
    end
end
StaticNumber{X}(x::Number) where {X} = x==X ? StaticNumber{X}() : throw(InexactError(:StaticNumber, StaticReal{X}, x))
Base.@pure StaticNumber(x::Number) = StaticNumber{x}()

"""
`Static{X}` is short-hand for the `Union` of `StaticInteger{X}`, `StaticReal{X}`
and `StaticNumber{X}`.
"""
const Static{X} = Union{StaticInteger{X}, StaticReal{X}, StaticNumber{X}}

# We'll define this constructor, but not recommend it.
Static{X}() where X = static(X)

# This is the recommended constructor.
"""
`static(X)` is shorthand for `StaticInteger{X}()`, `StaticReal{X}()` or `StaticNumber{X}()`,
depending on the type of `X`.
"""
static(x::StaticInteger) = x
static(x::StaticReal) = x
static(x::StaticNumber) = x
Base.@pure static(x::Integer) = StaticInteger(x)
Base.@pure static(x::Real) = StaticReal(x)
Base.@pure static(x::Number) = StaticNumber(x)
static(x::Irrational) = x # These are already defined by their type.

# Functions that take only `Int` may be too restrictive.
# The StaticOrInt type union is often a better choice.
const StaticOrInt = Union{StaticInteger, Int}

# Promotion
Base.promote_rule(::Type{<:Static{X}}, ::Type{<:Static{X}}) where {X} =
    typeof(X)
Base.promote_rule(::Type{<:AbstractIrrational}, ::Type{<:Static{X}}) where {X} =
    promote_type(Float64, typeof(X))

# We need to override promote and promote_typeof because they don't even call
# promote_rule for all-same types.
for T in (StaticInteger, StaticReal, StaticNumber)
    @eval Base.promote(::$T{X}, ys::$T{X}...) where {X} = ntuple(i->X, 1+length(ys))
    @eval Base.promote_type(::Type{$T{X}}, ::Type{$T{X}}) where {X} = typeof(X)
    @eval Base.promote_typeof(::$T{X}, ::$T{X}...) where {X} = typeof(X)
    # To avoid infinite recursion, we need this:
    @eval Base.promote_type(::Type{$T{X}}, S::Type...) where {X} = promote_type(typeof(X), promote_type(S...))
end

Base.promote_rule(::Type{<:Static{X}}, ::Type{<:Static{Y}}) where {X,Y} =
    promote_type(typeof(X),typeof(Y))

Base.promote_rule(::Type{<:Static{X}}, ::Type{T}) where {X,T<:Number} =
    promote_type(typeof(X), T)

# Bool has a special rule that we need to override?
#Base.promote_rule(::Type{Bool}, ::Type{StaticInteger{X}}) where {X} = promote_type(Bool, typeof(X))

function Base.convert(T::Type{<:Static{X}}, y::Number) where {X}
    X == y || throw(InexactError(:convert, T, y))
    return T()
end

Base.convert(::Type{T}, ::Static{X}) where {T<:Number,X} = convert(T, X)

"ofstatictype(x,y) - like oftype(x,y), but return a `Static` `x` is a `Static`."
ofstatictype(::Static{X}, y) where {X} = static(oftype(X, y))
ofstatictype(x, y) = oftype(x, y)

# TODO: Constructors to avoid Static{Static}

# Some of the more common constructors that do not default to `convert`
# TODO:  We should have a (::Type{T})(x::Static) where {T<:Number} constructor
# instead of all of these.
# Need to figure out a way to avoid ambiguities.
for T in (:Bool, :Int32, :UInt32, :Int64, :UInt64, :Int128, :Integer)
    @eval Base.$T(::StaticInteger{X}) where X = $T(X)
end
(::Type{T})(x::Union{StaticReal{X}, StaticInteger{X}}) where {T<:AbstractFloat, X} = T(X)
(::Type{T})(x::Static{X}) where {T<:Complex, X} = T(X)
for T in (:ComplexF32, :ComplexF64, :Complex)
    @eval Base.$T(::Static{X}) where X = $T(X)
end
Rational{T}(::Union{StaticInteger{X}, StaticReal{X}}) where {T,X} = Rational{T}(X)
Complex{T}(::Static{X}) where {T,X} = Complex{T}(X)
# big(x) still defaults to convert.

# Single-argument functions that do not already work.
# Note: We're not attempting to support every function in Base.
# TODO: Should have a macro for easily extending support.
for fun in (:-, :zero, :one, :oneunit, :trailing_zeros, :widen, :decompose)
    @eval Base.$fun(::Static{X}) where X = Base.$fun(X)
end
for fun in (:trunc, :floor, :ceil, :round)
    @eval Base.$fun(::Union{StaticReal{X}, StaticNumber{X}}) where {X} = Base.$fun(X)
end
for fun in (:zero, :one, :oneunit)
    @eval Base.$fun(::Type{<:Static{X}}) where {X} = Base.$fun(typeof(X))
end

# It's a pity there's no AbstractBool supertype.
const StaticBool = Union{StaticInteger{false}, StaticInteger{true}}
const StaticOrBool = Union{StaticBool, Bool}
Base.:!(x::StaticBool) = !Bool(x)

# For complex-valued inputs, there's no auto-convert to floating-point.
# We only support a limited subset of functions, which the user can extend
# as needed.
# TODO: Should have a macro for making functions accept Static input.
for fun in (:abs, :cos, :sin, :exp, :log, :isinf, :isfinite, :isnan)
    @eval Base.$fun(::StaticNumber{X}) where {X} = Base.$fun(X)
end

# Other functions that do not already work
Base.:(<<)(::StaticInteger{X}, y::UInt64) where {X} = X << y
Base.:(>>)(::StaticInteger{X}, y::UInt64) where {X} = X >> y

# Two-argument functions that have methods in promotion.jl that give no_op_err:
for f in (:+, :-, :*, :/, :^)
    @eval Base.$f(::Static{X}, ::Static{X}) where {X} = $f(X,X)
end
# ...where simplifications are possible:
Base.:&(::StaticInteger{X}, ::StaticInteger{X}) where {X} = X
Base.:|(::StaticInteger{X}, ::StaticInteger{X}) where {X} = X
Base.:xor(::StaticInteger{X}, ::StaticInteger{X}) where {X} = zero(X)
Base.:<(::Static{X}, ::Static{X}) where {X} = false
Base.:<=(::Static{X}, ::Static{X}) where {X} = true
Base.:rem(::Static{X}, ::Static{X}) where {X} = zero(X)
Base.:mod(::Static{X}, ::Static{X}) where {X} = zero(X)

# Three-argument function that gives no_op_err
fma(x::Static{X}, y::Static{X}, z::Static{X}) where {X} = fma(X,X,X)

# Static powers using Base.literal_pow.
# This avoids DomainError in some cases?
for T in (Int32, Int64, Float32, Float64, ComplexF32, ComplexF64)
    Base.:^(x::T, ::StaticInteger{p}) where {p} = Base.literal_pow(^, x, Val(p))
end
Base.:^(x::Static{X}, ::StaticInteger{p}) where {X,p} = Base.literal_pow(^, X, Val(p))
Base.:^(x::Static{X}, ::StaticInteger{X}) where {X} = Base.literal_pow(^, X, Val(X)) #disambig

# ntuple accepts Val, so it should also accept static
@inline Base.ntuple(f::F, ::StaticInteger{N}) where {F,N} = Base.ntuple(f, Val(N))

# For brevity, all `Static` numbers are displayed as `static(X)`, rather than, for
# example, `StaticInteger{X}()`. It is possible to discern between the different
# types of `Static` by looking at `X`.
# To get the default behaviour back, run:
#   methods(Base.show, (IO, Static{X} where X)) |> first |> Base.delete_method
function Base.show(io::IO, x::Static{X}) where X
    print(io, "static(")
    show(io, X)
    print(io, ")")
end

include("macros.jl")

include("LengthRanges.jl")

include("trystatic.jl")

end # module


text:
using PtFEM, Test

data = Dict(
  # Frame(nels, nn, ndim, nst, nip, finite_element(nod, nodof))
  :struc_el => Frame(3, 4, 3, 1, 1, Line(2, 3)),
  :properties => [
    4.0e6 1.0e6 0.3e6 0.3e6],
  :x_coords => [0.0, 5.0, 5.0, 5.0],
  :y_coords => [5.0, 5.0, 5.0, 0.0],
  :z_coords => [5.0, 5.0, 0.0, 0.0],
  :gamma => [0.0, 0.0, 90.0],
  :g_num => [
    1 3 4;
    2 2 3],
  :support => [
    (1, [0 0 0 0 0 0]),
    (4, [0 0 0 0 0 0])
    ],
  :loaded_nodes => [
    (2, [0.0 -100.0 0.0 0.0 0.0 0.0])]
)

@time m, dis_df, fm_df = p44(data)

@test round.(m.actions[:, 3], digits=3) == [
  -2.431, 36.287, 27.539, 24.027, 9.49, 80.616,
  2.431, -36.287, -27.539, 113.666, -9.49, -68.461
]


text:
using PiGPIO

red_pin = 18
p=Pi()
set_mode(p, red_pin, OUTPUT)
try
	for i in 1:10
		write(p, red_pin, HIGH)
		sleep(0.5)
		write(p, red_pin, LOW)
		sleep(0.5)
	end
finally
	println("Cleaning up!")
	set_mode(p, red_pin, INPUT)
end


text:
@testset "construction" begin
    @test DateTime(nd) == daytime
    @test Date(nd) == dayt
    @test Time(nd) == tyme
end

@testset "conversion" begin
    @test convert(DateTime, nd) == daytime
    @test convert(Date, nd) == dayt
    @test convert(Time, nd) == tyme
end


text:
export
    ConstraintSet,
	add_constraint!,
	findmax_violation

############################################################################################
#  								CONSTRAINT SETS 										   #
############################################################################################

""" $(TYPEDEF) Set of all constraints for a trajectory optimization problem
Holds a vector of [`ConstraintVals`](@ref) that specify where in the trajectory each constraint
is applied. The `ConstraintSet` efficiently dispatches functions to all of the constraints.

# Constructors:
```julia
ConstraintSet(n,m,N)
ConstraintSet(n,m,Vector{<:ConstraintVals},N)
```

"""
struct ConstraintSet{T}
	n::Int
	m::Int
	constraints::Vector{ConstraintVals}
	p::Vector{Int}
	c_max::Vector{T}
end

function ConstraintSet(n::Int,m::Int,constraints, N)
	@assert !isempty(constraints)
	p = zeros(Int,N)
	c_max = zeros(length(constraints))
	for con in constraints
		for k = 1:N
			p[k] += length(con, k)
		end
	end
	cons = ConstraintVals[]
	append!(cons, constraints)
	ConstraintSet(n,m, cons, p, c_max)
end

function ConstraintSet(n,m,N)
	constraints = Vector{ConstraintVals}()
	p = zeros(Int,N)
	c_max = zeros(0)
	ConstraintSet(n,m,constraints,p,c_max)
end

"Get size of state and control dimensions"
Base.size(conSet::ConstraintSet) = conSet.n, conSet.m, length(conSet.constraints)
"Get number of separate constraints (i.e. ConstraintVals) in the set"
Base.length(conSet::ConstraintSet) = length(conSet.constraints)
Base.copy(conSet::ConstraintSet) = ConstraintSet(conSet.n, conSet.m,
	copy(conSet.constraints), copy(conSet.p), copy(conSet.c_max))

"Re-calculate the number of constraints in the constraint set"
function num_constraints!(conSet::ConstraintSet)
	p = conSet.p
	p .*= 0
	for con in conSet.constraints
		for k = 1:length(p)
			p[k] += length(con, k)
		end
	end
end

"""```julia
num_constraints(::ConstraintSet)
num_constraints(::AbstractSolver)
num_constraints(::Problem)
```
Get the total number of constraints at each time step
"""
@inline num_constraints(conSet::ConstraintSet) = conSet.p

"""```julia
add_constraint!(conSet, conVal::ConstraintVals, idx=-1)
add_constraint!(conSet, con::AbstractConstraint, inds::UnitRange, idx=-1)
```
Add a constraint to the constraint set. You can directly add a `ConstraintVals` type,
but when adding a normal `AbstractConstraint` you must specify the range of knotpoints
to which the constraint applies.

The optional `idx` argument allows you to specify exactly where in the vector of constraints
you'd like the constraint to be added. The ordering will effect the order in which the constraint
appears in concatenated vectors of constraint values (that show up in direct methods), potentially
effecting the band structure of the resultant Jacobian. See solvers for more details.
"""
function add_constraint!(conSet::ConstraintSet, conVal::ConstraintVals, idx=-1)
	if idx == -1
		push!(conSet.constraints, conVal)
		push!(conSet.c_max, 0)
	else
		insert!(conSet.constraints, idx, conVal)
		insert!(conSet.c_max, idx, 0)
	end
	num_constraints!(conSet)
end

function add_constraint!(conSet::ConstraintSet, con::AbstractConstraint,
		inds::UnitRange, idx=-1)
	conVal = ConstraintVals(con, inds)
	add_constraint!(conSet, conVal, idx)
end

"""```julia
max_violation(conSet::ConstraintSet)
max_violation(conSet::ConstraintSet, Z::Traj)
max_violation(prob::Problem, Z=prob.Z)
max_violation(solver::AbstractSolver)
max_violation(solver::AbstractSolver, Z)
```
Calculate the maximum constraint violation for the entire constraint set.
	If the a trajectory is not passed in, the violation is computed from the currently
	stored constraint values; otherwise, the constraints are re-computed using the
	trajectory passed in.
"""
function max_violation(conSet::ConstraintSet)
	max_violation!(conSet)
	maximum(conSet.c_max)
end

function max_violation(conSet::ConstraintSet, Z::Traj)
	evaluate!(conSet, Z)
	max_violation(conSet)
end

function max_violation!(conSet::ConstraintSet{T}) where T
	for i in eachindex(conSet.constraints)
		con = conSet.constraints[i]
		max_violation!(con)
		conSet.c_max[i] = maximum(con.c_max::Vector{T})
	end
end

function findmax_violation(conSet::ConstraintSet)
	max_violation!(conSet)
	c_max0, j_con = findmax(conSet.c_max) # which constraint
	if c_max0 < 1e-16
		return "No constraints violated"
	end
	con = conSet.constraints[j_con]
	i_con = findmax(con.c_max)[2]  # whicn index
	k_con = con.inds[i_con] # time step
	con_sense = sense(con)
	viol = violation(con_sense, con.vals[i_con])
	c_max, i_max = findmax(viol)  # index into constraint
	@assert c_max == c_max0
	con_name = string(typeof(con.con).name)
	return con_name * " at time step $k_con at " * con_label(con.con, i_max)
end

"Calculate the maximum penalty parameter across all constraints"
function max_penalty(conSet::ConstraintSet)
	max_penalty!(conSet)
	maximum(conSet.c_max)
end

function max_penalty!(conSet::ConstraintSet{T}) where T
	for (i,con) in enumerate(conSet.constraints)
		max_penalty!(con)
		conSet.c_max[i] = maximum(con.c_max::Vector{T})
	end
end

function shift_fill!(conSet::ConstraintSet)
	for con in conSet.constraints
		shift_fill!(con)
	end
end

"""```julia
evaluate!(conSet::ConstraintSet, Z::Traj)
```
Compute constraint values for all constraints for the entire trajectory
"""
function evaluate!(conSet::ConstraintSet, Z::Traj)
	for con in conSet.constraints
		evaluate!(con, Z)
	end
end

"""```julia
jacobian!(conSet::ConstraintSet, Z::Traj)
```
Compute constraint Jacobians for all constraints for the entire trajectory
"""
function jacobian!(conSet::ConstraintSet, Z::Traj)
	for con in conSet.constraints
		jacobian!(con, Z)
	end
end

"""```julia
update_active_set!(conSet::ConstraintSet, Z::Traj, ::Val{tol})
```
Compute the active set for the current constraint values, with tolerance tol.
	Uses a value type to avoid an allocation down the line.
"""
function update_active_set!(conSet::ConstraintSet, Z::Traj, val::Val{tol}=Val(0.0)) where tol
	for con in conSet.constraints
		update_active_set!(con, val)
	end
end

Base.iterate(conSet::ConstraintSet) = @inbounds (conSet.constraints[1], 1)
Base.iterate(conSet::ConstraintSet, i) = @inbounds i >= length(conSet.constraints) ? nothing : (conSet.constraints[i+1], i+1)
@inline Base.getindex(conSet::ConstraintSet, i) = conSet.constraints[i]

"Reset all the Lagrange multipliers and constraint values to zero and
	penalties the their initial value"
function reset!(conSet::ConstraintSet)
	for con in conSet.constraints
		reset!(con)
	end
end

function reset_duals!(conSet::ConstraintSet)
	for con in conSet.constraints
		reset_duals!(con)
	end
end

function reset_penalties!(conSet::ConstraintSet)
	for con in conSet.constraints
		reset_penalties!(con)
	end
end

""" ```julia
change_dimension(conSet::ConstraintSet, n, m)
```
Change the dimensionality of the constraint set to one that is strictly larger
than the current dimensions. Useful for state and control augmentation. If the dimension
change affects a particular constraint, it will be wrapped in an [`IndexedConstraint`](@ref)
that simply passes the original sizes to the original constraints. Right now, this assumes
that the original state and controls are first (i.e. the new states or controls are appended
to the end of the state and control vectors).
"""
function change_dimension(conSet::ConstraintSet, n, m)
	cons = map(conSet.constraints) do con

		# Check if the new dimensions match the old ones
		if check_dims(con,n,m)
			return con
		end

		con_idx = IndexedConstraint(n,m,con.con)
		con_val = ConstraintVals(con_idx, con.inds)
	end
	ConstraintSet(n,m,cons, length(conSet.p))
end


text:
using Documenter
using Millboard

makedocs(
    build = joinpath(@__DIR__, "local" in ARGS ? "build_local" : "build"),
    modules = [Millboard],
    clean = false,
    format = Documenter.HTML(
        prettyurls = !("local" in ARGS),
        assets = ["assets/custom.css"],
    ),
    sitename = "Millboard.jl",
    authors = "WooKyoung Noh",
    pages = Any[
        "Home" => "index.md",
    ],
)


text:
module RocketListTest

using Test
using Rocket

import Rocket: List, pushnode!, remove

@testset "List" begin

    println("Testing: List")

    @testset begin
        list = List(Any)

        push!(list, 1)
        push!(list, 2.0)
        push!(list, "Hello")

        @test collect(list) == [ 1, 2.0, "Hello" ]
    end

    @testset begin
        list = List(Int)

        @test_throws MethodError push!(list, 2.0)
        @test_throws MethodError push!(list, "Hello")
    end

    @testset begin
        list = List(Int)

        push!(list, 1)
        push!(list, 2)
        push!(list, 3)

        @test collect(list) == 1:3

        i = 1
        for item in list
            @test item === i
            i += 1
        end

        empty!(list)

        @test collect(list) == []

        empty!(list)

        @test collect(list) == []
    end

    @testset begin
        list = List(Int)

        n1 = pushnode!(list, 1)
        n2 = pushnode!(list, 2)
        n3 = pushnode!(list, 3)

        nodes = [ n1, n2, n3 ]

        i = 1
        for item in list
            @test item === i
            @test collect(list) == i:3
            remove(nodes[i])
            @test collect(list) == (i+1):3
            i += 1
        end

        @test isempty(list)
    end

    @testset begin
        list = List(Int)

        n1 = pushnode!(list, 1)
        n2 = pushnode!(list, 2)
        n3 = pushnode!(list, 3)

        @test collect(list) == [ 1, 2, 3 ]

        remove(n1)

        @test collect(list) == [ 2, 3 ]

        remove(n2)

        @test collect(list) == [ 3 ]

        remove(n3)

        @test collect(list) == [ ]
    end

    @testset begin
        list = List(Int)

        n1 = pushnode!(list, 1)
        n2 = pushnode!(list, 2)
        n3 = pushnode!(list, 3)

        @test collect(list) == [ 1, 2, 3 ]

        remove(n3)

        @test collect(list) == [ 1, 2 ]

        remove(n2)

        @test collect(list) == [ 1 ]

        remove(n1)

        @test collect(list) == [ ]
    end

    @testset begin
        list = List(Int)

        n1 = pushnode!(list, 1)
        n2 = pushnode!(list, 2)
        n3 = pushnode!(list, 3)

        @test collect(list) == [ 1, 2, 3 ]

        remove(n2)

        @test collect(list) == [ 1, 3 ]

        remove(n1)

        @test collect(list) == [ 3 ]

        remove(n3)

        @test collect(list) == [ ]
    end

    @testset begin
        list = List(Int)

        n1 = pushnode!(list, 1)
        n2 = pushnode!(list, 2)
        n3 = pushnode!(list, 3)

        @test collect(list) == [ 1, 2, 3 ]

        remove(n2)

        @test collect(list) == [ 1, 3 ]

        remove(n3)

        @test collect(list) == [ 1 ]

        remove(n1)

        @test collect(list) == [ ]
    end

end

end


text:
import Pkg

using Plots
using VisualRegressionTests
using SpikingNN
using Test
using Distributions
using HypothesisTests
using Random

# environment settings
isci = "CI" ∈ keys(ENV)
datadir = joinpath(@__DIR__, "data")
if !isci
    Pkg.add("Gtk")
    using Gtk
end

# test files to include
testfiles = [
    "excite-test.jl",
    "lif-test.jl",
    "population-test.jl",
    "srm0-test.jl",
    "srm0-threshold-test.jl",
    "stdp-test.jl"
]

@testset "Visual Regression Tests" begin
    ENV["GKSwstype"] = "100"
    for testfile in testfiles
        include(testfile)
    end
end

@testset "Input Models" begin
    include("input-test.jl")
end

@testset "Threshold" begin
    include("threshold-test.jl")
end


text:
export SerreGreenNaghdi

"""
    SerreGreenNaghdi(param;kwargs)

Define an object of type `AbstractModel` in view of solving the initial-value problem for
the Serre-Green-Naghdi model

# Argument
`param` is of type `NamedTuple` and must contain
- dimensionless parameters `ϵ` (nonlinearity) and `μ` (dispersion);
- numerical parameters to construct the mesh of collocation points as `mesh = Mesh(param)`.

## Optional keyword arguments
- `iterative`: solve the elliptic problem through GMRES if `true`, LU decomposition if `false` (default is `true`);
- `precond`: use a (left) preconditioner for GMRES if `true` (default), choose `precond` as the preconditioner if provided;
- `gtol`: relative tolerance of the GMRES algorithm (default is `1e-14`);
- `restart`: the corresponding option of the GMRES algorithm (default is `100`);
- `maxiter`: the corresponding option of GMRES (default is `nothing`);
- `ktol`: tolerance of the Krasny filter (default is `0`, i.e. no filtering);
- `dealias`: dealiasing with Orlicz rule `1-dealias/(dealias+2)` (default is `0`, i.e. no dealiasing);
- `label`: a label for future references (default is `"Green-Naghdi"`);
- `verbose`: prints information if `true` (default is `true`).

# Return values
Generate necessary ingredients for solving an initial-value problem via `solve!`:
1. a function `SerreGreenNaghdi.f!` to be called in explicit time-integration solvers;
2. a function `SerreGreenNaghdi.mapto` which from `(η,v)` of type `InitialData` provides the raw data matrix on which computations are to be executed;
3. a function `SerreGreenNaghdi.mapfro` which from such data matrix returns the Tuple of real vectors `(η,v)`, where
    - `η` is the surface deformation;
    - `v` is the derivative of the trace of the velocity potential;
4. additionally, a handy function `SerreGreenNaghdi.mapfrofull` which from data matrix returns the Tuple of real vectors `(η,v,u)`, where
    - `u` corresponds to the layer-averaged velocity.

"""
mutable struct SerreGreenNaghdi <: AbstractModel

	label   :: String
	f!		:: Function
	mapto	:: Function
	mapfro	:: Function
	mapfrofull	:: Function

    function SerreGreenNaghdi(param::NamedTuple;
							dealias = 0,
							ktol	= 0,
							iterate	= true,
							gtol	= 1e-14,
							precond	= true,
							restart	= nothing,
							maxiter	= nothing,
							label	= "Green-Naghdi",
							verbose	= true)

		m=WhithamGreenNaghdi(param;
							SGN=true,
							dealias=dealias,
							ktol=ktol,
							iterate=iterate,
							gtol=gtol,
							precond=precond,
							restart=restart,
							maxiter=maxiter,
							label=label,
							verbose=verbose)

		new(m.label, m.f!, m.mapto, m.mapfro, m.mapfrofull)
    end
end


text:
"""
Phenotype
"""

mutable struct Network
    nodes::Dict{Int, Node}
    input_id::Array
    hidden_id::Array
    output_id::Array

    function Network(config)
        println("build")

        nodes = Dict{Int, Node}()

        num_inputs = config["num_inputs"]
        num_outputs = config["num_outputs"]

        net = new(nodes, [], [], [])

        for i = 1:num_inputs
            add_node!(net, Node(i, INPUT::NodeType))
            for j = num_inputs+1:num_outputs+num_inputs
                if i == 1
                    add_node!(net, Node(j, OUTPUT::NodeType))
                end
                add_connection!(net, Connection(i, j))
            end
        end
        return net
    end

end

function add_node!(net::Network, node::Node)
    @assert !haskey(net.nodes, node.id) "ID already used"
    net.nodes[node.id] = node

    if node.dtype == INPUT
        push!(net.input_id, node.id)
    elseif node.dtype == OUTPUT
        push!(net.output_id, node.id)
    else
        push!(net.hidden_id, node.id)
    end
end

function add_connection!(net::Network, conn::Connection)
    push!(net.nodes[conn.input].connections["in"], conn)
    push!(net.nodes[conn.output].connections["out"], conn)
end


text:
function sinh(x::DoubleFloat{T}) where {T<:AbstractFloat}
    x < 0 && return -sinh(-x)
    iszero(x) && return zero(x)
    !isfinite(x) && return nan(typeof(x))

    y = exp(x) - exp(-x)
    z = Double64(y.hi*0.5, y.lo*0.5)
    return z
end

function cosh(x::DoubleFloat{T}) where {T<:AbstractFloat}
    x < 0 && return cosh(-x)
    iszero(x) && return zero(x)
    !isfinite(x) && return nan(typeof(x))

    y = exp(x) + exp(-x)
    z = DoubleFloat{T}(y.hi*0.5, y.lo*0.5)
    return z
end

function tanh(x::DoubleFloat{T}) where {T<:AbstractFloat}
    x < 0 && return -tanh(-x)
    iszero(x) && return zero(x)
    !isfinite(x) && return nan(typeof(x))
    epos = exp(x)
    eneg = exp(-x)
    z = (epos - eneg) / (epos + eneg)
    return z
end

csch(x::DoubleFloat{T}) where {T<:AbstractFloat} =
    inv(sinh(x))

sech(x::DoubleFloat{T}) where {T<:AbstractFloat} =
    inv(cosh(x))

coth(x::DoubleFloat{T}) where {T<:AbstractFloat} =
    inv(tanh(x))

#=
function sinh(x::DoubleFloat{T}) where {T<:AbstractFloat}
    x < 0 && return -sinh(-x)
    iszero(x) && return zero(x)
    !isfinite(x) && return nan(typeof(x))

    if abs(x.hi) > 0.05
        ex = exp(x)
        invex = inv(ex)
        y = ex - invex
        z = DoubleFloat{T}(y.hi/2, y.lo/2)
    else
        s = x
        t = x
        r = sqrt(t)
        m = 1.0
        thresh = x.hi * eps(eps(x.hi))

        while (abs(t) > thresh)
            m += 2.0
            t *= r
            t = t / ((m-1) * m)
            s += t
        end
        s
    end
end

function cosh(x::DoubleFloat{T}) where {T<:AbstractFloat}
    x < 0 && return cosh(-x)
    iszero(x) && return zero(x)
    !isfinite(x) && return nan(typeof(x))

    ex = exp(x)
    invex = inv(ex)
    y = ex - invex
    z = DoubleFloat{T}(y.hi/2, y.lo/2)

    return z
end

function tanh(x::DoubleFloat{T}) where {T<:AbstractFloat}
    x < 0 && return -tanh(-x)
    iszero(x) && return zero(x)
    !isfinite(x) && return nan(typeof(x))

    if (abs(x.hi) > 0.05)
        ex = exp(x)
        invex = inv(ex)
        y1 = ex - invex
        y2 = ea + invea
        z = y1 / y2
    else
        s = sinh(x)
        c = sqrt(1.0 + square(s))
        z = s / c
    end

    return z
end

csch(x::DoubleFloat{T}) where {T<:AbstractFloat} =
    inv(sinh(x))

sech(x::DoubleFloat{T}) where {T<:AbstractFloat} =
    inv(cosh(x))

coth(x::DoubleFloat{T}) where {T<:AbstractFloat} =
    inv(tanh(x))
=#


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

mutable struct IoK8sApiCoreV1VsphereVirtualDiskVolumeSource <: SwaggerModel
    fsType::Any # spec type: Union{ Nothing, String } # spec name: fsType
    storagePolicyID::Any # spec type: Union{ Nothing, String } # spec name: storagePolicyID
    storagePolicyName::Any # spec type: Union{ Nothing, String } # spec name: storagePolicyName
    volumePath::Any # spec type: Union{ Nothing, String } # spec name: volumePath

    function IoK8sApiCoreV1VsphereVirtualDiskVolumeSource(;fsType=nothing, storagePolicyID=nothing, storagePolicyName=nothing, volumePath=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1VsphereVirtualDiskVolumeSource, Symbol("fsType"), fsType)
        setfield!(o, Symbol("fsType"), fsType)
        validate_property(IoK8sApiCoreV1VsphereVirtualDiskVolumeSource, Symbol("storagePolicyID"), storagePolicyID)
        setfield!(o, Symbol("storagePolicyID"), storagePolicyID)
        validate_property(IoK8sApiCoreV1VsphereVirtualDiskVolumeSource, Symbol("storagePolicyName"), storagePolicyName)
        setfield!(o, Symbol("storagePolicyName"), storagePolicyName)
        validate_property(IoK8sApiCoreV1VsphereVirtualDiskVolumeSource, Symbol("volumePath"), volumePath)
        setfield!(o, Symbol("volumePath"), volumePath)
        o
    end
end # type IoK8sApiCoreV1VsphereVirtualDiskVolumeSource

const _property_map_IoK8sApiCoreV1VsphereVirtualDiskVolumeSource = Dict{Symbol,Symbol}(Symbol("fsType")=>Symbol("fsType"), Symbol("storagePolicyID")=>Symbol("storagePolicyID"), Symbol("storagePolicyName")=>Symbol("storagePolicyName"), Symbol("volumePath")=>Symbol("volumePath"))
const _property_types_IoK8sApiCoreV1VsphereVirtualDiskVolumeSource = Dict{Symbol,String}(Symbol("fsType")=>"String", Symbol("storagePolicyID")=>"String", Symbol("storagePolicyName")=>"String", Symbol("volumePath")=>"String")
Base.propertynames(::Type{ IoK8sApiCoreV1VsphereVirtualDiskVolumeSource }) = collect(keys(_property_map_IoK8sApiCoreV1VsphereVirtualDiskVolumeSource))
Swagger.property_type(::Type{ IoK8sApiCoreV1VsphereVirtualDiskVolumeSource }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_IoK8sApiCoreV1VsphereVirtualDiskVolumeSource[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1VsphereVirtualDiskVolumeSource }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1VsphereVirtualDiskVolumeSource[property_name]

function check_required(o::IoK8sApiCoreV1VsphereVirtualDiskVolumeSource)
    (getproperty(o, Symbol("volumePath")) === nothing) && (return false)
    true
end

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


text:
module TestGeneticVariation

using Test

import BioCore.Testing:
    get_bio_fmt_specimens,
    random_seq,
    random_interval
import BioCore.Exceptions.MissingFieldException
using BioSequences, GeneticVariation
import BufferedStreams: BufferedInputStream
import IntervalTrees: IntervalValue
import YAML

function random_seq(::Type{A}, n::Integer) where A <: Alphabet
    nts = alphabet(A)
    probs = Vector{Float64}(undef, length(nts))
    fill!(probs, 1 / length(nts))
    return BioSequence{A}(random_seq(n, nts, probs))
end

fmtdir = get_bio_fmt_specimens()

include("vcf.jl")
include("bcf.jl")
include("site_counting.jl")
include("minhash.jl")
include("allele_freq.jl")
include("diversity_measures.jl")
include("seg_sites.jl")

end # Module TestGeneticVariation


text:
using IntervalArithmetic, IntervalOptimisation

G(X) = 1 + sum(abs2, X) / 4000 - prod( cos(X[i] / √i) for i in 1:length(X) )

@time minimise(G, -600..600)
@time minimise(G, -600..600)

for N in (10, 30, 50)
    @show N
    @time minimise(G, IntervalBox(-600..600, N))
    @time minimise(G, IntervalBox(-600..600, N))
end


text:
export symmetric_cauchy_kernel, symmetric_blob_kernel

# w(ξ) = ∑_J ΓJ k_{symmetric, cauchy}(ξ - z)
symmetric_cauchy_kernel(ξ, z) = (ξ - z) != zero(z) ? imag(z)/(π*(ξ - z)*(ξ - conj(z))) : zero(z)

# w(ξ) = ∑_J ΓJ k_{symmetric, blob}(ξ - z, δ)
symmetric_blob_kernel(ξ, z, δ::Float64) = -im/(2*π)*(conj(ξ - z)/(abs2(ξ - z) + δ^2) - (conj(ξ) - z)/(abs2(ξ - conj(z)) + δ^2))


text:
using ScatteringTransform
sourceFolder = "/home/dsweber/Dropbox/matlab/sonarProject/data/sf"
destFolder = "/home/dsweber/Dropbox/matlab/sonarProject/scatteredData/sf"
loadThis = loadSyntheticMatFile
nonlinear = abs; subsam = bspline; stType = "full"
(fullMatrix, rando) = loadSyntheticMatFile("$(sourceFolder)/1000/data00.mat")
layers = stParallel(3, fullMatrix[1,:], [8, 4, 2, 1])
transformFolder(sourceFolder, destFolder, layers, false)
# using JLD
# results = load("/home/dsweber/Dropbox/matlab/sonarProject/scatteredData/sf/1000/data110.mat.jld")
# awe = results["result"]
# awe


text:
# This file is a part of TypeDBClient.  License is MIT: https://github.com/Humans-of-Julia/TypeDBClient.jl/blob/main/LICENSE

# Create a suitable object according to the type specified in the ProtoBuf object `p`.
# Return value will be a subtype of `AbstractType`.

function instantiate(t::Proto._Type)
    label = hasproperty(t, :scope) ? Label(t.scope, t.label) : Label(t.label)
    if t.encoding == Proto.Type_Encoding.ROLE_TYPE
        return RoleType(label, t.root)
    elseif t.encoding == Proto.Type_Encoding.RELATION_TYPE
        return RelationType(label, t.root)
    elseif t.encoding == Proto.Type_Encoding.ATTRIBUTE_TYPE
        return AttributeType(label, t.root, t.value_type)
    elseif t.encoding == Proto.Type_Encoding.ENTITY_TYPE
        return EntityType(label, t.root)
    elseif t.encoding == Proto.Type_Encoding.THING_TYPE
        return ThingType(label, t.root)
    else
        throw(TypeDBClientException(CONCEPT_BAD_ENCODING, t.encoding))
    end
end

# secures the equality of two types by comparing the two labels and the
# same type.
function Base.:(==)(type1::AbstractType, type2::AbstractType)
    if type1 === type2
        return true
    elseif typeof(type1) == typeof(type2) && type1.label == type2.label
        return true
    else
        return false
    end
end

# Return the ProtoBuf encoding value for the given concept type `t`.
function encoding end

encoding(t::AbstractEntityType) = Proto.Type_Encoding.ENTITY_TYPE
encoding(t::AbstractRelationType) = Proto.Type_Encoding.RELATION_TYPE
encoding(t::AbstractAttributeType) = Proto.Type_Encoding.ATTRIBUTE_TYPE
encoding(t::AbstractRoleType) = Proto.Type_Encoding.ROLE_TYPE
encoding(t::AbstractThingType) = Proto.Type_Encoding.THING_TYPE

# TODO The UNRECOGNIZED enum is missing; need to regenerate from proto?
# encoding(t::AbstractType) = Proto.Type_Encoding.UNRECOGNIZED

# I don't think this is unnecessary since the user can broadcast themselves
# proto(types::AbstractVector{T}) where {T <: AbstractType} = proto.(types)


text:
module Autotuning

using Test
using Dates
using DocStringExtensions
using DataFrames
using Random
using Distributions

import Base: getproperty, getindex, keys, values, pairs, push!

export Configuration, Parameter, PowerOfTwo, Measurement,
    perturb, configuration_table, getproperty, getindex,
    keys, values, pairs, push!

include("configuration.jl")
include("database.jl")

end