text:
export pseudo_basis, basis_pmatrix

add_verbose_scope(:NfRelOrd)

################################################################################
#
#  Is maximal order known
#
################################################################################

function ismaximal_order_known(K::T) where T <: Union{NfRel, NfRelNS}
  try
    _get_maximal_order_of_nf_rel(K)
    return true
  catch e
    if !isa(e, AccessorNotSetError)
      rethrow(e)
    end
    return false
  end
end

################################################################################
#
#  Basic field access
#
################################################################################

parent(O::NfRelOrd) = O.parent

base_ring(O::NfRelOrd) = order(pseudo_basis(O, copy = false)[1][2])

elem_type(::NfRelOrd{T, S, U}) where {T, S, U} = NfRelOrdElem{T, U}

elem_type(::Type{NfRelOrd{T, S, U}}) where {T, S, U} = NfRelOrdElem{T, U}

ideal_type(::NfRelOrd{T, S, U}) where {T, S, U} = NfRelOrdIdl{T, S, U}

ideal_type(::Type{NfRelOrd{T, S, U}}) where {T, S, U} = NfRelOrdIdl{T, S, U}

fractional_ideal_type(::NfRelOrd{T, S, U}) where {T, S, U} = NfRelOrdFracIdl{T, S, U}

fractional_ideal_type(::Type{NfRelOrd{T, S, U}}) where {T, S, U} = NfRelOrdFracIdl{T, S, U}

################################################################################
#
#  Modulus
#
################################################################################

# Return a small OK-ideal I such that I * OK^n \subseteq OL (with respect to
# the basis of OK)
# This is required for the modular computations with ideals.
function _modulus(O::NfRelOrd)
  D = get_attribute(O, :modulus)
  if D isa Nothing
    D = reduce(lcm, coefficient_ideals(basis_pmatrix(O)))
    # Let D = N/d
    # In particular (N/d) * OK^n \subseteq OL
    # N*OK \subseteq N/d * OK \subseteq OL
    # So we take N
    M = numerator(simplify(D))
    set_attribute!(O, :modulus => M)
    return M
  else
    return D::ideal_type(base_ring(O))
  end
end

################################################################################
#
#  "Assure" functions for fields
#
################################################################################

function assure_has_basis_pmatrix(O::NfRelOrd{T, S, U}) where {T, S, U}
  if isdefined(O, :basis_pmatrix)
    return nothing
  end
  if !isdefined(O, :pseudo_basis)
    error("No pseudo_basis and no basis_pmatrix defined.")
  end
  pb = pseudo_basis(O, copy = false)
  L = nf(O)
  M = zero_matrix(base_field(L), degree(O), degree(O))
  C = Vector{S}()
  for i = 1:degree(O)
    elem_to_mat_row!(M, i, pb[i][1])
    push!(C, deepcopy(pb[i][2]))
  end
  O.basis_pmatrix = PseudoMatrix(M, C)
  return nothing
end

function assure_has_pseudo_basis(O::NfRelOrd{T, S, U}) where {T, S, U}
  if isdefined(O, :pseudo_basis)
    return nothing
  end
  if !isdefined(O, :basis_pmatrix)
    error("No pseudo_basis and no basis_pmatrix defined.")
  end
  P = basis_pmatrix(O, copy = false)
  L = nf(O)
  pseudo_basis = Vector{Tuple{U, S}}()
  for i = 1:degree(O)
    a = elem_from_mat_row(L, P.matrix, i)
    push!(pseudo_basis, (a, deepcopy(P.coeffs[i])))
  end
  O.pseudo_basis = pseudo_basis
  return nothing
end

function assure_has_basis_nf(O::NfRelOrd{T, S, U}) where {T, S, U}
  if isdefined(O, :basis_nf)
    return nothing
  end
  L = nf(O)
  pb = pseudo_basis(O)
  basis_nf = Vector{U}()
  for i = 1:degree(O)
    push!(basis_nf, pb[i][1])
  end
  O.basis_nf = basis_nf
  return nothing
end

function assure_has_basis_matrix(O::NfRelOrd)
  if isdefined(O, :basis_matrix)
    return nothing
  end
  O.basis_matrix = basis_pmatrix(O).matrix
  return nothing
end

function assure_has_basis_mat_inv(O::NfRelOrd)
  if isdefined(O, :basis_mat_inv)
    return nothing
  end
  O.basis_mat_inv = inv(basis_matrix(O, copy = false))
  return nothing
end

function assure_has_inv_coeff_ideals(O::NfRelOrd)
  if isdefined(O, :inv_coeff_ideals)
    return nothing
  end
  pb = pseudo_basis(O, copy = false)
  O.inv_coeff_ideals = [ inv(pb[i][2]) for i in 1:degree(O) ]
  return nothing
end

################################################################################
#
#  Pseudo basis / basis pseudo-matrix
#
################################################################################

@doc Markdown.doc"""
      pseudo_basis(O::NfRelOrd{T, S}) -> Vector{Tuple{NumFieldElem{T}, S}}

Returns the pseudo-basis of $\mathcal O$.
"""
function pseudo_basis(O::NfRelOrd; copy::Bool = true)
  assure_has_pseudo_basis(O)
  if copy
    return deepcopy(O.pseudo_basis)
  else
    return O.pseudo_basis
  end
end

@doc Markdown.doc"""
      basis_pmatrix(O::NfRelOrd) -> PMat

Returns the basis pseudo-matrix of $\mathcal O$ with respect to the power basis
of the ambient number field.
"""
function basis_pmatrix(O::NfRelOrd; copy::Bool = true)
  assure_has_basis_pmatrix(O)
  if copy
    return deepcopy(O.basis_pmatrix)
  else
    return O.basis_pmatrix
  end
end

@doc Markdown.doc"""
      inv_coeff_ideals(O::NfRelOrd{T, S}) -> Vector{S}

Returns the inverses of the coefficient ideals of the pseudo basis of $O$.
"""
function inv_coeff_ideals(O::NfRelOrd; copy::Bool = true)
  assure_has_inv_coeff_ideals(O)
  if copy
    return deepcopy(O.inv_coeff_ideals)
  else
    return O.inv_coeff_ideals
  end
end

################################################################################
#
#  Basis / (inverse) basis matrix
#
################################################################################

@doc Markdown.doc"""
      basis_nf(O::NfRelOrd) -> Vector{NumFieldElem}

Returns the elements of the pseudo-basis of $\mathcal O$ as elements of the
ambient number field.
"""
function basis_nf(O::NfRelOrd{S, T, U}; copy::Bool = true) where {S, T, U}
  assure_has_basis_nf(O)
  if copy
    return deepcopy(O.basis_nf)::Vector{U}
  else
    return O.basis_nf::Vector{U}
  end
end

function basis_matrix(O::NfRelOrd; copy::Bool = true)
  assure_has_basis_matrix(O)
  if copy
    return deepcopy(O.basis_matrix)
  else
    return O.basis_matrix
  end
end

function basis_mat_inv(O::NfRelOrd; copy::Bool = true)
  assure_has_basis_mat_inv(O)
  if copy
    return deepcopy(O.basis_mat_inv)
  else
    return O.basis_mat_inv
  end
end

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

function show(io::IO, S::NfRelOrdSet)
  print(io, "Set of orders of the number field ")
  print(io, S.nf)
end

function show(io::IO, O::NfRelOrd)
  compact = get(io, :compact, false)
  if compact
    if ismaximal_known_and_maximal(O)
      print(io, "Relative maximal order with pseudo-basis ")
    else
      print(io, "Relative order with pseudo-basis ")
    end
    pb = pseudo_basis(O, copy = false)
    for i = 1:degree(O)
      print(io, "($(pb[i][1])) * ")
      show(IOContext(io, :compact => true), pb[i][2])
      if i != degree(O)
        print(io, ", ")
      end
    end
  else
    if ismaximal_known_and_maximal(O)
      print(io, "Relative maximal order of ")
    else
      print(io, "Relative order of ")
    end
    println(io, nf(O))
    print(io, "with pseudo-basis ")
    pb = pseudo_basis(O, copy = false)
    for i = 1:degree(O)
      print(io, "\n(")
      print(io, pb[i][1])
      print(io, ", ")
      show(IOContext(io, :compact => true), pb[i][2])
      print(io, ")")
    end
  end
end

################################################################################
#
#  Discriminant
#
################################################################################

function assure_has_discriminant(O::NfRelOrd{nf_elem, NfOrdFracIdl, NfRelElem{nf_elem}})
  if isdefined(O, :disc_abs)
    return nothing
  end
  if isequation_order(O)
    K = nf(O)
    F = base_field(K)
    OF = maximal_order(F)
    d = OF(discriminant(K.pol))
    O.disc_abs = ideal(OF, d)
    return nothing
  end
  d = det(trace_matrix(O, copy = false))
  pb = pseudo_basis(O, copy = false)
  a = pb[1][2]^2
  for i = 2:degree(O)
    a *= pb[i][2]^2
  end
  disc = d*a
  simplify(disc)
  O.disc_abs = numerator(disc)
  return nothing
end

function assure_has_discriminant(O::NfRelOrd{nf_elem, NfOrdFracIdl, NfRelNSElem{nf_elem}})
  if isdefined(O, :disc_abs)
    return nothing
  end
  if isequation_order(O)
    K = nf(O)
    F = base_field(K)
    OF = maximal_order(F)
    pols = K.pol
    Fx, _ = PolynomialRing(F, "x", cached = false)
    pol = to_univariate(Fx, pols[1])
    d = OF(discriminant(pol))^(div(degree(K), degree(pol)))
    for i = 2:length(pols)
      pol = to_univariate(Fx, pols[i])
      d *= OF(discriminant(pol))^(div(degree(K), degree(pol)))
    end
    O.disc_abs = ideal(OF, d)
    return nothing
  end
  d = det(trace_matrix(O, copy = false))
  pb = pseudo_basis(O, copy = false)
  a = pb[1][2]^2
  for i = 2:degree(O)
    a *= pb[i][2]^2
  end
  disc = d*a
  simplify(disc)
  O.disc_abs = numerator(disc)
  return nothing
end

function assure_has_discriminant(O::NfRelOrd{T, S, U}) where {T, S, U}
  if isdefined(O, :disc_rel)
    return nothing
  end
  d = det(trace_matrix(O, copy = false))
  pb = pseudo_basis(O, copy = false)
  a = pb[1][2]^2
  for i = 2:degree(O)
    a *= pb[i][2]^2
  end
  disc = d*a
  simplify(disc)
  O.disc_rel = numerator(disc)
  return nothing
end

function discriminant(O::NfRelOrd{nf_elem, S}) where S
  assure_has_discriminant(O)
  return deepcopy(O.disc_abs)
end

function discriminant(O::NfRelOrd{T, S}) where {T <: NumFieldElem{U} where U, S}
  assure_has_discriminant(O)
  return deepcopy(O.disc_rel)::ideal_type(order_type(base_field(nf(O))))
end

function absolute_discriminant(O::NfRelOrd{T, S}) where {T, S}
  d = norm(discriminant(O))
  d1 = absolute_discriminant(base_ring(O))
  di = abs(d1^(degree(O))*d)

  # Compute the sign using Brill's theorem
  _, s = signature(nf(O))
  if mod(s, 2) == 0
    return di
  else
    return -di
  end
end

function absolute_discriminant(O::NfAbsOrd)
  return discriminant(O)
end

################################################################################
#
#  Codifferent
#
################################################################################

function codifferent(O::NfRelOrd)
  T = trace_matrix(O, copy = false)
  R = base_ring(O)
  K = nf(R)
  pb = pseudo_basis(O)
  pm = pseudo_matrix(inv(change_base_ring(K, T)),
                     fractional_ideal_type(base_ring(O))[inv(y) for (_, y) in pb])
  return fractional_ideal(O, pm)
end

function different(O::NfRelOrd)
  if !ismaximal_known_and_maximal(O)
    throw(error("Order not known to be maximal"))
  end
  return ideal(O, basis_pmatrix(inv(codifferent(O))))
end

################################################################################
#
#  Degree
#
################################################################################

degree(O::NfRelOrd) = degree(nf(O))

################################################################################
#
#  Deepcopy
#
################################################################################

function Base.deepcopy_internal(O::NfRelOrd{T, S, U}, dict::IdDict) where {T, S, U}
  z = NfRelOrd{T, S, U}(O.nf)
  for x in fieldnames(typeof(O))
    if x != :nf && x != :parent && isdefined(O, x)
      setfield!(z, x, Base.deepcopy_internal(getfield(O, x), dict))
    end
  end
  z.nf = O.nf
  z.parent = O.parent
  return z
end

################################################################################
#
#  Inclusion of number field elements
#
################################################################################

function _check_elem_in_order(a::U, O::NfRelOrd{T, S, U}, short::Type{Val{V}} = Val{false}) where {T, S, U, V}
  b_pmat = basis_pmatrix(O, copy = false)
  t = zero_matrix(base_field(nf(O)), 1, degree(O))
  elem_to_mat_row!(t, 1, a)
  t = t*basis_mat_inv(O, copy = false)
  if short == Val{true}
    for i = 1:degree(O)
      if !(t[1, i] in b_pmat.coeffs[i])
        return false
      end
    end
    return true
  else
    for i = 1:degree(O)
      if !(t[1, i] in b_pmat.coeffs[i])
        return false, Vector{T}()
      end
    end
    v = Vector{T}(undef, degree(O))
    for i in 1:degree(O)
      v[i] = deepcopy(t[1, i])
    end
    return true, v
  end
end

function in(a::U, O::NfRelOrd{T, S, U}) where {T, S, U}
  return _check_elem_in_order(a, O, Val{true})
end

################################################################################
#
#  Construction
#
################################################################################

function Order(L::NfRel{nf_elem}, M::Generic.Mat{nf_elem})
  return NfRelOrd{nf_elem, NfOrdFracIdl, NfRelElem{nf_elem}}(L, deepcopy(M))
end

function Order(L::NfRelNS{nf_elem}, M::Generic.Mat{nf_elem})
  return NfRelOrd{nf_elem, NfOrdFracIdl, NfRelNSElem{nf_elem}}(L, deepcopy(M))
end

function Order(L::NumField{S}, M::Generic.Mat{S}) where S <: NumFieldElem{T} where T
  return NfRelOrd{S, NfRelOrdFracIdl{T}, elem_type(L)}(L, deepcopy(M))
end

#=
@doc Markdown.doc"""
      Order(K::NumField, M::PMat) -> NfRelOrd

Returns the order which has basis pseudo-matrix $M$ with respect to the power basis
of $K$.
"""
=#
function Order(L::NumField{T}, M::PMat{T, S}) where {T, S}
  # checks
  return NfRelOrd{T, S, elem_type(L)}(L, deepcopy(M))
end

function EquationOrder(L::NumField)
  M = identity_matrix(base_field(L), degree(L))
  PM = PseudoMatrix(M)
  O = Order(L, PM)
  O.basis_mat_inv = M
  O.isequation_order = true
  return O
end

function EquationOrder(L::NfRel{nf_elem})
  M = identity_matrix(base_field(L), degree(L))
  PM = PseudoMatrix(M)
  O = Order(L, PM)
  O.basis_mat_inv = M
  O.isequation_order = true
  O.index = ideal(maximal_order(base_field(L)), 1)
  return O
end

equation_order(L::NumField) = EquationOrder(L)

function MaximalOrder(L::NumField)
  try
    O = _get_maximal_order_of_nf_rel(L)
    return O::order_type(L)
  catch e
    if !isa(e, AccessorNotSetError)
      rethrow(e)
    end
    O = MaximalOrder(EquationOrder(L))
    O.ismaximal = 1
    _set_maximal_order_of_nf_rel(L, O)
    return O::order_type(L)
  end
end

function maximal_order_via_absolute(L::NfRel)
  Labs, LabsToL = absolute_simple_field(L)
  Oabs = maximal_order(Labs)
  return relative_order(Oabs, LabsToL)
end

function maximal_order_via_absolute(m::NfToNfRel)
  Oabs = maximal_order(domain(m))
  return relative_order(Oabs, m)
end

function maximal_order_via_simple(L::NfRelNS)
  Ls, m = simple_extension(L)
  Os = maximal_order(Ls)
  return non_simple_order(Os, m)
end

function maximal_order_via_simple(m::NumFieldMor{<:NfRel, <:NfRelNS})
  Os = maximal_order(domain(m))
  return non_simple_order(Os, m)
end

function maximal_order_via_relative(K::AnticNumberField, m::NfToNfRel)
  try
    O = _get_maximal_order(K)
    return O
  catch e
    if !isa(e, AccessorNotSetError)
      rethrow(e)
    end
  end
  L = codomain(m)
  OL = maximal_order(L)
  B = absolute_basis(OL, L)
  OK = Order(K, [ m\b for b in B ], check = false, isbasis = true)
  OK.ismaximal = 1
  _set_maximal_order(K, OK)
  return OK
end

function absolute_basis(O::NfRelOrd{T, S, U}, K::NumField{T}) where {T, S, U}
  pb = pseudo_basis(O, copy = false)
  res = Vector{elem_type(K)}()
  for i = 1:degree(O)
    for b in absolute_basis(pb[i][2])
      push!(res, b*pb[i][1])
    end
  end
  return res
end

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

function ==(R::NfRelOrd, S::NfRelOrd)
  nf(R) != nf(S) && return false
  return basis_pmatrix(R, copy = false) == basis_pmatrix(S, copy = false)
end

################################################################################
#
#  Trace matrix
#
################################################################################

#=
@doc Markdown.doc"""
    trace_matrix(O::NfRelOrd{T, S}) -> Generic.Mat{T}

Returns the trace matrix of $\mathcal O$ with respect to $(\omega_1,\dotsc,\omega_d)$
where $(\mathfrak c_i, \omega_i)$ is the pseudo-basis of $\mathcal O$.
"""
=#
function trace_matrix(O::NfRelOrd; copy::Bool = true)
  if isdefined(O, :trace_mat)
    if copy
      return deepcopy(O.trace_mat)
    else
      return O.trace_mat
    end
  end
  L = nf(O)
  K = base_field(L)
  b = basis_nf(O, copy = false)
  d = degree(L)
  g = zero_matrix(K, d, d)
  for i = 1:d
    t = tr(b[i]*b[i])
    g[i, i] = t
    for j = (i + 1):d
      t = tr(b[i]*b[j])
      g[i, j] = t
      g[j, i] = t
    end
  end
  O.trace_mat = g
  if copy
    return deepcopy(O.trace_mat)
  else
    return O.trace_mat
  end
end

################################################################################
#
#  Dedekind criterion
#
################################################################################

function fq_nmod_poly_to_nf_elem_poly(R::Generic.PolyRing{nf_elem}, m::InverseMap, f::fq_nmod_poly)
  @assert codomain(m) == base_ring(R)
  @assert domain(m) == base_ring(parent(f))

  g = zero(R)
  for i = 0:degree(f)
    setcoeff!(g, i, m(coeff(f, i)))
  end
  return g
end

function fq_poly_to_nf_elem_poly(R::Generic.PolyRing{T}, m::InverseMap, f::fq_poly) where {T <: Union{nf_elem, NfRelElem}}
  return map_coefficients(m, f, parent = R)
end

# Algorithm IV.6. in "Berechnung relativer Ganzheitsbasen mit dem
# Round-2-Algorithmus" by C. Friedrichs.
function dedekind_test(O::NfRelOrd{U1, V, Z}, p::Union{NfAbsOrdIdl, NfRelOrdIdl}, compute_order::Type{Val{S}} = Val{true}) where {S, U1, V, Z <: NfRelElem}
  !isequation_order(O) && error("Order must be an equation order")

  L = nf(O)
  K = base_field(L)
  T = L.pol
  Kx = parent(T)
  OK = maximal_order(K)
  F, mF = ResidueField(OK, p)
  mmF = extend(mF, K)
  immF = pseudo_inv(mmF)
  Fy, y = PolynomialRing(F,"y", cached=false)

  Tmodp = map_coefficients(mmF, T, parent = Fy)
  fac = factor(Tmodp)
  g = Kx(1)
  for (t, e) in fac
    mul!(g, g, fq_poly_to_nf_elem_poly(Kx, immF, t))
  end
  gmodp = map_coefficients(mmF, g, parent = Fy)
  hmodp = divexact(Tmodp, gmodp)
  h = fq_poly_to_nf_elem_poly(Kx, immF, hmodp)
  a = anti_uniformizer(p)
  f = a*(g*h - T)
  fmodp = map_coefficients(mmF, f, parent = Fy)

  d = gcd(fmodp, gcd(gmodp, hmodp))

  if compute_order == Val{false}
    return isone(d)
  else
    if isone(d)
      return true, O
    end

    Umodp = divexact(Tmodp, d)
    U = fq_poly_to_nf_elem_poly(Kx, immF, Umodp)
    PM = PseudoMatrix(representation_matrix(a*U(gen(L))), [ K(1)*OK for i = 1:degree(O) ])
    PN = vcat(basis_pmatrix(O), PM)
    PN = sub(pseudo_hnf_full_rank_with_modulus(PN, p, :lowerleft), degree(O) + 1:2*degree(O), 1:degree(O))
    OO = typeof(O)(L, PN)
    OO.isequation_order = false
    return false, OO
  end
end

dedekind_ispmaximal(O::NfRelOrd, p::Union{NfAbsOrdIdl, NfRelOrdIdl}) = dedekind_test(O, p, Val{false})

dedekind_poverorder(O::NfRelOrd, p::Union{NfAbsOrdIdl, NfRelOrdIdl}) = dedekind_test(O, p)[2]

################################################################################
#
#  p-overorder
#
################################################################################

#=
@doc Markdown.doc"""
      poverorder(O::NfRelOrd, p::Union{NfOrdIdl, NfRelOrdIdl}) -> NfRelOrd

This function tries to find an order that is locally larger than $\mathcal O$
at the prime $p$.
"""
=#
function poverorder(O::NfRelOrd, p::Union{NfAbsOrdIdl, NfRelOrdIdl})
  if isequation_order(O) && issimple(O)
    @vprint :NfRelOrd 3 "Applying Dedekind criterion\n"
    return dedekind_poverorder(O, p)
  else
    @vprint :NfRelOrd 3 "Computing pradical\n"
    @vtime :NfRelOrd 4 Ip = pradical(O, p)
    @vprint :NfRelOrd 3 "Computing ring of multipliers\n"
    @vtime :NfRelOrd 4 Op = ring_of_multipliers(Ip)
    return Op
  end
end

function poverorder(O::NfRelOrd{S, T, NfRelElem{nf_elem}}, p::NfOrdIdl) where {S, T}
  if isequation_order(O)
    return overorder_polygons(O, p)
  end
  @vprint :NfRelOrd 3 "Computing pradical\n"
  @vtime :NfRelOrd 4 Ip = pradical(O, p)
  @vprint :NfRelOrd 3 "Computing ring of multipliers\n"
  @vtime :NfRelOrd 4 Op = ring_of_multipliers(Ip)
  return Op
end

################################################################################
#
#  p-maximal overorder
#
################################################################################
#=
@doc Markdown.doc"""
      pmaximal_overorder(O::NfRelOrd, p::Union{NfOrdIdl, NfRelOrdIdl}) -> NfRelOrd

This function finds a $p$-maximal order $R$ containing $\mathcal O$.
"""
=#
function pmaximal_overorder(O::NfRelOrd, p::Union{NfAbsOrdIdl, NfRelOrdIdl})
  d = discriminant(O)
  if valuation(d, p) < 2
    return O
  end
  OO = poverorder(O, p)
  dd = discriminant(OO)
  while d != dd
    if valuation(dd, p) < 2
      break
    end
    d = dd
    OO = poverorder(OO, p)
    dd = discriminant(OO)
  end
  return OO
end

function MaximalOrder(O::NfRelOrd)
  disc = discriminant(O)
  fac = factor(disc)
  OO = O
  for (p, e) in fac
    if e == 1
      continue
    end
    OO = sum_as_OK_modules(OO, pmaximal_overorder(O, p))
  end
  OO.ismaximal = 1
  return OO
end

function MaximalOrder(O::NfRelOrd{S, T, U}) where {S, T, U <: NonSimpleNumFieldElem}
  L = nf(O)
  K = base_field(L)
  Obase_K = maximal_order(K)
  fields = Vector{Tuple{NfRel{S}, morphism_type(NfRel{S}, NfRelNS{S})}}(undef, length(L.pol))
  for i = 1:length(L.pol)
    fields[i] = component(L, i)
  end
  #Now, bring the maximal order of every component in L
  discs = Vector{ideal_type(Obase_K)}(undef, ngens(L))
  B = Vector{Vector{Tuple{elem_type(L), fractional_ideal_type(Obase_K)}}}(undef, length(fields))
  for i = 1:length(fields)
    OK = maximal_order(fields[i][1])
    discs[i] = discriminant(OK)
    BOK = pseudo_basis(OK, copy = false)
    BK = Vector{Tuple{elem_type(L), fractional_ideal_type(Obase_K)}}(undef, degree(OK))
    for j = 1:length(BK)
      BK[j] = (fields[i][2](BOK[j][1]), BOK[j][2])
    end
    B[i] = BK
  end
  Bp = product_pseudobasis(B)
  MOstart = PseudoMatrix(basis_matrix(U[x[1] for x in Bp]), fractional_ideal_type(Obase_K)[x[2] for x in Bp])
  Ostart = Order(L, MOstart)
  lp = ideal_type(Obase_K)[]
  for i = 1:length(fields)
    for j = i+1:length(fields)
      push!(lp, gcd(discs[i], discs[j]))
    end
  end
  cp = coprime_base(lp)
  for I in cp
    lP = factor(I)
    for P in keys(lP)
      Ostart = pmaximal_overorder(Ostart, P)
    end
  end
  return Ostart
end

function product_pseudobasis(l::Vector{Vector{Tuple{T, S}}}) where {S, T <: Union{NfAbsOrdElem, NfRelOrdElem, NumFieldElem}}
  nelems = 1
  for i = 1:length(l)
    nelems *= length(l[i])
  end
  B = typeof(l[1])(undef, nelems)
  ind = length(l[1])
  for i = 1:ind
    B[i] = l[1][i]
  end
  for jj = 2:length(l)
    new_deg = length(l[jj])
    for i = 2:new_deg
      for j = 1:ind
        B[(i-1)* ind + j] = (B[j][1] * l[jj][i][1], B[j][2] * l[jj][i][2])
      end
    end
    ind *= new_deg
  end
  return B
end

################################################################################
#
#  Addition of orders
#
################################################################################
#TODO: this is wrong unless the indices are coprime
function +(a::NfRelOrd{T, S, U}, b::NfRelOrd{T, S, U}) where {T, S, U}
  # checks
  @assert nf(a) == nf(b)
  d = degree(a)
  aB = basis_pmatrix(a, copy = false)
  bB = basis_pmatrix(b, copy = false)
  M = vcat(aB, bB)
  PM = sub(pseudo_hnf(M, :lowerleft, true), d + 1:2*d, 1:d)
  return NfRelOrd{T, S, U}(nf(a), PM)
end

function +(a::NfRelOrd{T, S, U}, b::NfRelOrd{T, S, U}) where {T, S, U <: Union{NfRelElem{nf_elem}, NfRelNSElem{nf_elem}}}
  # checks
  @assert nf(a) == nf(b)
  K = base_field(nf(a))
  d = degree(a)
  aB = basis_pmatrix(a)
  bB = basis_pmatrix(b)
  M = vcat(aB, bB)
  z = _make_integral!(M)
  d1 = numerator(det(sub(M, 1:d, 1:d)))
  d2 = numerator(det(sub(M, d+1:2*d, 1:d)))
  m = d1 + d2
  M = sub(pseudo_hnf_mod(M, m, :lowerleft), d + 1:2*d, 1:d)
  for i in 1:nrows(M)
    M.coeffs[i] = M.coeffs[i]*inv(K(z))
    simplify(M.coeffs[i])
  end
  return NfRelOrd{T, S, U}(nf(a), M)
end

function sum_as_OK_modules(a::NfRelOrd{T, S, U}, b::NfRelOrd{T, S, U}) where {T, S, U <: NfRelElem{nf_elem}}
  if !isdefined(a, :index) || !isdefined(b, :index)
    return a+b
  end
  K = base_field(nf(a))
  d = degree(a)
  aB = basis_pmatrix(a)
  bB = basis_pmatrix(b)
  M = vcat(aB, bB)
  new_index = a.index * b.index
  M = sub(pseudo_hnf_full_rank_with_modulus(M, new_index, :lowerleft), d+1:2*d, 1:d)
  res = NfRelOrd{T, S, U}(nf(a), M)
  res.index = new_index
  res.disc_abs = simplify(gcd(discriminant(a), discriminant(b)))
  return res
end

function sum_as_OK_modules(a::NfRelOrd{T, S, U}, b::NfRelOrd{T, S, U}) where {T, S, U}
  return a+b
end

################################################################################
#
#  Absolute to relative
#
################################################################################

function relative_order(O::NfOrd, m::NfToNfRel)
  L = codomain(m)
  Labs = domain(m)
  @assert nf(O) == Labs
  K = base_field(L)
  OK = maximal_order(K)
  B = basis(O, Labs, copy = false)
  if ismaximal_known_and_maximal(O)
    E = EquationOrder(L)
    els = elem_type(L)[m(x) for x in B]
    return add_to_order(E, els)
  else
    return _order(elem_type(L)[mp(x) for x in B])
  end
end

################################################################################
#
#  Simple to non-simple
#
################################################################################

function non_simple_order(O::NfRelOrd, m::NumFieldMor{<:NfRel, <:NfRelNS})
  L = domain(m)
  L_ns = codomain(m)
  @assert nf(O) == L
  K = base_field(L)
  B = basis_nf(O, copy = false)
  d = degree(L)
  M = zero_matrix(K, d, d)
  for i = 1:d
    elem_to_mat_row!(M, i, m(L(B[i])))
  end
  PM = pseudo_hnf(PseudoMatrix(M, Hecke.basis_pmatrix(O).coeffs), :lowerleft, true)
  return NfRelOrd{typeof(PM.matrix[1, 1]), typeof(PM.coeffs[1]), elem_type(L_ns)}(L_ns, PM)
end

################################################################################
#
#  Denominator in an order
#
################################################################################

@doc Markdown.doc"""
    denominator(a::NumFieldElem, O::NfOrd) -> fmpz

Returns the smallest positive integer $k$ such that $k \cdot a$ is contained in
$\mathcal O$.
"""
function denominator(a::NumFieldElem, O::NfRelOrd)
  t = zero_matrix(base_field(nf(O)), 1, degree(O))
  elem_to_mat_row!(t, 1, a)
  t = t*basis_mat_inv(O, copy = false)
  d = fmpz(1)
  icv = inv_coeff_ideals(O, copy = false)
  for i = 1:degree(O)
    tt = icv[i]*t[1, i]
    tt = simplify(tt)
    d = lcm(d, denominator(tt))
  end
  return d
end

################################################################################
#
#  Random elements
#
################################################################################

RandomExtensions.maketype(O::NfRelOrd, ::Int) = elem_type(O)

function rand(rng::AbstractRNG, sp::SamplerTrivial{<:Make2{<:NfRelOrdElem,<:NfRelOrd,Int}})
  O, B = sp[][1:2]
  pb = pseudo_basis(O, copy = false)
  z = nf(O)()
  for i = 1:degree(O)
    t = rand(rng, pb[i][2], B)
    z += t*pb[i][1]
  end
  return O(z)
end

rand(O::NfRelOrd, B::Int) = rand(GLOBAL_RNG, O, B)
rand(rng::AbstractRNG, O::NfRelOrd, B::Int) = rand(rng, make(O, B))

################################################################################
#
#  Order generated by a set of elements
#
################################################################################

Order(elt::Vector{S}; check::Bool = false, cached::Bool = false) where {S <: Union{NfRelElem, NfRelNSElem}} = _order(elt, check = check)

function _order(elt::Vector{S}; check::Bool = false) where {S <: Union{NfRelElem, NfRelNSElem}}
  @assert length(elt) > 0
  K = parent(elt[1])
  n = degree(K)
  L = base_field(K)
  OL = maximal_order(L)
  bas = S[one(K)]

  for e = elt
    if check
      fe = minpoly(e)
      for i = 1:degree(fe)
        if !(coeff(fe, i-1) in OL)
          error("data does not define an order, $e is non-integral")
        end
      end
      df = degree(fe)-1
    else
      df = n-1
    end
    f = one(K)
    for i=1:df
      f *= e
      b = S[e*x for x in bas]
      append!(bas, b)
      if length(bas) >= n
        BK = basis_matrix(bas)
        B = pseudo_hnf(PseudoMatrix(BK), :lowerleft)
        rk = nrows(BK) - n + 1
        while iszero_row(BK, rk)
          rk += 1
        end
        B = sub(B, rk:nrows(B), 1:n)
        bas = _get_gens(B)
      end
    end
  end
  if nrows(B) != degree(K)
    error("Data does not define an order")
  end

  # Make an explicit check
  @hassert :NfRelOrd 1 defines_order(K, B)[1]
  return Order(K, B)
end

function _get_gens(M::PMat)
  mat = M.matrix
  ids = M.coeffs
  gens = Vector{NfRelElem{nf_elem}}()
  for i = 1:nrows(M)
    el = elem_from_mat_row(K, B.matrix, i)
    if isone(ids[i].num)
      push!(gens, divexact(el, ids[i].den))
    else
      push!(gens, divexact(el* ids[i].num.gen_one, ids[i].den))
      push!(gens, divexact(el*ids[i].num.gen_two.elem_in_nf, ids[i].den))
    end
  end
  return gens
end

function _get_gens(O::NfRelOrd)
  B = pseudo_basis(O)
  gens = Vector{NfRelElem{nf_elem}}()
  for i = 1:length(B)
    el = B[i][1]
    I = B[i][2]
    if isone(I.num)
      push!(gens, divexact(el, I.den))
    else
      push!(gens, divexact(I.num.gen_one*el, I.den))
      push!(gens, el*divexact(I.num.gen_two.elem_in_nf, I.den))
    end
  end
  return gens
end

function add_to_order(O::NfRelOrd, elt::Vector{T}; check::Bool = false) where T
  K = nf(O)
  k = base_field(K)
  Ok = maximal_order(k)
  B = basis_pmatrix(O)
  n = degree(O)
  lelt = length(elt)
  count = 0
  for e = elt
    count += 1
    @vprint :NfRelOrd 1 "Element $count / $lelt \n"
    if e in O
      continue
    end
    if check
      fe = minpoly(e)
      for i = 1:degree(fe)
        if !(coeff(fe, i-1) in Ok)
          error("data does not define an order, $e is non-integral")
        end
      end
      df = degree(fe)-1
    else
      df = n-1
    end
    f = one(K)

    for i=1:df
      f *= e
      if f in O
        break
      end
      bas = _get_gens(O)
      els_to_add = Vector{NfRelElem{nf_elem}}()
      for x in bas
        el = e*x
        if el in O
          continue
        end
        push!(els_to_add, el)
      end
      if isempty(els_to_add)
        break
      end
      BK = PseudoMatrix(basis_matrix(els_to_add))
      BK = vcat(B, BK)
      B = pseudo_hnf(BK, :lowerleft, true)
      rk = nrows(BK) - n + 1
      while iszero_row(BK.matrix, rk)
        rk += 1
      end
      B = sub(B, rk:nrows(B), 1:n)
      O = Order(K, B)
    end
  end
  return O
end

################################################################################
#
#  Dedekind composite
#
################################################################################

function dedekind_test_composite(O::NfRelOrd{U1, V, Z}, P::Union{NfRelOrdIdl, NfOrdIdl}) where {U1, V, Z <: NfRelElem}
  !isequation_order(O) && error("Order must be an equation order")

  L = nf(O)
  K = base_field(L)
  T = L.pol
  Kx = parent(T)
  OK = maximal_order(K)
  F, mF = quo(OK, P)
  Fy, y = PolynomialRing(F,"y", cached=false)

  t = map_coefficients(mF, map_coefficients(OK, T), parent = Fy)
  fail, g = gcd_with_failure(t, derivative(t))
  if !isone(fail)
    return K(fail.elem), O
  end
  h = divrem(t, g)[1]

  G = map_coefficients(K, map_coefficients(x -> x.elem, g), parent = Kx)::typeof(T)
  H = map_coefficients(K, map_coefficients(x -> x.elem, h), parent = Kx)::typeof(T)
  assure_2_normal(P)
  pi = anti_uniformizer(P)
  F = pi*(G*H - T)
  f = map_coefficients(mF, map_coefficients(OK, F), parent = Fy)

  fail, dd = gcd_with_failure(g, h)
  if !isone(fail)
    return K(fail.elem), O
  end
  fail, d = gcd_with_failure(f, dd)
  if !isone(fail)
    return K(fail.elem), O
  end

  if isone(d)
    return one(K), O
  end

  u = divrem(t, d)[1]
  U = map_coefficients(K, map_coefficients(x -> x.elem, u), parent = Kx)
  M = representation_matrix(pi*L(U))
  PM = PseudoMatrix(representation_matrix(pi*L(U)), [ K(1)*OK for i = 1:degree(O) ])
  BM = basis_pmatrix(O)
  PN = vcat(BM, PM)
  PN = sub(pseudo_hnf_full_rank_with_modulus(PN, P, :lowerleft), degree(O) + 1:2*degree(O), 1:degree(O))
  OO = typeof(O)(L, PN)
  OO.index = P
  OO.isequation_order = false
  return one(K), OO
end

function prefactorization_discriminant(K::NfRel, d::Union{NfRelOrdIdl, NfAbsOrdIdl})
  OK = order(d)
  @assert nf(OK) == base_field(K)
  f = K.pol
  factors = typeof(d)[]
  moduli = prefactorization(d)
  while !isempty(moduli)
    I = pop!(moduli)
    I = ispower(I)[2]
    if isprime(absolute_minimum(I))
      push!(factors, I)
      continue
    end
    Q, mQ = quo(OK, I)
    Qx = PolynomialRing(Q, cached = false)[1]
    fQ = map_coefficients(mQ, map_coefficients(OK, f) , parent = Qx)
    fQ1 = derivative(fQ)
    fail = gcd_with_failure(fQ, fQ1)[1]
    if isone(fail)
      push!(factors, I)
      continue
    end
    J = ideal(OK, fail.elem)
    cp = coprime_base(typeof(d)[J, I])
    append!(moduli, typeof(I)[Inew for Inew in cp if !iscoprime(I, Inew)])
  end
  return factors
end

function prefactorization(I::NfRelOrdIdl)
  OK = order(I)
  m = absolute_minimum(I)
  ideals = typeof(I)[]
  pp, r = Hecke._factors_trial_division(m)
  for p in pp
    push!(ideals, I + ideal(OK, p))
  end
  r = ispower(r)[2]
  if !isone(r)
    push!(ideals, I + ideal(OK, r))
  end
  return ideals
end

function maximal_order(O::NfRelOrd{S, T, U}) where {S, T, U <: NfRelElem}
  K = nf(O)
  L = base_field(K)
  OL = maximal_order(L)
  d = discriminant(O)
  facts = prefactorization_discriminant(K, d)
  sort!(facts, by = x -> absolute_minimum(x), rev = true)
  @vprint :NfRelOrd 1 "Factors of the discriminant lying over $([minimum(x) for x in facts]) \n"
  E = EquationOrder(K)
  OO = O
  while !isempty(facts)
    p = pop!(facts)
    pm = absolute_minimum(p)
    if isprime_power(pm)
      @vprint :NfRelOrd 1 "Factoring ideal over $(pm)\n"
      @vtime :NfRelOrd 1 lf = factor(p)
      for q in keys(lf)
        @vprint :NfRelOrd 1 "Computing pmaximal order for $(q)\n"
        @vtime :NfRelOrd 1 Oq = pmaximal_overorder(O, q)
        @vtime :NfRelOrd 1 OO = sum_as_OK_modules(OO, Oq)
      end
    else
      @vprint :NfRelOrd 1 "Dedekind test for ideal lying over $(pm)\n"
      @vtime :NfRelOrd 1 fail, E1 = Hecke.dedekind_test_composite(E, p)
      if !isone(fail)
        J = ideal(OL, OL(fail))
        cp = coprime_base(typeof(p)[J, p])
        for q in cp
          if !iscoprime(q, p)
            push!(facts, q)
          end
        end
        continue
      end
      if discriminant(E) != discriminant(E1)
        @vtime :NfRelOrd 1 OO = sum_as_OK_modules(OO, E1)
      end
      g = gcd(discriminant(OO), p)
      if !isone(g)
        @vprint :NfRelOrd 1 "Factoring ideal over $(absolute_minimum(g))\n"
        @vtime :NfRelOrd 1 lf = factor(g)
        for q in keys(lf)
          @vprint :NfRelOrd 1 "Computing pmaximal order for $(q)\n"
          @vtime :NfRelOrd 1 Oq = pmaximal_overorder(O, q)
          @vtime :NfRelOrd 1 OO = sum_as_OK_modules(OO, Oq)
        end
      end
    end
  end
  OO.ismaximal = 1
  return OO
end

function overorder_polygons(O::NfRelOrd{S, T, NfRelElem{nf_elem}}, p::NfOrdIdl) where {S, T}
  @assert isequation_order(O)
  K = nf(O)
  f = K.pol
  k = base_field(K)
  kt = parent(f)
  Ok = maximal_order(k)
  F, mF = ResidueField(Ok, p)
  mF1 = extend(mF, k)
  f1 = map_coefficients(mF1, f)
  sqf = factor_squarefree(f1)
  l = powers(gen(K), degree(K)-1)
  maxv = 0
  regular = true
  vdisc = 0
  for (gg, m) in sqf
    isone(m) && continue
    fac = factor(gg)
    for (g, m1) in fac
      phi = map_coefficients(pseudo_inv(mF1), g, parent = kt)
      dev, quos = phi_development_with_quos(f, phi)
      N = _newton_polygon(dev, p)
      for i = 1:m
        v = _floor_newton_polygon(N, i)
        maxv = max(v, maxv)
        if v > 0
          vdisc += v*degree(phi)
          pow_anti = anti_uniformizer(p)^v
          for j = 1:degree(phi)
            q1 = shift_left(quos[i], j-1)
            push!(l, mod(K(q1)*pow_anti, minimum(p, copy = false)))
          end
        end
      end
    end
  end
  B = basis_matrix(l)
  M = pseudo_matrix(B)
  index = p^maxv
  M = sub(pseudo_hnf_full_rank_with_modulus(M, index, :lowerleft), length(l)-degree(K)+1:length(l), 1:degree(K))
  O1 = typeof(O)(K, M)
  O1.index = index
  O1.disc_abs = divexact(O.disc_abs, p^(2*vdisc))
  return O1
end

################################################################################
#
#   Is ramified
#
################################################################################

function isramified(R::NfRelOrd, p::T) where T <: Union{NfAbsOrdIdl, NfRelOrdIdl, fmpz, Int}
  @assert isprime(p)
  D = prime_decomposition(R, p)
  for (_, e) in D
    if e > 1
      return true
    end
  end
  return false
end


text:
############################################################
# Groups
############################################################
"""
    Base.keys(g::NCDatasets.Groups)

Return the names of all subgroubs of the group `g`.
"""
function Base.keys(g::Groups)
    return String[nc_inq_grpname(ncid)
                  for ncid in nc_inq_grps(g.ncid)]
end

"""
    group = getindex(g::NCDatasets.Groups,groupname::AbstractString)

Return the NetCDF `group` with the name `groupname`.
For example:

```julia-repl
julia> ds = NCDataset("results.nc", "r");
julia> forecast_group = ds.group["forecast"]
julia> forecast_temp = forecast_group["temperature"]
```

"""
function Base.getindex(g::Groups,groupname::AbstractString)
    grp_ncid = nc_inq_grp_ncid(g.ncid,groupname)
    return NCDataset(grp_ncid,g.isdefmode)
end

"""
    defGroup(ds::NCDataset,groupname, attrib = []))

Create the group with the name `groupname` in the dataset `ds`.
`attrib` is a list of attribute name and attribute value pairs (see `NCDataset`).
"""
function defGroup(ds::NCDataset,groupname; attrib = [])
    grp_ncid = nc_def_grp(ds.ncid,groupname)
    ds = NCDataset(grp_ncid,ds.isdefmode)

    # set global attributes for group
    for (attname,attval) in attrib
        ds.attrib[attname] = attval
    end

    return ds
end
export defGroup

group(ds::AbstractDataset,groupname) = ds.group[groupname]

"""
    groupname(ds::NCDataset)
Return the group name of the NCDataset `ds`
"""
groupname(ds::NCDataset) = nc_inq_grpname(ds.ncid)
export groupname

ext: jl


text:
"""
    accumulate(collection, operation)
Takes in a collection and applies the operation on that collection then returns
a new collection.
"""
function accumulate(collection, operation)

    new_collection = operation.(collection, collection)

end


text:
module RemoteS

using GMT, Printf, Statistics, Requires, Dates

const SCENE_HALFW = Dict("AQUA" => 1163479, "TERRA" => 1163479, "LANDSAT8" => 92500)	# half widths

if isdefined(Base, :Experimental) && isdefined(Base.Experimental, Symbol("@optlevel"))
	@eval Base.Experimental.@optlevel 1
end

export
	cutcube, subcube, dn2temperature, dn2radiance, dn2reflectance, reflectance_surf, grid_at_sensor, truecolor,
	clg, clre, evi, evi2, gndvi, mndwi, mtci, mcari, msavi, mbri, ndvi, ndwi, ndwi2, ndrei1,
	ndrei2, satvi, savi, slavi,
	clip_orbits, findscenes, sat_scenes, sat_tracks, reportbands

include("grid_at_sensor.jl")
include("spectral_indices.jl")
include("utils.jl")

function __init__()
	@require SatelliteToolbox="6ac157d9-b43d-51bb-8fab-48bf53814f4a" include("sat_tracks.jl")
end

end # module

max_stars_repo_name: GenericMappingTools/RemoteS

max_issues_repo_name: GenericMappingTools/RemoteS

max_forks_repo_name: GenericMappingTools/RemoteS


text:
using PackageCompiler

# FIXME this probably would compile all functions in these packages? If that's
# slow, try identify the key functions to compile
create_sysimage([:CuArrays, :Zygote,
                 :Luxor,
                 :JSON,
                 :ProgressMeter,
                 :ProfileView,
                 :Flux
                 ],
                # sysimage_path="myimage.so",
                replace_default=true)

## run this file to create the image
# julia --project precompile.jl

## to use the image
# julia --sysimage myimage.so --project main.jl
# restore_default_sysimage()

max_stars_repo_name: lihebi/bhdl

max_issues_repo_name: lihebi/bhdl

max_forks_repo_name: lihebi/bhdl


text:
module YasolSolver
using Revise

import MathOptInterface
const MOI = MathOptInterface
using EzXML

include("MOI_wrapper/MOI_wrapper.jl")

include("JuMP.jl")

# change parameter in Yasol.ini file
function setInitialParameter(yasolDir::String, parameter::String, value::Int64)
    # create temp copy
    mv(joinpath(yasolDir, "Yasol.ini"), joinpath(yasolDir, "Yasol_temp.ini"))

    # clear original file
    open(joinpath(yasolDir, "Yasol.ini"), "w") do f
        write(f, "")
    end

    # copy updated values
    open(joinpath(yasolDir, "Yasol_temp.ini"), "r") do f
        open(joinpath(yasolDir, "Yasol.ini"), "a") do i
            while !eof(f)
              x = readline(f)
              par = rsplit(x, "=")
              if par[1] == parameter
                  # overwrite parameter value
                  write(i, parameter * "=" * string(value) * "\n")
              else
                  write(i, x * "\n")
              end
            end
        end
    end

    # delete temp file
    rm(joinpath(yasolDir, "Yasol_temp.ini"))
end

# return all initial parameter
function getInitialParameter(yasolDir::String)
    result = []
    open(joinpath(yasolDir, "Yasol.ini")) do file
        for ln in eachline(file)
            if ln != "END"
                push!(result, ln)
            end
        end
    end
    return result
end

# import and return solution
function importSolution(solPath::String)
    doc = readxml(solPath)

    objective_value = 0.0
    runtime = 0.0
    solutionStatus = ""
    gap = 0.0
    values = Dict{String,Float64}()

    for node in eachelement(doc.root)
        if node.name == "header"
            objective_value = parse(Float64, node["ObjectiveValue"])
            runtime = parse(Float64, rsplit(node["Runtime"], " ")[1])
        elseif node.name == "quality"
            solutionStatus = node["SolutionStatus"]
            gap = parse(Float64, node["Gap"])
        elseif node.name == "variables"
            for var in eachelement(node)
                values[var["name"]] = parse(Float64, var["value"])
            end
        end
    end

    res = _Results(objective_value, runtime, solutionStatus, gap, values)

    return res
end

# print solution
function printSolution(res::YasolSolver._Results)
    println("---Solution---")
    println("Objective value: " * string(res.objective_value))
    println("Runtime: " * string(res.runtime))
    println("Solution status: " * string(res.solutionStatus))
    println("Gap: " * string(res.gap))
    println("Variable values: ")
    for (key, value) in res.values
        println(key * ": " * string(value))
    end
    println("---End---")
end

end


text:
function getparticledata( dataobject::InfoType,
                        Nvarp::Int,
                        nvarp_corr::Array{Int,1},
                        stars::Bool,
                        lmax::Real,
                        ranges::Array{Float64,1},
                        print_filenames::Bool,
                        verbose::Bool,
                        read_cpu::Bool )

    output = dataobject.output
    lmin = dataobject.levelmin
    path = dataobject.path
    ndim = dataobject.ndim
    boxlen = dataobject.boxlen

    omega_m = dataobject.omega_m
    omega_l = dataobject.omega_l
    omega_k = dataobject.omega_k
    h0  = dataobject.H0
    aexp = dataobject.aexp

    xmin, xmax, ymin, ymax, zmin, zmax = ranges # 6 elements

    idom = zeros(Int32, 8)
    jdom = zeros(Int32, 8)
    kdom = zeros(Int32, 8)

    bounding_min = zeros(Float64, 8)
    bounding_max = zeros(Float64, 8)

    cpu_min = zeros(Int32, 8)
    cpu_max = zeros(Int32, 8)

    dmax= maximum([xmax-xmin,ymax-ymin,zmax-zmin])
    ilevel = 1 # define ilevel before loop
    for il=1:lmax
        ilevel=il
        dx=0.5^ilevel
        if dx < dmax break end
    end

    bit_length=ilevel-1
    maxdom=2^bit_length

    # Todo: floor = fotran int ?
    if bit_length > 0
        imin=floor(Int32, xmin * maxdom)
        imax=imin+1
        jmin=floor(Int32, ymin * maxdom)
        jmax=jmin+1
        kmin=floor(Int32, zmin * maxdom)
        kmax=kmin+1
    else
        imin=0
        imax=0
        jmin=0
        jmax=0
        kmin=0
        kmax=0
    end

    dkey=(2^(dataobject.grid_info.nlevelmax+1)/maxdom)^dataobject.ndim

    if bit_length>0 ndom=8 else ndom=1 end

    idom = [imin, imax, imin, imax, imin, imax, imin, imax]
    jdom = [jmin, jmin, jmax, jmax, jmin, jmin, jmax, jmax]
    kdom = [kmin, kmin, kmin, kmin, kmax, kmax, kmax, kmax]

    order_min=0.0e0 #todo: predeclare
    for i=1:ndom
        if bit_length > 0
            order_min = hilbert3d(idom[i],jdom[i],kdom[i],bit_length,1)
        else
            order_min=0.0e0
        end

        bounding_min[i]=(order_min)*dkey
        bounding_max[i]=(order_min+1.)*dkey
    end

    for impi=1:dataobject.ncpu
        for i=1:ndom
            if (dataobject.grid_info.bound_key[impi] <= bounding_min[i] &&
                dataobject.grid_info.bound_key[impi+1] > bounding_min[i])
                cpu_min[i]=impi
            end
            if (dataobject.grid_info.bound_key[impi] < bounding_max[i] &&
                dataobject.grid_info.bound_key[impi+1] >= bounding_max[i])
                cpu_max[i]=impi
            end
        end
    end

    cpu_read = copy(dataobject.grid_info.cpu_read)
    cpu_list = zeros(Int32, dataobject.ncpu)
    ncpu_read=Int32(0)
    for i=1:ndom #Todo
        for j=(cpu_min[i]):(cpu_max[i])
            if cpu_read[j]==false
                ncpu_read=ncpu_read+1
                cpu_list[ncpu_read]=j
                cpu_read[j]=true
            end
        end
    end

    if read_cpu
        if dataobject.descriptor.pversion == 0
            pos_1D, vars_1D, cpus_1D, identity_1D, levels_1D = readpart( dataobject,
                                     Nvarp=Nvarp, nvarp_corr=nvarp_corr,
                                     lmax=lmax,
                                     ranges=ranges,
                                     cpu_list=cpu_list,
                                     ncpu_read=ncpu_read,
                                     stars=stars,
                                     read_cpu=read_cpu,
                                     verbose=verbose,
                                     print_filenames=print_filenames )

            return pos_1D, vars_1D, cpus_1D, identity_1D, levels_1D

        elseif dataobject.descriptor.pversion > 0
            pos_1D, vars_1D, cpus_1D, family_1D, tag_1D, levels_1D = readpart( dataobject,
                                     Nvarp=Nvarp, nvarp_corr=nvarp_corr,
                                     lmax=lmax,
                                     ranges=ranges,
                                     cpu_list=cpu_list,
                                     ncpu_read=ncpu_read,
                                     stars=stars,
                                     read_cpu=read_cpu,
                                     verbose=verbose,
                                     print_filenames=print_filenames )

            return  pos_1D, vars_1D, cpus_1D, identity_1D, family_1D, tag_1D, levels_1D
        end # if pversion
    else # if read_cpu
        if dataobject.descriptor.pversion == 0
            pos_1D, vars_1D, identity_1D, levels_1D = readpart( dataobject,
                                     Nvarp=Nvarp, nvarp_corr=nvarp_corr,
                                     lmax=lmax,
                                     ranges=ranges,
                                     cpu_list=cpu_list,
                                     ncpu_read=ncpu_read,
                                     stars=stars,
                                     read_cpu=read_cpu,
                                     verbose=verbose,
                                     print_filenames=print_filenames )

            return pos_1D, vars_1D, identity_1D, levels_1D

        elseif dataobject.descriptor.pversion > 0
            pos_1D, vars_1D, family_1D, identity_1D, tag_1D, levels_1D = readpart( dataobject,
                                     Nvarp=Nvarp, nvarp_corr=nvarp_corr,
                                     lmax=lmax,
                                     ranges=ranges,
                                     cpu_list=cpu_list,
                                     ncpu_read=ncpu_read,
                                     stars=stars,
                                     read_cpu=read_cpu,
                                     verbose=verbose,
                                     print_filenames=print_filenames )

            return pos_1D, vars_1D, identity_1D, family_1D, tag_1D, levels_1D
        end # if pversion
    end
end

function readpart(dataobject::InfoType;
                            Nvarp::Int,
                            nvarp_corr::Array{Int,1},
                            lmax::Int=dataobject.levelmax,
                            ranges::Array{Float64,1}=[0.,1.],
                            cpu_list::Array{Int32,1}=[1],
                            ncpu_read::Int=0,
                            stars::Bool=true,
                            read_cpu::Bool=false,
                            verbose::Bool=verbose_mode,
                            print_filenames::Bool=false)

    boxlen = dataobject.boxlen
    path = dataobject.path
    ndim = dataobject.ndim

    fnames = createpath(dataobject.output, path)

    vars_1D = ElasticArray{Float64}(undef, Nvarp, 0)
    r1, r2, r3, r4, r5, r6 = ranges .* boxlen

    pos_1D = ElasticArray{Float64}(undef, 3, 0)
    if read_cpu cpus_1D = Array{Int}(undef, 0) end #zeros(Int, npart)
    identity_1D  = Array{Int32}(undef, 0) #zeros(Int32, npart)
    #if dataobject.descriptor.pversion == 0

    if dataobject.descriptor.pversion > 0
        #identity_1D  = Array{Int32}(undef, 0) #zeros(Int32, npart)
        family_1D  = Array{Int8}(undef, 0) #zeros(Int32, npart)
        tag_1D  = Array{Int8}(undef, 0) #zeros(Int32, npart)
    end
    levels_1D =  Array{Int32}(undef, 0) #zeros(Int32, npart) #manu

    ndim2=Int32(0)
    #ngrida=0
    parti=Int32(0) # particle iterator

    @showprogress 1 "Reading data..." for k=1:ncpu_read

       icpu=cpu_list[k]

       partpath = getproc2string(fnames.particles, icpu)
       f_part = FortranFile(partpath)

       skiplines(f_part, 1) #ncpu2
       ndim2 = read(f_part, Int32) #, (Int32, ndim2))
       npart2 = read(f_part, Int32) #, (Int32, npart2))
       skiplines(f_part, 1)
       nstar = read(f_part, Int32)
       skiplines(f_part, 3)

       if npart2 != 0
           pos_1D_buffer = zeros(Float64, 3, npart2)

           vars_1D_buffer = zeros(Float64, Nvarp, npart2 )

           if dataobject.descriptor.pversion > 0

               family_1D_buffer = zeros(Int8, npart2)
               tag_1D_buffer = zeros(Int8, npart2)
           end
           identity_1D_buffer = zeros(Int32, npart2)
           levels_1D_buffer   = zeros(Int32, npart2)

           # Read position
           for i=1:ndim
               #x[:,i] = read(f_part, (Float64, npart2)
               pos_1D_buffer[i, 1:npart2] = read(f_part, (Float64, npart2) )
           end

           pos_selected = (pos_1D_buffer[1,:] .>= r1) .& (pos_1D_buffer[1,:] .<= r2) .&
                           (pos_1D_buffer[2,:] .>= r3) .& (pos_1D_buffer[2,:] .<= r4) .&
                           (pos_1D_buffer[3,:] .>= r5) .& (pos_1D_buffer[3,:] .<= r6)

            ls = Int32(length( pos_1D_buffer[1, pos_selected] ) ) #selected length #todo int32

            if ls != 0

                append!(pos_1D, pos_1D_buffer[:, pos_selected] )

                # Read velocity
                for i=1:ndim
                    if nvarp_corr[i] != 0
                    #vel[:,i] = read(f_part, (Float64, npart2)
                        vars_1D_buffer[nvarp_corr[i], 1:npart2] = read(f_part, (Float64, npart2) )
                    else
                        skiplines(f_part, 1)
                    end
                    #read(f_part)
                end

                # Read mass
                if nvarp_corr[4] != 0
                    vars_1D_buffer[nvarp_corr[4], 1:npart2] = read(f_part, (Float64, npart2) )
                else
                    skiplines(f_part, 1)
                end

                # Read identity
                #if dataobject.descriptor.pversion == 0
                identity_1D_buffer[1:npart2] = read(f_part, (Int32, npart2) ) # identity
                append!(identity_1D, identity_1D_buffer[pos_selected]) # identity
                #end

                # Read level
                levels_1D_buffer[1:npart2] = read(f_part, (Int32, npart2) ) # level
                append!(levels_1D, levels_1D_buffer[pos_selected]) # level

                # Read family, tag
                if dataobject.descriptor.pversion > 0
                    #skiplines(f_part, 1) # skip identity
                    family_1D_buffer[1:npart2] = read(f_part, (Int8, npart2) ) # family
                    tag_1D_buffer[1:npart2] = read(f_part, (Int8, npart2) ) # tag

                    append!(family_1D, family_1D_buffer[pos_selected]) # family
                    append!(tag_1D, tag_1D_buffer[pos_selected]) # tag

                    # read birth
                    if nstar>0 && nvarp_corr[7] != 0
                       #birth = read(f_part, (Float64, npart2)
                       #skiplines(f_part, 1)
                       vars_1D_buffer[nvarp_corr[7], 1:npart2] = read(f_part, (Float64, npart2) ) #age (birth)
                    elseif nstar>0 && nvarp_corr[7] == 0
                        skiplines(f_part, 1)
                    end

                elseif dataobject.descriptor.pversion == 0
                    # read birth
                    if nstar>0 && nvarp_corr[5] != 0
                       #birth = read(f_part, (Float64, npart2)
                       #skiplines(f_part, 1)
                       vars_1D_buffer[nvarp_corr[5], 1:npart2] = read(f_part, (Float64, npart2) ) #age (birth)
                    elseif nstar>0 && nvarp_corr[5] == 0
                        skiplines(f_part, 1)
                    end
                end

                # read additional variables
                if Nvarp>7
                    for iN = 8:Nvarp
                        if nvarp_corr[iN] != 0
                            vars_1D_buffer[nvarp_corr[iN], 1:npart2] = read(f_part, (Float64, npart2) )
                        end
                    end
                end

                append!(vars_1D, vars_1D_buffer[:, pos_selected] )
                if read_cpu append!(cpus_1D, fill(k,ls) ) end
                parti = parti + ls
            end
        end
        close(f_part)

    end #for

    if read_cpu
        if dataobject.descriptor.pversion == 0
            return pos_1D, vars_1D, cpus_1D, identity_1D, levels_1D
        elseif dataobject.descriptor.pversion > 0
            return pos_1D, vars_1D, cpus_1D, identity_1D, family_1D, tag_1D, levels_1D
        end
    else
        if dataobject.descriptor.pversion == 0
            return pos_1D, vars_1D, identity_1D, levels_1D
        elseif dataobject.descriptor.pversion > 0
            return pos_1D, vars_1D, identity_1D, family_1D, tag_1D, levels_1D
        end
    end
end


text:
export is_uniform

"""
`NV(H::SimpleHypergraph)` is the number of vertices in `H`.
"""
NV(H::SimpleHypergraph) = length(H.V)

"""
`NE(H::SimpleHypergraph)` is the number of (hyper)edges in `H`.
"""
NE(H::SimpleHypergraph) = length(H.E)

function show(io::IO, H::SimpleHypergraph{T}) where {T}
    suffix = " (n=$(NV(H)), m=$(NE(H)))"
    print(io, "SimpleHypergraph{$T}" * suffix)
end

"""
`has(H::SimpleHypergraph{T}, v::T)` tests if `v` is a vertex of `H`.

`has(H::SimpleHypergraph{T}, e::Set{T})` tests if `e` is an edge of `H`.
Also works if `e` is a `Vector{T}` or a list of two or more arguments.
"""
function has(H::SimpleHypergraph{T}, v::T)::Bool where {T}
    return in(v, H.V)
end

function has(H::SimpleHypergraph{T}, e::Set{T})::Bool where {T}
    return in(e, H.E)
end

function has(H::SimpleHypergraph{T}, e::Vector{T})::Bool where {T}
    return has(H, Set(e))
end

function has(H::SimpleHypergraph{T}, e::T...)::Bool where {T}
    return has(H, set(e))
end

"""
`H[v]` for a hypergraph `H` and vertex `v` returns the set of edges
that contain vertex `v`. This throws an error if `v` is not a vertex
of `H`.
"""
function getindex(H::SimpleHypergraph{T}, v::T)::Set{Set{T}} where {T}
    if !has(H, v)
        error("No such vertex: $v")
    end
    return H.VE[v]
end

function deg(H::SimpleHypergraph{T}, v::T)::Int where {T}
    return length(H[v])
end

function eltype(H::SimpleHypergraph{T}) where {T}
    return T
end

"""
`vlist(H::SimpleHypergraph)` returns a list of the vertices
in `H`.
"""
function vlist(H::SimpleHypergraph)
    return collect(H.V)
end

"""
`elist(H::SimpleHypergraph)` returns a list of the edges
in `H`.
"""
function elist(H::SimpleHypergraph)
    return collect(H.E)
end

function is_uniform(H::SimpleHypergraph{T})::Bool where {T}
    if NE(H) < 2
        return true
    end
    klist = length.(collect(H.E))  # get sizes of all edges
    n = length(unique(klist))      # throw away duplicates
    return n == 1
end


text:
using Distributed
addprocs(2)
@everywhere using ArrayChannels

function main()
    local iterations, payload
    if length(ARGS) == 0
        iterations, payload = 1000, 1000000
    elseif length(ARGS) == 2
        iterations, payload = map(x -> parse(Int,x), ARGS)
    else
        println("Usage: ping_pong_array_channels <# iterations> <vector_length>")
        exit(1)
    end

    # Allocate the array channel
    vector_channel = ArrayChannel(Float64, workers(), payload)

    @sync for proc in workers()
        @spawnat proc work(vector_channel, iterations, payload)
    end
end

@everywhere mean(l) = sum(l) / length(l)

@everywhere function work(vector_channel, iterations, payload)
    other_worker = myid() == 2 ? 3 : 2
    local t0, t1
    for k in 1 : 2 * iterations
        if k == iterations + 1
            t0 = time_ns()
        end

        if k % 2 == myid() - 2
            j = k % payload + 1
            i = if (j != 0) j else payload end
            vector_channel[i] += Float64(k)
            put!(vector_channel, other_worker)
        else
            take!(vector_channel, other_worker)
        end
    end
    t1 = time_ns()

    throughput = iterations * payload * 8.0 * 1e-6 / ((t1-t0)*1e-9)

    if myid() == 2
        println("Rate (MB/s): $throughput")
    end
end

main()


text:
# Modelo de tanque de armazenamento com quatro alimentacoes
type Tanque4
	Tanque4()=begin
		PP2=outers.PP2
		new(
			DanaPlugin(Dict{Symbol,Any}(
				:Brief=>"Steam tables"
			)),
			Corrente(Dict{Symbol,Any}(
				:Symbol=>"_{in1}"
			)),
			Corrente(Dict{Symbol,Any}(
				:Symbol=>"_{in2}"
			)),
			Corrente(Dict{Symbol,Any}(
				:Symbol=>"_{in3}"
			)),
			Corrente(Dict{Symbol,Any}(
				:Symbol=>"_{in4}"
			)),
			Corrente(Dict{Symbol,Any}(
				:Symbol=>"_{out}"
			)),
			[
				:(Fout.F = Fin1.F + Fin2.F + Fin3.F + Fin4.F),
				:(Fout.F * Fout.H = Fin1.F * Fin1.H + Fin2.F * Fin2.H + Fin3.F * Fin3.H + Fin4.F * Fin4.H),
				:([Fout.S,Fout.T] = PP2.propPH(Fout.P,Fout.H)),
			],
			[
				"","","",
			],
			[:PP2,],
			[:Fin1,:Fin2,:Fin3,:Fin4,:Fout,]
		)
	end
	PP2::DanaPlugin
	Fin1::Corrente
	Fin2::Corrente
	Fin3::Corrente
	Fin4::Corrente
	Fout::Corrente
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	parameters::Array{Symbol,1}
	variables::Array{Symbol,1}
end
export Tanque4
function setEquationFlow(in::Tanque4)
	addEquation(1)
	addEquation(2)
	addEquation(3)
end


text:
abstract type AbstractFeature end
abstract type AbstractPartOfSpeech <: AbstractFeature end
abstract type AbstractNominal <: AbstractPartOfSpeech end
abstract type AbstractState <: AbstractNominal end
abstract type AbstractCase <: AbstractNominal end
abstract type AbstractPreposition <: AbstractPartOfSpeech end
abstract type AbstractParticle <: AbstractPartOfSpeech end
abstract type AbstractDisLetters <: AbstractPartOfSpeech end
abstract type AbstractConjunction <: AbstractPartOfSpeech end
abstract type AbstractPerson <: AbstractPartOfSpeech end
abstract type AbstractGender <: AbstractPartOfSpeech end
abstract type AbstractNumber <: AbstractPartOfSpeech end
abstract type AbstractNoun <: AbstractPartOfSpeech end
abstract type AbstractDerivedNominal <: AbstractPartOfSpeech end
abstract type AbstractPronoun <: AbstractPartOfSpeech end
abstract type AbstractAdverb <: AbstractPartOfSpeech end
abstract type AbstractVerb <: AbstractPartOfSpeech end
abstract type AbstractAspect <: AbstractVerb end
abstract type AbstractMood <: AbstractVerb end
abstract type AbstractVoice <: AbstractVerb end
abstract type AbstractForm <: AbstractVerb end
abstract type AbstractPrefix <: AbstractPartOfSpeech end
abstract type AbstractSuffix <: AbstractPartOfSpeech end
abstract type AbstractDerivedNoun <: AbstractPartOfSpeech end

macro partofspeech(name, parent, data, desc, ar_label)
    esc(quote
        struct $name <: $parent
            data::Symbol
            desc::String
            ar_label::String
        end

        $name() = $name($data, $desc, $ar_label)
    end)
end

@partofspeech Noun AbstractNoun Symbol("N") "Noun" "اسم"
@partofspeech ProperNoun AbstractNoun Symbol("PN") "Proper noun" "اسم علم"

@partofspeech Adjective AbstractDerivedNominal Symbol("ADJ") "Adjective" "صفة"
@partofspeech ImperativeVerbalNoun AbstractDerivedNominal Symbol("IMPN") "Imperative verbal noun" "اسم فعل أمر"

@partofspeech Personal AbstractPronoun Symbol("PRON") "Personal pronoun" "ضمير"
@partofspeech Demonstrative AbstractPronoun Symbol("DEM") "Demonstrative pronoun" "اسم اشارة"
@partofspeech Relative AbstractPronoun Symbol("REL") "Relative pronoun" "اسم موصول"

@partofspeech Time AbstractAdverb Symbol("T") "Time adverb" "ظرف زمان"
@partofspeech Location AbstractAdverb Symbol("LOC") "Location adverb" "ظرف مكان"

@partofspeech Preposition AbstractPreposition Symbol("P") "Preposition" "حرف جر"

@partofspeech EmphaticLam AbstractPrefix Symbol("EMPH") "Emphatic lam prefix" "لام التوكيد"
@partofspeech ImperativeLam AbstractPrefix Symbol("IMPV") "Imperative lam prefix" "لام الامر"
@partofspeech PurposeLam AbstractPrefix Symbol("PRP") "Purpose lam prefix" "لام التعليل"
@partofspeech EmphaticNun AbstractPrefix Symbol("+n:EMPH") "Emphatic lam prefix" "لام التوكيد"

@partofspeech Coordinating AbstractConjunction Symbol("CONJ") "Coordinating conjunction" "حرف عطف"
@partofspeech Subordinating AbstractConjunction Symbol("SUB") "Subordinating particle" "حرف مصدري"

@partofspeech Accusative AbstractParticle Symbol("ACC") "Accusative particle" "حرف نصب"
@partofspeech Amendment AbstractParticle Symbol("AMD") "Amendment particle" "حرف استدراك"
@partofspeech Answer AbstractParticle Symbol("ANS") "Answer particle" "حرف جواب"
@partofspeech Aversion AbstractParticle Symbol("AVR") "Aversion particle" "حرف ردع"
@partofspeech Cause AbstractParticle Symbol("CAUS") "Particle of cause" "حرف سببية"
@partofspeech Certainty AbstractParticle Symbol("CERT") "Particle of certainty" "حرف تحقيق"
@partofspeech Circumstantial AbstractParticle Symbol("CIRC") "Circumstantial particle" "حرف حال"
@partofspeech Comitative AbstractParticle Symbol("COM") "Comitative particle" "واو المعية"
@partofspeech Conditional AbstractParticle Symbol("COND") "Conditional particle" "حرف شرط"
@partofspeech Equalization AbstractParticle Symbol("EQ") "Equalization particle" "حرف تسوية"
@partofspeech Exhortation AbstractParticle Symbol("EXH") "Exhortation particle" "حرف تحضيض"
@partofspeech Explanation AbstractParticle Symbol("EXL") "Explanation particle" "حرف تفصيل"
@partofspeech Exceptive AbstractParticle Symbol("EXP") "Exceptive particle" "أداة استثناء"
@partofspeech Future AbstractParticle Symbol("FUT") "Future particle" "حرف استقبال"
@partofspeech Inceptive AbstractParticle Symbol("INC") "Inceptive particle" "حرف ابتداء"
@partofspeech Interpretation AbstractParticle Symbol("INT") "Inceptive particle" "حرف تفسير"
@partofspeech Interogative AbstractParticle Symbol("INTG") "Interogative particle" "حرف استفهام"
@partofspeech Negative AbstractParticle Symbol("NEG") "Negative particle" "حرف نفي"
@partofspeech Preventive AbstractParticle Symbol("PREV") "Preventive particle" "حرف كاف"
@partofspeech Prohibition AbstractParticle Symbol("PRO") "Prohibition particle" "حرف نهي"
@partofspeech Resumption AbstractParticle Symbol("REM") "Resumption particle" "حرف استئنافية"
@partofspeech Restriction AbstractParticle Symbol("RES") "Restriction particle" "أداة حصر"
@partofspeech Retraction AbstractParticle Symbol("RET") "Retraction particle" "حرف اضراب"
@partofspeech Result AbstractParticle Symbol("RSLT") "Result particle" "حرف واقع في جواب الشرط"
@partofspeech Supplemental AbstractParticle Symbol("SUP") "Suplemental particle" "حرف زائد"
@partofspeech Surprise AbstractParticle Symbol("SUR") "Surprise particle" "حرف فجاءة"
@partofspeech Vocative AbstractParticle Symbol("VOC") "Vocative particle" "حرف نداء"

@partofspeech DisconnectedLetters AbstractDisLetters Symbol("INL") "Quranic initials" "حروف مقطعة"

@partofspeech FirstPerson AbstractPerson Symbol("1") "First person" "الاسناد"
@partofspeech SecondPerson AbstractPerson Symbol("2") "Second person" "الاسناد"
@partofspeech ThirdPerson AbstractPerson Symbol("3") "Third person" "الاسناد"

@partofspeech Masculine AbstractGender Symbol("M") "Masculine" "الجنس"
@partofspeech Feminine AbstractGender Symbol("F") "Feminine" "الجنس"

@partofspeech Singular AbstractNumber Symbol("S") "Singular" "العدد"
@partofspeech Dual AbstractNumber Symbol("D") "Dual" "العدد"
@partofspeech Plural AbstractNumber Symbol("P") "Plural" "العدد"

@partofspeech Verb AbstractPartOfSpeech Symbol("V") "Verb" "فعل"
@partofspeech Perfect AbstractAspect Symbol("PERF") "Perfect verb" "فعل ماض"
@partofspeech Imperfect AbstractAspect Symbol("IMPF") "Imperfect verb" "فعل مضارع"
@partofspeech Imperative AbstractAspect Symbol("IMPV") "Imperative verb" "فعل أمر"
@partofspeech Indicative AbstractMood Symbol("IND") "Indicative mood (default)" "مرفوع"
@partofspeech Subjunctive AbstractMood Symbol("SUBJ") "Subjunctive mood" "منصوب"
@partofspeech Jussive AbstractMood Symbol("JUS") "Jussive mood" "مجزوم"
@partofspeech Active AbstractVoice Symbol("ACT") "Active voice (default)" "مبني للمعلوم"
@partofspeech Passive AbstractVoice Symbol("PASS") "Passive voice" "مبني للمجهول"

@partofspeech VerbFormI AbstractForm Symbol("I") "First verb form (default)" "فعل"
@partofspeech VerbFormII AbstractForm Symbol("II") "Second verb form" "فعل"
@partofspeech VerbFormIII AbstractForm Symbol("III") "Third verb form" "فعل"
@partofspeech VerbFormIV AbstractForm Symbol("IV") "Fourth verb form" "فعل"
@partofspeech VerbFormV AbstractForm Symbol("V") "Fifth verb form" "فعل"
@partofspeech VerbFormVI AbstractForm Symbol("VI") "Sixth verb form" "فعل"
@partofspeech VerbFormVII AbstractForm Symbol("VII") "Seventh verb form" "فعل"
@partofspeech VerbFormVIII AbstractForm Symbol("VIII") "Eighth verb form" "فعل"
@partofspeech VerbFormIX AbstractForm Symbol("IX") "Ninth verb form" "فعل"
@partofspeech VerbFormX AbstractForm Symbol("X") "Tenth verb form" "فعل"
@partofspeech VerbFormXI AbstractForm Symbol("XI") "Eleventh verb form" "فعل"
@partofspeech VerbFormXII AbstractForm Symbol("XII") "Twelfth verb form" "فعل"

@partofspeech ActiveParticle AbstractDerivedNoun Symbol("ACT PCPL") "Active particle" "اسم فاعل"
@partofspeech PassiveParticle AbstractDerivedNoun Symbol("PASS PCPL") "Passive particle" "اسم مفعول"
@partofspeech VerbalNoun AbstractDerivedNoun Symbol("VN") "Verbal noun" "مصدر"

@partofspeech Definite AbstractState Symbol("DEF") "Definite state" "معرفة"
@partofspeech Indefinite AbstractState Symbol("INDEF") "Indefinite state" "نكرة"
@partofspeech Nominative AbstractCase Symbol("NOM") "Nominative case" "مرفوع"
@partofspeech Genetive AbstractCase Symbol("GEN") "Genetive case" "مجرور"


text:
include("physics.jl")
include("graphics.jl")

function engine(contraption::Contraption;
                settings::PhysicsSettings = PhysicsSettings(),
                Δt::Float64 = 0.01,
                t_total::Float64 = 10.,
                render::Bool = false,
                logging::Bool = true,
                realtime::Bool = true)

    @assert Δt > 0
    @assert t_total > 0
    @assert render | logging

    n = length(contraption.mass)
    n_steps = Int(div(t_total, Δt))

    # Avoid making changes to initial object passed in
    contraption = deepcopy(contraption)

    # Check that contraption is initialized in bounds
    for i in 1:n
    in_bounds(s) = (settings.bounds.x[1] <= s[1] <= settings.bounds.x[2]) &
                   (settings.bounds.y[1] <= s[2] <= settings.bounds.y[2])

    @assert in_bounds(contraption.position[:, i])
    end

    # Store position and velocity evolution as an array
    # for display, debugging and evaluation
    # Axes are: time step, x vs y, point number
    if logging
        positions = Array{Float64}(undef, n_steps, 2, n)
        velocities = Array{Float64}(undef, n_steps, 2, n)
    end

    for t in 1:n_steps
        if logging
            positions[t, :, :] = contraption.position
            velocities[t, :, :] = contraption.velocity
        end

        if render
            IJulia.clear_output(true)

            scene = build(contraption, settings.bounds)
            draw(SVGJS(), scene)
            if realtime
                sleep(Δt) # FIXME: not true constant framerate
            end
        end

        # Update state
        contraption.position, contraption.velocity = stepped_dynamics(contraption, settings, Δt)
    end

    if logging
        return (position = positions, velocity = velocities)
    else
        return nothing
    end
end


text:
using BenchmarkTools, TaylorModels

SUITE = BenchmarkGroup()

include("arithmetic.jl")
include("daisy/daisy.jl")


text:
using Test
using FlightMechanics.Models
using LinearAlgebra

# Constructor
rs = RigidSolid(15.0, [1.0, 1.0, 1.0], Matrix(1.0I, 3, 3))
@test isapprox(15.0, rs.mass)
@test isapprox([1.0, 1.0, 1.0], rs.cg)
@test isapprox(Matrix(1.0I, 3, 3), rs.inertia)

# Constructor
pm = PointMass(15.0, [1.0, 1.0, 1.0])
@test isapprox(15.0, pm.mass)
@test isapprox([1.0, 1.0, 1.0], pm.cg)
@test isapprox(zeros(3, 3), pm.inertia)

# Getters
@test isapprox(rs.mass, get_mass(rs))
@test isapprox(rs.cg, get_cg(rs))
@test isapprox(rs.inertia, get_inertia(rs))

# get_inertia in other point (Steiner's theorem)
mass = 10
point0 = [0, 0, 0]
inertia = [1 0 0;
           0 2 0;
           0 0 3]

rs = RigidSolid(mass, point0, inertia)
point1 = [10, 0, 0]
inertia_p1 = get_inertia(rs, point1)
exp_inertia = [0       0     0;
               0    1000     0;
               0       0  1000] + inertia
@test isapprox(inertia_p1, exp_inertia)

# composition
m1, m2 = 10.0, 100.0
p1 = [1.0, 2.0, 3.0]
inertia = Matrix(1.0I, 3, 3)

rs1 = RigidSolid(m1, p1, inertia)
rs2 = RigidSolid(m2, p1, inertia)

rs = rs1 + rs2
@test isapprox(m1 + m2, get_mass(rs))
@test isapprox(p1, get_cg(rs))
@test isapprox(2*inertia, get_inertia(rs))

p2 = [-1.0, -2.0, -3.0]
rs2 = RigidSolid(m2, p2, inertia)

rs = rs1 + rs2
@test isapprox(m1 + m2, get_mass(rs))
@test isapprox((p1*m1 + p2*m2) / (m1 + m2), get_cg(rs))
@test isapprox(get_inertia(rs1, get_cg(rs)) + get_inertia(rs2, get_cg(rs)),
               get_inertia(rs)
               )


text:
"""
Default implementation of set_initial_condition_value
"""
function set_ic_quantity!(
    ic::InitialCondition{T, PJ.ParameterRef},
    var_value::Float64,
) where {T <: InitialConditionType}
    @assert isfinite(var_value) ic
    PJ.set_value(ic.value, var_value)
    return
end

"""
Default implementation of set_initial_condition_value
"""
function set_ic_quantity!(
    ic::InitialCondition{T, Float64},
    var_value::Float64,
) where {T <: InitialConditionType}
    @assert isfinite(var_value) ic
    @debug "Initial condition value set with Float64. Won't update the model until rebuild" _group =
        LOG_GROUP_BUILD_INITIAL_CONDITIONS
    ic.value = var_value
    return
end


text:
##################################################
# Setup
##################################################

bnd = MCBN.BayesNetDAI(5)
bnd2 = MCBN.BayesNetDAI(5)
bnd3 = MCBN.BayesNetDAI(5)

##################################################
# Entropy and KLD Tests
##################################################

@test_approx_eq MCBN.entropy(bnd) 5.0
@test_approx_eq MCBN.naive_entropy(bnd) 5.0
@test_approx_eq MCBN.kld(bnd,bnd2) 0.0

function test_net(bnd)
    MCBN.check_bnd(bnd)
    @test_approx_eq MCBN.naive_entropy(bnd) MCBN.entropy(bnd)
    @test_approx_eq MCBN.kld(bnd,bnd) 0.0
end

function test_kld(bnd1, bnd2)
    kldval = MCBN.kld(bnd1, bnd2)
    kldvalrev = MCBN.kld(bnd2, bnd1)
    @test kldval >= 0
    @test kldvalrev >= 0
end

test_net(bnd)
test_kld(bnd, bnd2)

MCBN.add_edge!(bnd3, 1, 2)
#MCBN.check_bnd(bnd3)

MCBN.set_factor!(bnd3, 2, [1.0, 0.5, 0.0, 0.5], false)
MCBN.set_factor!(bnd3, 1, [0.8, 0.2])
MCBN.check_bnd(bnd3)

@test_approx_eq MCBN.entropy(bnd3) 3.9219280948873623
@test_approx_eq MCBN.entropy(bnd3) MCBN.naive_entropy(bnd3)

test_kld(bnd, bnd3)

MCBN.restoreFactors!(bnd3.fg)
bnd3.dirty = true

@test_approx_eq MCBN.entropy(bnd3) MCBN.naive_entropy(bnd3)
@test_approx_eq MCBN.entropy(bnd3) 5.0

@test_approx_eq MCBN.kld(bnd,bnd3) 0.0
@test_approx_eq MCBN.kld(bnd3,bnd) 0.0

test_net(bnd3)
test_kld(bnd, bnd3)

##################################################
# Entropy and KLD Tests
##################################################

bns = MCBN.BayesNetSampler(3, rand(1:2, 5,3))
origE = MCBN.energy(bns)
origfvalue = copy(bns.fvalue)

@time for i=1:1000
    scheme = MCBN.propose!(bns)
    MCBN.check_bns(bns)
    MCBN.reject!(bns)
    MCBN.check_bns(bns)
    E = MCBN.energy(bns)
    @test_approx_eq bns.fvalue origfvalue
    @test_approx_eq E origE
end

@time for i=1:1000
    s = MCBN.propose!(bns)
    MCBN.check_bns(bns)
    E1 = MCBN.energy(bns)
    f1 = copy(bns.fvalue)
    bns.changelist = [1:3;]
    E2 = MCBN.energy(bns)
    f2 = copy(bns.fvalue)
    MCBN.save!(bns)
    @test_approx_eq f1 f2
    @test_approx_eq E1 E2
end


text:
## sgetypes.jl.  Iterators for use in conjunction with Sun Grid Engine
## (c) 2014 David A. van Leeuwen
##

## The inner constructor returns the iterator itself if not running in
## an SGE environment.  I don't know if that is the Julian way...
type SGEArray{T}
    first::Int
    step::Int
    iterator::T
    state::Any
    next::Any
    function SGEArray(it::T)
        if !issgearray()
            new(0, 1, it)
        else
            first = parse(Int, ENV["SGE_TASK_ID"]) - 1
            step = parse(Int, ENV["SGE_TASK_LAST"])
            if step < first
                error("Incorrect SGE array parameters")
            end
            new(first, step, it)
        end
    end
end
## There is something about inner constructors, they need this:
SGEArray{T}(it::T) = SGEArray{T}(it)


text:
using IcanProj
using NetCDF
using DataFrames
using ProgressMeter
using Statistics
using CSV

function forcings_table(path)

    df_all = CSV.File(joinpath(dirname(pathof(IcanProj)), "..", "data", "df_links.csv")) |> DataFrame

    variables = ["ilwr", "iswr", "pres", "rainf", "rhum", "snowf", "tair", "wind"]

    for variable in variables

        @show variable

        file = joinpath(path, "$(variable)_1km.nc")

        idname = Symbol(ncgetatt(file, "dim_space", "id"))

        dataname = Symbol("$(variable)_mean")

        ids = convert(Array{Int64}, ncread(file, "id"))

        ntime = length(ncread(file, "dim_time"))

        nspace = length(ncread(file, "dim_space"))

        dataave = fill(0.0, nspace)

        for itime in 1:ntime

            dataave += ncread(file, variable, start = [1,itime], count=[-1,1])/ntime

        end

        df_nc = DataFrame()

        df_nc[idname] = ids

        df_nc[dataname] = dataave[:]

        df_all = join(df_all, df_nc, on = idname)

    end

    return df_all

end

path = "/data04/jmg/fsm_simulations/netcdf/forcings_st"

df_all = forcings_table(path)

file = joinpath(dirname(pathof(IcanProj)), "..", "data", "forcings_summary.txt")

df_all |> CSV.write(file)


text:
export Consensus

struct Consensus{T <: Real}
    index::T
end

Consensus() = Consensus(0.0)

using Statistics
function project!(s::Consensus{<:Real}, v::Matrix{T}, y::Matrix{T}) where {T <: Real}
  v = repeat(Statistics.mean(y, dims = 2), inner = [1, size(y, 2)])
end

function project(s::Consensus{<:Real}, y::Matrix{<:Real})
    project!(s, similar(y), y)
end


text:
# The "iris" dataset contains 150 measurements on 3 species
# of Iris flowers, namely I. setosa, I. versicolor and I. virginica.
# The columns 1 to 4 respectively contain the:
# sepal length - sepal width -petal length - petal width.
# The species is specified in the last column.
# The dataset can be found in the file iris.txt

# 1- extracting the length measurements from the dataset in the array:
lengths = Array(Float64, 1, 4) # create an empty (initialized with arbitrary values) array
# with 1 row, 4 columns
fname = "iris.txt"
open(fname) do file
    for line in eachline(file)
        # The columns 1 to 16 respectively contain the lengths of 1 measurement:
        mlengths = line[1:16] # measurement lengths
        arr_temp = split(mlengths) # split on spaces to an array of strings
        # convert to an array of Float64 values (a column vector 4 x 1)
        arr_mlength = [parse(Float64, str) for str in arr_temp]
        # transpose arr_mlength to obtain a 1 x 4 array with the ' operator
        # and append it vertically to lengths with vcat
        global lengths = vcat(lengths, arr_mlength')
        # we must indicate that we mean the variable lengths
        # defined in global scope by prepending it with global
    end
end
lengths = lengths[2:end, :] # cut off the first row which contained arbitrary values
# println(lengths) #
show(lengths) #
println()
show(typeof(lengths)) # Array{Float64,2}
println()
show(size(lengths, 1)) # 150 rows
println()
show(size(lengths, 2)) # 4 columns
println()

# 2- calculations on the array lengths containing all measurements:
sum(lengths) # 2078.7
sum(lengths, 1) # sum of each column
# 1x4 Array{Float64,2}:
#  876.5  458.6  563.7  179.9
sum(lengths, 2) # sum of each row
# 150x1 Array{Float64,2}:
#  10.2
#   9.5
#   9.4
#   9.4
#   ...

lengths_setosa = lengths[1:50, :]
# 50x4 Array{Float64,2}:
#  5.1  3.5  1.4  0.2
#  4.9  3.0  1.4  0.2
#  4.7  3.2  1.3  0.2
sum(lengths_setosa, 1)
# 1x4 Array{Float64,2}:
#  250.3  171.4  73.1  12.3

mean(lengths_setosa, 1)
# 1x4 Array{Float64,2}:
#  5.006  3.428  1.462  0.246

mean(lengths, 1)
# 1x4 Array{Float64,2}:
#  5.84333  3.05733  3.758  1.19933

# smallest value:
sort(lengths_setosa, 1)[1, :]
# 1x4 Array{Float64,2}:
#  4.3  2.3  1.0  0.1

# largest value:
sort(lengths_setosa, 1)[end, :]
# 1x4 Array{Float64,2}:
#  5.8  4.4  1.9  0.6

sortrows(lengths_setosa) # sorts array by the first element in each row
# 50x4 Array{Float64,2}:
#  4.3  3.0  1.1  0.1
#  4.4  2.9  1.4  0.2

max_stars_repo_name: lytemar/Julia-for-Data-Science-Video

max_issues_repo_name: lytemar/Julia-for-Data-Science-Video

max_forks_repo_name: lytemar/Julia-for-Data-Science-Video


text:
export Client,
    me,
    enable_cache!,
    disable_cache!,
    add_handler!,
    delete_handler!,
    DEFAULT_HANDLER_TAG

"""
Tag assigned to default handlers, which you can use to delete them.
"""
const DEFAULT_HANDLER_TAG = :DJL_DEFAULT

# Messages are created regularly, and lose relevance quickly.
const DEFAULT_TTLS = TTLDict(
    Guild          => nothing,
    DiscordChannel => nothing,
    User           => nothing,
    Member         => nothing,
    Presence       => nothing,
    Message        => Hour(6),
)

mutable struct Handler
    f::Function
    expiry::Union{Int, DateTime, Nothing}
end
Handler(f::Function, expiry::Period) = Handler(f, now(UTC) + expiry)

function isexpired(h::Handler)
    return if h.expiry === nothing
        false
    elseif h.expiry isa Int
        h.expiry <= 0
    else
        now(UTC) > h.expiry
    end
end

struct Conn
    io
    v::Int
end

"""
    Client(
        token::String;
        presence::Union{Dict, NamedTuple}=Dict(),
        ttls::$TTLDict=Dict(),
        version::Int=$API_VERSION,
    ) -> Client

A Discord bot. `Client`s can connect to the gateway, respond to events, and make REST API
calls to perform actions such as sending/deleting messages, kicking/banning users, etc.

# Keywords
- `presence::Union{Dict, NamedTuple}=Dict()`: Client's presence set upon connection.
  The schema [here](https://discordapp.com/developers/docs/topics/gateway#update-status-gateway-status-update-structure)
  must be followed.
- `ttls::$TTLDict=Dict()`: Cache lifetime overrides. Values of `nothing` indicate no
  expiry. Keys can be any of the following: [`Guild`](@ref), [`DiscordChannel`](@ref),
  [`Message`](@ref), [`User`](@ref), [`Member`](@ref), or [`Presence`](@ref). For most
  workloads, the defaults are sufficient.
- `version::Int=$API_VERSION`: Version of the Discord API to use. Using anything but
  $API_VERSION is not officially supported by the Discord.jl developers.
"""
mutable struct Client
    token::String       # Bot token, always with a leading "Bot ".
    hb_interval::Int    # Milliseconds between heartbeats.
    hb_seq::Union{Int, Nothing}  # Sequence value sent by Discord for resuming.
    last_hb::DateTime   # Last heartbeat send.
    last_ack::DateTime  # Last heartbeat ack.
    ttls::TTLDict       # Cache lifetimes.
    version::Int        # Discord API version.
    state::State        # Client state, cached data, etc.
    shards::Int         # Number of shards in use.
    shard::Int          # Client's shard index.
    limiter::Limiter    # Rate limiter.
    handlers::Dict{Type{<:AbstractEvent}, Dict{Symbol, Handler}}  # Event handlers.
    ready::Bool         # Client is connected and authenticated.
    use_cache::Bool     # Whether or not to use the cache for REST ops.
    conn::Conn          # WebSocket connection.

    function Client(
        token::String;
        presence::Union{Dict, NamedTuple}=Dict(),
        ttls::TTLDict=TTLDict(),
        version::Int=API_VERSION,
    )
        token = startswith(token, "Bot ") ? token : "Bot $token"
        ttls = merge(DEFAULT_TTLS, ttls)
        state = State(presence, ttls)

        c = new(
            token,        # token
            0,            # hb_interval
            nothing,      # hb_seq
            DateTime(0),  # last_hb
            DateTime(0),  # last_ack
            ttls,         # ttls
            version,      # version
            state,        # state
            nprocs(),     # shards
            myid() - 1,   # shard
            Limiter(),    # limiter
            Dict(),       # handlers
            false,        # ready
            true,         # use_cache
            # conn left undef, it gets assigned in open.
        )

        add_handler!(c, Defaults; tag=DEFAULT_HANDLER_TAG)
        return c
    end
end

function Base.show(io::IO, c::Client)
    print(io, "Discord.Client(shard=$(c.shard + 1)/$(c.shards), api=$(c.version), ")
    isopen(c) || print(io, "not ")
    print(io, "logged in)")
end

"""
    me(c::Client) -> Union{User, Nothing}

Get the [`Client`](@ref)'s bot user.
"""
me(c::Client) = c.state.user

"""
    enable_cache!(c::Client)
    enable_cache!(f::Function c::Client)

Enable the cache. `do` syntax is also supported.
"""
enable_cache!(c::Client) = c.use_cache = true
enable_cache!(f::Function, c::Client) = set_cache(f, c, true)

"""
    disable_cache!(c::Client)
    disable_cache!(f::Function, c::Client)

Disable the cache. `do` syntax is also supported.
"""
disable_cache!(c::Client) = c.use_cache = false
disable_cache!(f::Function, c::Client) = set_cache(f, c, false)

"""
    add_handler!(
        c::Client,
        T::Type{<:AbstractEvent},
        func::Function;
        tag::Symbol=gensym(),
        expiry::Union{Int, Period, Nothing}=nothing,
    )
    add_handler!(
        func::Function;
        c::Client,
        T::Type{<:AbstractEvent},
        tag::Symbol=gensym(),
        expiry::Union{Int, Period, Nothing}=nothing,
    )

Add an event handler. The handler should be a function which takes two arguments: A
[`Client`](@ref) and an [`AbstractEvent`](@ref) (or a subtype). The handler is appended to
the event's current handlers. You can also define a single handler for multiple event types
by using a `Union`. `do` syntax is also accepted.

# Keywords
- `tag::Symbol=gensym()`: A label for the handler, which can be used to remove it with
  [`delete_handler!`](@ref).
- `expiry::Union{Int, Period, Nothing}=nothing`: The handler's expiry. If an `Int` is given,
  the handler will run that many times before expiring. If a `Period` is given, the handler
  will expire after it elapsed. The default of `nothing` indicates no expiry.

!!! note
    There is no guarantee on the order in which handlers run, except that catch-all
    ([`AbstractEvent`](@ref)) handlers run before specific ones.
"""
function add_handler!(
    c::Client,
    T::Type{<:AbstractEvent},
    func::Function;
    tag::Symbol=gensym(),
    expiry::Union{Int, Period, Nothing}=nothing,
)
    if T isa Union
        add_handler!(c, T.a, func; tag=tag, expiry=expiry)
        add_handler!(c, T.b, func; tag=tag, expiry=expiry)
        return
    end

    if !hasmethod(func, (Client, T))
        throw(ArgumentError("Handler function must accept (::Client, ::$T)"))
    end

    h = Handler(func, expiry)
    if isexpired(h)
        throw(ArgumentError("Can't add a handler that will never run"))
    end

    if haskey(c.handlers, T)
        c.handlers[T][tag] = h
    else
        c.handlers[T] = Dict(tag => h)
    end
end

function add_handler!(
    func::Function,
    c::Client,
    T::Type{<:AbstractEvent};
    tag::Symbol=gensym(),
    expiry::Union{Int, Period, Nothing}=nothing,
)
    return add_handler!(c, T, func; tag=tag, expiry=expiry)
end

"""
    add_handler!(
        c::Client,
        m::Module;
        tag::Symbol=gensym(),
        expiry::Union{Int, Period, Nothing}=nothing,
    )

Add all of the event handlers defined in a module. Any function you wish to use as a
handler must be exported. Only functions with correct type signatures (see above) are used.

!!! note
    If you specify a `tag` and/or `expiry`, it's applied to all of the handlers in the
    module. That means if you add two handlers for the same event type, one of them will be
    immediately overwritten.
"""
function add_handler!(
    c::Client,
    m::Module;
    tag::Symbol=gensym(),
    expiry::Union{Int, Period, Nothing}=nothing,
)
    for f in filter(f -> f isa Function, map(n -> getfield(m, n), names(m)))
        for m in methods(f).ms
            ts = m.sig.types[2:end]
            length(m.sig.types) == 3 || continue
            if m.sig.types[2] === Client && m.sig.types[3] <: AbstractEvent
                add_handler!(c, m.sig.types[3], f; tag=tag, expiry=expiry)
            end
        end
    end
end

"""
    delete_handler!(c::Client, T::Type{<:AbstractEvent})
    delete_handler!(c::Client, T::Type{<:AbstractEvent}, tag::Symbol)

Delete event handlers. If no `tag` is supplied, all handlers for the event are deleted.
Using the tagless method is generally not recommended because it also clears default
handlers which maintain the client state. If you do want to delete a default handler, use
[`DEFAULT_HANDLER_TAG`](@ref).
"""
delete_handler!(c::Client, T::Type{<:AbstractEvent}) = delete!(c.handlers, T)

function delete_handler!(c::Client, T::Type{<:AbstractEvent}, tag::Symbol)
    delete!(get(c.handlers, T, Dict()), tag)
end

function handlers(c::Client, T::Type{<:AbstractEvent})
    return collect(filter!(p -> !isexpired(p.second), get(c.handlers, T, Dict())))
end

function allhandlers(c::Client, T::Type{<:AbstractEvent})
    catchalls = T === AbstractEvent ? Handler[] : handlers(c, AbstractEvent)
    specifics = handlers(c, T)
    fallbacks = T === FallbackEvent ? Handler[] : handlers(c, FallbackEvent)

    return if isempty(catchalls) && isempty(specifics)
        fallbacks
    elseif isempty(catchalls) && all(h -> h.first === DEFAULT_HANDLER_TAG, specifics)
        [specifics; fallbacks]
    else
        [catchalls; specifics]
    end
end

function hasdefault(c::Client, T::Type{<:AbstractEvent})
    return haskey(get(c.handlers, T, Dict()), DEFAULT_HANDLER_TAG)
end

@enum LogLevel DEBUG INFO WARN ERROR

function logmsg(c::Client, level::LogLevel, msg::AbstractString; kwargs...)
    msg = c.shards > 1 ? "[Shard $(c.shard)] $msg" : msg
    msg = "$(now()) $msg"

    if level === DEBUG
        @debug msg kwargs...
    elseif level === INFO
        @info msg kwargs...
    elseif level === WARN
        @warn msg kwargs...
    elseif level === ERROR
        @error msg kwargs...
    else
        error("Unknown log level $level")
    end
end

function Base.tryparse(c::Client, T::Type, data)
    return try
        T <: Vector ? eltype(T).(data) : T(data), nothing
    catch e
        logmsg(c, ERROR, catchmsg(e))
        push!(c.state.errors, data)
        nothing, e
    end
end

function set_cache(f::Function, c::Client, use_cache::Bool)
    old = c.use_cache
    c.use_cache = use_cache
    try
        f()
    finally
        # Usually the above function is going to be calling REST endpoints. The cache flag
        # is checked asynchronously, so by the time it happens there's a good chance we've
        # already returned and set the cache flag back to its original value.
        sleep(Millisecond(1))
        c.use_cache = old
    end
end


text:
"""
    LegendreTransform(; modes)

Constructs a legendre transform with modes `modes` that
operates on the 2nd, 3rd, ...N-modes+1-st dimension of the data.

# Example
```
julia> LegendreTransform(modes = (12, 3, 4, 12))
```
"""
Base.@kwdef struct LegendreTransform{N, T} <: AbstractTransform
    modes::NTuple{N, T}
end

function forward(::LegendreTransform, x)
    L = legendre_transform_forward_matrix.(size(x)[2:(end - 1)])
    return apply_legendre(L, x)
end

function inverse(::LegendreTransform, x)
    L⁻¹ = legendre_transform_inverse_matrix.(size(x)[2:(end - 1)])
    return apply_legendre(L⁻¹, x)
end

function truncate_modes(tr::LegendreTransform{N}, coeff) where {N}
    # Returns a low-pass filtered version of coeff assuming
    # that coeff is a tensor of spectral weights.
    # Ex.: tr.modes = (3,)
    # [0, 1, 2, 3, 4, 5] -> [0, 1, 2]
    # [a, b, c, d, e, f] -> [a, b, c]

    # indices for the spectral coefficients that we need to retain
    inds = [collect(1:m) for (s, m) in zip(size(coeff)[2:(end - 1)], tr.modes)]
    coeff_truncated = OperatorFlux.mview(coeff, inds, Val(N))

    return coeff_truncated
end

function pad_modes(
    ::LegendreTransform,
    coeff,
    size_pad::NTuple{N, T},
) where {N, T}
    # return a padded-with-zeros version of coeff assuming
    # that coeff is a tensor of spectral weights, thereby inflating coeff.
    size_space = size(coeff)[2:(end - 1)] # sizes of space-like dimensions of coeff

    # generate a zero array and indices that point to the filled-in
    # locations
    coeff_padded =
        zeros(eltype(coeff), (size(coeff)[1], size_pad..., size(coeff)[end]))
    inds = [collect(1:m) for (s, m) in zip(size_pad, size(coeff)[2:(end - 1)])]

    coeff_padded_view = OperatorFlux.mview(coeff_padded, inds, Val(N))
    coeff_padded_view .= coeff

    return coeff_padded
end

# utils
function vandermonde(x, N)
    # create view to assign values
    P = zeros(length(x), N + 1)
    P⁰ = view(P, :, 0 + 1)
    @. P⁰ = 1

    # explicitly compute second coefficient
    if N == 0
        return P
    end

    P¹ = view(P, :, 1 + 1)
    @. P¹ = x

    if N == 1
        return P
    end

    for n in 1:(N - 1)
        # get views for ith, i-1th, and i-2th columns
        Pⁿ⁺¹ = view(P, :, n + 1 + 1)
        Pⁿ = view(P, :, n + 0 + 1)
        Pⁿ⁻¹ = view(P, :, n - 1 + 1)

        # compute coefficients for ith column
        @. Pⁿ⁺¹ = ((2n + 1) * x * Pⁿ - n * Pⁿ⁻¹) / (n + 1)
    end

    return P
end

function legendre_transform_inverse_matrix(N)
    # get the legendre points to construct transformation
    # matrix. For small transforms this is performant.
    x, _ = GaussQuadrature.legendre(N, GaussQuadrature.both)
    return vandermonde(x, N - 1)
end

function legendre_transform_forward_matrix(N)
    return legendre_transform_inverse_matrix(N) \ I
end

function apply_legendre(w::Tuple{S}, x) where {S}
    F1, = w
    @tullio c[s, i, b] := F1[i, ii] * x[s, ii, b]
end

function apply_legendre(w::Tuple{S, T}, x) where {S, T}
    F1, F2 = w
    @tullio c[s, i, j, b] := F1[i, ii] * F2[j, jj] * x[s, ii, jj, b]
end

function apply_legendre(w::Tuple{S, T, V}, x) where {S, T, V}
    F1, F2, F3 = w
    @tullio c[s, i, j, k, b] :=
        F1[i, ii] * F2[j, jj] * F3[k, kk] * x[s, ii, jj, kk, b]
end

# Base extensions
Base.ndims(tr::LegendreTransform) = length(tr.modes)
Base.eltype(::LegendreTransform) = Float32
Base.size(tr::LegendreTransform) = tr.modes

function Base.show(io::IO, tr::LegendreTransform)
    print(io, "LegendreTransform(modes = $(tr.modes))")
end


text:
using Documenter
push!(LOAD_PATH,"../src/")
using ClusteredLowRankSolver
using DocumenterCitations

bib = CitationBibliography(joinpath(@__DIR__, "references.bib"), sorting=:nyt)

makedocs(bib,
    sitename = "ClusteredLowRankSolver.jl Documentation",
    format = Documenter.HTML(),
    modules = [ClusteredLowRankSolver],
    pages = [
        "Home" => "index.md",
        "Manual" => [
            "Interface" => "manual/interface.md",
            "manual/solver.md",
        ],
        "Examples" => [
            "Delsarte LP bound" =>"examples/delsarte.md",
            "Sphere packing" => "examples/sphere_packing.md",
            "Symmetric polynomial optimization" =>"examples/poly_opt.md",
        ],
        "references.md",
    ]
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
deploydocs(
    repo = "github.com/nanleij/ClusteredLowRankSolver.jl.git",
    devbranch="main"
)


text:
module Evaluate

import Base.start, Base.next, Base.done, Base.strip
export parseit, exec_cmd

# A special module in which a documents code is executed.
module WeaveSandbox
    # Replace OUTPUT_STREAM references so we can capture output.
    OUTPUT_STREAM = IOBuffer()
    print(x) = Base.print(OUTPUT_STREAM, x)
    println(x) = Base.println(OUTPUT_STREAM, x)
end

# An iterator for the parse function: parsit(source) will iterate over the
# (cmd, expressions) in a string.
type ParseIt
    value::String
end

parseit(value::String) = ParseIt(value)
start(it::ParseIt) = 1
function next(it::ParseIt, pos)
    (ex,newpos) = parse(it.value, pos)
    ((it.value[pos:(newpos-1)], ex), newpos)
end
done(it::ParseIt, pos) = pos > length(it.value)

## Execute a julia command returning (result, output_stream) tuple
function exec_cmd(expr)
    result = eval(WeaveSandbox, expr)

    seek(WeaveSandbox.OUTPUT_STREAM, 0)
    output = takebuf_array(WeaveSandbox.OUTPUT_STREAM)
    truncate(WeaveSandbox.OUTPUT_STREAM, 0)
    (result, output)
end

end