text:
################################################################################
#
#  Lattices
#
################################################################################

mutable struct ModAlgAssLat{S, T, U}
  base_ring::S
  V::T
  basis::U
  basis_inv::U

  function ModAlgAssLat{S, T, U}(base_ring::S, V::T, basis::U) where {S, T, U}
    z = new{S, T, U}()
    z.base_ring = base_ring
    z.V = V
    z.basis = basis
    return z
  end
end

#ModAlgAssLat(base_ring::S, V::T, basis::U) where {S, T, U} = ModAlgAssLat{S, T, U}(base_ring, V, basis)

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

function Base.show(io::IO, L::ModAlgAssLat)
  print(io, "Lattice of rank ", rank(L), " over ", base_ring(L.base_ring))
end

function _defines_lattice(V::ModAlgAss, O, B)
  Binv = inv(B)
  for g in basis(O)
    T = action(V, elem_in_algebra(g))
    BB = B * T * Binv
    if !isone(abs(denominator(BB)))
      return false
    end
  end
  return true
end

################################################################################
#
#  Lattice arithmetic
#
################################################################################

function _hnf_nonzero(a::fmpq_mat)
  b = fmpq_mat(hnf(FakeFmpqMat(a)))
  i = 1
  while is_zero_row(b, i)
    i += 1
  end
  return b[i:nrows(b), 1:ncols(b)]
end

function *(m::Int, L::ModAlgAssLat)
  return lattice(L.V, L.base_ring, _hnf_nonzero(m * basis_matrix(L)))
end

function +(L::T, M::T) where {T <: ModAlgAssLat}
  @req L.base_ring === M.base_ring "Lattices must be defined over the same order"
  @req L.V === M.V "Lattices must have same ambient module"
  return lattice(L.V, L.base_ring, _hnf_nonzero(vcat(basis_matrix(L),
                                                     basis_matrix(M))))
end

function intersect(L::T, M::T) where {T <: ModAlgAssLat}
  @req L.base_ring === M.base_ring "Lattices must be defined over the same order"
  @req L.V === M.V "Lattices must have same ambient module"
  BM = basis_matrix(M)
  BN = basis_matrix(N)
  dM = denominator(BM)
  dN = denominator(BN)
  d = lcm(dM, dN)
  BMint = change_base_ring(FlintZZ, d * BM)
  BNint = change_base_ring(FlintZZ, d * BN)
  H = vcat(BMint, BNint)
  k, K = left_kernel(H)
  BI = divexact(change_base_ring(FlintQQ, hnf(view(K, 1:k, 1:nrows(BM)) * BMint)), d)
  return lattice(L.V, L.base_ring, BI)
end

################################################################################
#
#  Z-lattices
#
################################################################################

Markdown.doc"""
    lattice(V::ModAlgAss, O::AlgAssAbsOrd, B::MatElem)

Given a module with matrix action over a $\mathbf{Q}$-algebra $A$, a
$\mathbf{Z}$-order of $A$, return the lattice with $O$-lattice basis matrix $B$.
"""
function lattice(V::ModAlgAss{FlintRationalField}, O::AlgAssAbsOrd, B::MatElem; check::Bool = true)
  if B isa fmpq_mat
    return _lattice(V, O, B, check = check)
  else
    return _lattice(V, O, change_base_ring(QQ, B)::fmpq_mat, check = check)
  end
end

# internal function to construct lattice
function _lattice(V::ModAlgAss{FlintRationalField}, O::AlgAssAbsOrd, B::fmpq_mat; check::Bool = true, is_hnf::Bool = false)
  if check
    fl = _defines_lattice(V, O, B)
    @req fl "Z-lattice with this basis matrix is not invariant under order"
  end

  @hassert :ModLattice _defines_lattice(V, O, B)

  if !is_hnf
    BB = fmpq_mat(hnf!(FakeFmpqMat(B), :upperright))
  else
    BB = B
  end
  return ModAlgAssLat{typeof(O), typeof(V), typeof(BB)}(O, V, BB)
end

rank(L::ModAlgAssLat) = nrows(basis_matrix(L))

################################################################################
#
#  Basis matrix
#
################################################################################

function basis_matrix(L::ModAlgAssLat)
  @req base_ring(L.base_ring) isa FlintIntegerRing "Order of lattice must be a Z-order"
  @req has_matrix_action(L.V) "Action on module must be given by matrices"
  return L.basis
end

function basis_matrix_inverse(L::ModAlgAssLat)
  @req base_ring(L.base_ring) isa FlintIntegerRing "Order of lattice must be a Z-order"
  @req has_matrix_action(L.V) "Action on module must be given by matrices"
  if isdefined(L, :basis_inv)
    return L.basis_inv
  else
    M = inv(basis_matrix(L))
    L.basis_inv = M
    return M
  end
end

################################################################################
#
#  Reduction
#
################################################################################

@doc Markdown.doc"""
    reduction(L::ModAlgAssLat, p::IntegerUnion) -> ModAlgAss

Given an $L$ over an $\mathbf{Z}$-order and a prime $p$, return the module $L/pL$
over the field $\mathbf{F}_p$.

Note that the module will be defined without algebra and the action will be
given by $\rank(O)$ many generators. To obtain a module for the
$\mathbf{F}_p$-algebra $O/pO$, the algebra first has to be constructed using
`reduction(O, p)`.

See also `change_coefficient_ring`.
"""
function reduction(L::ModAlgAssLat, p::IntegerUnion)
  @req base_ring((L.base_ring)) isa FlintIntegerRing "Order must be a Z-order"
  F = GF(p, cached = false)
  a = action_of_basis(L)
  amodp = map(m -> change_base_ring(F, m), a)
  return Module(amodp)
end

@doc Markdown.doc"""
    change_coefficient_ring(R::Ring, L::ModAlgAssLat{FlintIntegerRing}) -> ModAlgAss

Given a lattice $L$ over an $\mathbf{Z}$-order $L$, return the $L \otimes R$
over the ring $R$.

Note that the module will be defined without algebra and the action will be
given by $\rank(O)$ many generators.
"""
function change_coefficient_ring(R::Ring, L::ModAlgAssLat)
  @req base_ring((L.base_ring)) isa FlintIntegerRing "Order must be a Z-order"
  a = action_of_basis(L)
  aR = map(m -> change_base_ring(R, m), a)
  return Module(aR)
end

@doc Markdown.doc"""
    action(L::ModAlgAssLat, x)

Given a lattice $L$ over an order $O$ and an element $x$ of $O$, return
the matrix with which $x$ is acting on $L$.
"""
function action(L::ModAlgAssLat, x)
  T = basis_matrix(L)
  Tinv = basis_matrix_inverse(L)
  M = T * action(L.V, elem_in_algebra(x))
  mul!(M, M, Tinv)
  return M
end

function action_of_basis(L::ModAlgAssLat)
  A = action_of_order_basis(L.V, L.base_ring)
  T = basis_matrix(L)
  Tinv = basis_matrix_inverse(L)
  res = Vector{eltype(A)}(undef, length(A))
  for i in 1:length(res)
    M = T * A[i]
    res[i] = mul!(M, M, Tinv)
  end
  return res
end

#function algebra(M::Vector{T}) where {T <: MatElem}
#  @assert length(M) > 0
#  A = M[1]
#  n = nrows(A)
#  n2 = n^2
#  @assert n == ncols(A)
#  K = base_ring(A)
#  Mprod = M
#  Morig = copy(M)
#
#  current_dim = -1
#
#  B = zero_matrix(K, length(Mprod) + 1, n2)
#
#  l = 0
#  while true
#    if l != 0
#      B = zero_matrix(K, length(Mprod), n2)
#    end
#    for k in 1:length(Mprod)
#      for i in 1:n
#        for j in 1:n
#          B[k, (i - 1)* n  + j] = Mprod[k][i, j]
#        end
#      end
#    end
#    # Add the identity
#    if l == 0
#      for i in 1:n
#        B[length(M) + 1, (i - 1)*n + i] = one(K)
#      end
#    end
#    new_dim = rref!(B)
#    if new_dim == current_dim
#      break
#    end
#    current_dim = new_dim
#    M = [ matrix(K, n, n, [B[k, (i - 1)*n + j] for i in 1:n for j in 1:n]) for k in 1:new_dim]
#    Mprod = [ M[i] * M[j] for i in 1:length(M) for j in 1:length(M) ]
#    l = l + 1
#  end
#
#  dim = current_dim
#  B = sub(B, 1:dim, 1:ncols(B))
#
#  basis = [ matrix(K, n, n, [B[k, (i - 1)*n + j] for i in 1:n for j in 1:n]) for k in 1:dim]
#
#  @assert isone(basis[1])
#
#  v = zero_matrix(K, 1, n2)
#
#  structure = Array{elem_type(K), 3}(dim, dim, dim)
#
#  for k in 1:dim
#    for l in 1:dim
#      N = basis[k] * basis[l]
#      for i in 1:n
#        for j in 1:n
#          v[1, (i - 1)* n  + j] = N[i, j]
#        end
#      end
#      b, u = can_solve_with_solution(B, v, side = :left)
#      error("NOT HERE!")
#      @assert b
#      @assert N == sum(u[i]*basis[i] for i in 1:dim)
#      for m in 1:dim
#        structure[k, l, m] = u[m, 1]
#      end
#    end
#  end
#
#  A = AlgAss(K, structure)
#
#  gens = Vector{AlgAssElem{elem_type(K)}}(length(Morig))
#
#  for l in 1:length(Morig)
#    N = Morig[l]
#    for i in 1:n
#      for j in 1:n
#        v[1, (i - 1)* n  + j] = N[i, j]
#      end
#    end
#    b, u = can_solve_with_solution(B, v, side = :left)
#    gens[l] =  A([u[1, m] for m in 1:dim])
#  end
#
#  A.gens = gens
#
#  return A
#end

#function gens(A::AlgAss{T}) where {T}
#  #return A.gens::Vector{AlgAssElem{T}}
#end

##

function natural_lattice(O::AlgAssAbsOrd{<:AlgMat{fmpq, fmpq_mat}})
  A = algebra(O)
  if all(x -> isone(denominator(matrix(elem_in_algebra(x)))), basis(O, copy = false))
    return lattice(Module(A, matrix.(basis(A))), O, identity_matrix(QQ, degree(algebra(O))))
  else
    error("Order is not contained in M_n(Z)")
  end
end

################################################################################
#
#  Local isomorphism
#
################################################################################

function is_locally_isomorphic_with_isomophism(L::ModAlgAssLat, M::ModAlgAssLat, p::fmpz)
  @req L.base_ring === M.base_ring "Orders of lattices must agree"
  @req base_ring(L.base_ring) isa FlintIntegerRing "Order must be a Z-order"
  return _islocally_isomorphic_with_isomophism(L, M, p, Val{true})
end

function is_locally_isomorphic(L::ModAlgAssLat, M::ModAlgAssLat, p::fmpz)
  @req L.base_ring === M.base_ring "Orders of lattices must agree"
  @req base_ring(L.base_ring) isa FlintIntegerRing "Order must be a Z-order"
  return _islocally_isomorphic_with_isomophism(L, M, p, Val{false})
end

function _islocally_isomorphic_with_isomophism(L::ModAlgAssLat, M::ModAlgAssLat, p::fmpz, with_isomorphism::Type{Val{S}} = Val{true}) where {S}
  @req L.base_ring === M.base_ring "Orders of lattices must agree"
  @req base_ring(L.base_ring) isa FlintIntegerRing "Order must be a Z-order"
  # We are assuming that L.V === M.V is absolutely irreducible
  # I will not even check this.
  # TODO: Relax this, Tommy knows how to do this
  @assert L.V === M.V
  T = basis_matrix(L) * basis_matrix_inverse(M)
  d = denominator(T)
  T = d * T
  if with_isomorphism === Val{true}
    fl = iszero(valuation(det(T), p))
    if fl
      error("Tell the developers to finally do it!")
    else
      return false, T
    end
  else
    return iszero(valuation(det(T), p))
  end
end

function _lift_to_Q(K::gfp_mat)
  z = zero_matrix(QQ, nrows(K), ncols(K))
  for i in 1:nrows(K)
    for j in 1:ncols(K)
      u = K[i, j].data
      if iszero(u)
        continue
      else
        z[i, j] = u
      end
    end
  end
  return z
end

function _lift_to_Q!(z, K::gfp_mat)
  for i in 1:nrows(K)
    for j in 1:ncols(K)
      u = K[i, j].data
      if iszero(u)
        continue
      else
        z[i, j] = u
      end
    end
  end
  return z
end

function addmul!(A::gfp_mat, B::gfp_mat, C::gfp_elem, D::gfp_mat)
  ccall((:nmod_mat_scalar_addmul_ui, libflint), Cvoid, (Ref{gfp_mat}, Ref{gfp_mat}, Ref{gfp_mat}, UInt), A, B, D, C.data)
  return A
end

function mul!(A::gfp_mat, B::gfp_elem, D::gfp_mat)
  ccall((:nmod_mat_scalar_mul_ui, libflint), Cvoid, (Ref{gfp_mat}, Ref{gfp_mat}, UInt), A, C, B.data)
end

function pmaximal_sublattices(L::ModAlgAssLat, p::Int; filter = nothing, composition_factors = nothing)
  res = typeof(L)[]
  if composition_factors === nothing
    V = reduction(L, p)
    F = V.base_ring
    comp_fac = Hecke.composition_factors(V)
  else
    comp_fac = composition_factors
    F = coefficient_ring(comp_fac[1])
    V = change_coefficient_ring(F, L)
  end
  for C in comp_fac
    H = basis_of_hom(V, C)
    if length(H) == 0
      continue
    end
    # We have to loop over all morphisms, but we can discard multiples
    # So we do a projective space thingy
    pivs = Int[]
    Kl = zero_matrix(QQ, rank(L), rank(L))
    for v in enumerate_lines(F, length(H))
      zero!(Kl)
      empty!(pivs)
      #k = sum(H[i] * v[i] for i in 1:length(H))
      if length(H) == 1
        k = H[1]
      else
        k = v[1] * H[1]
        for i in 2:length(H)
          if iszero(v[i])
            continue
          end
          addmul!(k, k, v[i], H[i])
          #k = add!(k ,k, v[i] * H[i])
        end
      end
      r, K = left_kernel(k)
      _, K = rref(K)
      # The lattice we are looking is generated by lift(K) + pZ^n
      # The HNF of this lattice is/would be the lift(rref(K) and p's added in
      # the zero-rows such that the rank is full.
      # We scale the rows and check the pivots
      l = 1
      m = r + 1
      for i in 1:r
        while iszero(K[i, l])
          l += 1
        end
        # The rows are not normalized
        push!(pivs, l)
      end
      # We lift directly to Q
      Kl = _lift_to_Q!(Kl, K)
      # Set the zero rows correctly
      for i in 1:nrows(K)
        if !(i in pivs)
          Kl[m, i] = p
          m += 1
        end
      end
      # Kl has the same Z-span as fmpq_mat(hnf_modular_eldiv(lift(K), fmpz(p)))
      # We need the basis matrix with respect to
      _bmat = mul!(Kl, Kl, L.basis)
      LL = lattice(L.V, L.base_ring, _bmat, check = false)
      if any(LLL -> LLL.basis == LL.basis, res)
        continue
      end
      if filter === nothing ||
         (filter == :local_isomorphism && all(LLL -> !is_locally_isomorphic(LLL, LL, fmpz(p)), res))
        push!(res, LL)
      end
    end
  end
  return res
end

################################################################################
#
#  Centering algorithm
#
################################################################################

function sublattice_classes(L::ModAlgAssLat, p::Int)
  res = typeof(L)[L]
  to_check = typeof(L)[L]
  while !isempty(to_check)
    M = pop!(to_check)
    X = pmaximal_sublattices(M, p, filter = :local_isomorphism)
    for N in X
      if any(LLL -> is_locally_isomorphic(LLL, N, fmpz(p))[1], res)
        continue
      else
        push!(res, N)
        push!(to_check, N)
      end
    end
    #@show length(res)
  end
  return res
end

function is_sublattice(L::ModAlgAssLat, M::ModAlgAssLat)
  return isone(denominator(basis_matrix(L) * basis_matrix_inverse(M)))
end

function _issublattice(L::ModAlgAssLat, M::ModAlgAssLat, tmp)
  tmp = mul!(tmp, basis_matrix(L), basis_matrix_inverse(M))
  return isone(denominator(tmp))
end

function sublattices(L::ModAlgAssLat, p::Int, level = inf)
  res = typeof(L)[L]
  to_check = typeof(L)[L]
  pL = p*L
  temp = typeof(L)[]
  i = 0
  F = GF(p)
  comp_fac = composition_factors(change_coefficient_ring(F, L))
  tmp_ = zero_matrix(QQ, rank(L), rank(L))
  while !isempty(to_check)
    if i >= level
      break
    end
    i += 1
    empty!(temp)

    for M in to_check
      pmax = pmaximal_sublattices(M, p, composition_factors = comp_fac)
      for N in pmax
        if _issublattice(N, pL, tmp_)
          continue
        end
        if any(LLL -> LLL.basis == N.basis, temp) || any(LLL -> LLL.basis == N.basis, res)
          continue
        else
          push!(temp, N)
        end
      end
    end
    empty!(to_check)
    append!(to_check, temp)
    append!(res, temp)
  end
  return res
end


text:
# SBP Node data #
#N=1
Q_GLegendre_N1 = 6;
Q_GLegendre_rq_N1 = [-0.577350269189626,0.577350269189626,-1.000000000000000,-1.000000000000000,-0.577350269189626,0.577350269189626];
Q_GLegendre_sq_N1 = [-1.000000000000000,-1.000000000000000,-0.577350269189626,0.577350269189626,0.577350269189626,-0.577350269189626];
Q_GLegendre_wq_N1 = [0.333333333333333,0.333333333333333,0.333333333333333,0.333333333333333,0.333333333333333,0.333333333333333];
GQ_N1 = (Q_GLegendre_rq_N1, Q_GLegendre_sq_N1, Q_GLegendre_wq_N1)
#N=2
Q_GLegendre_N2 = 10;
Q_GLegendre_rq_N2 = [-0.774596669241483,0.000000000000000,0.774596669241483,-1.000000000000000,-1.000000000000000,-1.000000000000000,-0.774596669241483,0.000000000000000,0.774596669241483,-0.333333333333333];
Q_GLegendre_sq_N2 = [-1.000000000000000,-1.000000000000000,-1.000000000000000,-0.774596669241483,0.000000000000000,0.774596669241483,0.774596669241483,0.000000000000000,-0.774596669241483,-0.333333333333333];
Q_GLegendre_wq_N2 = [0.083333333333333,0.200000000000000,0.083333333333333,0.083333333333333,0.200000000000000,0.083333333333333,0.083333333333333,0.200000000000000,0.083333333333333,0.900000000000000];
GQ_N2 = (Q_GLegendre_rq_N2, Q_GLegendre_sq_N2, Q_GLegendre_wq_N2)
#N=3
Q_GLegendre_N3 = 18;
Q_GLegendre_rq_N3 = [-0.861136311594054,-0.339981043584857,0.339981043584856,0.861136311594053,-1.000000000000000,-1.000000000000000,-1.000000000000000,-1.000000000000000,-0.861136311594054,-0.339981043584857,0.339981043584856,0.861136311594053,0.168314227951316,-0.542461986333870,0.168314227951316,-0.625852241617446,-0.542461986333870,-0.625852241617446];
Q_GLegendre_sq_N3 = [-1.000000000000000,-1.000000000000000,-1.000000000000000,-1.000000000000000,-0.861136311594054,-0.339981043584857,0.339981043584856,0.861136311594053,0.861136311594054,0.339981043584857,-0.339981043584856,-0.861136311594053,-0.542461986333870,0.168314227951316,-0.625852241617446,0.168314227951316,-0.625852241617446,-0.542461986333870];
Q_GLegendre_wq_N3 = [0.030198029745131,0.080913081365980,0.080913081365980,0.030198029745131,0.030198029745131,0.080913081365980,0.080913081365980,0.030198029745131,0.030198029745131,0.080913081365980,0.080913081365980,0.030198029745131,0.222222222222222,0.222222222222222,0.222222222222222,0.222222222222222,0.222222222222222,0.222222222222222];
GQ_N3 = (Q_GLegendre_rq_N3, Q_GLegendre_sq_N3, Q_GLegendre_wq_N3)

#N=4
Q_GLegendre_N4 = 22;
Q_GLegendre_rq_N4 = [-0.906179845938664,-0.538469310105683,0.000000000000000,0.538469310105683,0.906179845938664,-1.000000000000000,-1.000000000000000,-1.000000000000000,-1.000000000000000,-1.000000000000000,-0.906179845938664,-0.538469310105683,0.000000000000000,0.538469310105683,0.906179845938664,-0.333333333333333,-0.721132537169100,-0.721132537169100,0.442265074338200,-0.123152095118386,-0.123152095118386,-0.753695809763228];
Q_GLegendre_sq_N4 = [-1.000000000000000,-1.000000000000000,-1.000000000000000,-1.000000000000000,-1.000000000000000,-0.906179845938664,-0.538469310105683,0.000000000000000,0.538469310105683,0.906179845938664,0.906179845938664,0.538469310105683,0.000000000000000,-0.538469310105683,-0.906179845938664,-0.333333333333333,-0.721132537169100,0.442265074338200,-0.721132537169100,-0.123152095118386,-0.753695809763228,-0.123152095118386];
Q_GLegendre_wq_N4 = [0.013202630162003,0.041060919360860,0.037074169667907,0.041060919360860,0.013202630162003,0.013202630162003,0.041060919360860,0.037074169667907,0.041060919360860,0.013202630162003,0.013202630162003,0.041060919360860,0.037074169667907,0.041060919360860,0.013202630162003,0.182199822395302,0.210858659241677,0.210858659241677,0.210858659241677,0.249473464579591,0.249473464579591,0.249473464579591];
GQ_N4 = (Q_GLegendre_rq_N4, Q_GLegendre_sq_N4, Q_GLegendre_wq_N4)

#N=5
Q_GLegendre_N5 = 34;
Q_GLegendre_rq_N5 = [-0.932469514203153,-0.661209386466264,-0.238619186083197,0.238619186083197,0.661209386466263,0.932469514203152,-1.000000000000000,-1.000000000000000,-1.000000000000000,-1.000000000000000,-1.000000000000000,-1.000000000000000,-0.932469514203153,-0.661209386466264,-0.238619186083197,0.238619186083197,0.661209386466263,0.932469514203152,-0.333333333333333,-0.550932095205921,-0.550932095205921,0.101864190411842,0.611382635695146,-0.744884752682278,0.611382635695146,-0.866497883012868,-0.744884752682278,-0.866497883012868,0.154088161163644,-0.813213547833420,0.154088161163644,-0.340874613330224,-0.813213547833420,-0.340874613330224];
Q_GLegendre_sq_N5 = [-1.000000000000000,-1.000000000000000,-1.000000000000000,-1.000000000000000,-1.000000000000000,-1.000000000000000,-0.932469514203153,-0.661209386466264,-0.238619186083197,0.238619186083197,0.661209386466263,0.932469514203152,0.932469514203153,0.661209386466264,0.238619186083197,-0.238619186083197,-0.661209386466263,-0.932469514203152,-0.333333333333333,-0.550932095205921,0.101864190411842,-0.550932095205921,-0.744884752682278,0.611382635695146,-0.866497883012868,0.611382635695146,-0.866497883012868,-0.744884752682278,-0.813213547833420,0.154088161163644,-0.340874613330224,0.154088161163644,-0.340874613330224,-0.813213547833420];
Q_GLegendre_wq_N5 = [0.006658503841347,0.016313739546177,0.025165807590991,0.025165807590991,0.016313739546177,0.006658503841347,0.006658503841347,0.016313739546177,0.025165807590991,0.025165807590991,0.016313739546177,0.006658503841347,0.006658503841347,0.016313739546177,0.025165807590991,0.025165807590991,0.016313739546177,0.006658503841347,0.209054393645172,0.095225751929022,0.095225751929022,0.095225751929022,0.060381677758004,0.060381677758004,0.060381677758004,0.060381677758004,0.060381677758004,0.060381677758004,0.142358329691442,0.142358329691442,0.142358329691442,0.142358329691442,0.142358329691442,0.142358329691442];
GQ_N5 = (Q_GLegendre_rq_N5, Q_GLegendre_sq_N5, Q_GLegendre_wq_N5)

#N=6
Q_GLegendre_N6 = 42;
Q_GLegendre_rq_N6 = [-0.949107912342759,-0.741531185599395,-0.405845151377397,0.000000000000000,0.405845151377397,0.741531185599396,0.949107912342757,-1.000000000000000,-1.000000000000000,-1.000000000000000,-1.000000000000000,-1.000000000000000,-1.000000000000000,-1.000000000000000,-0.949107912342759,-0.741531185599395,-0.405845151377397,0.000000000000000,0.405845151377397,0.741531185599396,0.949107912342757,-0.896694477209711,-0.896694477209711,0.793388954419422,-0.210196994443282,-0.210196994443282,-0.579606011113436,-0.593408612745712,-0.593408612745712,0.186817225491424,-0.273069434548136,-0.867319561224946,-0.273069434548136,0.140388995773082,-0.867319561224946,0.140388995773082,-0.874214255458302,0.514759945734612,-0.874214255458302,-0.640545690276310,0.514759945734612,-0.640545690276310];
Q_GLegendre_sq_N6 = [-1.000000000000000,-1.000000000000000,-1.000000000000000,-1.000000000000000,-1.000000000000000,-1.000000000000000,-1.000000000000000,-0.949107912342759,-0.741531185599395,-0.405845151377397,0.000000000000000,0.405845151377397,0.741531185599396,0.949107912342757,0.949107912342759,0.741531185599395,0.405845151377397,0.000000000000000,-0.405845151377397,-0.741531185599396,-0.949107912342757,-0.896694477209711,0.793388954419422,-0.896694477209711,-0.210196994443282,-0.579606011113436,-0.210196994443282,-0.593408612745712,0.186817225491424,-0.593408612745712,-0.867319561224946,-0.273069434548136,0.140388995773082,-0.273069434548136,0.140388995773082,-0.867319561224946,0.514759945734612,-0.874214255458302,-0.640545690276310,-0.874214255458302,-0.640545690276310,0.514759945734612];
Q_GLegendre_wq_N6 = [0.002716104484117,0.009406550665867,0.013917356931396,0.015633693986814,0.013917356931396,0.009406550665867,0.002716104484117,0.002716104484117,0.009406550665867,0.013917356931396,0.015633693986814,0.013917356931396,0.009406550665867,0.002716104484117,0.002716104484117,0.009406550665867,0.013917356931396,0.015633693986814,0.013917356931396,0.009406550665867,0.002716104484117,0.032730621299045,0.032730621299045,0.032730621299045,0.140317170156143,0.140317170156143,0.140317170156143,0.117404101260644,0.117404101260644,0.117404101260644,0.087979462909756,0.087979462909756,0.087979462909756,0.087979462909756,0.087979462909756,0.087979462909756,0.066271064990874,0.066271064990874,0.066271064990874,0.066271064990874,0.066271064990874,0.066271064990874];
GQ_N6 = (Q_GLegendre_rq_N6, Q_GLegendre_sq_N6, Q_GLegendre_wq_N6)

#N=1
Q_GLobatto_N1 = 6;
Q_GLobatto_rq_N1 = [-1.000000000000000,0.000000000000000,1.000000000000000,0.000000000000000,-1.000000000000000,-1.000000000000000];
Q_GLobatto_sq_N1 = [-1.000000000000000,-1.000000000000000,-1.000000000000000,0.000000000000000,1.000000000000000,0.000000000000000];
Q_GLobatto_wq_N1 = [0.333333333333333,0.333333333333333,0.333333333333333,0.333333333333333,0.333333333333333,0.333333333333333];
GL_N1 = (Q_GLobatto_rq_N1, Q_GLobatto_sq_N1, Q_GLobatto_wq_N1)

#N=2
Q_GLobatto_N2 = 10;
Q_GLobatto_rq_N2 = [-1.000000000000000,-0.447213595499958,0.447213595499958,1.000000000000000,0.447213595499958,-0.447213595499958,-1.000000000000000,-1.000000000000000,-1.000000000000000,-0.333333333333333];
Q_GLobatto_sq_N2 = [-1.000000000000000,-1.000000000000000,-1.000000000000000,-1.000000000000000,-0.447213595499958,0.447213595499958,1.000000000000000,0.447213595499958,-0.447213595499958,-0.333333333333333];
Q_GLobatto_wq_N2 = [0.033333333333333,0.166666666666667,0.166666666666667,0.033333333333333,0.166666666666667,0.166666666666667,0.033333333333333,0.166666666666667,0.166666666666667,0.900000000000000];
GL_N2 = (Q_GLobatto_rq_N2, Q_GLobatto_sq_N2, Q_GLobatto_wq_N2)

#N=3
Q_GLobatto_N3 = 15;
Q_GLobatto_rq_N3 = [-1.000000000000000,-0.654653670707977,0.000000000000000,0.654653670707977,1.000000000000000,0.654653670707977,0.000000000000000,-0.654653670707977,-1.000000000000000,-1.000000000000000,-1.000000000000000,-1.000000000000000,-0.585309648672816,-0.585309648672816,0.170619297345632];
Q_GLobatto_sq_N3 = [-1.000000000000000,-1.000000000000000,-1.000000000000000,-1.000000000000000,-1.000000000000000,-0.654653670707977,0.000000000000000,0.654653670707977,1.000000000000000,0.654653670707977,0.000000000000000,-0.654653670707977,-0.585309648672816,0.170619297345632,-0.585309648672816];
Q_GLobatto_wq_N3 = [0.014698618394802,0.056443964486202,0.097526003618641,0.056443964486202,0.014698618394802,0.056443964486202,0.097526003618641,0.056443964486202,0.014698618394802,0.056443964486202,0.097526003618641,0.056443964486202,0.441554115680819,0.441554115680819,0.441554115680819];
GL_N3 = (Q_GLobatto_rq_N3, Q_GLobatto_sq_N3, Q_GLobatto_wq_N3)

#N=4
Q_GLobatto_N4 = 24;
Q_GLobatto_rq_N4 = [-1.000000000000000,-0.765055323929465,-0.285231516480645,0.285231516480645,0.765055323929466,1.000000000000000,0.765055323929465,0.285231516480645,-0.285231516480645,-0.765055323929466,-1.000000000000000,-1.000000000000000,-1.000000000000000,-1.000000000000000,-1.000000000000000,-0.153862961347394,-0.153862961347394,-0.692274077305212,0.462899656119068,-0.801694541056790,0.462899656119068,-0.661205115062278,-0.801694541056790,-0.661205115062278];
Q_GLobatto_sq_N4 = [-1.000000000000000,-1.000000000000000,-1.000000000000000,-1.000000000000000,-1.000000000000000,-1.000000000000000,-0.765055323929465,-0.285231516480645,0.285231516480645,0.765055323929466,1.000000000000000,0.765055323929465,0.285231516480645,-0.285231516480645,-0.765055323929466,-0.153862961347394,-0.692274077305212,-0.153862961347394,-0.801694541056790,0.462899656119068,-0.661205115062278,0.462899656119068,-0.661205115062278,-0.801694541056790];
Q_GLobatto_wq_N4 = [0.007531613457659,0.021086271203965,0.045504525339300,0.045504525339300,0.021086271203965,0.007531613457659,0.021086271203965,0.045504525339300,0.045504525339300,0.021086271203965,0.007531613457659,0.021086271203965,0.045504525339300,0.045504525339300,0.021086271203965,0.304805308340987,0.304805308340987,0.304805308340987,0.110574075890745,0.110574075890745,0.110574075890745,0.110574075890745,0.110574075890745,0.110574075890745];
GL_N4 = (Q_GLobatto_rq_N4, Q_GLobatto_sq_N4, Q_GLobatto_wq_N4)

#N=5
Q_GLobatto_N5 = 33;
Q_GLobatto_rq_N5 = [-1.000000000000000,-0.830223896278567,-0.468848793470714,0.000000000000000,0.468848793470714,0.830223896278567,1.000000000000000,0.830223896278567,0.468848793470714,0.000000000000000,-0.468848793470714,-0.830223896278567,-1.000000000000000,-1.000000000000000,-1.000000000000000,-1.000000000000000,-1.000000000000000,-1.000000000000000,-0.481290025456663,-0.481290025456663,-0.037419949086674,-0.085394918585138,-0.085394918585138,-0.829210162829724,-0.826342421973033,-0.826342421973033,0.652684843946066,-0.834845895799092,0.320885646218972,-0.834845895799092,-0.486039750419880,0.320885646218972,-0.486039750419880];
Q_GLobatto_sq_N5 = [-1.000000000000000,-1.000000000000000,-1.000000000000000,-1.000000000000000,-1.000000000000000,-1.000000000000000,-1.000000000000000,-0.830223896278567,-0.468848793470714,0.000000000000000,0.468848793470714,0.830223896278567,1.000000000000000,0.830223896278567,0.468848793470714,0.000000000000000,-0.468848793470714,-0.830223896278567,-0.481290025456663,-0.037419949086674,-0.481290025456663,-0.085394918585138,-0.829210162829724,-0.085394918585138,-0.826342421973033,0.652684843946066,-0.826342421973033,0.320885646218972,-0.834845895799092,-0.486039750419880,-0.834845895799092,-0.486039750419880,0.320885646218972];
Q_GLobatto_wq_N5 = [0.002314778534840,0.013700973144398,0.020117879303463,0.023426958130206,0.020117879303463,0.013700973144398,0.002314778534840,0.013700973144398,0.020117879303463,0.023426958130206,0.020117879303463,0.013700973144398,0.002314778534840,0.013700973144398,0.020117879303463,0.023426958130206,0.020117879303463,0.013700973144398,0.176274129326032,0.176274129326032,0.176274129326032,0.113088210164246,0.113088210164246,0.113088210164246,0.076346374759077,0.076346374759077,0.076346374759077,0.103789255428272,0.103789255428272,0.103789255428272,0.103789255428272,0.103789255428272,0.103789255428272];
GL_N5 = (Q_GLobatto_rq_N5, Q_GLobatto_sq_N5, Q_GLobatto_wq_N5)

#N=6
Q_GLobatto_N6 = 40;
Q_GLobatto_rq_N6 = [-1.000000000000000,-0.871740148509606,-0.591700181433142,-0.209299217902479,0.209299217902479,0.591700181433143,0.871740148509605,1.000000000000000,0.871740148509606,0.591700181433142,0.209299217902479,-0.209299217902479,-0.591700181433143,-0.871740148509605,-1.000000000000000,-1.000000000000000,-1.000000000000000,-1.000000000000000,-1.000000000000000,-1.000000000000000,-1.000000000000000,-0.333333333333333,-0.832518933286299,-0.832518933286299,0.665037866572598,-0.018986264273055,-0.018986264273055,-0.962027471453890,-0.130967625384817,-0.130967625384817,-0.738064749230366,-0.591781036586590,-0.591781036586590,0.183562073173180,-0.465390307260470,0.362915057330052,-0.465390307260470,-0.897524750069582,0.362915057330052,-0.897524750069582];
Q_GLobatto_sq_N6 = [-1.000000000000000,-1.000000000000000,-1.000000000000000,-1.000000000000000,-1.000000000000000,-1.000000000000000,-1.000000000000000,-1.000000000000000,-0.871740148509606,-0.591700181433142,-0.209299217902479,0.209299217902479,0.591700181433143,0.871740148509605,1.000000000000000,0.871740148509606,0.591700181433142,0.209299217902479,-0.209299217902479,-0.591700181433143,-0.871740148509605,-0.333333333333333,-0.832518933286299,0.665037866572598,-0.832518933286299,-0.018986264273055,-0.962027471453890,-0.018986264273055,-0.130967625384817,-0.738064749230366,-0.130967625384817,-0.591781036586590,0.183562073173180,-0.591781036586590,0.362915057330052,-0.465390307260470,-0.897524750069582,-0.465390307260470,-0.897524750069582,0.362915057330052];
Q_GLobatto_wq_N6 = [0.000886028699507,0.010480045161785,0.011623548297153,0.001211835242253,0.001211835242253,0.011623548297153,0.010480045161785,0.000886028699507,0.010480045161785,0.011623548297153,0.001211835242253,0.001211835242253,0.011623548297153,0.010480045161785,0.000886028699507,0.010480045161785,0.011623548297153,0.001211835242253,0.001211835242253,0.011623548297153,0.010480045161785,0.171525496248007,0.074048529130209,0.074048529130209,0.074048529130209,0.049562169655280,0.049562169655280,0.049562169655280,0.136384241605594,0.136384241605594,0.136384241605594,0.139937548771021,0.139937548771021,0.139937548771021,0.081021062993335,0.081021062993335,0.081021062993335,0.081021062993335,0.081021062993335,0.081021062993335];
GL_N6 = (Q_GLobatto_rq_N6, Q_GLobatto_sq_N6, Q_GLobatto_wq_N6)

#N=1
 Tri_diage_N1 = 7;
 Tri_diage_rq_N1 = [0.0, 1.0, 0.0, 0.5, 0.5, 0.0, 0.3333333333333333];
 Tri_diage_sq_N1 = [0.0, 0.0, 1.0, 0.0, 0.5, 0.5, 0.3333333333333333];
 Tri_diage_wq_N1 = [0.024999999999999998, 0.024999999999999998, 0.024999999999999998, 0.06666666666666667, 0.06666666666666667, 0.06666666666666667, 0.22500000000000006];
 Tri_diage_N1 = (2*Tri_diage_rq_N1 .-1, 2*Tri_diage_sq_N1 .-1, 4*Tri_diage_wq_N1)
#N=2
 Tri_diage_N2 = 12;
 Tri_diage_rq_N2 = [0.0, 1.0, 0.0, 0.21285435711180825, 0.21285435711180825, 0.5742912857763836, 0.27639320225002106, 0.7236067977499789, 0.7236067977499789, 0.27639320225002106, 0.0, 0.0];
 Tri_diage_sq_N2 = [0.0, 0.0, 1.0, 0.5742912857763836, 0.21285435711180825, 0.21285435711180825, 0.0, 0.0, 0.27639320225002106, 0.7236067977499789, 0.7236067977499789, 0.27639320225002106];
 Tri_diage_wq_N2 = [0.006261126504899741, 0.006261126504899741, 0.006261126504899741, 0.10675793966098839, 0.10675793966098839, 0.10675793966098839, 0.026823800250389242, 0.026823800250389242, 0.026823800250389242, 0.026823800250389242, 0.026823800250389242, 0.026823800250389242];
 Tri_diage_N2 = (2*Tri_diage_rq_N2 .-1, 2*Tri_diage_sq_N2 .-1, 4*Tri_diage_wq_N2)
#N=3
 Tri_diage_N3 = 18;
 Tri_diage_rq_N3 = [0.0, 1.0, 0.0, 0.5, 0.5, 0.0, 0.4243860251718814, 0.4243860251718814, 0.1512279496562372, 0.14200508409677795, 0.14200508409677795, 0.7159898318064442, 0.17267316464601146, 0.8273268353539885, 0.8273268353539885, 0.17267316464601146, 0.0, 0.0];
 Tri_diage_sq_N3 = [0.0, 0.0, 1.0, 0.0, 0.5, 0.5, 0.1512279496562372, 0.4243860251718814, 0.4243860251718814, 0.7159898318064442, 0.14200508409677795, 0.14200508409677795, 0.0, 0.0, 0.17267316464601146, 0.8273268353539885, 0.8273268353539885, 0.17267316464601146];
 Tri_diage_wq_N3 = [0.0022825661430496253, 0.0022825661430496253, 0.0022825661430496253, 0.015504052643022513, 0.015504052643022513, 0.015504052643022513, 0.07467669469983994, 0.07467669469983994, 0.07467669469983994, 0.051518167995569394, 0.051518167995569394, 0.051518167995569394, 0.011342592592592586, 0.011342592592592586, 0.011342592592592586, 0.011342592592592586, 0.011342592592592586, 0.011342592592592586];
 Tri_diage_N3 = (2*Tri_diage_rq_N3 .-1, 2*Tri_diage_sq_N3 .-1, 4*Tri_diage_wq_N3)

#N=4
 Tri_diage_N4 = 27;
 Tri_diage_rq_N4 = [0.0, 1.0, 0.0, 0.10367750814280517, 0.10367750814280517, 0.7926449837143896, 0.2653313804842097, 0.2653313804842097, 0.46933723903158064, 0.35738424175967753, 0.6426157582403225, 0.6426157582403225, 0.35738424175967753, 0.0, 0.0, 0.11747233803526758, 0.8825276619647324, 0.8825276619647324, 0.11747233803526758, 0.0, 0.0, 0.5870855671333673, 0.3246404722650515, 0.3246404722650515, 0.5870855671333673, 0.0882739606015811, 0.0882739606015811];
 Tri_diage_sq_N4 = [0.0, 0.0, 1.0, 0.7926449837143896, 0.10367750814280517, 0.10367750814280517, 0.46933723903158064, 0.2653313804842097, 0.2653313804842097, 0.0, 0.0, 0.35738424175967753, 0.6426157582403225, 0.6426157582403225, 0.35738424175967753, 0.0, 0.0, 0.11747233803526758, 0.8825276619647324, 0.8825276619647324, 0.11747233803526758, 0.0882739606015811, 0.0882739606015811, 0.5870855671333673, 0.3246404722650515, 0.3246404722650515, 0.5870855671333673];
 Tri_diage_wq_N4 = [0.001090393904993471, 0.001090393904993471, 0.001090393904993471, 0.02839719066391149, 0.02839719066391149, 0.02839719066391149, 0.03996004802785181, 0.03996004802785181, 0.03996004802785181, 0.0069669428714637, 0.0069669428714637, 0.0069669428714637, 0.0069669428714637, 0.0069669428714637, 0.0069669428714637, 0.005519747637357106, 0.005519747637357106, 0.005519747637357106, 0.005519747637357106, 0.005519747637357106, 0.005519747637357106, 0.03612282652613417, 0.03612282652613417, 0.03612282652613417, 0.03612282652613417, 0.03612282652613417, 0.03612282652613417];
 Tri_diage_N4 = (2*Tri_diage_rq_N4 .-1, 2*Tri_diage_sq_N4 .-1, 4*Tri_diage_wq_N4)

GQ_SBP = (GQ_N1, GQ_N2, GQ_N3, GQ_N4, GQ_N5, GQ_N6);
GL_SBP = (GL_N1, GL_N2, GL_N3, GL_N4, GL_N5, GL_N6);
Tri_diage = (Tri_diage_N1, Tri_diage_N2, Tri_diage_N3, Tri_diage_N4);

ext: jl

max_stars_repo_name: yiminllin/ESDG-PosLimit

max_issues_repo_name: yiminllin/ESDG-PosLimit

max_forks_repo_name: yiminllin/ESDG-PosLimit


text:
function Header(io::IOStream)
    # check magic number
    (read(io, UInt8) == 0x6c && read(io, UInt8) == 0x1b) || throw(ArgumentError("wrong magic number in the input file"))
    # 2 bytes so far
    storage_mode = read(io, UInt8) # 3 bytes
    storage_mode == 0x10 || throw(ArgumentError("Only the standard format supported for now"))
    n_variants = read(io, UInt32) # 7 bytes
    n_samples = read(io, UInt32) # 11 bytes

    if n_samples <= 2 ^ 8
        bytes_per_sample_id = 1
    elseif n_samples <= 2 ^ 16
        bytes_per_sample_id = 2
    elseif n_samples <= 2 ^ 24
        bytes_per_sample_id = 3
    else
        bytes_per_sample_id = 4
    end

    format_num = read(io, UInt8) # 12 bytes

    # figure out bits_per_variant_type and bytes_per_record_length
    format_num_0_3 = format_num & 0x0f

    (bits_per_variant_type, bytes_per_record_length) = variant_type_lengths[format_num_0_3]

    bytes_allele_counts = (format_num & 0x30) >> 4

    # provisional reference
    provisional_reference = (format_num & 0xc0) >> 6

    n_blocks = ceil_int(n_variants, 2 ^ 16) #Int(ceil(n_variants / 2 ^ 16))
    variant_block_offsets = Vector{UInt64}(undef, n_blocks)
    read!(io, variant_block_offsets) # 12 + 8n_blocks
    #variant_block_offsets = reinterpret(UInt64, data[13:(12 + n_blocks * 8)])

    offset = convert(UInt64, 12 + 8n_blocks)

    t_variant_types = nothing
    sectors_variant_types = []

    dtype_variant_length = bytes_to_UInt[bytes_per_record_length]
    t_variant_sizes = nothing
    sectors_variant_lengths = []

    if bytes_allele_counts != 0
        dtype_allele_counts = bytes_to_UInt[bytes_allele_counts]
        t_allele_counts = nothing
        sectors_allele_counts = []
    else
        sectors_allele_counts = nothing
    end
    if provisional_reference == 0x03
        t_provisional_reference_flags = nothing
        sectors_provisional_reference = []
    else
        sectors_provisional_reference = nothing
    end
    for i in 1:n_blocks
        block_size = i < n_blocks ? convert(UInt, 2 ^ 16) : begin
            rem = n_variants % convert(UInt, 2 ^ 16)
            rem > 0 ? rem : convert(UInt, 2 ^ 16)
        end

        # read variant type track
        data_per_byte = 8 ÷ bits_per_variant_type
        size_variant_types = ceil_int(block_size, data_per_byte)
        #convert(UInt, ceil(block_size / (8 ÷ bits_per_variant_type)))
        arr = read(io, size_variant_types) #data[offset + 1 : offset + size_variant_types]
        push!(sectors_variant_types, Ref(arr))
        t_variant_types = typeof(arr)
        offset += size_variant_types

        # read variant record length track
        arr = read(io, block_size * bytes_per_record_length) #data[offset + 1 : offset + block_size * bytes_per_record_length]
        reinterpreted = reinterpret(dtype_variant_length, arr)
        t_variant_sizes = typeof(reinterpreted)
        push!(sectors_variant_lengths, Ref(reinterpreted))
        offset += bytes_per_record_length * block_size

        # read allele counts track
        if sectors_allele_counts !== nothing
            arr = read(io, block_size * bytes_allele_counts)
            #arr = data[offset + 1:offset + block_size * bytes_allele_counts]
            reinterpreted = reinterpret(dtype_allele_counts, arr)
            t_allele_counts = typeof(reinterpreted)
            push!(sectors_allele_counts, Ref(reinterpreted))
            offset += bytes_allele_counts * block_size
        end

        # read provisional reference flags track
        if sectors_provisional_reference !== nothing
            size_pr = ceil_int(block_size, 8)#convert(UInt, ceil(block_size / 8))
            arr = read(io, size_pr)
            #arr = data[offset + 1 : offset + size_pr]
            t_provisional_reference_flags = typeof(arr)
            push!(sectors_provisional_reference, Ref(arr))
            offset += size_pr
        end
    end
    if bits_per_variant_type == 8
        VTT = ScatteredVector{UInt8, t_variant_types}
        variant_types = ScatteredVector{UInt8, t_variant_types}(
            sectors_variant_types, n_variants, n_blocks)
    else
        VTT = ScatteredBitsVector{t_variant_types}
        variant_types = ScatteredBitsVector{t_variant_types}(
            sectors_variant_types, bits_per_variant_type,
            n_variants, n_blocks)
    end
    VLT = ScatteredVector{dtype_variant_length,t_variant_sizes}
    variant_lengths = ScatteredVector{dtype_variant_length,t_variant_sizes}(
        sectors_variant_lengths,
        n_variants, n_blocks)
    if sectors_allele_counts !== nothing
        ACT = ScatteredVector{dtype_allele_counts,t_allele_counts}
        allele_counts = ScatteredVector{dtype_allele_counts,t_allele_counts}(
            sectors_allele_counts,
            n_variants, n_blocks)
    else
        ACT = Nothing
        allele_counts = nothing
    end
    if sectors_provisional_reference !== nothing
        PRFT = ScatteredBitsVector{t_provisional_reference_flags}
        provisional_reference_flags = ScatteredBitsVector{t_provisional_reference_flags}(
            sectors_provisional_reference, 1,
            n_variants, n_blocks)
    else
        PRFT = Nothing
        provisional_reference_flags = nothing
    end

    # Store the most recent non-LD-compressed variant for each LD-compressed variant
    # Useful for random-access of variants
    most_recent_non_ld = Dict{UInt, Variant}()
    buf = nothing
    offset = variant_block_offsets[1]
    for j in 1:n_variants
        vt = variant_types[j]
        vl = variant_lengths[j]
        if !((vt & 0x07 == 0x02) || (vt & 0x07 == 0x03))
            buf = Variant(j, offset, vt, vl)
        else
            most_recent_non_ld[j] = buf
        end
        offset += vl
    end

    Header{VTT, VLT, ACT, PRFT}(storage_mode, n_variants, n_samples,
        bits_per_variant_type, bytes_per_record_length,
        bytes_allele_counts, bytes_per_sample_id, provisional_reference,
        n_blocks, variant_block_offsets,
        variant_types, variant_lengths, allele_counts, provisional_reference_flags,
        most_recent_non_ld)
end

function Header(filename::String)
    io = open(filename)
    Header(io)
end


text:
using OSF
using Documenter

DocMeta.setdocmeta!(OSF, :DocTestSetup, :(using OSF); recursive=true)

makedocs(;
    modules=[OSF],
    authors="Zaki A <zaki@live.ca> and contributors",
    repo="https://github.com/notZaki/OSF.jl/blob/{commit}{path}#{line}",
    sitename="OSF.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://notZaki.github.io/OSF.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/notZaki/OSF.jl",
)


text:
using RegressionBenchmarks, Distributions, MLDataUtils, DataFrames, Gurobi,
    CPLEX

srand(1)

# The scale of data we care about
nrange = collect(100:20:500)

# What form will the correlation in X take
Xcorr = NoCorrelation()
# True sparsity
sparsity = 10

noisedist = Normal(0.0, 1.0)

# What distribution will X come from
for Xdist in [Uniform, MvNormal]
    for wdist in [BinChoice(), Uniform(), Normal(0.0, 1.0)]
        for SNR in [20.0, 5.0, 1.0, 0.5]
              for d in [1000, 5000, 10000]

                  # Create an object for our X data
                  Xdata = XData(Xdist, Xcorr, d)

                  # Create our data object
                  bd = BenchmarkData(Xdata = Xdata,
                                    wdist = wdist,
                                    noisedist = noisedist,
                                    SNR = SNR,
                                    n = nrange,
                                    nfeatures = d,
                                    sparsity = sparsity)
                  # Choose model
                  for m in [RelaxDualSubgradient(0.0, PolyakStepping(2.0, 30), 100);
                            RelaxDualCutting(0.0, GurobiSolver(OutputFlag=0))]
                      # Get results
                      results_table, val_results = benchmark(bd, m)
                      # Save
                      resdir = joinpath(Pkg.dir("RegressionBenchmarks"), "results")
                      !isdir(resdir) && mkdir(resdir)
                      tstamp = round(Int64, time() * 1000)
                      datadir = joinpath(resdir, data2str(bd), method2str(m), "$(tstamp)")
                      !ispath(datadir) && mkpath(datadir)
                      df = convert(DataFrame, results_table)
                      colnames = [:accuracy, :false_detection, :train_r2, :test_r2, :time, :gamma]
                      names!(df, colnames)
                      insert!(df, 1, nrange, :nrange)
                      writetable(joinpath(datadir, "sparsity$(sparsity)_results.csv"), df)
                      # Archive validation results also
                      open(joinpath(datadir, "sparsity$(sparsity)_validation.csv"), "w") do io
                          RegressionBenchmarks.valid2io(io, val_results, nrange)
                      end
                      # Make plots
                      for field in ["accuracy", "false_detection", "test_r2", "time"]
                          makeplot(datadir, field, df, sparsity)
                      end
                  end # method
              end # d
        end  # SNR
    end  # wdist
end # xdist


text:
using SDDiP, GLPKMathProgInterface, JuMP, SDDP, Base.Test
# using Ipopt  # uncomment for the Level Method

struct AllBlacksData
  T::Int                  # Number of time periods
  N::Int                  # Number of seats
  R::Array{Float64,2}     # R_ij = evenue from selling seat i at time j
  offer::Array{Float64,2} # offer_ij = whether an offer for seat i will come at time j
end

data = AllBlacksData(3, 2, [3 3 6; 3 3 6], [1 1 0; 1 0 1])

m=SDDPModel(stages=data.T, objective_bound=100.0, sense=:Max, solver=GLPKSolverMIP()) do sp, stage

    # Seat remaining?
    @binarystate(sp, x[i=1:data.N], x0==1, Bin)

    # Action: accept offer, or don't accept offer
    @variable(sp, accept_offer, Bin)

    @constraints(sp, begin
        # Balance on seats
        dynamics[i=1:data.N], x[i] == x0[i] - data.offer[i, stage] * accept_offer
    end)

    @stageobjective(sp, sum(data.R[i, stage] * data.offer[i, stage] * accept_offer for i=1:data.N))

    # Level method example:
    # setSDDiPsolver!(sp, method=LevelMethod(quadsolver=IpoptSolver(print_level=0)))
    # Subgradient descent example:
    setSDDiPsolver!(sp, method=SubgradientMethod())
    # setSDDiPsolver!(sp, method=KelleyMethod())

end

srand(11111)
solution = solve(m, max_iterations=8)
@test isapprox(getbound(m), 9.0)


text:
using InfiniteOpt, Ipopt, Plots
#using KNITRO

exponential_relaxation(z, δ) = exp(one(typeof(z)) - z / δ) - one(typeof(z)) - log(δ)
function relaxed_log_barrier(z; δ=0.3f0)
    return max(z > δ ? -log(z) : exponential_relaxation(z, δ), zero(typeof(z)))
end
# Set the SEIR parameters
γ = 0.303
β = 0.727
N = 1e4
ξ = 0.2

# Set the domain information
i_max = 0.001
t0 = 0
tf = 50

# Set the intial condition values
e0 = 1 / N
i0 = 0
r0 = 0
s0 = 1 - 1 / N

model = InfiniteModel(Ipopt.Optimizer)

@infinite_parameter(model, t ∈ [t0, tf], num_supports = 201,
    derivative_method = OrthogonalCollocation(2))

@variable(model, s ≥ 0, Infinite(t))
@variable(model, e ≥ 0, Infinite(t))
@variable(model, i ≥ 0, Infinite(t))
@variable(model, r ≥ 0, Infinite(t))
@variable(model, 0.2 ≤ u ≤ 0.8, Infinite(t), start = 0.2)

@objective(model, Min, ∫(u^2, t))

# Define the initial conditions
@constraint(model, s(0) == s0)
@constraint(model, e(0) == e0)
@constraint(model, i(0) == i0)
@constraint(model, r(0) == r0)

# Define the SEIR equations
@constraint(model, s_constr, ∂(s, t) == -(1 - u) * β * s * i)
@constraint(model, e_constr, ∂(e, t) == (1 - u) * β * s * i - ξ * e)
@constraint(model, i_constr, ∂(i, t) == ξ * e - γ * i)
@constraint(model, r_constr, ∂(r, t) == γ * i)

# Define the infection rate limit
@constraint(model, imax_constr, i ≤ i_max)

print(model)

optimize!(model)

r_opt = value(r, ndarray=true) * 100 # make the population fractions into percentages
s_opt = value(s, ndarray=true) * 100
i_opt = value(i, ndarray=true) * 100
e_opt = value(e, ndarray=true) * 100
u_opt = value(u)
@show obj_opt = objective_value(model)
ts = value(t)

p = plot(ts, i_opt, label="i(t)")
plot!(ts, e_opt, label="e(t)")
ylabel!("Pop. (%)")
xlabel!("Time (Days)")
display(p)

p1 = plot(ts, u_opt, label="u(t)")
xlabel!("Time (Days)")
ylabel!("Distancing Ratio")
display(p1)

using DiffEqFlux, DifferentialEquations
const solveeq = DifferentialEquations.solve

tspan = (0.0f0, 50.0f0)
ann = FastChain(FastDense(1, 32, relu), FastDense(32, 1), (x, p) -> 0.6 * sigmoid.(x) .+ 0.2)
θ = initial_params(ann)
function sir_nn!(du, u, p, t)
    S, E, I, R, z = u
    du[1] = -(1 - ann([t], p)[1]) * β * S * I
    du[2] = (1 - ann([t], p)[1]) * β * S * I - ξ * E
    du[3] = ξ * E - γ * I
    du[4] = γ * I
    du[5] = (ann([t], p)[1])^2
    du[6] = (relu(I - i_max))^2
end
u0 = [s0, e0, i0, r0, 0, 0]
ts = Float32.(collect(0.0:1:tspan[2]))
prob = ODEProblem(sir_nn, u0, tspan, θ)
sol_init = solveeq(prob, Vern9())
plot(sol_init)
function predict_adjoint(θ)
    Array(solveeq(prob, Vern9(), p=θ, saveat=ts, sensealg=InterpolatingAdjoint(autojacvec=ReverseDiffVJP(true))))
end

# function quadratic_relaxation(z, δ)
#     return one(typeof(z)) / 2 * (((z - 2δ) / δ)^2 - one(typeof(z))) - log(δ)
# end
#
λ = 10000
for ii = 1:2
    λ = 10* λ
    function loss_adjoint(θ)
        x = predict_adjoint(θ)
        # I_barrier_loss = map(y -> relaxed_log_barrier(i_max - y; δ = 100 / N), x[3, :])
        λ * x[6, end] + x[5, end]
    end
    @show l = loss_adjoint(θ)
    cb = function (θ, l)
        println(l)
        p = plot(Array(solveeq(remake(prob, p=θ), Tsit5(), saveat=0.5))[3, :], lw=3)
        p1 = plot(ts, [first(ann([t], θ)) for t in ts], label="u(t)", lw=3)
        display(p)
        display(p1)
        return false
    end
    # Display the ODE with the current parameter values.
    #cb(θ, l)
    #loss1 = loss_adjoint(θ)
    res1 = DiffEqFlux.sciml_train(loss_adjoint, θ, ADAM(0.05), cb=cb, maxiters=100)
    θ = res1.u
    res2 = DiffEqFlux.sciml_train(loss_adjoint, θ,
        BFGS(initial_stepnorm=0.01), cb=cb, maxiters=100,
        allow_f_increases=false)
    θ = res2.u
end
p = plot(Array(solveeq(remake(prob, p=θ), Tsit5(), saveat=0.5))[3, :], lw=3)
# Part Three: Saving the results
display(p)
using DataFrames
using CSV
p_value = res1.u
results_saving_opt = DataFrame()
results_saving_opt.t = tspan
results_saving_opt.I = I_value
results_saving_opt.u = u_value
results_saving_opt.I_nn = Array(solveeq(remake(prob, p=p_value), Tsit5(), saveat=0.2))[3, :]
results_saving_opt.u_nn = [first(ann([t], p_value)) for t in ts]
CSV.write("results_saving_SEIRconstraintmodel.csv", results_saving_opt)

plot(results_saving_opt.t, results_saving_opt.I, lw=2, label="I(t) by direct collection")
plot!(results_saving_opt.t, results_saving_opt.I_nn, lw=2, label="I(t) by deep learning")
xlabel!("t(day)")
ylabel!("I(t)")
savefig("SEIRconstraint_I.png")

plot(results_saving_opt.t[2:end-1], results_saving_opt.u[2:end-1], lw=2, label="u(t) by direct collection")
plot!(results_saving_opt.t[2:end-1], results_saving_opt.u_nn[2:end-1], lw=2, label="u(t) by deep learning")
xlabel!("t(day)")
ylabel!("u(t)")
savefig("SEIRconstraint_control.png")

# Loss saving

opt_method_loss = 13.90369069230782
nn_method_loss = 13.711496454176164


text:
# RUN: julia -e "import Brutus; Brutus.lit(:emit_lowered)" --startup-file=no %s 2>&1 | FileCheck %s

index(A, i) = A[i]
emit(index, Array{Int64, 3}, Int64)
#

# CHECK: Core.MethodMatch(Tuple{typeof(Main.Main.index), Array{Int64, 3}, Int64}, svec(), index(A, i) in Main.Main at /home/mccoy/Dev/brutus/test/Codegen/lower/index_linear_3D.jl:3, true)after translating to MLIR in JLIR dialect:module  {
# CHECK-NEXT:   func nested @"Tuple{typeof(Main.index), Array{Int64, 3}, Int64}"(%arg0: !jlir<"typeof(Main.index)">, %arg1: !jlir<"Array{Int64, 3}">, %arg2: !jlir.Int64) -> !jlir.Int64 attributes {llvm.emit_c_interface} {
# CHECK-NEXT:     "jlir.goto"()[^bb1] : () -> ()
# CHECK-NEXT:   ^bb1:  // pred: ^bb0
# CHECK-NEXT:     %0 = "jlir.constant"() {value = #jlir.Core.arrayref} : () -> !jlir<"typeof(Core.arrayref)">
# CHECK-NEXT:     %1 = "jlir.constant"() {value = #jlir.true} : () -> !jlir.Bool
# CHECK-NEXT:     %2 = "jlir.call"(%0, %1, %arg1, %arg2) : (!jlir<"typeof(Core.arrayref)">, !jlir.Bool, !jlir<"Array{Int64, 3}">, !jlir.Int64) -> !jlir.Int64
# CHECK-NEXT:     "jlir.return"(%2) : (!jlir.Int64) -> ()
# CHECK-NEXT:   }
# CHECK-NEXT: }

# CHECK: module  {
# CHECK-NEXT:   func nested @"Tuple{typeof(Main.index), Array{Int64, 3}, Int64}"(%arg0: !jlir<"typeof(Main.index)">, %arg1: memref<?x?x?xi64>, %arg2: i64) -> i64 attributes {llvm.emit_c_interface} {
# CHECK-NEXT:     %c0 = constant 0 : index
# CHECK-NEXT:     %c1 = constant 1 : index
# CHECK-NEXT:     %0 = "jlir.convertstd"(%arg2) : (i64) -> index
# CHECK-NEXT:     %1 = subi %0, %c1 : index
# CHECK-NEXT:     %2 = load %arg1[%c0, %c0, %1] : memref<?x?x?xi64>
# CHECK-NEXT:     return %2 : i64
# CHECK-NEXT:   }
# CHECK-NEXT: }

# CHECK:   llvm.func @"Tuple{typeof(Main.index), Array{Int64, 3}, Int64}"(%arg0: !llvm.ptr<struct<"struct_jl_value_type", opaque>>, %arg1: !llvm.ptr<i64>, %arg2: !llvm.ptr<i64>, %arg3: i64, %arg4: i64, %arg5: i64, %arg6: i64, %arg7: i64, %arg8: i64, %arg9: i64, %arg10: i64) -> i64 attributes {llvm.emit_c_interface, sym_visibility = "nested"} {
# CHECK-NEXT:     %0 = llvm.mlir.undef : !llvm.struct<(ptr<i64>, ptr<i64>, i64, array<3 x i64>, array<3 x i64>)>
# CHECK-NEXT:     %1 = llvm.insertvalue %arg1, %0[0] : !llvm.struct<(ptr<i64>, ptr<i64>, i64, array<3 x i64>, array<3 x i64>)>
# CHECK-NEXT:     %2 = llvm.insertvalue %arg2, %1[1] : !llvm.struct<(ptr<i64>, ptr<i64>, i64, array<3 x i64>, array<3 x i64>)>
# CHECK-NEXT:     %3 = llvm.insertvalue %arg3, %2[2] : !llvm.struct<(ptr<i64>, ptr<i64>, i64, array<3 x i64>, array<3 x i64>)>
# CHECK-NEXT:     %4 = llvm.insertvalue %arg4, %3[3, 0] : !llvm.struct<(ptr<i64>, ptr<i64>, i64, array<3 x i64>, array<3 x i64>)>
# CHECK-NEXT:     %5 = llvm.insertvalue %arg7, %4[4, 0] : !llvm.struct<(ptr<i64>, ptr<i64>, i64, array<3 x i64>, array<3 x i64>)>
# CHECK-NEXT:     %6 = llvm.insertvalue %arg5, %5[3, 1] : !llvm.struct<(ptr<i64>, ptr<i64>, i64, array<3 x i64>, array<3 x i64>)>
# CHECK-NEXT:     %7 = llvm.insertvalue %arg8, %6[4, 1] : !llvm.struct<(ptr<i64>, ptr<i64>, i64, array<3 x i64>, array<3 x i64>)>
# CHECK-NEXT:     %8 = llvm.insertvalue %arg6, %7[3, 2] : !llvm.struct<(ptr<i64>, ptr<i64>, i64, array<3 x i64>, array<3 x i64>)>
# CHECK-NEXT:     %9 = llvm.insertvalue %arg9, %8[4, 2] : !llvm.struct<(ptr<i64>, ptr<i64>, i64, array<3 x i64>, array<3 x i64>)>
# CHECK-NEXT:     %10 = llvm.mlir.constant(0 : index) : i64
# CHECK-NEXT:     %11 = llvm.mlir.constant(1 : index) : i64
# CHECK-NEXT:     %12 = llvm.sub %arg10, %11  : i64
# CHECK-NEXT:     %13 = llvm.extractvalue %9[1] : !llvm.struct<(ptr<i64>, ptr<i64>, i64, array<3 x i64>, array<3 x i64>)>
# CHECK-NEXT:     %14 = llvm.extractvalue %9[4, 0] : !llvm.struct<(ptr<i64>, ptr<i64>, i64, array<3 x i64>, array<3 x i64>)>
# CHECK-NEXT:     %15 = llvm.mul %10, %14  : i64
# CHECK-NEXT:     %16 = llvm.extractvalue %9[4, 1] : !llvm.struct<(ptr<i64>, ptr<i64>, i64, array<3 x i64>, array<3 x i64>)>
# CHECK-NEXT:     %17 = llvm.mul %10, %16  : i64
# CHECK-NEXT:     %18 = llvm.add %15, %17  : i64
# CHECK-NEXT:     %19 = llvm.add %18, %12  : i64
# CHECK-NEXT:     %20 = llvm.getelementptr %13[%19] : (!llvm.ptr<i64>, i64) -> !llvm.ptr<i64>
# CHECK-NEXT:     %21 = llvm.load %20 : !llvm.ptr<i64>
# CHECK-NEXT:     llvm.return %21 : i64
# CHECK-NEXT:   }
# CHECK-NEXT:   llvm.func @"_mlir_ciface_Tuple{typeof(Main.index), Array{Int64, 3}, Int64}"(%arg0: !llvm.ptr<struct<"struct_jl_value_type", opaque>>, %arg1: !llvm.ptr<struct<(ptr<i64>, ptr<i64>, i64, array<3 x i64>, array<3 x i64>)>>, %arg2: i64) -> i64 attributes {llvm.emit_c_interface, sym_visibility = "nested"} {
# CHECK-NEXT:     %0 = llvm.load %arg1 : !llvm.ptr<struct<(ptr<i64>, ptr<i64>, i64, array<3 x i64>, array<3 x i64>)>>
# CHECK-NEXT:     %1 = llvm.extractvalue %0[0] : !llvm.struct<(ptr<i64>, ptr<i64>, i64, array<3 x i64>, array<3 x i64>)>
# CHECK-NEXT:     %2 = llvm.extractvalue %0[1] : !llvm.struct<(ptr<i64>, ptr<i64>, i64, array<3 x i64>, array<3 x i64>)>
# CHECK-NEXT:     %3 = llvm.extractvalue %0[2] : !llvm.struct<(ptr<i64>, ptr<i64>, i64, array<3 x i64>, array<3 x i64>)>
# CHECK-NEXT:     %4 = llvm.extractvalue %0[3, 0] : !llvm.struct<(ptr<i64>, ptr<i64>, i64, array<3 x i64>, array<3 x i64>)>
# CHECK-NEXT:     %5 = llvm.extractvalue %0[3, 1] : !llvm.struct<(ptr<i64>, ptr<i64>, i64, array<3 x i64>, array<3 x i64>)>
# CHECK-NEXT:     %6 = llvm.extractvalue %0[3, 2] : !llvm.struct<(ptr<i64>, ptr<i64>, i64, array<3 x i64>, array<3 x i64>)>
# CHECK-NEXT:     %7 = llvm.extractvalue %0[4, 0] : !llvm.struct<(ptr<i64>, ptr<i64>, i64, array<3 x i64>, array<3 x i64>)>
# CHECK-NEXT:     %8 = llvm.extractvalue %0[4, 1] : !llvm.struct<(ptr<i64>, ptr<i64>, i64, array<3 x i64>, array<3 x i64>)>
# CHECK-NEXT:     %9 = llvm.extractvalue %0[4, 2] : !llvm.struct<(ptr<i64>, ptr<i64>, i64, array<3 x i64>, array<3 x i64>)>
# CHECK-NEXT:     %10 = llvm.call @"Tuple{typeof(Main.index), Array{Int64, 3}, Int64}"(%arg0, %1, %2, %3, %4, %5, %6, %7, %8, %9, %arg2) : (!llvm.ptr<struct<"struct_jl_value_type", opaque>>, !llvm.ptr<i64>, !llvm.ptr<i64>, i64, i64, i64, i64, i64, i64, i64, i64) -> i64
# CHECK-NEXT:     llvm.return %10 : i64
# CHECK-NEXT:   }
# CHECK-NEXT: }

max_stars_repo_name: JuliaLabs/brutus

max_issues_repo_name: JuliaLabs/brutus

max_forks_repo_name: JuliaLabs/brutus


text:
# utilities.jl
#
# A few simple bindings because I am lazy.

function linspace(a::Real, b::Real, n::Int64)
  return LinRange(a, b, n);
end

function logspace(a::Real, b::Real, n::Int64)
  return 10.0.^LinRange(a, b, n);
end

max_stars_repo_path: Machine Learning Summer School 2019 (London, UK)/tutorials/optimization/utilities.jl


text:
using LightGraphs, SimpleWeightedGraphs
d = map.(parse, matchall.(r"\d+", readlines("input.txt")))
g = SimpleWeightedGraph(length(d))
foreach(x -> foreach(y -> add_edge!(g, x[1]+1, y+1), x), d)
comps = connected_components(g)
println(length(comps[1]))
println(length(comps))


text:
module ImmersedBoundaries

export ImmersedBoundaryGrid, GridFittedBoundary, GridFittedBottom, ImmersedBoundaryCondition

using Adapt

using Oceananigans.Grids
using Oceananigans.Operators
using Oceananigans.Fields
using Oceananigans.Utils

using Oceananigans.TurbulenceClosures: AbstractTurbulenceClosure, time_discretization
using Oceananigans.Grids: size_summary, inactive_node, peripheral_node

using Oceananigans.TurbulenceClosures:
    viscous_flux_ux,
    viscous_flux_uy,
    viscous_flux_uz,
    viscous_flux_vx,
    viscous_flux_vy,
    viscous_flux_vz,
    viscous_flux_wx,
    viscous_flux_wy,
    viscous_flux_wz,
    diffusive_flux_x,
    diffusive_flux_y,
    diffusive_flux_z

using Oceananigans.Advection:
    advective_momentum_flux_Uu,
    advective_momentum_flux_Uv,
    advective_momentum_flux_Uw,
    advective_momentum_flux_Vu,
    advective_momentum_flux_Vv,
    advective_momentum_flux_Vw,
    advective_momentum_flux_Wu,
    advective_momentum_flux_Wv,
    advective_momentum_flux_Ww,
    advective_tracer_flux_x,
    advective_tracer_flux_y,
    advective_tracer_flux_z,
    WENOVectorInvariant

import Base: show, summary
import Oceananigans.Utils: cell_advection_timescale

import Oceananigans.Grids:
        cpu_face_constructor_x,
        cpu_face_constructor_y,
        cpu_face_constructor_z

import Oceananigans.Grids: architecture, on_architecture, with_halo
import Oceananigans.Grids: xnode, ynode, znode, all_x_nodes, all_y_nodes, all_z_nodes
import Oceananigans.Grids: inactive_cell
import Oceananigans.Coriolis: φᶠᶠᵃ

import Oceananigans.Advection:
    _advective_momentum_flux_Uu,
    _advective_momentum_flux_Uv,
    _advective_momentum_flux_Uw,
    _advective_momentum_flux_Vu,
    _advective_momentum_flux_Vv,
    _advective_momentum_flux_Vw,
    _advective_momentum_flux_Wu,
    _advective_momentum_flux_Wv,
    _advective_momentum_flux_Ww,
    _advective_tracer_flux_x,
    _advective_tracer_flux_y,
    _advective_tracer_flux_z

import Oceananigans.TurbulenceClosures:
    _viscous_flux_ux,
    _viscous_flux_uy,
    _viscous_flux_uz,
    _viscous_flux_vx,
    _viscous_flux_vy,
    _viscous_flux_vz,
    _viscous_flux_wx,
    _viscous_flux_wy,
    _viscous_flux_wz,
    _diffusive_flux_x,
    _diffusive_flux_y,
    _diffusive_flux_z,
    κᶠᶜᶜ,
    κᶜᶠᶜ,
    κᶜᶜᶠ,
    νᶜᶜᶜ,
    νᶠᶠᶜ,
    νᶜᶠᶠ,
    νᶠᶜᶠ

"""
    abstract type AbstractImmersedBoundary

Abstract supertype for immersed boundary grids.
"""
abstract type AbstractImmersedBoundary end

#####
##### ImmersedBoundaryGrid
#####

struct ImmersedBoundaryGrid{FT, TX, TY, TZ, G, I, Arch} <: AbstractGrid{FT, TX, TY, TZ, Arch}
    architecture :: Arch
    underlying_grid :: G
    immersed_boundary :: I

    # Internal interface
    function ImmersedBoundaryGrid{TX, TY, TZ}(grid::G, ib::I) where {TX, TY, TZ, G <: AbstractUnderlyingGrid, I}
        FT = eltype(grid)
        arch = architecture(grid)
        Arch = typeof(arch)
        return new{FT, TX, TY, TZ, G, I, Arch}(arch, grid, ib)
    end
end

const IBG = ImmersedBoundaryGrid

@inline Base.getproperty(ibg::IBG, property::Symbol) = get_ibg_property(ibg, Val(property))
@inline get_ibg_property(ibg::IBG, ::Val{property}) where property = getfield(getfield(ibg, :underlying_grid), property)
@inline get_ibg_property(ibg::IBG, ::Val{:immersed_boundary}) = getfield(ibg, :immersed_boundary)
@inline get_ibg_property(ibg::IBG, ::Val{:underlying_grid}) = getfield(ibg, :underlying_grid)

@inline architecture(ibg::IBG) = architecture(ibg.underlying_grid)

Adapt.adapt_structure(to, ibg::IBG{FT, TX, TY, TZ}) where {FT, TX, TY, TZ} =
    ImmersedBoundaryGrid{TX, TY, TZ}(adapt(to, ibg.underlying_grid), adapt(to, ibg.immersed_boundary))

with_halo(halo, ibg::ImmersedBoundaryGrid) = ImmersedBoundaryGrid(with_halo(halo, ibg.underlying_grid), ibg.immersed_boundary)

function Base.summary(grid::ImmersedBoundaryGrid)
    FT = eltype(grid)
    TX, TY, TZ = topology(grid)

    return string(size_summary(size(grid)),
                  " ImmersedBoundaryGrid{$FT, $TX, $TY, $TZ} on ", summary(architecture(grid)),
                  " with ", size_summary(halo_size(grid)), " halo")
end

function show(io::IO, ibg::ImmersedBoundaryGrid)
    print(io, summary(ibg), ":", '\n',
              "├── immersed_boundary: ", summary(ibg.immersed_boundary), '\n',
              "├── underlying_grid: ", summary(ibg.underlying_grid), '\n')

    return show(io, ibg.underlying_grid, false)
end

#####
##### Interface for immersed_boundary
#####

"""
    immersed_cell(i, j, k, grid)

Return true if a `cell` is "completely" immersed, and thus
is not part of the prognostic state.
"""
@inline immersed_cell(i, j, k, grid) = false

# Unpack to make defining new immersed boundaries more convenient
@inline immersed_cell(i, j, k, grid::ImmersedBoundaryGrid) =
    immersed_cell(i, j, k, grid.underlying_grid, grid.immersed_boundary)

"""
    inactive_cell(i, j, k, grid::ImmersedBoundaryGrid)

Return `true` if the tracer cell at `i, j, k` either (i) lies outside the `Bounded` domain
or (ii) lies within the immersed region of `ImmersedBoundaryGrid`.

Example
=======

Consider the configuration

```
   Immersed      Fluid
  =========== ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅

       c           c
      i-1          i

 | ========= |           |
 × === ∘ === ×     ∘     ×
 | ========= |           |

i-1          i
 f           f           f
```

We then have

    * `inactive_node(f, c, c, i, 1, 1, grid) = false`

As well as

    * `inactive_node(c, c, c, i,   1, 1, grid) = false`
    * `inactive_node(c, c, c, i-1, 1, 1, grid) = true`
    * `inactive_node(f, c, c, i-1, 1, 1, grid) = true`
"""
@inline inactive_cell(i, j, k, ibg::IBG) = immersed_cell(i, j, k, ibg) | inactive_cell(i, j, k, ibg.underlying_grid)

# Isolate periphery of the immersed boundary
@inline immersed_peripheral_node(LX, LY, LZ, i, j, k, ibg::IBG) =  peripheral_node(LX, LY, LZ, i, j, k, ibg) &
                                                                  !peripheral_node(LX, LY, LZ, i, j, k, ibg.underlying_grid)

#####
##### Utilities
#####

const IBG = ImmersedBoundaryGrid
const c = Center()
const f = Face()

@inline Base.zero(ibg::IBG) = zero(ibg.underlying_grid)
@inline cell_advection_timescale(u, v, w, ibg::IBG) = cell_advection_timescale(u, v, w, ibg.underlying_grid)
@inline φᶠᶠᵃ(i, j, k, ibg::IBG) = φᶠᶠᵃ(i, j, k, ibg.underlying_grid)

@inline xnode(LX, i, ibg::IBG) = xnode(LX, i, ibg.underlying_grid)
@inline ynode(LY, j, ibg::IBG) = ynode(LY, j, ibg.underlying_grid)
@inline znode(LZ, k, ibg::IBG) = znode(LZ, k, ibg.underlying_grid)

@inline xnode(LX, LY, LZ, i, j, k, ibg::IBG) = xnode(LX, LY, LZ, i, j, k, ibg.underlying_grid)
@inline ynode(LX, LY, LZ, i, j, k, ibg::IBG) = ynode(LX, LY, LZ, i, j, k, ibg.underlying_grid)
@inline znode(LX, LY, LZ, i, j, k, ibg::IBG) = znode(LX, LY, LZ, i, j, k, ibg.underlying_grid)

all_x_nodes(loc, ibg::IBG) = all_x_nodes(loc, ibg.underlying_grid)
all_y_nodes(loc, ibg::IBG) = all_y_nodes(loc, ibg.underlying_grid)
all_z_nodes(loc, ibg::IBG) = all_z_nodes(loc, ibg.underlying_grid)

@inline cpu_face_constructor_x(ibg::IBG) = cpu_face_constructor_x(ibg.underlying_grid)
@inline cpu_face_constructor_y(ibg::IBG) = cpu_face_constructor_y(ibg.underlying_grid)
@inline cpu_face_constructor_z(ibg::IBG) = cpu_face_constructor_z(ibg.underlying_grid)

function on_architecture(arch, ibg::IBG)
    underlying_grid   = on_architecture(arch, ibg.underlying_grid)
    immersed_boundary = on_architecture(arch, ibg.immersed_boundary)
    return ImmersedBoundaryGrid(underlying_grid, immersed_boundary)
end

isrectilinear(ibg::IBG) = isrectilinear(ibg.underlying_grid)

#####
##### Diffusivities (for VerticallyImplicit)
##### (the diffusivities on the immersed boundaries are kept)
#####

for (locate_coeff, loc) in ((:κᶠᶜᶜ, (f, c, c)),
                            (:κᶜᶠᶜ, (c, f, c)),
                            (:κᶜᶜᶠ, (c, c, f)),
                            (:νᶜᶜᶜ, (c, c, c)),
                            (:νᶠᶠᶜ, (f, f, c)),
                            (:νᶠᶜᶠ, (f, c, f)),
                            (:νᶜᶠᶠ, (c, f, f)))

    @eval begin
        @inline $locate_coeff(i, j, k, ibg::IBG{FT}, coeff) where FT =
            ifelse(inactive_node(loc..., i, j, k, ibg), $locate_coeff(i, j, k, ibg.underlying_grid, coeff), zero(FT))
    end
end

include("immersed_grid_metrics.jl")
include("grid_fitted_immersed_boundaries.jl")
include("conditional_fluxes.jl")
include("immersed_boundary_condition.jl")
include("conditional_derivatives.jl")
include("mask_immersed_field.jl")
include("immersed_reductions.jl")

end # module


text:
struct PoseData{T}
    configuration::TreeJointSegmentedVector{T}
    marker_positions::Dict{RigidBody{T}, Vector{Point3DS{T}}}
end

mutable struct CalibrationProblem{T}
    mechanism::Mechanism{T}
    calibration_param_bounds::Dict{<:Joint{T}, Vector{Tuple{T, T}}}
    free_joint_configuration_bounds::Dict{<:Joint{T}, Vector{Tuple{T, T}}}
    marker_location_bounds::Dict{RigidBody{T}, Vector{Tuple{Point3DS{T}, Point3DS{T}}}}
    pose_data::Vector{PoseData{T}}
    body_weights::Dict{RigidBody{T}, T}
    ordered_marker_bodies::Vector{RigidBody{T}}

    function CalibrationProblem(
            mechanism::Mechanism{T},
            calibration_param_bounds::Dict{<:Joint{T}, Vector{Tuple{T, T}}},
            free_joint_configuration_bounds::Dict{<:Joint{T}, Vector{Tuple{T, T}}},
            marker_location_bounds::Dict{RigidBody{T}, Vector{Tuple{Point3DS{T}, Point3DS{T}}}},
            pose_data::Vector{PoseData{T}},
            body_weights::Dict{RigidBody{T}, T} = Dict(b => one(T) for b in keys(marker_location_bounds))) where {T}
        @assert isempty(symdiff(keys(marker_location_bounds), keys(body_weights)))
        canonicalize!(marker_location_bounds)
        ordered_marker_bodies = intersect(bodies(mechanism), keys(body_weights))
        new{T}(mechanism, calibration_param_bounds, free_joint_configuration_bounds,
            marker_location_bounds, pose_data, body_weights, ordered_marker_bodies)
    end
end

num_poses(problem::CalibrationProblem) = length(problem.pose_data)
num_calibration_params(problem::CalibrationProblem) = sum(length, values(problem.calibration_param_bounds))
num_calibration_params(problem::CalibrationProblem, joint::Joint) = length(problem.calibration_param_bounds[joint])
num_markers(problem::CalibrationProblem) = sum(length, values(problem.marker_location_bounds))
num_markers(problem::CalibrationProblem, body::RigidBody) = length(problem.marker_location_bounds[body])
RigidBodyDynamics.num_bodies(problem::CalibrationProblem) = length(problem.body_weights)
calibration_joints(problem::CalibrationProblem) = keys(problem.calibration_param_bounds)
free_joints(problem::CalibrationProblem) = keys(problem.free_joint_configuration_bounds)

function Base.show(io::IO, problem::CalibrationProblem{T}) where {T}
    msg = """CalibrationProblem{$T} with:
    * $(num_poses(problem)) poses
    * $(num_calibration_params(problem)) calibration parameters
    * $(num_markers(problem)) markers attached to $(num_bodies(problem)) bodies
    """
    println(io, msg)
end


text:
# AVX-512 (64 bytes)
const shuffle_48_bytes_to_64 = ( shuffle_24_bytes_to_32..., shuffle_24_bytes_to_32 .+ 24... )
const shuffle_48_bytes_to_64_val = Val{shuffle_48_bytes_to_64}()
const shuffle_48_bytes_to_64_m8_val = Val{shuffle_48_bytes_to_64 .+ 8}()
const shuffle_48_bytes_to_64_m16_val = Val{shuffle_48_bytes_to_64 .+ 16}()
const choose_even_odd_ints_32 = Val( ntuple(i->isodd(i-1)*32 + i -1, 32) )

"""
    unpack_uint12_to_uint16_512bits(A::Vector{UInt8}, out::Vector{UInt16}, [ i ])

    Unpack 12-bit integers into 16-bit integers. Two 12-bit integers are packed consecutively into three bytes.
"""
@inline function unpack_uint12_to_uint16_512bits(
    A::SIMD.FastContiguousArray{UInt8,1},
    out::SIMD.FastContiguousArray{UInt16,1},
    i,
    offset = 0,
    shuffle_48_bytes_to_64::Val{T} = shuffle_48_bytes_to_64_val,
    ) where T
    # Load 32 bytes (we only use the first 24)
    a = @inbounds vload(Vec{64,UInt8}, A, i-offset)

    a16 = reinterpret( Vec{32, UInt16}, shufflevector(a, shuffle_48_bytes_to_64) )

    a16 = shufflevector(a16 & 0xfff, a16 >> 4, choose_even_odd_ints_32)

    @inbounds vstore(a16, out, 1 + 2( ( i-1 ) ÷ 3) )
end
@inline unpack_uint12_to_uint16_512bits_m8(A, out, i) = unpack_uint12_to_uint16_512bits(A, out, i, 8, shuffle_48_bytes_to_64_m8_val)
@inline unpack_uint12_to_uint16_512bits_m16(A, out, i) = unpack_uint12_to_uint16_512bits(A, out, i, 16, shuffle_48_bytes_to_64_m16_val)

"""
    unpack_uint12_to_uint16_512bits(A::SIMD.FastContiguousArray{UInt8,1}, out::Vector{UInt16})

    Unpack entire array with a preallocated out buffer for compressed data
"""
function unpack_uint12_to_uint16_512bits(A::SIMD.FastContiguousArray{UInt8,1}, out::SIMD.FastContiguousArray{UInt16,1})
    unpack_uint12_to_uint16_512bits(A, out, 1)
    idx = 49:48:length(A)
    @inbounds for i = idx[1:end-1]
       unpack_uint12_to_uint16_512bits_m8(A, out, i)
    end
    unpack_uint12_to_uint16_512bits_m16(A, out, length(A)-47)
    out
end

"""
    unpack_uint12_to_uint16_512bits(A::SIMD.FastContiguousArray{UInt8,1})

    Unpack entire array, allocates and returns output buffer
"""
function unpack_uint12_to_uint16_512bits(A::SIMD.FastContiguousArray{UInt8,1})
    out = Vector{UInt16}(undef, length(A) ÷ 3 * 2)
    unpack_uint12_to_uint16_512bits(A, out)
end


text:
#=
Lee Model Reference: Simon Saunders, “Antennas and Propagation for Wireless Communication Systems”, John Wiley & Sons,1999
Levis, C., Johnson, J. T., & Teixeira, F. L. (2010). Radiowave propagation: physics and applications. John Wiley & Sons.
=#

@doc """
 Lee Point-to-Point Model
 freq: signal frequency
 txH: base station's height
 rxH: mobile's height
 leeArea : area type
 FreeSpace       n=2.0 Po= 5.0
 OpenArea        n=4.35 Po=49.0
 SubUrban        n=3.84 Po=61.7
 Philadelphia    n=3.68 Po=70.0
 Newark=4        n=4.31 Po=64.0
 Tokyo=5         n=3.05 Po=84.0
 NewYorkCity =6  n=4.80  Po=77.0
"""
type LeeModel<:Model
  freq::Float64
  txH::Float64
  rxH::Float64
  leeArea::Tuple{Float64,Float64}
  LeeModel() = new(900,50,1.5,LeeArea.SubUrban)
end

"Calculate Path loss for distance dist (km) with Lee Point-to-Point Model"
function pathloss_op(model::LeeModel, dist::AbstractFloat)
  #FreeSpace = 0   # n=2.0 Po=5.0
  #OpenArea = 1    # n=4.35 Po=49.0
  #SubUrban= 2     # n=3.84 Po=61.7
  #Philadelphia=3  # n=3.68 Po=70.0
  #Newark=4        # n=4.31 Po=64.0
  #Tokyo=5         # n=3.05 Po=84.0
  #NewYorkCity =6  # n=4.8 Po=77.0
  hm::Float64 = model.rxH
  hb::Float64 = model.txH
  d::Float64 = dist
  f::Float64 = model.freq
  no::Float64 = model.leeArea[1]
  po::Float64 = model.leeArea[2]
   nf::Float64 = f>= 850.0 ? 3:2
   L1 = -20 *log10(hb/30)
   X2 = hm >= 3 ? 2:1
   L2 = -10*X2*log10(hm/3)
   Lo = 50.3 + po - 10*no*log10(1.61)-10*nf*log10(900)
   L = Lo + 10*no*log10(d)+10*nf*log10(f)+L1+L2
   return L
end # pathloss_op end


text:
infstep(x0, d) = convert(eltype(d), Inf)
stopbygradient(tol::Real) = (x, xpre, y, ypre, g) -> norm(g, 2) ≤ tol

#=
A generic way to create a mutable N×N matrix from a length-N vector.
Behaves like `similar(vec)` but returns `MMatrix` if the
argument is a `StaticVector` (`similar(staticvec, (n,n))` creates a
`Base.Matrix` which is less performant, and I don't want an explicit
dependency on `StaticArrays.jl` - V.P.)
=#
"""
    sqmatr(vec::AbstractVector, [element_type = eltype(vec)])

Create an uninitialized mutable N×N matrix with the given element type, given a
vector of length N.
"""
function sqmatr(vec::AbstractVector{T}, F::DataType = T) where {T}
    sz = axes(vec)
    return similar(vec, F, (sz..., sz...))
end

"""
Return the maximum absolute values of diagonal and off-diagonal elements of `A`.
"""
function γξ(A::AbstractMatrix)
    n = LinearAlgebra.checksquare(A)
    γ = ξ = zero(eltype(A))
    for j in 1:n
        for i in 1:j-1
            ξ = max(ξ, abs(A[i,j]))
        end
        γ = max(γ, abs(A[j,j]))
    end
    return γ, ξ
end

"""
    mcholesky!(A::AbstractMatrix)

Perform an in-place modified Cholesky decomposition on matrix `A`
(Gill, Murray, Wright, Practical optimization (1981), p.111)
"""
function mcholesky!(A::AbstractMatrix{T}; δ=eps(T)^(2/3)) where {T}
    n = LinearAlgebra.checksquare(A)
    γ, ξ = γξ(A)
    ν = max(1, sqrt(n^2 - 1))
    β² = max(γ, ξ / ν, eps(T))
    u = UpperTriangular(A)
    δA = δ * γ # similar to the choice of Schnabel-Eskow (1999)
    @inbounds for j in 1:n
        c_jj = u[j,j]
        θ = zero(c_jj)
        for k in 1:j-1
            u[k,j] /= u[k,k]
        end
        for i in j+1:n
            c_ij = u[j,i]
            for k in 1:j-1
                c_ij -= u[k,j] * u[k,i] / u[k,k]
            end
            θ = max(abs(c_ij), θ)
            u[j,i] = c_ij
        end
        d_j = max(abs(c_jj), θ^2 / β², δA)
        u_jj = sqrt(d_j)
        u[j,j] = u_jj
        for i in j+1:n
            u[i,i] -= (u[j,i] / u_jj)^2
        end
    end
    return Cholesky(A, 'U', 0)
end

#####LOGGING#####

function metafmt_noprefix_nosuffix(level::LogLevel, _module, group, id, file, line)
    @nospecialize
    if level < Info
        color = default_logcolor(level)
        prefix = string(
            level == OptLogLevel ? "Optimizer" :
            level == LSLogLevel ? "Linesearch" :
            string(level),
            ':'
        )
        suffix::String = ""
        return color, prefix, suffix
    else
        return default_metafmt(level, _module, group, id, file, line)
    end
end


text:
@testset "work items" begin
    @on_device get_work_dim()

    @on_device get_global_size(0)
    @on_device get_global_id(0)

    @on_device get_local_size(0)
    @on_device get_enqueued_local_size(0)
    @on_device get_local_id(0)

    @on_device get_num_groups(0)
    @on_device get_group_id(0)

    @on_device get_global_offset(0)

    @on_device get_global_linear_id()
    @on_device get_local_linear_id()
end

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

@testset "math" begin
    @testset "log10" begin
        @test testf(a->log10.(a), Float32[100])
    end

    for op in (exp, exp2, exp10, expm1)
        @testset "$op" begin
            for T in (Float32, Float64)
                @test testf(x->op.(x), rand(T, 1))
                @test testf(x->op.(x), -rand(T, 1))
            end

        end
    end
end

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

endline = Sys.iswindows() ? "\r\n" : "\n"

@testset "formatted output" begin
    _, out = @grab_output @on_device oneAPI.@printf("")
    @test out == ""

    _, out = @grab_output @on_device oneAPI.@printf("Testing...\n")
    @test out == "Testing...$endline"

    # narrow integer
    _, out = @grab_output @on_device oneAPI.@printf("Testing %d %d...\n", Int32(1), Int32(2))
    @test out == "Testing 1 2...$endline"

    # wide integer
    _, out = @grab_output if Sys.iswindows()
        @on_device oneAPI.@printf("Testing %lld %lld...\n", Int64(1), Int64(2))
    else
        @on_device oneAPI.@printf("Testing %ld %ld...\n", Int64(1), Int64(2))
    end
    @test out == "Testing 1 2...$endline"

    _, out = @grab_output @on_device begin
        oneAPI.@printf("foo")
        oneAPI.@printf("bar\n")
    end
    @test out == "foobar$endline"

    # c argument promotions
    function kernel(A)
        oneAPI.@printf("%f %f\n", A[1], A[1])
        return
    end
    x = oneArray(ones(2, 2))
    _, out = @grab_output begin
        @oneapi kernel(x)
        synchronize()
    end
    @test out == "1.000000 1.000000$endline"
end

@testset "@print" begin
    # basic @print/@println

    _, out = @grab_output @on_device oneAPI.@print("Hello, World\n")
    @test out == "Hello, World$endline"

    _, out = @grab_output @on_device oneAPI.@println("Hello, World")
    @test out == "Hello, World$endline"

    # argument interpolation (by the macro, so can use literals)

    _, out = @grab_output @on_device oneAPI.@print("foobar")
    @test out == "foobar"

    _, out = @grab_output @on_device oneAPI.@print(:foobar)
    @test out == "foobar"

    _, out = @grab_output @on_device oneAPI.@print("foo", "bar")
    @test out == "foobar"

    _, out = @grab_output @on_device oneAPI.@print("foobar ", 42)
    @test out == "foobar 42"

    _, out = @grab_output @on_device oneAPI.@print("foobar $(42)")
    @test out == "foobar 42"

    _, out = @grab_output @on_device oneAPI.@print("foobar $(4)", 2)
    @test out == "foobar 42"

    _, out = @grab_output @on_device oneAPI.@print("foobar ", 4, "$(2)")
    @test out == "foobar 42"

    _, out = @grab_output @on_device oneAPI.@print(42)
    @test out == "42"

    _, out = @grab_output @on_device oneAPI.@print(4, 2)
    @test out == "42"

    # bug: @println failed to invokce @print with endline in the case of interpolation
    _, out = @grab_output @on_device oneAPI.@println("foobar $(42)")
    @test out == "foobar 42$endline"

    # argument types

    # we're testing the generated functions now, so can't use literals
    function test_output(val, str)
        canary = rand(Int32) # if we mess up the main arg, this one will print wrong
        _, out = @grab_output @on_device oneAPI.@print(val, " (", canary, ")")
        @test out == "$(str) ($(Int(canary)))"
    end

    for typ in (Int16, Int32, Int64, UInt16, UInt32, UInt64)
        test_output(typ(42), "42")
    end

    for typ in (Float32, Float64)
        test_output(typ(42), "42.000000")
    end

    test_output(Cchar('c'), "c")

    for typ in (Ptr{Cvoid}, Ptr{Int})
        ptr = convert(typ, Int(0x12345))
        test_output(ptr, Sys.iswindows() ? "0000000000012345" : "0x12345")
    end

    test_output(true, "1")
    test_output(false, "0")

    # escaping

    kernel1(val) = (oneAPI.@print(val); nothing)
    _, out = @grab_output @on_device kernel1(42)
    @test out == "42"

    kernel2(val) = (oneAPI.@println(val); nothing)
    _, out = @grab_output @on_device kernel2(42)
    @test out == "42$endline"
end

@testset "@show" begin
    function kernel()
        seven_i32 = Int32(7)
        three_f64 = Float64(3)
        oneAPI.@show seven_i32
        oneAPI.@show three_f64
        oneAPI.@show 1f0 + 4f0
        return nothing
    end

    _, out = @grab_output @on_device kernel()
    @test out == "seven_i32 = 7$(endline)three_f64 = 3.000000$(endline)1.0f0 + 4.0f0 = 5.000000$(endline)"
end

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

# a composite type to test for more complex element types
@eval struct RGB{T}
    r::T
    g::T
    b::T
end

@testset "local memory" begin

n = 256

@testset "constructors" begin
    # static
    @on_device oneLocalArray(Float32, 1)
    @on_device oneLocalArray(Float32, (1,2))
    @on_device oneLocalArray(Tuple{Float32, Float32}, 1)
    @on_device oneLocalArray(Tuple{Float32, Float32}, (1,2))
    @on_device oneLocalArray(Tuple{RGB{Float32}, UInt32}, 1)
    @on_device oneLocalArray(Tuple{RGB{Float32}, UInt32}, (1,2))
end

@testset "static" begin

@testset "statically typed" begin
    function kernel(d, n)
        t = get_local_id(0)
        tr = n-t+1

        s = oneLocalArray(Float32, 1024)
        s2 = oneLocalArray(Float32, 1024)  # catch aliasing

        s[t] = d[t]
        s2[t] = 2*d[t]
        barrier()
        d[t] = s[tr]

        return
    end

    a = rand(Float32, n)
    d_a = oneArray(a)

    @oneapi items=n kernel(d_a, n)
    @test reverse(a) == Array(d_a)
end

@testset "parametrically typed" begin
    @testset for typ in [Int32, Int64, Float32, Float64]
        function kernel(d::oneDeviceArray{T}, n) where {T}
            t = get_local_id(0)
            tr = n-t+1

            s = oneLocalArray(T, 1024)
            s2 = oneLocalArray(T, 1024)  # catch aliasing

            s[t] = d[t]
            s2[t] = d[t]
            barrier()
            d[t] = s[tr]

            return
        end

        a = rand(typ, n)
        d_a = oneArray(a)

        @oneapi items=n kernel(d_a, n)
        @test reverse(a) == Array(d_a)
    end
end

end

end

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

@testset "atomics (low level)" begin

@testset "atomic_add($T)" for T in [Int32, UInt32]
    a = oneArray([zero(T)])

    function kernel(a, b)
        oneAPI.atomic_add!(pointer(a), b)
        return
    end

    @oneapi items=256 kernel(a, one(T))
    @test Array(a)[1] == T(256)
end

@testset "atomic_sub($T)" for T in [Int32, UInt32]
    a = oneArray([T(256)])

    function kernel(a, b)
        oneAPI.atomic_sub!(pointer(a), b)
        return
    end

    @oneapi items=256 kernel(a, one(T))
    @test Array(a)[1] == T(0)
end

@testset "atomic_inc($T)" for T in [Int32, UInt32]
    a = oneArray([zero(T)])

    function kernel(a)
        oneAPI.atomic_inc!(pointer(a))
        return
    end

    @oneapi items=256 kernel(a)
    @test Array(a)[1] == T(256)
end

@testset "atomic_dec($T)" for T in [Int32, UInt32]
    a = oneArray([T(256)])

    function kernel(a)
        oneAPI.atomic_dec!(pointer(a))
        return
    end

    @oneapi items=256 kernel(a)
    @test Array(a)[1] == T(0)
end

@testset "atomic_min($T)" for T in [Int32, UInt32]
    a = oneArray([T(256)])

    function kernel(a, T)
        i = get_global_id()
        oneAPI.atomic_min!(pointer(a), i%T)
        return
    end

    @oneapi items=256 kernel(a, T)
    @test Array(a)[1] == one(T)
end

@testset "atomic_max($T)" for T in [Int32, UInt32]
    a = oneArray([zero(T)])

    function kernel(a, T)
        i = get_global_id()
        oneAPI.atomic_max!(pointer(a), i%T)
        return
    end

    @oneapi items=256 kernel(a, T)
    @test Array(a)[1] == T(256)
end

@testset "atomic_and($T)" for T in [Int32, UInt32]
    a = oneArray([T(1023)])

    function kernel(a, T)
        i = get_global_id() - 1
        k = 1
        for i = 1:i
            k *= 2
        end
        b = 1023 - k  # 1023 - 2^i
        oneAPI.atomic_and!(pointer(a), T(b))
        return
    end

    @oneapi items=10 kernel(a, T)
    @test Array(a)[1] == zero(T)
end

@testset "atomic_or($T)" for T in [Int32, UInt32]
    a = oneArray([zero(T)])

    function kernel(a, T)
        i = get_global_id()
        b = 1  # 2^(i-1)
        for i = 1:i
            b *= 2
        end
        b /= 2
        oneAPI.atomic_or!(pointer(a), T(b))
        return
    end

    @oneapi items=10 kernel(a, T)
    @test Array(a)[1] == T(1023)
end

@testset "atomic_xor($T)" for T in [Int32, UInt32]
    a = oneArray([T(1023)])

    function kernel(a, T)
        i = get_global_id()
        b = 1  # 2^(i-1)
        for i = 1:i
            b *= 2
        end
        b /= 2
        oneAPI.atomic_xor!(pointer(a), T(b))
        return
    end

    @oneapi items=10 kernel(a, T)
    @test Array(a)[1] == zero(T)
end

@testset "atomic_xchg($T)" for T in [Int32, UInt32, Float32]
    a = oneArray([zero(T)])

    function kernel(a, b)
        oneAPI.atomic_xchg!(pointer(a), b)
        return
    end

    @oneapi items=256 kernel(a, one(T))
    @test Array(a)[1] == one(T)
end

end

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

@testset "atomics (high-level)" begin

@testset "add" begin
    @testset for T in [Int32, UInt32, Float32]
        a = oneArray([zero(T)])

        function kernel(T, a)
            oneAPI.@atomic a[1] = a[1] + 1
            oneAPI.@atomic a[1] += 1
            return
        end

        @oneapi items=256 kernel(T, a)
        @test Array(a)[1] == 512
    end
end

@testset "sub" begin
    @testset for T in [Int32, UInt32, Float32]
        a = oneArray(T[1024])

        function kernel(T, a)
            oneAPI.@atomic a[1] = a[1] - 1
            oneAPI.@atomic a[1] -= 1
            return
        end

        @oneapi items=256 kernel(T, a)
        @test Array(a)[1] == 512
    end
end

@testset "and" begin
    @testset for T in [Int32, UInt32]
        a = oneArray([~zero(T), ~zero(T)])

        function kernel(T, a)
            i = get_local_id()
            mask = ~(T(1) << (i-1))
            oneAPI.@atomic a[1] = a[1] & mask
            oneAPI.@atomic a[2] &= mask
            return
        end

        @oneapi items=8*sizeof(T) kernel(T, a)
        @test Array(a)[1] == zero(T)
        @test Array(a)[2] == zero(T)
    end
end

@testset "or" begin
    @testset for T in [Int32, UInt32]
        a = oneArray([zero(T), zero(T)])

        function kernel(T, a)
            i = get_local_id()
            mask = T(1) << (i-1)
            oneAPI.@atomic a[1] = a[1] | mask
            oneAPI.@atomic a[2] |= mask
            return
        end

        @oneapi items=8*sizeof(T) kernel(T, a)
        @test Array(a)[1] == ~zero(T)
        @test Array(a)[2] == ~zero(T)
    end
end

@testset "xor" begin
    @testset for T in [Int32, UInt32]
        a = oneArray([zero(T), zero(T)])

        function kernel(T, a)
            i = get_local_id()
            mask = T(1) << ((i-1)%(8*sizeof(T)))
            oneAPI.@atomic a[1] = a[1] ⊻ mask
            oneAPI.@atomic a[2] ⊻= mask
            return
        end

        nb = 4
        @oneapi items=(8*sizeof(T)+nb) kernel(T, a)
        @test Array(a)[1] == ~zero(T) & ~((one(T) << nb) - one(T))
        @test Array(a)[2] == ~zero(T) & ~((one(T) << nb) - one(T))
    end
end

@testset "max" begin
    @testset for T in [Int32, UInt32, Float32]
        a = oneArray([zero(T)])

        function kernel(T, a)
            i = get_local_id()
            oneAPI.@atomic a[1] = max(a[1], i)
            return
        end

        @oneapi items=32 kernel(T, a)
        @test Array(a)[1] == 32
    end
end

@testset "min" begin
    @testset for T in [Int32, UInt32, Float32]
        a = oneArray([typemax(T)])

        function kernel(T, a)
            i = get_local_id()
            oneAPI.@atomic a[1] = min(a[1], i)
            return
        end

        @oneapi items=32 kernel(T, a)
        @test Array(a)[1] == 1
    end
end

@testset "mul" begin
    @testset for T in [Int32, UInt32, Float32]
        a = oneArray(T[1])

        function kernel(T, a)
            oneAPI.@atomic a[1] = a[1] * 2
            oneAPI.@atomic a[1] *= 2
            return
        end

        @oneapi items=8 kernel(T, a)
        @test Array(a)[1] == 65536
    end
end

@testset "div" begin
    @testset for T in [Int32, UInt32, Float32]
        a = oneArray(T[65536])

        function kernel(T, a)
            oneAPI.@atomic a[1] = a[1] / 2
            oneAPI.@atomic a[1] /= 2
            return
        end

        @oneapi items=8 kernel(T, a)
        @test Array(a)[1] == 1
    end
end

@testset "macro" begin
    using oneAPI: AtomicError

    @test_throws_macro AtomicError("right-hand side of an @atomic assignment should be a call") @macroexpand begin
        oneAPI.@atomic a[1] = 1
    end
    @test_throws_macro AtomicError("right-hand side of an @atomic assignment should be a call") @macroexpand begin
        oneAPI.@atomic a[1] = b ? 1 : 2
    end

    @test_throws_macro AtomicError("right-hand side of a non-inplace @atomic assignment should reference the left-hand side") @macroexpand begin
        oneAPI.@atomic a[1] = a[2] + 1
    end

    @test_throws_macro AtomicError("unknown @atomic expression") @macroexpand begin
        oneAPI.@atomic wat(a[1])
    end

    @test_throws_macro AtomicError("@atomic should be applied to an array reference expression") @macroexpand begin
        oneAPI.@atomic a = a + 1
    end
end

end


text:
metre = 1e2
meter = metre
cm = metre / 1e2
mm = metre / 1e3
um = metre / 1e6
nm = metre / 1e9

second = 1e3
ms = second / 1e3
Hz = 1 / second

voltage = 1e3
mV = voltage / 1e3

ampere = 1e6
mA = ampere / 1e3
uA = ampere / 1e6
nA = ampere / 1e9

farad=1e6
uF = farad/1e6
ufarad = uF

siemens = 1e3
mS = siemens / 1e3
msiemens = mS
nS = siemens / 1e9
nsiemens = nS

max_stars_repo_name: russelljjarvis/juliaSNN

max_issues_repo_name: russelljjarvis/juliaSNN

max_forks_repo_name: russelljjarvis/juliaSNN


text:
type Bag{T}
	v::Vector{T}
	size::Int
	Bag(maxsize::Int) = new(Array(T, maxsize), 0)
end

Base.length(b::Bag) = b.size
Base.isempty(b::Bag) = b.size == 0
Base.start(b::Bag) = 1
Base.next(b::Bag, state) = b[state], state + 1
Base.done(b::Bag, state) = state > length(b)

function Base.getindex(b::Bag, i::Int)
	@assert 1 <= i <= b.size
	b.v[i]
end

function Base.push!{T}(b::Bag{T}, item::T)
	@assert b.size < size(b.v, 1)
	b.size += 1
	b.v[b.size] = item
	b
end

function Base.deleteat!(b::Bag, i::Int)
	@assert 1 <= i <= b.size
	item = b.v[i]
	b.v[i] = b.v[b.size]
	b.size -= 1
	b
end

max_stars_repo_name: yurivish/Julia-Playground

max_issues_repo_name: yurivish/Julia-Playground

max_forks_repo_name: yurivish/Julia-Playground


text:
isdefined(Base, :__precompile__) && __precompile__()

module SHA

using Compat

# Export convenience functions, context types, update!() and digest!() functions
export sha1, SHA1_CTX, update!, digest!
export sha224, sha256, sha384, sha512
export sha2_224, sha2_256, sha2_384, sha2_512
export sha3_224, sha3_256, sha3_384, sha3_512
export SHA224_CTX, SHA256_CTX, SHA384_CTX, SHA512_CTX
export SHA2_224_CTX, SHA2_256_CTX, SHA2_384_CTX, SHA2_512_CTX
export SHA3_224_CTX, SHA3_256_CTX, SHA3_384_CTX, SHA3_512_CTX

include("constants.jl")
include("types.jl")
include("base_functions.jl")
include("sha1.jl")
include("sha2.jl")
include("sha3.jl")
include("common.jl")

# Create data types and convenience functions for each hash implemented
for (f, ctx) in [(:sha1, :SHA1_CTX),
                 (:sha224, :SHA224_CTX),
                 (:sha256, :SHA256_CTX),
                 (:sha384, :SHA384_CTX),
                 (:sha512, :SHA512_CTX),
                 (:sha2_224, :SHA2_224_CTX),
                 (:sha2_256, :SHA2_256_CTX),
                 (:sha2_384, :SHA2_384_CTX),
                 (:sha2_512, :SHA2_512_CTX),
                 (:sha3_224, :SHA3_224_CTX),
                 (:sha3_256, :SHA3_256_CTX),
                 (:sha3_384, :SHA3_384_CTX),
                 (:sha3_512, :SHA3_512_CTX),]
    @eval begin
        # Our basic function is to process arrays of bytes
        function $f(data::Array{UInt8,1})
            ctx = $ctx()
            update!(ctx, data)
            return digest!(ctx)
        end

        # AbstractStrings are a pretty handy thing to be able to crunch through
        $f(str::AbstractString) = $f(convert(Array{UInt8,1}, str))

        # Convenience function for IO devices, allows for things like:
        # open("test.txt") do f
        #     sha256(f)
        # done
        $f(io::IO) = $f(read(io))
    end
end

end #module SHA


text:
import Base.copy

Base.@kwdef mutable struct NestedInterval{T₁ <: AbstractRange} <: AbstractModel
    start::Float64
    stop::Float64
    interval::T₁
end

function get_start(ni::NestedInterval)
    ni.start
end

function set_start!(ni::NestedInterval, val::Number)
    ni.start = val > first(ni.interval) ? val : first(ni.interval)
end

function get_stop(ni::NestedInterval)
    ni.stop
end

function set_stop!(ni::NestedInterval, val::Number)
    ni.stop = val < last(ni.interval) ? val : last(ni.interval)
end

function get_interval(ni::NestedInterval)
    ni.interval
end

function set_interval!(ni::NestedInterval, interval::AbstractRange)
    ni.interval = interval
end

function copy(ni::NestedInterval)
    NestedInterval(get_start(ni), get_stop(ni), get_interval(ni))
end


text:
using Pkg
Pkg.activate(pwd())

# # Arithmetic operators

1 + 2
2*3
4/3

#+

x = 1;
y = 3;

(x + 2)/(y - 1) - 4*(x - 2)^2

#+

2(3 + 4) # equivalent to 2*(3 + 4)

# ### Exercise:
# Determine the value and type of $y$ given by the following expression
# $$
# y = \frac{(x + 2)^2 - 4}{(x - 2)^{p - 2}},
# $$
# where $x = 4$ and $p = 5$.
#
# ---
# ### Solution:

# ---
#
# ## Promotion system

x = 1.0 # Float64
y = 2 # Int64

xp, yp = promote(x, y)

typeof(xp)
typeof(yp)

#+

promote(1, 2f0, true, 4.5, Int32(1))
promote_type(Float64, Int64, Bool, Int32)

#+

x = 1 # Int64
y = 2f0 # Float32
z = x + y

typeof(z)

# ### Exercise:
# All of these values represent number ``1``. Determine the smallest type which can
# represent them.

x = 1
y = 1f0
z = true
w = Int32(1)

# ---
# ### Solution:

# ---
#
# ## Updating operators

x = 1
x += 3 # x = x + 3
x *= 4 # x = x * 4
x /= 2 # x = x / 2
x \= 16 # x = x \ 16 = 16 / x

# ### Exercise:
# Compute the value of $y$ given by the following expression
# $$
# y = \frac{(x + 4)^{\frac{3}{2}}}{(x + 1)^{p - 1}},
# $$
# where $x = 5$ and $p = 3$. Then multiply the result by $8$, add $3$, divide by $3$, and
# subtract $1$. What are all the intermediate results and the final result?
#
# ---
# ### Solution:

# ---
#
# ## Numeric comparison

1 == 1
1 == 1.0

#+

-1 <= 1
-1 ≥ 1

#+

3 > 2 > 1
3 > 2 & 2 > 1

#+

1 < 2 <= 2 < 3 == 3 > 2 >= 1 == 1 < 3 != 5

#+

NaN == NaN
NaN != NaN
NaN < NaN

#+

isequal(NaN, NaN)
!isequal(NaN, NaN)

#+

!true
!false

# ## Rounding functions

x = 3141.5926

round(x)
floor(x)
ceil(x)

#+

round(Int64, x)
floor(Int32, x)
ceil(Int16, x)

#+

round(x; digits = 3)
round(x; sigdigits = 3)

# ### Exercise:
# Use rounding functions to solve the following tasks:
# - Round `1252.1518` to the nearest larger integer and convert the resulting value to
#   `Int64`.
# - Round `1252.1518` to the nearest smaller integer and convert the resulting value to
#   `Int16`.
# - Round `1252.1518` to `2` digits after the decimal point.
# - Round `1252.1518` to `3` significant digits.
#
# ---
# ### Solution:

# ---
#
# ## Numerical Conversions

convert(Float32, 1.234)
Float32(1.234)

#+

convert(Int64, 1.0)
Int64(1.0)

#+

convert(Int64, 1.234)
Int64(1.234)

# ### Exercise:
# Use the proper numeric conversion to get the correct result (not approximate) of summing
# the following two numbers
#
# ```julia
# x = 1//3
# y = 0.5
# ```
#
# **Hint:** rational numbers can be summed without approximation.
#
# ---
# ### Solution: