text:
export conductor, isconductor, norm_group, maximal_abelian_subfield, genus_field, content_ideal, subfields, isnormal, iscentral, normal_closure

########################################################################################
#
#  Tools for conductor
#
########################################################################################

function _norm_group_gens_small(C::ClassField)

  mp=C.mq

  mR = C.rayclassgroupmap
  mS = C.quotientmap

  R=domain(mR)
  cond, inf_plc1 = defining_modulus(mR)
  O = order(cond)
  E=order(domain(mp))
  expo=Int(exponent(domain(mp)))
  K=O.nf

  mS=pseudo_inv(mS)
  dom=domain(mS)
  M=zero_matrix(FlintZZ,ngens(dom), ngens(codomain(mS)))
  for i=1:ngens(dom)
    elem=mS(dom[i]).coeff
    for j=1:ngens(codomain(mS))
      M[i,j]=elem[1,j]
    end
  end
  S1=Hecke.GrpAbFinGenMap(domain(mS),codomain(mS),M)
  T,mT=Hecke.kernel(S1)

  Sgens=find_gens_sub(mR,mT)

  return Sgens

end

#
#  Find small primes generating a subgroup of the ray class group
#

function find_gens_sub(mR::MapRayClassGrp, mT::GrpAbFinGenMap)

  O = order(codomain(mR))
  R = domain(mR)
  T = domain(mT)
  m = Hecke._modulus(mR)
  l = minimum(m)
  lp = NfOrdIdl[]
  sR = GrpAbFinGenElem[]

  if isdefined(mR, :prime_ideal_cache)
    S = mR.prime_ideal_cache
  else
    S = prime_ideals_up_to(O, 1000)
    mR.prime_ideal_cache = S
  end
  q, mq = quo(T, sR, false)
  for (i,P) in enumerate(S)
    if divisible(l,P.minimum)
      continue
    end
    if haskey(mR.prime_ideal_preimage_cache, P)
      f = mR.prime_ideal_preimage_cache[P]
    else
      f = mR\P
      mR.prime_ideal_preimage_cache[P] = f
    end
    bool, pre = haspreimage(mT, f)
    if !bool
      continue
    end
    if iszero(mq(pre))
      continue
    end
    push!(sR, pre)
    push!(lp, P)
    q, mq = quo(T, sR, false)
    if order(q) == 1
      break
    end
  end
  if order(q) == 1
    return lp
  else
    error("Not enough primes")
  end
end

#
#  This functions constructs generators for 1+p^u/1+p^u+1
#

function _1pluspk_1pluspk1(K::AnticNumberField, p::NfOrdIdl, pk::NfOrdIdl, pv::NfOrdIdl, powers::Vector{Tuple{NfOrdIdl, NfOrdIdl}}, a::Union{Int, fmpz}, n::Int)

  O = maximal_order(K)
  b = basis(pk, copy = false)
  N = basis_matrix(pv, copy = false)*basis_mat_inv(pk, copy = false)
  G = abelian_group(N.num)
  S, mS = snf(G)
  #Generators
  gens = Vector{NfOrdElem}(undef, ngens(S))
  for i=1:ngens(S)
    gens[i] = one(O)
    for j = 1:ngens(G)
      mult = mod(mS.map[i,j], S.snf[end])
      if !iszero(mult)
        add!(gens[i], gens[i], mult*b[j])
      end
    end
  end
  if length(powers) > 1
    i = findfirst(x -> x[1] == p, powers)
    q = powers[i][2]
    i_without_p = prod([powers[j][2] for j = 1:length(powers) if j != i])
    alpha, beta = idempotents(q, i_without_p)
    for i in 1:length(gens)
      mul!(gens[i], gens[i], beta)
      add!(gens[i], gens[i], alpha)
    end
  end
  if mod(n,2)==0
    for i=1:length(gens)
      gens[i] = make_positive(gens[i], fmpz(a))
    end
  end
  return gens
end

#######################################################################################
#
#  Signature
#
#######################################################################################

@doc Markdown.doc"""
    signature(C::ClassField) -> Int, Int

Return the signature of the number field defined by $C$.
"""
function signature(C::ClassField)
  mR = C.rayclassgroupmap
  mS = C.quotientmap
  inf_plc = defining_modulus(mR)[2]
  K = base_field(C)
  rK, sK = signature(K)
  if isempty(inf_plc)
    r = degree(C)*rK
    s = div(degree(K)*degree(C) - r, 2)
    return r, s
  end
  if degree(K) == 1
    OK = base_ring(C)
    el = mR\(ideal(OK, 1-minimum(defining_modulus(mR)[1])))
    if iszero(mS(el))
      return degree(C), 0
    else
      return 0, div(degree(C), 2)
    end
  end
  D = mR.disc_log_inf_plc
  r = rK - length(D)
  for (P, el) in D
    if iszero(mS(el))
      r += 1
    end
  end
  r *= degree(C)
  s = div(degree(K)*degree(C) - r, 2)
  return r, s
end

function signature(C::ClassField{MapClassGrp, GrpAbFinGenMap})
  K = base_field(C)
  rK, sK = signature(K)
  r = degree(C)*rK
  s = div(degree(K)*degree(C) - r, 2)
  return r, s
end

#######################################################################################
#
#  Conductor functions
#
#######################################################################################

@doc Markdown.doc"""
    conductor(C::ClassField) -> NfOrdIdl, Vector{InfPlc}

Return the conductor of the abelian extension corresponding to $C$.
"""
function conductor(C::T) where T <:Union{ClassField, ClassField_pp}

  if isdefined(C,:conductor)
    return C.conductor
  end
  mR = C.rayclassgroupmap
  mS = C.quotientmap
  mp = pseudo_inv(mS)* mR
  G = domain(mp)
  #
  #  First, we need to find the subgroup
  #

  cond, inf_plc = defining_modulus(C)
  O = order(cond)
  if isone(cond) && isempty(inf_plc)
    return ideal(O,1), InfPlc[]
  end
  E = order(G)
  expo = Int(exponent(G))
  K = O.nf

  #
  #  Some of the factors of the modulus are unnecessary for order reasons:
  #
  L = Dict{NfOrdIdl, Int}()
  for (p, vp) in mR.fact_mod
    if !divisible(E, minimum(p, copy = false))
      if !iscoprime(E, norm(p)-1)
        L[p] = 1
      end
    else
      if !isone(vp)
        L[p] = vp
      end
    end
  end

  #Finite part of the modulus
  mult_grps = mR.groups_and_maps
  powers = mR.powers
  for i = 1:length(mult_grps)
    mG = mult_grps[i][2]
    p = powers[i][1]
    if !haskey(L, p)
      continue
    end
    v = L[p]
    if isone(v)
      tmg = mG.tame[p]
      if iszero(mS(tmg.disc_log))
        Base.delete!(L, p)
      end
    else
      k1 = v-1
      k2 = v
      gens = GrpAbFinGenElem[]
      Q = abelian_group(Int[])
      while k1 >= 1
        multg = _1pluspk_1pluspk1(K, p, p^k1, p^k2, powers, minimum(cond), expo)
        for i = 1:length(multg)
          push!(gens, preimage(mp, ideal(O, multg[i])))
        end
        Q, mQ = quo(G,gens,false)
        if order(Q) != E
          L[p] = k2
          break
        end
        k1 -= 1
        k2 -= 1
      end
      if k2 == 1 && order(Q) == E
        tmgD = mG.tame
        if haskey(tmgD, p)
          push!(gens, mS(tmgD[p].disc_log))
          Q,mQ = quo(G, gens, false)
          if order(Q) == E
            delete!(L, p)
          else
            L[p] = 1
          end
        else
          delete!(L,p)
        end
      end
    end
  end
  cond = ideal(O,1)
  C.factored_conductor = L
  for (p,vp) in L
    cond *= p^vp
  end

  #Infinite part of the modulus
  cond_inf = InfPlc[]
  if !isempty(inf_plc)
    D = mR.disc_log_inf_plc
    for (Pl, el) in D
      if !iszero(mS(el))
        push!(cond_inf, Pl)
      end
    end
  end

  return cond, cond_inf

end

###############################################################################
#
#  isconductor function
#
###############################################################################

@doc Markdown.doc"""
    isconductor(C::Hecke.ClassField, m::NfOrdIdl, inf_plc::Vector{InfPlc}=InfPlc[]; check) -> NfOrdIdl, Vector{InfPlc}

Checks if (m, inf_plc) is the conductor of the abelian extension corresponding to $C$. If `check` is `false`, it assumes that the
given modulus is a multiple of the conductor.
This is usually faster than computing the conductor.
"""
function isconductor(C::Hecke.ClassField, m::NfOrdIdl, inf_plc::Vector{InfPlc} = InfPlc[]; check::Bool=true)
  if isdefined(C, :conductor)
    real_cond = C.conductor
    return real_cond[1] == m && Set(real_cond[2]) == Set(inf_plc)
  end
  mR = C.rayclassgroupmap
  mS = C.quotientmap
  G = codomain(mS)
  mp = pseudo_inv(mS) * mR

  R = domain(mR)
  cond, inf_plc2 = defining_modulus(mR)
  O = order(cond)
  E = order(G)
  expo = Int(exponent(G))
  K = O.nf

  if check
    mS1 = pseudo_inv(mS)
    dom = domain(mS1)
    M = zero_matrix(FlintZZ,ngens(dom), ngens(codomain(mS1)))
    for i = 1:ngens(dom)
      elem = mS1(dom[i]).coeff
      for j = 1:ngens(codomain(mS1))
        M[i, j] = elem[1,j]
      end
    end
    S1=Hecke.GrpAbFinGenMap(domain(mS1), codomain(mS1), M)
    T,mT = Hecke.kernel(S1)

    Sgens = find_gens_sub(mR, mT)

    r,mr = ray_class_group(m, inf_plc, n_quo = expo)
    quot = GrpAbFinGenElem[mr\s for s in Sgens]
    s,ms = quo(r, quot, false)
    if order(s) != E
      return false
    end
  end

  #  Some of the factors of the modulus are unnecessary for order reasons:
  L = factor(m)
  for (p,vp) in L
    if !haskey(mR.fact_mod, p) || vp>mR.fact_mod[p]
      return false
    end
    if !divisible(E,minimum(p))
      if gcd(E, norm(p)-1)==1
        return false
      elseif vp>1
        return false
      end
    elseif vp==1
      return false
    end
  end

  #Infinite part of the modulus
  if isodd(E) && !isempty(inf_plc)
    return false
  end
  for pl in inf_plc
    if !(pl in inf_plc2)
      return false
    end
  end

  if !isempty(inf_plc2)
    D = mR.disc_log_inf_plc
    for (Pl, el) in D
      Q, mQ = quo(G, mS(el), false)
      if order(Q) == E
        return false
      end
    end
  end

  #Finite part of the modulus
  powers = mR.powers
  g_and_maps = mR.groups_and_maps
  fact_def_mod = mR.fact_mod
  for i = 1:length(powers)
    P = powers[i][1]
    if !haskey(L, P) || fact_def_mod[P] < L[P]
      return false
    end
    v = L[P]
    if v == 1
      mG = g_and_maps[i][2]
      tmg = mG.tame[P]
      Q, mQ = quo(G, [mS(tmg.disc_log)], false)
      if order(Q) == E
        return false
      end
    else
      multg = _1pluspk_1pluspk1(K, P, P^(v-1), P^v, powers, cond.gen_one, expo)
      gens = Vector{GrpAbFinGenElem}(undef, length(multg))
      for i = 1:length(multg)
        gens[i] = preimage(mp, ideal(O, multg[i]))
      end
      Q,mQ = quo(G, gens, false)
      if order(Q) == E
        return false
      end
    end
  end
  C.conductor = (m, inf_plc)
  return true
end

####################################################################################
#
#  Discriminant function
#
####################################################################################

function discriminant_sign(C::ClassField)
# Compute the sign using Brill's theorem
  _, s = signature(C)
  if mod(s, 2) == 0
    return 1
  else
    return -1
  end
end

function absolute_discriminant(C::ClassField)
  OK = base_ring(C)
  return discriminant_sign(C) * norm(discriminant(C))*discriminant(OK)^degree(C)
end

function discriminant(C::ClassField, ::FlintRationalField)
  return absolute_discriminant(C)
end

@doc Markdown.doc"""
    discriminant(C::ClassField) -> NfOrdIdl

Using the conductor-discriminant formula, compute the (relative) discriminant of $C$.
This does not use the defining equations.
"""
function discriminant(C::ClassField)
  if isdefined(C,:conductor)
    m = C.conductor[1]
    inf_plc = C.conductor[2]
  else
    C.conductor = conductor(C)
    m = C.conductor[1]
    inf_plc = C.conductor[2]
  end
  O = order(m)
  if isdefined(C, :relative_discriminant)
    if isempty(C.relative_discriminant)
      return ideal(O, 1)
    else
      return prod([P^v for (P, v) in C.relative_discriminant])
    end
  end

  @assert typeof(m) == NfOrdIdl

  mR = C.rayclassgroupmap
  mS = C.quotientmap
  mp = pseudo_inv(mS) * mR
  R = domain(mp)
  n = order(R)
  relative_disc = Dict{NfOrdIdl,Int}()
  lp = factor(m)

  if isprime(n)
    for (p, v) in lp
      ap = n*v - v
      relative_disc[p] = ap
      continue
    end
    C.relative_discriminant = relative_disc
    if isempty(C.relative_discriminant)
      return ideal(O, 1)
    else
      return prod([P^v for (P, v) in C.relative_discriminant])
    end
  end

  expo = Int(exponent(R))
  K = O.nf
  a = minimum(m)
  g_and_maps = mR.groups_and_maps
  powers = mR.powers
  for i = 1:length(powers)
    p = powers[i][1]
    if !haskey(lp, p)
      continue
    end
    v = lp[p]
    mG = g_and_maps[i][2]
    if isone(v)
      tmg = mG.tame[p]
      el = mS(tmg.disc_log)
      Q, mQ = quo(R, GrpAbFinGenElem[el], false)
      relative_disc[p] = n - order(Q)
      continue
    end
    s = v
    ap = v*degree(C)
    @hassert :AbExt 1 s>=2
    els = GrpAbFinGenElem[]
    for k = 2:v
      s = s-1
      pk = p^s
      pv = pk*p
      gens = _1pluspk_1pluspk1(K, p, pk, pv, powers, a, expo)
      for i=1:length(gens)
        push!(els, mp\ideal(O, gens[i]))
      end
      ap -= Int(order(quo(R, els, false)[1]))
      @hassert :AbExt 1 ap>0
    end
    if haskey(mG.tame, p)
      push!(els, mS(mG.tame[p].disc_log[1]))
    end
    ap -= Int(order(quo(R, els, false)[1]))
    @hassert :AbExt 1 ap>0
    relative_disc[p] = ap
  end
  C.relative_discriminant = relative_disc
  if isempty(relative_disc)
    return ideal(O, 1)
  else
    return prod([P^v for (P, v) in relative_disc])
  end
end

##############################################################################
#
#  Is Abelian function
#
##############################################################################

function isabelian(K::NfRel)
  k = base_field(K)
  Ok = maximal_order(k)
  d = ideal(Ok, Ok(discriminant(K.pol)))
  r, mr = ray_class_group(d, real_places(k), n_quo = degree(K))
  s, ms = norm_group(K.pol, mr, false)
  deg = divexact(order(r), order(s))
  return deg == degree(K)
end

function isabelian(K::NfRelNS)
  k = base_field(K)
  kx, _ = PolynomialRing(k, "x", cached = false)
  Ok = maximal_order(k)
  pols = [to_univariate(kx, x) for x in K.pol]
  d = ideal(Ok, Ok(discriminant(pols[1])))
  for i = 2:length(pols)
    d = lcm(d, ideal(Ok, Ok(discriminant(pols[i]))))
  end
  r, mr = ray_class_group(d, real_places(k), n_quo = degree(K))
  s, ms = norm_group(pols, mr, false, cached = false)
  deg = divexact(order(r), order(s))
  return deg == degree(K)
end

function isabelian(K::AnticNumberField)
  c = get_attribute(K, :isabelian)
  if c !== nothing
    return c
  end
  if !isnormal_easy(K)
    return false
  end
  if ismaximal_order_known(K)
    d1 = discriminant(maximal_order(K))
  else
    d = discriminant(K)
    den = denominator(d)
    if !isone(den)
      d *= den^degree(K)
    end
    d1 = numerator(d)
  end
  KQ = rationals_as_number_field()[1]
  ZKQ = maximal_order(KQ)
  r, mr = ray_class_group(ideal(ZKQ, d1), real_places(KQ), n_quo = degree(K))
  s, ms = norm_group(map_coefficients(KQ, K.pol), mr, false, cached = false)
  deg = divexact(order(r), order(s))
  if deg == degree(K)
    set_attribute!(K, :isabelian => true)
    return true
  else
    set_attribute!(K, :isabelian => false)
    return false
  end
end

################################################################################
#
#  Norm group function
#
################################################################################

@doc Markdown.doc"""
    norm_group(K::NfRel{nf_elem}, mR::Hecke.MapRayClassGrp) -> Hecke.FinGenGrpAb, Hecke.FinGenGrpAbMap

    norm_group(K::NfRelNS{nf_elem}, mR::Hecke.MapRayClassGrp) -> Hecke.FinGenGrpAb, Hecke.FinGenGrpAbMap

Computes the subgroup of the Ray Class Group $R$ given by the norm of the extension.
"""
function norm_group(K::NfRel{nf_elem}, mR::T, isabelian::Bool = true; of_closure::Bool = false) where T <: Union{MapClassGrp, MapRayClassGrp}
  base_field(K) == nf(order(codomain(mR))) || error("field has to be over the same field as the ray class group")
  return norm_group(K.pol, mR, isabelian, of_closure = of_closure)
end
function norm_group(K::NfRelNS{nf_elem}, mR::T, isabelian::Bool = true; of_closure::Bool = false) where T <: Union{MapClassGrp, MapRayClassGrp}
  base_field(K) == nf(order(codomain(mR))) || error("field has to be over the same field as the ray class group")
  kx, = PolynomialRing(base_field(K), "x", cached = false)
  return norm_group([to_univariate(kx, x) for x = K.pol], mR, isabelian, of_closure = of_closure)
end

@doc Markdown.doc"""
    norm_group(f::Nemo.PolyElem, mR::Hecke.MapRayClassGrp, isabelian::Bool = true; of_closure::Bool = false) -> Hecke.FinGenGrpAb, Hecke.FinGenGrpAbMap

    norm_group(f::Array{PolyElem{nf_elem}}, mR::Hecke.MapRayClassGrp, isabelian::Bool = true; of_closure::Bool = false) -> Hecke.FinGenGrpAb, Hecke.FinGenGrpAbMap

Computes the subgroup of the Ray Class Group $R$ given by the norm of the extension generated by a/the roots of $f$. If `isabelian` is set to true, then the code assumes the field to be
abelian, hence the algorithm stops when the quotient by the norm group has the correct order.
Even though the algorithm is probabilistic by nature, in this case the result is guaranteed.
If `of_closure` is given, then the norm group of the splitting field of the polynomial(s)
is computed.
It is the callers responsibility to ensure that the ray class group passed in is large enough.
"""
function norm_group(f::Nemo.PolyElem, mR::T, isabelian::Bool = true; of_closure::Bool = false, cached::Bool = true, check::Bool = false)  where T <: Union{MapClassGrp, MapRayClassGrp}
  return norm_group(typeof(f)[f], mR, isabelian, of_closure = of_closure, cached = cached, check = check)
end

function norm_group(l_pols::Vector{T}, mR::U, isabelian::Bool = true; of_closure::Bool = false, cached::Bool = true, check::Bool = false) where {T <: PolyElem{nf_elem}, U <: Union{MapClassGrp, MapRayClassGrp}}

  R = domain(mR)
  O = order(codomain(mR))
  K = nf(O)
  @assert all(x->base_ring(x) == K, l_pols) "Polynomials must be over the same field"
  if check
    @assert all(x -> isirreducible(x), l_pols) "Input polynomials must be irreducible"
  end
  N1 = minimum(defining_modulus(mR)[1])
  #I don't want to compute the discriminant of the polynomials.
  #The degree might be large and so difficult to compute.
  #Thus I will check for every prime if the projection has discriminant 0

  n = lcm(Int[degree(x) for x = l_pols])
  if of_closure
    #we cannot work in the quotient, it "could" be lcm(factorial(degree(x)) for x = f)
    Q, mQ = quo(R, GrpAbFinGenElem[])
  else
    Q, mQ = quo(R, n, false)
  end

  p = maximum(degree(x)+1 for x = l_pols)

  listprimes = GrpAbFinGenElem[]

  # Adding small primes until it stabilizes
  B = prod(Int[degree(x) for x in l_pols])
  max_stable = 50*n
  stable = max_stable
  denom = lcm([denominator(coeff(x, i)) for x in l_pols for i = 0:degree(x)])
  indexO = index(O)
  while true
    @vprint :ClassField 3 "Order of Q: $(order(Q))\n"
    if isabelian && order(Q) == B
      break
    end
    if !isabelian && order(Q) <= B && stable <= 0
      break
    end
    p = next_prime(p)
    @vprint :ClassField 3 "Using prime $(p)\n"
    if divisible(N1, p) || divisible(denom, p)
      continue
    end
    if divides(indexO, fmpz(p))[1]
      continue
    end
    found = false
    L = prime_decomposition(O, p, 1)
    for i = 1:length(L)
      candidate = mR\L[i][1]
      if iszero(mQ(candidate))
        continue
      end
      F, mF = ResidueFieldSmallDegree1(O, L[i][1])
      mFp = extend_easy(mF, K)
      all_deg = Vector{Int}[]
      #=
        the idea, taking 2 polys:
          f splits in d_i
          g splits in e_i
        Then, over an extensions of degree d_i an irreducible of degree e_i
        splits into factors of degree e_i/gcd(d_i, e_i) so there are gcd() many
        (but this is not used). The total degree over base field is then
        d_i * e_i/gcd() = lcm(d_i, e_i)
        This can be iterated...
              =#
      fl = true
      for x = l_pols
        g = map_coefficients(mFp, x)
        if degree(g) != degree(x) || iszero(discriminant(g))
          fl = false
          break
        end
        D = factor_distinct_deg(g)
        push!(all_deg, Int[x[1] for x in D])
      end
      if !fl
        continue
      end
      if !of_closure
        all_f = Set{Int}()
        for d = Iterators.product(all_deg...)
          push!(all_f, lcm(collect(d)))
        end
        else
        all_f = Set(lcm([x for x = all_f]))
      end
      E = gcd(collect(all_f))
      candidate = E*candidate
      if !iszero(mQ(candidate))
        push!(listprimes, candidate)
        Q, mQ = quo(R, listprimes, false)
        found = true
        stable = max_stable
      end
    end
    if !found
      stable -= 1
    end
  end

  if !of_closure
    for i=1:ngens(R)
       push!(listprimes, n*R[i])
    end
  end
  return sub(R, listprimes, cached)
end

function defining_modulus(mC::MapClassGrp)
  OK = order(codomain(mC))
  I = ideal(OK, 1)
  lp = Vector{InfPlc}()
  return I, lp
end

function norm_group(mL::NfToNfMor, mR::Union{MapRayClassGrp, MapClassGrp}, expected_index::Int = 1)

  K = domain(mL)
  L = codomain(mL)
  R = domain(mR)
  O = order(codomain(mR))
  @assert nf(O) == K
  if iscoprime(exponent(R), divexact(degree(L), degree(K)))
    return sub(R, gens(R), true)
  end

  N = minimum(defining_modulus(mR)[1])

  els = GrpAbFinGenElem[]

  #  Adding small primes until it stabilizes
  n = divexact(degree(L), degree(K))
  max_stable = 20*n
  GRH_bound = (4*log(abs(discriminant(maximal_order(L)))) + 2.5*expected_index + 5)^2
  stable = max_stable
  p = 0
  Q, mQ = quo(R, els, false)
  while true
    if order(Q) == expected_index || (order(Q) <= n && stable <= 0)
      break
    end
    p = next_prime(p)
    if p > GRH_bound
      break
    end
    if !iscoprime(N, p)
      continue
    end
    found = false
    lP = prime_decomposition(O, p)
    for (P, e) in lP
      lQ = prime_decomposition_type(mL, P)
      s = gcd(Int[x[1] for x in lQ])
      candidate = s*(mR\P)
      if !iszero(mQ(candidate))
        push!(els, candidate)
        Q, mQ = quo(R, els, false)
        stable = max_stable
        found = true
      end
    end
    if !found
      stable -= 1
    end
  end
  return sub(R, els, !false)
end

function norm_group(KK::KummerExt, mp::NfToNfMor, mR::Union{MapRayClassGrp, MapClassGrp})
  k = domain(mp)
  K = codomain(mp)
  ZK = maximal_order(K)
  zk = order(codomain(mR))
  # disc(ZK/Q) = N(disc(ZK/zk)) * disc(zk)^deg
  # we need the disc ZK/k, well a conductor.

  n = degree(KK)
  els = GrpAbFinGenElem[]
  stable = 0
  max_stable = 15*n*degree(k)
  R = domain(mR)
  expo = exponent(R)
  Q, mQ = quo(R, els, false)
  modu = lcm(minimum(defining_modulus(mR)[1]), index(ZK))
  prev = length(els)
  #S = PrimesSet(minimum(defining_modulus(mR)[1]), fmpz(-1), minimum(defining_modulus(mR)[1]), fmpz(1))
  S = PrimesSet(200, -1, exponent(KK), 1)
  for p in S
    if !iscoprime(p, modu)
      continue
    end
    lp = prime_decomposition(zk, p, 1)
    if isempty(lp)
      continue
    end
    D = Vector{Vector{gfp_poly}}(undef, length(KK.gen))
    for i = 1:length(D)
      D[i] = Vector{gfp_poly}(undef, length(KK.gen[i].fac))
    end
    first = false
    for i = 1:length(lp)
      P = lp[i][1]
      if iszero(mR\P)
        continue
      end
      lP = prime_decomposition(mp, P)
      local z::GrpAbFinGenElem
      try
        z = _canonical_frobenius_with_cache(lP[1][1], KK, first, D)
        @hassert :ClassField 1 z == canonical_frobenius(lP[1][1], KK)
        first = true
      catch e
        if !isa(e, BadPrime)
          rethrow(e)
        end
        continue
      end
      f = order(z)*divexact(degree(lP[1][1]), degree(P))
      if divisible(f, expo)
        stable += 1
        continue
      end
      el = f*(mR\P)
      if !iszero(mQ(el))
        push!(els, el)
        Q, mQ = quo(R, els, false)
      else
        stable += 1
      end
    end
    if stable >= max_stable
      break
    end
  end
  return sub(R, els, false)
end

@doc Markdown.doc"""
    maximal_abelian_subfield(::Type{ClassField}, K::AnticNumberField) -> ClassField

The maximal abelian subfield of $K$ as a class field, i.e. the norm group
is computed and the corresponding `ray_class_field` created.
"""
function maximal_abelian_subfield(::Type{ClassField}, K::AnticNumberField)
  Zx, x = PolynomialRing(FlintZZ, cached = false)
  QQ = rationals_as_number_field()[1]
  R, mR = ray_class_group(discriminant(maximal_order(K))*maximal_order(QQ), infinite_places(QQ), n_quo = degree(K))
  f = hom(QQ, K, K(1), check = false)
  N, mN = norm_group(f, mR)
  return ray_class_field(mR, quo(R, N)[2])
end

function norm_group_map(R::ClassField{S, T}, r::Vector{<:ClassField}, map = false) where {S, T}
  @assert map != false || all(x -> base_ring(R) == base_ring(x), r)
#  @assert map == false && all(x -> base_ring(R) == base_ring(x), r)

  mR = defining_modulus(R)[1]
  @assert map != false || all(x->mR+defining_modulus(x)[1] == defining_modulus(x)[1], r)

  fR = compose(pseudo_inv(R.quotientmap), R.rayclassgroupmap)
  lp, sR = find_gens(MapFromFunc(x->preimage(fR, x), IdealSet(base_ring(R)), domain(fR)),
                             PrimesSet(100, -1), minimum(mR))
  if map == false
    h = [hom(sR, [preimage(compose(pseudo_inv(x.quotientmap), x.rayclassgroupmap), p) for p = lp]) for x = r]
  else
    h = [hom(sR, [preimage(compose(pseudo_inv(x.quotientmap), x.rayclassgroupmap), map(p)) for p = lp]) for x = r]
  end
  return h
end

function norm_group_map(R::ClassField, r::ClassField, map = false)
  return norm_group_map(R, [r], map)[1]
end

@doc Markdown.doc"""
    maximal_abelian_subfield(K::NfRel{nf_elem}; of_closure::Bool = false) -> ClassField

Using a probabilistic algorithm for the norm group computation, determine the maximal
abelian subfield in $K$ over its base field. If `of_closure` is set to true, then
the algorithm is applied to the normal closure of $K$ (without computing it).
"""
function maximal_abelian_subfield(K::NfRel{nf_elem}; of_closure::Bool = false)
  zk = maximal_order(base_field(K))
  if has_attribute(K, :maximal_order)
    ZK = get_attribute(K, :maximal_order)
    d = discriminant(ZK)
    @assert order(d) == zk
    dd = d
  else
    d = ideal(zk, discriminant(K))
    dd = d.num
  end

  r1, r2 = signature(base_field(K))
  C, mC = ray_class_group(dd, infinite_places(base_field(K))[1:r1], n_quo = degree(K))
  N, iN = norm_group(K, mC, of_closure = of_closure)
  return ray_class_field(mC, quo(C, N)[2])
end

@doc Markdown.doc"""
    maximal_abelian_subfield(A::ClassField, k::AnticNumberField) -> ClassField

The maximal abelian extension of $k$ contained in $A$. $k$ must be a subfield of
the base field of $A$.
"""
function maximal_abelian_subfield(A::ClassField, k::AnticNumberField)
  K = base_field(A)
  fl, mp = issubfield(k, K)
  @assert fl
  return maximal_abelian_subfield(A, mp)
end

function maximal_abelian_subfield(A::ClassField, ::FlintRationalField)
  return maximal_abelian_subfield(A, Hecke.rationals_as_number_field()[1])
end

function factored_modulus(A::ClassField{MapRayClassGrp, T}) where T
  return A.rayclassgroupmap.fact_mod
end

function factored_modulus(A::ClassField{MapClassGrp, T}) where T
  return Dict{NfOrdIdl, Int}()
end

function factored_modulus(A::ClassField_pp{MapRayClassGrp, T}) where T
  return A.rayclassgroupmap.fact_mod
end

function factored_modulus(A::ClassField_pp{MapClassGrp, T}) where T
  return Dict{NfOrdIdl, Int}()
end

function maximal_abelian_subfield(A::ClassField, mp::NfToNfMor)
  k = domain(mp)
  K = codomain(mp)
  ZK = maximal_order(K)
  zk = maximal_order(k)
  # disc(ZK/Q) = N(disc(ZK/zk)) * disc(zk)^deg
  # we need the disc ZK/k, well a conductor.
  d = div(discriminant(ZK), discriminant(zk)^div(degree(K), degree(k)))
  deg = divexact(degree(K), degree(k))
  if isautomorphisms_known(K) && isnormal(K)
    G, mG = automorphism_group(K)
    deg = min(lcm([order(x) for x in G]), deg)
  end
  expo = Int(exponent(codomain(A.quotientmap)))
  mR1 = A.rayclassgroupmap
  mC = pseudo_inv(A.quotientmap)*mR1
  #First, I construct a suitable modulus for A/k
  f_m0 = Dict{NfOrdIdl, Int}()
  fact_mod = factored_modulus(A)
  for (P, e) in fact_mod
    p = intersect_prime(mp, P)
    if haskey(f_m0, p)
      if !iscoprime(minimum(P, copy = false), deg*expo)
        f_m0[p] += e
      end
    else
      if !iscoprime(minimum(P, copy = false), deg*expo)
        f_m0[p] = e
      else
        f_m0[p] = 1
      end
    end
  end
  lp = factor(ideal(zk, d))
  for (P, e) in lp
    if haskey(f_m0, P)
      if !iscoprime(minimum(P, copy = false), deg*expo)
        f_m0[P] += e
      end
    else
      if !iscoprime(minimum(P, copy = false), deg*expo)
        f_m0[P] = e
      else
        f_m0[P] = 1
      end
    end
  end

  #Now, I extend this modulus to K
  f_M0 = Dict{NfOrdIdl, Int}()
  for (p, v) in f_m0
    lp = prime_decomposition(mp, p)
    if iscoprime(minimum(p, copy = false), expo*deg)
      for (P, e) in lp
        f_M0[P] = 1
      end
    else
      for (P, e) in lp
        f_M0[P] = e*v
      end
    end
  end

  R, mR = Hecke.ray_class_group(ZK, f_M0, real_places(K), n_quo = expo * deg)
  r, mr = Hecke.ray_class_group(zk, f_m0, real_places(k), n_quo = expo * deg)
  lP, gS = Hecke.find_gens(mR, coprime_to = minimum(defining_modulus(mR1)[1]))
  listn = NfOrdIdl[norm(mp, x) for x in lP]
  # Create the map between R and r by taking norms
  proj = hom(gS, GrpAbFinGenElem[mr\x for x in listn])
  #compute the norm group of A in R
  proj1 = hom(gS, GrpAbFinGenElem[mC\x for x in lP])
  S, mS = kernel(proj1)
  mS1 = compose(mS, proj)
  G, mG = Hecke.cokernel(mS1)
  return ray_class_field(mr, mG)
end

@doc Markdown.doc"""
    ray_class_field(K::NfRel{nf_elem}) -> ClassField

For a (relative) abelian extension, compute an abstract representation
as a class field.
"""
function ray_class_field(K::NfRel{nf_elem})
  C = maximal_abelian_subfield(K)
  @assert degree(C) <= degree(K)
  if degree(C) != degree(K)
    error("field is not abelian")
  end
  return C
end

@doc Markdown.doc"""
    genus_field(A::ClassField, k::AnticNumberField) -> ClassField

The maximal extension contained in $A$ that is the compositum of $K$
with an abelian extension of $k$.
"""
function genus_field(A::ClassField, k::AnticNumberField)
  B = maximal_abelian_subfield(A, k)
  K = base_field(A)
  fl, mp = issubfield(k, K)
  @assert fl
  h = norm_group_map(A, B, x -> norm(mp, x))
  return ray_class_field(A.rayclassgroupmap, GrpAbFinGenMap(A.quotientmap * quo(domain(h), kernel(h)[1])[2]))
end

@doc Markdown.doc"""
    subfields(C::ClassField, d::Int) -> Vector{ClassField}

Find all subfields of $C$ of degree $d$ as class fields.
Note: this will not find all subfields over $Q$, but only the ones
sharing the same base field.
"""
function subfields(C::ClassField, d::Int)
  mR = C.rayclassgroupmap
  mQ = C.quotientmap

  return ClassField[ray_class_field(mR, GrpAbFinGenMap(mQ*x)) for x = subgroups(codomain(mQ), index = d, fun = (x,y) -> quo(x, y, false)[2])]
end

@doc Markdown.doc"""
    subfields(C::ClassField) -> Vector{ClassField}

Find all subfields of $C$ as class fields.
Note: this will not find all subfields over $Q$, but only the ones
sharing the same base field.
"""
function subfields(C::ClassField)
  mR = C.rayclassgroupmap
  mQ = C.quotientmap

  return ClassField[ray_class_field(mR, GrpAbFinGenMap(mQ*x)) for x = subgroups(codomain(mQ), fun = (x,y) -> quo(x, y, false)[2])]
end

@doc Markdown.doc"""
    normal_closure(C::ClassField) -> ClassField

For a ray class field $C$ extending a normal base field $k$, compute the
normal closure over $Q$.
"""
function normal_closure(C::ClassField)
  c = defining_modulus(C)
  k = base_field(C)
  if length(c[2]) > 0
    inf = real_places(k)
  else
    inf = InfPlc[]
  end
  aut = automorphisms(k)
  @assert length(aut) == degree(k)
  fin = lcm([induce_image(x, c[1]) for x = aut])

  D = ray_class_field(fin, inf, n_quo = Int(exponent(codomain(C.quotientmap))))
  h = norm_group_map(D, C)
  aut1 = small_generating_set(aut)
  act = Hecke.induce_action(D, aut1)
  k = kernel(h, true)[1]

  k = intersect([x(k)[1] for x = act])

  q, mq = quo(domain(h), k)
  return ray_class_field(D.rayclassgroupmap, GrpAbFinGenMap(D.quotientmap * mq))
end

function rewrite_with_conductor(C::ClassField)
  c, inf = conductor(C)
  if defining_modulus(C) == (c, inf)
    return C
  end
  E = ray_class_field(C.rayclassgroupmap)
  D = ray_class_field(c, inf, n_quo = Int(exponent(codomain(C.quotientmap))))
  h = norm_group_map(E, D)
  q, mq = quo(codomain(h), h(GrpAbFinGenMap(E.quotientmap)(kernel(GrpAbFinGenMap(C.quotientmap), true)[1])[1])[1])
  C = ray_class_field(D.rayclassgroupmap, GrpAbFinGenMap(D.quotientmap*mq))
  return C
end

function induce_action(C::ClassField, Aut::Vector{Hecke.NfToNfMor} = Hecke.NfToNfMor[])
  return induce_action(C.rayclassgroupmap, Aut, C.quotientmap)
end

@doc Markdown.doc"""
    isnormal(C::ClassField) -> Bool

For a class field $C$ defined over a normal base field $k$, decide
if $C$ is normal over $Q$.
"""
function isnormal(C::ClassField)

  K = base_field(C)
  aut = automorphisms(K)
  if length(aut) == degree(K)
    return isnormal_easy(C)
  else
    return isnormal_difficult(C)
  end

end

function isnormal_easy(C::ClassField)
  aut = automorphisms(base_field(C))
  c, inf = conductor(C)
  if any(x-> c != induce_image(x, c), aut)
    return false
  end
  s1 = Set(inf)
  if any(x -> s1 != Set(induce_image(x, y) for y = s1), aut)
    return false
  end
  C = rewrite_with_conductor(C)
  mR = C.rayclassgroupmap
  new_aut = small_generating_set(aut)
  act = induce_action(mR, new_aut)
  mk = kernel(GrpAbFinGenMap(C.quotientmap), true)[2]
  #normal iff kernel is invariant
  return isstable(act, mk)
end

function isnormal_difficult(C::ClassField)

  #First, I check that the norm group of the splitting field
  #of the base field contains C

  K = base_field(C)
  nK = degree(K)
  O = maximal_order(K)
  f = K.pol
  I = ideal(O, discriminant(O))
  r, mr = ray_class_group(I, real_places(K))
  Kt, t = PolynomialRing(K, "t", cached = false)
  g = divexact(evaluate(f, t), t - gen(K))
  G, mG = norm_group(g, mr, of_closure = true)
  k, mk = cokernel(mG)
  C1 = ray_class_field(mr, mk)
  if rem(degree(C), degree(C1))!= 0 || !issubfield(C1, C)
    return false
  end
  if degree(C1) == degree(C)
    return true
  end

  # Claus's Idea: I don't want to compute the extension, I want to test the stability of the modulus under the action of the
  # automorphisms, so only the totally split primes!
  # In other words, I need to check that given a totally split prime p, all the primes lying
  # over p are either all zero or all non-zero in the ray class field

  p = 1
  d = (discriminant(O)^degree(C1))*numerator(norm(evaluate(FacElem(discriminant(C1)))))
  ld = (ceil(fmpz, log(d)))
  n = degree(C1)*nK
  bound = (4*ld + 2*n +5)^2
  mp = pseudo_inv(C.quotientmap) * C.rayclassgroupmap
  while p < bound
    p = next_prime(p)
    if divisible(discriminant(O), p)
      continue
    end
    lp = prime_decomposition(O, p)
    if !all([q[1].splitting_type[2] == 1 for q in lp])
      continue
    end
    q = lp[1][1]
    fl = iszero(mp\q)
    for i = 2:nK
      if fl != iszero(mp\lp[i][1])
        return false
      end
    end
  end
  return true
end

@doc Markdown.doc"""
    iscentral(C::ClassField) -> Bool

For a class field $C$ defined over a normal base field $k$, decide
if $C$ is central over $Q$.
"""
function iscentral(C::ClassField)
  aut = automorphisms(base_field(C))
  c, inf = conductor(C)

  if any(x-> c != induce_image(x, c), aut)
    return false
  end
  s1 = Set(inf)
  if any(x -> s1 != Set(induce_image(x, y) for y = s1), aut)
    return false
  end
  C = rewrite_with_conductor(C)
  mR = C.rayclassgroupmap
  act = induce_action(mR, aut)
  k = kernel(GrpAbFinGenMap(C.quotientmap), true)
  #central iff action is trivial on the kernel
  g = [k[2](k[1][i]) for i = 1:ngens(k[1])]

  return all(x -> all(y -> x(y) == y, g), act)
end

function lcm(A::AbstractArray{<:NfAbsOrdIdl})
  a = first(A)
  a = ideal(order(a), 1)
  for b = A
    a = lcm(a, b)
  end
  return a
end

#TODO: should be done in Nemo/AbstractAlgebra s.w.
#      needed by ^ (the generic power in Base using square and multiply)
Base.copy(f::Generic.MPoly) = deepcopy(f)
Base.copy(f::Generic.Poly) = deepcopy(f)

@doc Markdown.doc"""
    lorenz_module(k::AnticNumberField, n::Int) -> NfOrdIdl

Finds an ideal $A$ s.th. for all positive units $e = 1 \bmod A$ we have that
$e$ is an $n$-th power. Uses Lorenz, number theory, 9.3.1.
If `containing` is set, it has to be an integral ideal. The resulting ideal will be
a multiple of this.
"""
function lorenz_module(k::AnticNumberField, n::Int; containing=false)
  lf = factor(n)
  return Base.reduce(lcm, [lorenz_module_pp(k, Int(p), l, containing = containing) for (p,l) = lf.fac])
end

#TODO: is this the right interface???
@doc Markdown.doc"""
    (::NfAbsOrdIdlSet)(m::Map, I::NfOrdIdl) -> NfOrdIdl

Given an embedding $m:k\to K$ of number fields and an ideal $I$ in $k$,
find the ideal above $I$ in $K$.
"""
function (I::NfAbsOrdIdlSet{Nemo.AnticNumberField,Nemo.nf_elem})(mp::Map, i::NfOrdIdl)
  assure_2_normal(i)
  return ideal(order(I), i.gen_one, order(I)(mp(i.gen_two.elem_in_nf)))
end

#TODO: write code (map?) to change polynomial rings other than evaluate

@doc Markdown.doc"""
    norm(m::T, a::nf_elem) where T <: Map{AnticNumberField, AnticNumberField} -> nf_elem

Given an embedding $m:k\to K$ of number fields and an element in $K$, find the norm
$N_{K/k}(a)$.
"""
function norm(m::T, a::nf_elem) where T <: Map{AnticNumberField, AnticNumberField}
  K = codomain(m)
  #= shamelessly from Trager:
           K  Then: K = Q(c) = k(c) = Q(b)(c)
           |        f(c) = 0 in Q[t]
      k    |        h(c) = 0 in k[t]. Trager: N(h) = f. eta in Q[t] s.th. m(b) = eta(c)
      |    |        h = gcd(b - eta, f)
      Q    Q  so N_K/k(a) = res(h, a)
  =#
  @assert K == parent(a)
  k = domain(m)
  kt, t = PolynomialRing(k, cached = false)
  Qt = parent(K.pol)
  h = gcd(gen(k) - evaluate(Qt(m(gen(k))), t), evaluate(K.pol, t))
  return resultant(h, mod(evaluate(Qt(a), t), h))
end

function norm(m::T, a::FacElem{nf_elem, AnticNumberField}) where T <: Map{AnticNumberField, AnticNumberField}
  K = codomain(m)
  @assert K == base_ring(a)
  k = domain(m)
  kt, t = PolynomialRing(k, cached = false)
  Qt = parent(K.pol)
  h = gcd(gen(k) - evaluate(Qt(m(gen(k))), t), evaluate(K.pol, t))
  d = Dict{nf_elem, fmpz}()
  for (e,v) = a.fac
    n = resultant(h, mod(change_base_ring(k, Qt(e), parent = kt), h))
    if haskey(d, n)
      d[n] += v
    else
      d[n] = v
    end
  end
  return FacElem(k, d)
end

#TODO: change order!!! this only works for maximal orders
function Base.intersect(I::NfAbsOrdIdl, R::NfAbsOrd)
  @assert ismaximal(R)
  if number_field(R) == number_field(order(I))
    return I
  end
  fl, m = issubfield(number_field(R), number_field(order(I)))
  @assert fl
  return minimum(m, I)
end

Base.intersect(R::NfAbsOrd, I::NfAbsOrdIdl) = intersect(I, R)

function Base.intersect(I::NfOrdFracIdl, R::NfAbsOrd)
  @assert ismaximal(R)
  n, d = integral_split(I)
  return intersect(n, R)
end

Base.intersect(R::NfAbsOrd, I::NfOrdFracIdl) = intersect(I, R)

@doc Markdown.doc"""
    content_ideal(f::PolyElem{nf_elem}, R::NfAbsOrd) -> NfAbsOrdIdl

The fractional $R$-ideal generated by the coefficients of $f$.
"""
function content_ideal(f::PolyElem{nf_elem}, R::NfAbsOrd)
  @assert number_field(R) == base_ring(f)
  i = sum(coeff(f, i)*R for i=0:degree(f) if !iszero(coeff(f, i)))
  return i
end

@doc Markdown.doc"""
    content_ideal(f::PolyElem{NfAbsOrdElem}) -> NfAbsOrdIdl

The ideal generated by the coefficients of $f$.
"""
function content_ideal(f::PolyElem{NfAbsOrdElem})
  R = base_ring(f)
  return sum(coeff(f, i)*R for i=0:degree(f) if !iszero(coeff(f, i)))
end

#TODO: check the math
# - I think in the p==2 l is too large in general
# - I probably only the p-part of c is needed
# - possibly even only the p-th cyclo field, although I really don't know
function lorenz_module_pp(k::AnticNumberField, p::Int, l::Int; containing=false)
  if p == 2
    l = max(l, lorenz_eta_level(k))
    l += 1
  end
  n = p^l
  C = cyclotomic_extension(k, n)
  Ka = C.Ka
  ZK = maximal_order(Ka)
  c, mc = class_group(Ka)
  lp = prime_decomposition(ZK, p)
  S = [P[1] for P = lp]
  s = [P[1] for P = prime_decomposition(maximal_order(k), p)]

  fc = false
  if containing != false
    @assert typeof(containing) == NfOrdIdl
    fc = factor(containing)
    s = union(s, collect(keys(fc)))
    fc = factor(parent(S[1])(C.mp[2], containing))
    S = union(S, collect(keys(fc)))
  end
  Q, mQ = quo(c, [mc\P for P = S])

  a, _ = find_gens(pseudo_inv(mc)*mQ, PrimesSet(degree(k), -1), p*numerator(discriminant(Ka)))
  S = Set(intersect_nonindex(C.mp[2], P) for P = a)
  union!(S, s)

  d = Dict{typeof(first(S)), Int}()
  for P = S
    # need x = 1 mod P^l -> x = y^n in k_P
    # Newton: x^n-1 has derivative nx^(n-1) and need l > 2*val(n, P)
    v = 2*valuation(p, P) + 1
    if containing != false
      v = max(v, valuation(containing, P))
    end
    d[P] = v
  end
  return numerator(evaluate(FacElem(d), coprime = true))
end

function lorenz_eta_level(k::AnticNumberField)
  # find max r s.th. eta_r in k, eta_(r+1) not in k
  # where eta_r = (zeta_(2^r) + 1/zeta_(2^r))
  r = 2
  x = PolynomialRing(FlintZZ, cached = false)[2]
  f = cos_minpoly(2^r, x)
  while hasroot(f, k)[1]
    r += 1
    f = cos_minpoly(2^r, x)
  end
  return r - 1
end


text:
using Pkg; for p in ("Knet",); haskey(Pkg.installed(),p) || Pkg.add(p); end

using Knet

struct Embed; w; end

Embed(vocab::Int,embed::Int)=Embed(param(embed,vocab))

(e::Embed)(x) = e.w[:,x]

struct Linear; w; b; end

Linear(input::Int, output::Int)=Linear(param(output,input), param0(output))

(l::Linear)(x) = l.w * mat(x,dims=1) .+ l.b

# Let's define a chain of layers
struct Chain
    layers
    Chain(layers...) = new(layers)
end
(c::Chain)(x) = (for l in c.layers; x = l(x); end; x)
(c::Chain)(x,y) = nll(c(x),y)
(c::Chain)(d::Knet.Data) = mean(c(x,y) for (x,y) in d)

CharLM(vocab::Int,input::Int,hidden::Int; o...) =
    Chain(Embed(vocab,input), RNN(input,hidden; h=0,c=0,o...), Linear(hidden,vocab))

# Sample from trained model

function generate(model,chars,n)
    reset!(m::Chain)=(for r in m.layers; r isa RNN && (r.c=r.h=0); end)
    function sample(y)
        p = Array(exp.(y)); r = rand()*sum(p)
        for j=1:length(p); (r -= p[j]) < 0 && return j; end
    end
    x = 1
    reset!(model)
    for i=1:n
        y = model([x])
        x = sample(y)
        print(chars[x])
    end
    println()
end

@info("Loading Shakespeare data")
include(Knet.dir("data","gutenberg.jl"))
trn,tst,shake_chars1 = shakespeare()
shake_text = String(shake_chars1[vcat(trn,tst)])

@info("Loading Shakespeare model")
isfile("shakespeare113.jld2") || download("http://people.csail.mit.edu/deniz/models/tutorial/shakespeare113.jld2","shakespeare113.jld2")
shake_model, shake_chars = Knet.load("shakespeare113.jld2","model","chars")

@info("Reading Julia files")
base = joinpath(Sys.BINDIR, Base.DATAROOTDIR, "julia")
julia_text = ""
for (root,dirs,files) in walkdir(base)
    for f in files
        global julia_text
        f[end-2:end] == ".jl" || continue
        julia_text *= read(joinpath(root,f),String)
    end
    # println((root,length(files),all(f->contains(f,".jl"),files)))
end

@info("Loading Julia model")
isfile("juliacharlm113.jld2") || download("http://people.csail.mit.edu/deniz/models/tutorial/juliacharlm113.jld2","juliacharlm113.jld2")
julia_model, julia_chars = Knet.load("juliacharlm113.jld2","model","chars")

nothing

ext: jl


text:
type View
    ptr::Ptr{Void}

    function View(ptr::Ptr{Void})
        v = new(ptr)
        # finalizer(v, destroy)
        # v
    end
end

function View()
    View(ccall((:sfView_create, libcsfml_graphics), Ptr{Void}, ()))
end

function View(rect::FloatRect)
    View(ccall((:sfView_createFromRect, libcsfml_graphics), Ptr{Void}, (FloatRect,), rect))
end

function View(center::Vector2f, size::Vector2f)
    v = View()
    set_center(v, center)
    set_size(v, size)
    v
end

function copy(view::View)
    return View(ccall((:sfView_copy, libcsfml_graphics), Ptr{Void}, (Ptr{Void},), view.ptr))
end

function destroy(view::View)
    println("destroyed view")
    ccall((:sfView_destroy, libcsfml_graphics), Void, (Ptr{Void},), view.ptr)
end

function set_center(view::View, pos::Vector2f)
    ccall((:sfView_setCenter, libcsfml_graphics), Void, (Ptr{Void}, Vector2f,), view.ptr, pos)
end

function set_size(view::View, size::Vector2f)
    ccall((:sfView_setSize, libcsfml_graphics), Void, (Ptr{Void}, Vector2f,), view.ptr, size)
end

function set_rotation(view::View, angle::Real)
    ccall((:sfView_setRotation, libcsfml_graphics), Void, (Ptr{Void}, Cfloat,), view.ptr, angle)
end

function set_viewport(view::View, rect::FloatRect)
    ccall((:sfView_setViewport, libcsfml_graphics), Void, (Ptr{Void}, FloatRect,), view.ptr, rect)
end

function reset(view::View, rect::FloatRect)
    ccall((:sfView_reset, libcsfml_graphics), Void, (Ptr{Void}, FloatRect,), view.ptr, rect)
end

function get_center(view::View)
    return ccall((:sfView_getCenter, libcsfml_graphics), Vector2f, (Ptr{Void},), view.ptr)
end

function get_size(view::View)
    return ccall((:sfView_getSize, libcsfml_graphics), Vector2f, (Ptr{Void},), view.ptr)
end

function get_rotation(view::View)
    return Real(ccall((:sfView_getRotation, libcsfml_graphics), Cfloat, (Ptr{Void},), view.ptr))
end

function get_viewport(view::View)
    return ccall((:sfView_getViewport, libcsfml_graphics), FloatRect, (Ptr{Void},), view.ptr)
end

function move(view::View, offset::Vector2f)
    ccall((:sfView_move, libcsfml_graphics), Void, (Ptr{Void}, Vector2f,), view.ptr, offset)
end

function rotate(view::View, angle::Real)
    ccall((:sfView_rotate, libcsfml_graphics), Void, (Ptr{Void}, Cfloat,), view.ptr, angle)
end

function zoom(view::View, factor::Real)
    ccall((:sfView_zoom, libcsfml_graphics), Void, (Ptr{Void}, Cfloat,), view.ptr, factor)
end


text:
@testset "Running Calculations" begin
    Random.seed!(SEED)
    @testset "Array" begin
        x = cumsum(randn(N))
        X = cumsum(randn(N, 4), dims=1)
        tmp = diffn(x)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
        @test sum(isnan.(tmp)) != N
        tmp = diffn(X)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == size(X,2)
        @test sum(isnan.(tmp)) != N
        tmp = runmean(x, cumulative=true)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
        @test sum(isnan.(tmp)) != N
        tmp = runmean(x, cumulative=false)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
        @test sum(isnan.(tmp)) != N
        tmp = runsum(x, cumulative=true)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
        @test sum(isnan.(tmp)) != N
        tmp = wilder_sum(x)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
        @test sum(isnan.(tmp)) != N
        tmp = runmad(x, cumulative=true)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
        @test sum(isnan.(tmp)) != N
        tmp = runmad(x, cumulative=false)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
        @test sum(isnan.(tmp)) != N
        tmp = runvar(x, cumulative=true)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
        @test sum(isnan.(tmp)) != N
        tmp = runvar(x, cumulative=false)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
        @test sum(isnan.(tmp)) != N
        tmp = runsd(x)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
        @test sum(isnan.(tmp)) != N
        tmp = runcov(x, x.*rand(N), cumulative=true)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
        @test sum(isnan.(tmp)) != N
        tmp = runcov(x, x.*rand(N), cumulative=false)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
        @test sum(isnan.(tmp)) != N
        tmp = runcor(x, x.*rand(N), cumulative=true)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
        @test sum(isnan.(tmp)) != N
        tmp = runcor(x, x.*rand(N), cumulative=false)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
        @test sum(isnan.(tmp)) != N
        tmp = runmin(x, cumulative=true, inclusive=true)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
        @test sum(isnan.(tmp)) != N
        tmp = runmin(x, cumulative=true, inclusive=false)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
        @test sum(isnan.(tmp)) != N
        tmp = runmin(x, cumulative=false, inclusive=true)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
        @test sum(isnan.(tmp)) != N
        tmp = runmin(x, cumulative=false, inclusive=false)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
        @test sum(isnan.(tmp)) != N
        tmp = runmax(x, cumulative=true, inclusive=true)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
        @test sum(isnan.(tmp)) != N
        tmp = runmax(x, cumulative=true, inclusive=false)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
        @test sum(isnan.(tmp)) != N
        tmp = runmax(x, cumulative=false, inclusive=true)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
        @test sum(isnan.(tmp)) != N
        tmp = runmax(x, cumulative=false, inclusive=false)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
        @test sum(isnan.(tmp)) != N
        tmp = mode(map(xi->round(xi), x))
        @test size(tmp, 1) == 1
        @test size(tmp, 2) == 1
        @test !isnan(tmp)
        tmp = runquantile(x, cumulative=true)
        @test !isnan(tmp[2]) && isnan(tmp[1])
        @test tmp[10] == quantile(x[1:10], 0.05)
        tmp = runquantile(x, cumulative=false)
        @test tmp[10] == quantile(x[1:10], 0.05)
        n = 20
        tmp = runacf(x, n=n, maxlag=15, cumulative=true)
        @test all(tmp[n:end,1] .== 1.0)
        tmp = runacf(x, n=n, maxlag=15, cumulative=false)
        @test all(tmp[n:end,1] .== 1.0)
    end
    @testset "Temporal" begin
        # running calculations
        x = TS(cumsum(randn(N)))
        X = TS(cumsum(randn(N, 4), dims=1))
        tmp = runmean(x, cumulative=true)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
        tmp = runmean(x, cumulative=false)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
        tmp = runsum(x, cumulative=true)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
        tmp = wilder_sum(x)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
        tmp = runmad(x, cumulative=true)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
        tmp = runmad(x, cumulative=false)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
        tmp = runvar(x, cumulative=true)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
        tmp = runvar(x, cumulative=false)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
        tmp = runsd(x)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
        tmp = runcov(X[:,1], X[:,4], cumulative=true)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
        tmp = runcov(X[:,1], X[:,4], cumulative=false)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
        tmp = runcor(X[:,1], X[:,4], cumulative=true)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
        tmp = runcor(X[:,1], X[:,4], cumulative=false)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
        tmp = runmin(x, cumulative=true, inclusive=true)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
        tmp = runmin(x, cumulative=true, inclusive=false)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
        tmp = runmin(x, cumulative=false, inclusive=true)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
        tmp = runmin(x, cumulative=false, inclusive=false)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
        tmp = runmax(x, cumulative=true, inclusive=true)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
        tmp = runmax(x, cumulative=true, inclusive=false)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
        tmp = runmax(x, cumulative=false, inclusive=true)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
        tmp = runmax(x, cumulative=false, inclusive=false)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
        tmp = runquantile(x, cumulative=true)
        @test !isnan(tmp.values[2,1]) && isnan(tmp.values[1,1])
        @test tmp.values[10,1] == quantile(x.values[1:10,1], 0.05)
        tmp = runquantile(x, cumulative=false)
        @test tmp.values[10,1] == quantile(x.values[1:10,1], 0.05)
        n = 20
        tmp = runacf(x, n=n, maxlag=15, cumulative=true)
        @test all(tmp.values[n:end,1] .== 1.0)
        tmp = runacf(x, n=n, maxlag=15, cumulative=false)
        @test all(tmp.values[n:end,1] .== 1.0)
    end
end


text:
# ## `items/image.jl`
#
# ### Item
#
# We first define the [`Image`](#) item. Since we need to keep
# track of the spatial bounds for projective transformations
# we add them as a field. By default, they will simply
# correspond to the image axes.

"""
    Image(image[, bounds])

Item representing an N-dimensional image with element type T.

## Examples

{cell=image}
```julia
using DataAugmentation, Images

imagedata = rand(RGB, 100, 100)
item = Image(imagedata)
showitems(item)
```

If `T` is not a color, the image will be interpreted as grayscale:

{cell=image}
```julia
imagedata = rand(Float32, 100, 100)
item = Image(imagedata)
showitems(item)
```

"""
struct Image{N,T} <: AbstractArrayItem{N,T}
    data::AbstractArray{T,N}
    bounds::Bounds{N}
end

Image(data) = Image(data, Bounds(axes(data)))

function Image(data::AbstractArray{T,N}, sz::NTuple{N,Int}) where {T,N}
    return Image(data, Bounds(sz))
end

Base.show(io::IO, item::Image{N,T}) where {N,T} =
    print(io, "Image{$N, $T}() with bounds $(item.bounds)")

function showitem!(img, image::Image{2, <:Colorant})
    showimage!(img, itemdata(image))
end

function showitem!(img, image::Image{2, <:AbstractFloat})
    return showimage!(img, colorview(Gray, itemdata(image)))
end

# ### Projective transformations
#
# To support projective transformations, we need to implement [`getbounds`](#)
# and [`project`](#).

getbounds(image::Image) = image.bounds

# For the projection, we use `warp` from
# [ImageTransformations.jl](https://github.com/JuliaImages/ImageTransformations.jl).
# We have to pass the inverse of the projection `P` as it uses backward
# mode warping.

function project(P, image::Image{N, T}, bounds::Bounds) where {N, T}
    # TODO: make interpolation scheme and boundary conditions configurable
    data_ = warp(
        itemdata(image),
        inv(P),
        bounds.rs,
        zero(T))
    return Image(data_, bounds)
end

# The inplace version `project!` is quite similar. Note `indices` are not needed
# as they are implicitly given by the buffer.

function project!(bufimage::Image, P, image::Image{N, T}, bounds::Bounds{N}) where {N, T}
    a = OffsetArray(parent(itemdata(bufimage)), bounds.rs)
    res = warp!(
        a,
        box_extrapolation(itemdata(image), zero(T)),
        inv(P),
    )
    return Image(res, bounds)
end


text:
struct Ginsu{N,JA,JB,JC,JD}
    A::JA
    B::JB
    C::JC
    D::JD
    rₒ::NTuple{N,Float64}
    δr::NTuple{N,Float64}
end

"""
    g = Ginsu(r0, dr, nr, sour, recr, padr, ndamp; dims=(:z,:y,:x), stencilhalfwidth=2, vector_width=8)

Create a Ginsu object that defines the model aperture for a given shot.  `g` is used to subset earth
models using the `sub` and `sub!` methods, and provides the inverse operation using the `super`,
`super!` and `super_accumulate!` methods.  The padding for the earth model subset is determined
from the source and receiver positions, along with the padding `padr` and damping region `ndamp`
parameters.

# required parameters[1,2]
* `r0::NTuple{N,Real}` model origin in each model dimension
* `dr::NTuple{N,Real}` model cell widths in each model dimension
* `nr::NTuple{N,Int}` model cell counts in each model dimension
* `sour::NTuple{N,Vector{Float}}` source locations in each model dimension
* `recr::NTuple{N,Vector{Float}}` receiver locations in each model dimension
* `padr::NTuple{N,NTuple{Real,Real}}` padding beyond the survey aperture in each model dimension[3,4]
* `ndamp::NTuple{N,NTuple{Int,Int}}` damping region for absorbing boundaies in each model dimension[3,4]

# named optional parameters
* `dim=(:z,:y,:x)` axis ordering with the default being `z` fast, and `x` slow
* `stencilhalfwidth=0` if there is a free-surface, set `stencilhalfwidth`, padr, and ndamp accordingly. This will add `stencilhalfwidth` cells above the free surface to allow copying the mirrored model to implement the free surface boundary condition.
* `vector_width=8` sets the width required for vectorization of the code, and ensures that the size(ginsu) divided by vector_width has remainder zero for each dimension.

# type specification
* `g.lextrng` is the logical (1-based) indices for the Ginsu model subset
* `g.lintrng` is the logical (1-based) indices for the Ginsu model subset interior
* `rₒ` is the physical origin of the Ginsu model subset
* `δr` is the grid cell sizes

# Notes
1. `N` is the number of model dimensions
2. `rₒ[i]` is the model origin along the ith model dimension
3. the model padding is the combination of `padr` and `ndamp`
4. `padr[i][1]` is the padding for the dimension start, and `padr[i][2]` for the end.  The same is true for `ndamp`.

## Example (free-surface)
```
o---------------c---x--------------x---c-------o
|               |   |              |   |       |
|               |   |              |   |       |
|               |   |              |   |       |
|               |   |              |   |       |
|               |   x--------------x   |       |
|               |                      |       |
o---------------c----------------------c-------o
```
* `o` denotes the original model
* `c` denotes the Ginsu module subset
* `x` denotes the Ginsu model subset interior
* the region falling in-between the Ginsu model subset and the Ginsu model subset interior is the absorbing boundary region.
"""
function Ginsu(
        r0::NTuple{N,Real},
        dr::NTuple{N,Real},
        nr::NTuple{N,Int64},
        sour::NTuple{N,AbstractArray{Float64,1}},
        recr::NTuple{N,AbstractArray{Float64,1}},
        padr::NTuple{N,Tuple{Real,Real}},
        ndamp::NTuple{N,Tuple{Int,Int}};
        T::DataType = Float32,
        dims::NTuple=(:z,:y,:x),
        stencilhalfwidth::Int=0,
        vector_width::Int = 8) where N
    nrec = length(recr[1])
    nsrc = length(sour[1])
    for idim = 2:N
        @assert length(recr[idim]) == nrec
        @assert length(sour[idim]) == nsrc
    end

    lextrng = Array{UnitRange{Int64}}(undef, N)
    lintrng = Array{UnitRange{Int64}}(undef, N)

    for idim = 1:N
        # source and receiver coords for this dimension:
        reci = recr[idim]
        soui = sour[idim]

        if dims[idim] == :z
            lb, ub = r0[idim] - padr[idim][1] - ndamp[idim][1]*dr[idim] - stencilhalfwidth*dr[idim], r0[idim] + (nr[idim]-1+ndamp[idim][2])*dr[idim] + padr[idim][2]
        else
            # midpoints for all sources:
            midi = Array{Float64}(undef, nsrc, nrec)
            for isrc = 1:nsrc
                midi[isrc,:] = (soui[isrc] .+ reci) / 2
            end

            # midpoint float range:
            lb, ub = minimum(midi) - padr[idim][1], maximum(midi) + padr[idim][2]

            # ensure that the range encloses the source:
            lb, ub = min(minimum(soui) - padr[idim][1], lb), max(maximum(soui) + padr[idim][2], ub)

            # ensure that the range encloses all receivers:
            lb, ub = min(minimum(reci) - padr[idim][1], lb), max(maximum(reci) + padr[idim][2], ub)

            # add damping region
            lb, ub = lb - ndamp[idim][1]*dr[idim], ub + ndamp[idim][2]*dr[idim]
        end

        # integer range:
        idx_lb, idx_ub = floor(Int, (lb - r0[idim]) / dr[idim]) + 1, ceil(Int, (ub - r0[idim]) / dr[idim]) + 1

        # ensure lengths are a scalar multiple of vector_width (for vector alignment)
        n = idx_ub - idx_lb + 1
        d,r = divrem(n,vector_width)
        if r != 0
            idx_ub = idx_lb + (d+1)*vector_width - 1
        end

        # interior/exterior ranges
        lextrng[idim] = idx_lb:idx_ub
        lintrng[idim] = (idx_lb+ndamp[idim][1]):(idx_ub-ndamp[idim][2])
    end

    A = JopPad(JetSpace(T,nr), ntuple(i->lintrng[i], N)..., extend=true)
    B = JopPad(JetSpace(T,nr), ntuple(i->lintrng[i], N)..., extend=false)
    C = JopPad(JetSpace(T,nr), ntuple(i->lextrng[i], N)..., extend=true)
    D = JopPad(JetSpace(T,nr), ntuple(i->lextrng[i], N)..., extend=false)

    Ginsu(A, B, C, D, ntuple(i->Float64(r0[i]), N), ntuple(i->Float64(dr[i]), N))
end

function Ginsu(
        r0::NTuple{N,Real},
        dr::NTuple{N,Real},
        nr::NTuple{N,Int64},
        sour::NTuple{N,AbstractArray{Array{Float64,1},1}},
        recr::NTuple{N,AbstractArray{Array{Float64,1},1}},
        padr::NTuple{N,Tuple{Real,Real}},
        ndamp::NTuple{N,Tuple{Int,Int}};
        T = Float32,
        dims::Tuple=(:z,:y,:x),
        stencilhalfwidth::Integer=0,
        vector_width::Int = 8) where N
    nshot = length(sour[1])

    ishot = 1
    g = Ginsu(r0, dr, nr, ntuple(idim->sour[idim][ishot], N), ntuple(idim->recr[idim][ishot], N), padr, ndamp, dims=dims, stencilhalfwidth=stencilhalfwidth, vector_width=vector_width)

    lextrng_beg = Int[lextents(g,interior=false)[idim][1] for idim=1:N]
    lextrng_end = Int[lextents(g,interior=false)[idim][end] for idim=1:N]
    lintrng_beg = Int[lextents(g,interior=true)[idim][1] for idim=1:N]
    lintrng_end = Int[lextents(g,interior=true)[idim][end] for idim=1:N]

    for ishot = 2:nshot
        g = Ginsu(r0, dr, nr, ntuple(idim->sour[idim][ishot], N), ntuple(idim->recr[idim][ishot], N), padr, ndamp, dims=dims, stencilhalfwidth=stencilhalfwidth, vector_width=vector_width)
        for idim = 1:N
            if lextents(g)[idim][1] < lextrng_beg[idim]
                lextrng_beg[idim] = lextents(g,interior=false)[idim][1]
                lintrng_beg[idim] = lextents(g,interior=true)[idim][1]
            end
            if lextents(g)[idim][end] > lextrng_end[idim]
                lextrng_end[idim] = lextents(g,interior=false)[idim][end]
                lintrng_end[idim] = lextents(g,interior=true)[idim][end]
            end
        end
    end

    A = JopPad(JetSpace(T,nr), ntuple(i->lintrng_beg[i]:lintrng_end[i], N)..., extend=true)
    B = JopPad(JetSpace(T,nr), ntuple(i->lintrng_beg[i]:lintrng_end[i], N)..., extend=false)
    C = JopPad(JetSpace(T,nr), ntuple(i->lextrng_beg[i]:lextrng_end[i], N)..., extend=true)
    D = JopPad(JetSpace(T,nr), ntuple(i->lextrng_beg[i]:lextrng_end[i], N)..., extend=false)

    Ginsu(A, B, C, D, r0, dr)
end

"""
    g = Ginsu(r0, dr, nr, aperturer, ndamp; T=Float32)

Create a Ginsu object from absolute aperture.

# required parameters[1,2]
* `aperturer::NTuple{Range}` aperture in each dimension
"""
function Ginsu(r0::NTuple{N,Real}, dr::NTuple{N,Real}, nr::NTuple{N,Int}, aperturer::NTuple{N,AbstractRange}, ndamp::NTuple{N,Tuple{Int,Int}}; T=Float32, stencilhalfwidth::Int=0, dims=(:z,:x,:y), vector_width::Int = 8) where N
    lextrng = Array{UnitRange{Int64}}(undef, N)
    lintrng = Array{UnitRange{Int64}}(undef, N)

    for idim = 1:N
        if dims[idim] == :z
            lb, ub = aperturer[idim][1] - (stencilhalfwidth+ndamp[idim][1])*dr[idim], aperturer[idim][end] + ndamp[idim][2]*dr[idim]
        else
            lb, ub = aperturer[idim][1] - ndamp[idim][1]*dr[idim], aperturer[idim][end] + ndamp[idim][2]*dr[idim]
        end

        # integer range:
        idx_lb, idx_ub = floor(Int, (lb - r0[idim]) / dr[idim]) + 1, ceil(Int, (ub - r0[idim]) / dr[idim]) + 1

        # ensure lengths are a scalar multiple of vector_width (for vector alignment)
        n = idx_ub - idx_lb + 1
        d,r = divrem(n,vector_width)
        if r != 0
            idx_ub = idx_lb + (d+1)*vector_width - 1
        end

        lextrng[idim] = idx_lb:idx_ub
        lintrng[idim] = (idx_lb+ndamp[idim][1]):(idx_ub-ndamp[idim][2])
    end

    A = JopPad(JetSpace(T,nr), ntuple(i->lintrng[i], N)..., extend=true)
    B = JopPad(JetSpace(T,nr), ntuple(i->lintrng[i], N)..., extend=false)
    C = JopPad(JetSpace(T,nr), ntuple(i->lextrng[i], N)..., extend=true)
    D = JopPad(JetSpace(T,nr), ntuple(i->lextrng[i], N)..., extend=false)

    Ginsu(A, B, C, D, r0, dr)
end

function _op(ginsu::Ginsu; extend=false, interior=false)
    interior  && extend  && return ginsu.A
    interior  && !extend && return ginsu.B
    !interior && extend  && return ginsu.C
    ginsu.D
end

Base.copy(ginsu::Ginsu) = Ginsu(ginsu.A, ginsu.B, ginsu.C, ginsu.D, ginsu.rₒ, ginsu.δr)

"""
    sub!(mg, ginsu, m; extend=true, interior=false)

Store the subset of `m` corresponding to `ginsu` in `mg`.
"""
sub!(prop_ginsu::AbstractArray, ginsu::Ginsu, prop::AbstractArray; extend=true, interior=false) = mul!(prop_ginsu, _op(ginsu;extend=extend,interior=interior), prop)

"""
    mg = sub(ginsu, m; extend=true, interior=false)

Store the subset of `m` corresponding to `ginsu` in `mg`.  New memory is allocated
for `mg`.  `interior=false` includes the sponge region.
"""
sub(ginsu::Ginsu, prop::AbstractArray; extend=true, interior=false) = _op(ginsu,extend=extend,interior=interior)*prop

"""
    super!(m, ginsu, mg; interior=false, accumulate=false)

Store the superset of `mg` corresponding to `ginsu` in `m`. `interior=false` includes the sponge
region.
"""
function super!(prop::AbstractArray, ginsu::Ginsu, prop_ginsu::AbstractArray; interior=false, accumulate=false)
    A = _op(ginsu, extend=false, interior=interior)
    state!(A, (accumulate=accumulate,))
    mul!(prop, A', prop_ginsu)
    prop
end

"""
    mi = interior(ginsu, mg)

Get a view of the Ginsu'd subset `mg` corresponding to its interior region which
excludes the sponge.
"""
function interior(ginsu::Ginsu{N}, prop_ginsu::AbstractArray{T,N}) where {T,N}
    rng = interior(ginsu)
    view(prop_ginsu, rng...)
end

"""
    I = interior(ginsu)

Get the index range corresponding to the interior region which excludes
the sponge.
"""
function interior(ginsu::Ginsu{N}) where {N}
    AI = _op(ginsu, extend=false, interior=true)
    AE = _op(ginsu, extend=false, interior=false)
    ntuple(i->begin
            strt = state(AI).pad[i][1] - state(AE).pad[i][1] + 1
            stop = strt + length(state(AI).pad[i]) - 1
            strt:stop
        end, Val{N}())
end

Base.size(gn::Ginsu, i::Int; interior=false) = length(state(_op(gn; interior=interior)).pad[i])
Base.size(gn::Ginsu{N}; interior=false) where {N} = ntuple(i->size(gn, i; interior=interior), Val{N}())

lextents(gn::Ginsu, i::Int; interior=false) = state(_op(gn; interior=interior)).pad[i]
lextents(gn::Ginsu{N}; interior=false) where {N} = ntuple(i->lextents(gn, i; interior=interior), Val{N}())

pextents(gn::Ginsu, i::Int; interior=false) = begin l=lextents(gn,i;interior=interior); (gn.rₒ[i] + (l[1]-1)*gn.δr[i]):(gn.δr[i]):(gn.rₒ[i] + (l[end]-1)*gn.δr[i]) end
pextents(gn::Ginsu{N}; interior=false) where {N} = ntuple(i->pextents(gn,i;interior=interior), Val{N}())

origin(gn::Ginsu, i::Int; interior=false) = gn.rₒ[i] + (first(lextents(gn,i;interior=interior))-1)*gn.δr[i]
origin(gn::Ginsu{N}; interior=false) where {N} = ntuple(i->origin(gn, i;interior=interior)[1], Val{N}())

export Ginsu, lextents, pextents, origin, interior, sub, sub!, super!

max_stars_repo_name: ChevronETC/JetPackWaveFD

max_issues_repo_name: ChevronETC/JetPackWaveFD

max_forks_repo_name: ChevronETC/JetPackWaveFD


text:
module FreeBSDPorts

include("./bindep.jl");

end # module


text:
@testset "compute_body_force_terms1" begin

    m = 10
    n = 10
    h = 1.0
    domain = example_static_domain1(m, n, h)
    fext = compute_body_force_terms1(domain, (x,y)->x^2+y^2+1)
    f = eval_f_on_gauss_pts((x,y)->x^2+y^2+1, m, n, h)
    pf_f = compute_fem_source_term1(f, m, n, h)
    @test fext ≈ pf_f

end

@testset "compute_stiffness_matrix_and_internal_force1" begin
    m = 10
    n = 10
    h = 1.0
    domain = example_static_domain1(m, n, h)
    K, fint = compute_stiffness_matrix_and_internal_force1(domain)
    pf_K = compute_fem_stiffness_matrix1(diagm(0=>ones(2)), m, n, h)
    @test pf_K≈K
end

@testset "compute_stiffness_matrix_and_internal_force1" begin
    m = 10
    n = 10
    h = 1.0
    domain = example_static_domain1(m, n, h)
    N = getNGauss(domain)
    k = zeros(N, 2, 2)
    for i = 1:N
        k[i,:,:] = diagm(0=>ones(2))
    end
    init_nnfem(domain)
    K = s_compute_stiffness_matrix1(k, domain)
    Kref, _ = compute_stiffness_matrix_and_internal_force1(domain)
    K = run(sess, K)
    @test Kref≈K
end


text:
using AbstractIterativeSolvers, LinearAlgebra
using Test

@testset "AbstractIterativeSolvers.jl: GMRES" begin
    n = 20
    A = randn(n,n)/(2*sqrt(n)) + I
    b = randn(n,1)
    out = GMRES(x -> A*x, b, (x,y) -> (y'*x)[1], 1e-10, n+1, x -> x)
    sol = sum([out[2][i]*out[1][i] for i=1:length(out[2]) ])
    @test norm(A*sol-b) < 1e-10
end


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENSE in the project root.
# ------------------------------------------------------------------

@recipe function f(partition::SpatialPartition)
  marker --> :auto
  for object in partition
    @series begin
      object
    end
  end
end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    Data

A domain together with a table of values for each element of the domain.
The i-th row of the table is a vector of features associated with the
i-th element of the domain. If the domain is a mesh, then the table
stores all the properties of the elements.
"""
abstract type Data end

"""
    domain(data)

Return underlying domain of the `data`.
"""
function domain end

"""
    values(data)

Return the values of `data` as a table.
"""
values(data::Data)

# ----------
# FALLBACKS
# ----------

# implement Domain traits for convenience
embeddim(data::Data) = embeddim(domain(data))
coordtype(data::Data) = coordtype(domain(data))
nelements(data::Data) = nelements(domain(data))
coordinates!(buff, data::Data, ind::Int) =
  coordinates!(buff, domain(data), ind)

==(data₁::Data, data₂::Data) =
  domain(data₁) == domain(data₂) && values(data₁) == values(data₂)

# -----------------
# TABLES INTERFACE
# -----------------

Tables.istable(::Type{<:Data}) = true
Tables.rowaccess(data::Data) = true
function Tables.rows(data::Data)
  rows = Tables.rows(values(data))
  elms = domain(data)
  ((row..., elms[i]) for (i, row) in Iterators.enumerate(rows))
end
function Tables.schema(data::Data)
  geomtype = eltype(domain(data))
  schema = Tables.schema(values(data))
  names, types = schema.names, schema.types
  Tables.Schema((names..., :geometry), (types..., geomtype))
end

# -------------------
# VARIABLE INTERFACE
# -------------------

"""
    getindex(data, var)

Returns the data for the variable `var` in `data` as a column vector.
"""
Base.getindex(data::Data, var::Symbol) =
  Tables.getcolumn(values(data), var)

# ----------
# UTILITIES
# ----------

"""
    asarray(data, var)

Returns the data for the variable `var` in `data` as a Julia array
with size equal to the size of the underlying domain if the size is
defined, otherwise returns a vector.
"""
function asarray(data::Data, var::Symbol)
  D = domain(data)
  hassize = hasmethod(size, (typeof(D),))
  hassize ? reshape(data[var], size(D)) : data[var]
end

# -----------
# IO METHODS
# -----------

function Base.show(io::IO, data::Data)
  name = nameof(typeof(data))
  nelm = nelements(domain(data))
  Dim  = embeddim(domain(data))
  T    = coordtype(domain(data))
  print(io, "$nelm $name{$Dim,$T}")
end

function Base.show(io::IO, ::MIME"text/plain", data::Data)
  𝒟 = domain(data)
  𝒯 = values(data)
  s = Tables.schema(𝒯)
  vars = zip(s.names, s.types)
  println(io, data)
  println(io, "  variables")
  varlines = ["    └─$var ($V)" for (var,V) in vars]
  println(io, join(sort(varlines), "\n"))
  print(  io, "  domain: ", 𝒟)
end


text:
export Geometry

module Geometry

include("point.jl")
include("mesh.jl")
include("line.jl")
include("aabb.jl")
include("rectangle.jl")
include("polygon.jl")

end

max_stars_repo_name: wdevore/Ranger-Julia-SDL

max_issues_repo_name: wdevore/Ranger-Julia-SDL

max_forks_repo_name: wdevore/Ranger-Julia-SDL


text:
# automatically generated by the FlatBuffers compiler, do not modify

# module: Example

import FlatBuffers

FlatBuffers.@STRUCT struct Vec3
    x::Float32
    y::Float32
    z::Float32
    test1::Float64
    test2::Color
    test3::Test
end
FlatBuffers.@ALIGN(Vec3, 16)

Vec3(buf::AbstractVector{UInt8}) = FlatBuffers.read(Vec3, buf)
Vec3(io::IO) = FlatBuffers.deserialize(io, Vec3)


text:
"""
    function most_frequent(f::Function, iterable)

Return the most frequent value of a given iterable.
"""
function most_frequent(f::Function, iterable)
    d = Dict{Int, Int}()
    for element ∈ iterable
        v = f(element)
        if haskey(d, v)
            d[v] += 1
        else
            d[v] = 1
        end
    end
    candidate = 0
    key = 0
    for (k, v) in d
        if v > candidate
            key = k
            candidate = v
        end
    end
    return key
end

"""
    nthbitset(n, a) = !Bool((a >> (n - 1)) & 1)

Return `true` if the n-th bit of `a` is set, `false` otherwise.
"""
nthbitset(n, a) = Bool((a >> n) & 1)

"""
$(METHODLIST)

Categorise the struct elements of a vector by a given field into a dictionary of
`T.field => Vector{T}`.
"""
function categorize end

"""
$(TYPEDSIGNATURES)

Examples
========

```
julia> using NeRCA

julia> struct PMT
         dom_id
         time
       end

julia> pmts = [PMT(2, 10.4), PMT(4, 23.5), PMT(2, 42.0)];

julia> categorize(:dom_id, pmts)
Dict{Any, Vector{PMT}} with 2 entries:
  4 => [PMT(4, 23.5)]
  2 => [PMT(2, 10.4), PMT(2, 42.0)]
```
"""
@inline function categorize(field::Symbol, elements::Vector)
    _categorize(Val{field}(), elements)
end

"""
$(TYPEDSIGNATURES)
"""
@inline function _categorize(field::Val{F}, elements::Vector{T}) where {T,F}
    out = Dict{fieldtype(T, F), Vector{T}}()
    for el ∈ elements
        key = getfield(el, F)
        if !haskey(out, key)
            out[key] = T[]
        end
        push!(out[key], el)
    end
    out
end


text:
@testset "nn" begin
    using LinearAlgebra
    k = NeuralNetworkKernel()
    v1 = rand(3); v2 = rand(3)
    @test k(v1,v2) ≈ asin(v1' * v2 / sqrt((1 + v1' * v1) * (1 + v2' * v2))) atol=1e-5

    # kernelmatrix tests
    m1 = rand(3,4)
    m2 = rand(3,4)
    @test kernelmatrix(k, m1, m1) ≈ kernelmatrix(k, m1) atol=1e-5
    @test_broken kernelmatrix(k, m1, m2) ≈ k(m1, m2) atol=1e-5

    x1 = rand()
    x2 = rand()
    @test kernelmatrix(k, x1*ones(1,1), x2*ones(1,1))[1] ≈ k(x1, x2) atol=1e-5

    @test k(v1, v2) ≈ k(v1, v2) atol=1e-5
    @test typeof(k(v1, v2)) <: Real

    @test_broken size(k(m1, m2)) == (4, 4)
    @test_broken size(k(m1)) == (4, 4)

    A1 = ones(4, 4)
    kernelmatrix!(A1, k, m1, m2)
    @test A1 ≈ kernelmatrix(k, m1, m2) atol=1e-5

    A2 = ones(4, 4)
    kernelmatrix!(A2, k, m1)
    @test A2 ≈ kernelmatrix(k, m1) atol=1e-5

    @test size(kerneldiagmatrix(k, m1)) == (4,)
    A3 = kernelmatrix(k, m1)
    @test kerneldiagmatrix(k, m1) ≈ [A3[i, i] for i in 1:LinearAlgebra.checksquare(A3)] atol=1e-5

    A4 = ones(4)
    kerneldiagmatrix!(A4, k, m1)
    @test kerneldiagmatrix(k, m1) ≈ A4 atol=1e-5

    A5 = ones(4,4)
    @test_throws AssertionError kernelmatrix!(A5, k, m1, m2, obsdim=3)
    @test_throws AssertionError kernelmatrix!(A5, k, m1, obsdim=3)
    @test_throws DimensionMismatch kernelmatrix!(A5, k, ones(4,3), ones(3,4))

    @test k([x1], [x2]) ≈ k(x1, x2) atol=1e-5
    test_ADs(NeuralNetworkKernel, ADs = [:ForwardDiff, :ReverseDiff])
    @test_broken "Zygote uncompatible with BaseKernel"
end


text:
using AbstractBijections
using Suppressor
using Test

@testset "AbstractBijections.jl" begin
  f = Bijection(["X1", "X2"], ["Y1", "Y2"])
  @test f["X1"] == "Y1"
  @test inv(f)["Y1"] == "X1"
  @test domain(f) == ["X1", "X2"]
  @test image(f) == ["Y1", "Y2"]

  f = bijection(["X1" => "Y1", "X2" => "Y2"])
  @test f["X1"] == "Y1"
  @test inv(f)["Y1"] == "X1"

  f = Bijection(["Y1", "Y2"])
  @test f[1] == "Y1"
  @test inv(f)["Y1"] == 1
end

@testset "examples.jl" begin
  examples_path = joinpath(pkgdir(AbstractBijections), "examples")
  files = readdir(examples_path; join=true)
  for file in files
    @suppress include(file)
  end
end


text:
# construct graph
const N = 5 # must be bigger than 4
F = Graph(N)
for n in 1:N - 2
    for _ in 1:n
        add_arc!(F, n, n + 1)
        add_arc!(F, n, n + 2)
    end
end

# ==
@test ArcIter(F, :, :) == arcs(F, :, :)
@test length(ArcIter(F, :, :)) == F.na
@test length(arcs(F, :, :)) == F.na

# arcs from
@test length(arcs(F, 1, :)) == 2
for i in 2:N - 2
    @test length(arcs(F, i, :)) == 2 * i
end

# arcs into
@test length(arcs(F, :, 1)) == 0
for i in 2:N - 2
    @test length(arcs(F, :, i)) == 2 * i - 3
end

# getindex and setindex!
for a in arcs(F)
    w = rand(1:5)
    F[a,:flow] = w
    @test F[a,:flow] == w
end

# random
test_ait = arcs(F, 1, :)
rand_test_ait = rand(test_ait)
@test rand_test_ait == (1, 2, 1) || rand_test_ait == (1, 3, 1)

# has_arc
for a in arcs(F)
    @test has_arc(F, a)
end

# add_arc rem_arc for arc type
while !isempty(arcs(F))
    a = ArcIter(F, :, :)

    recount = 0
    for n in keys(F.fadj), nbrs in keys(F.fadj[n]), key in keys(F.fadj[n][nbrs])
        recount += 1
    end

    @test recount == F.na

    @test add_arc!(F, a[1]) == a[1].key # nothing happens krn a[1] sudah ada
    @test add_arcs_from!(F, a[2:3]) == [a[i].key for i in 2:3] # nothing happens arc[2:3]
    @test rem_arcs_from!(F, a[1:3]) == length(1:3) # removed arcnya
    @test !rem_arc!(F, a[1]) # gaada sisanya
end

# source target key
@test source(Arc(1, 2, 1)) == 1
@test target(Arc(1, 2, 1)) == 2
@test key(Arc(1, 2, 1)) == 1

# ==
@test !has_arc(F, Arc(1, 2, 1))
@test Arc(1, 2, 1) == Arc((1, 2, 1))
add_arcs_from!(F,[(1, 2),(1, 3),(3, 1)])
@test F[1,2,1] == Dictionary()
@test (F[1,2,1,:flow] = 1) == 1
@test (F[1,2,3,:flow] = 2) == 2
@test collect(arcs(F, 1, 2)) == [Arc(1, 2, 1)]


text:
# sample code for magma binding

# MAGMA enum constants
# the whole file will be stored as enums.jl
#just like in JuliaGPU/MAGMA.jl

# MAGMA constants indicating the vectors status
# as input/output for some functions
# For example, the gesvd functions will use
# MagmaNoVec, MagmaSomeVec, MagmaAllVec and
# MagmaOverwriteVec to indicate the
# strategies that will be applied to the SVD
# U matrix and VT matrix in A = U Σ V**T
# (for MagmaOverwriteVec it is going to overwrite A)
const MagmaNoVec         = 301
const MagmaVec           = 302
const MagmaIVec          = 303
const MagmaAllVec        = 304
const MagmaSomeVec       = 305
const MagmaOverwriteVec  = 306
const MagmaBacktransVec  = 307

### some certain setting for this testing
using CUDAdrv
using CUDAapi
using CUDAnative
using CuArrays
# indicate the library position
const libmagma = "/usr/local/magma/lib/libmagma.so"


text:
using NiLang, NiLang.AD
using Test
using TensorOperations

@testset "HessianData" begin
    rings_init!()
    out! = beijingring!(6.0)
    a = beijingring!(2.0)
    b = beijingring!(3.0)

    @test nrings() == 3

    @test chfield(out!, Val(:x), 0.5).x == 0.5
    @test chfield(out!, value, 0.6).x == 0.6
    @test grad(chfield(a, grad, 0.5)) == 0.5

    @test hdata((out!, a)) == 0.0
    @instr hdata((out!, a)) ⊕ 0.5
    @instr hdata((a, out!)) ⊕ 0.5
    @test hdata((out!, a)) == 0.5
    @test hdata((a, out!)) == 0.5
end

@testset "hessian" begin
    h1 = (⊕(*)''(Loss(0.0), 2.0, 3.0); collect_hessian())
    h2 = nhessian(⊕(*), (Loss(0.0), 2.0, 3.0))
    @test h1 ≈ h2

    @i function test(a,b,c)
        a += b*c
        c += b^a
        ROT(a, b, c)
    end
    h1 = (test''(Loss(0.0), 2.0, 0.5); collect_hessian())
    h2 = nhessian(test, (Loss(0.0), 2.0, 0.5))
    @show h2
    @test isapprox(h1, h2, atol=1e-5)
end

function hessian_propagate(h, f, args; kwargs=())
    jac = jacobian(f, args; kwargs=kwargs)
    @tensor out[i,j,o] := jac[i,a]*h[a,b,o]*jac[j,b]
end

function hessian_propagate2(h, f, args; kwargs=())
    nargs = length(args)
    hes = zeros(nargs,nargs,nargs)
    @instr f(args...)
    for j=1:nargs
        # init rings
        rings_init!()
        largs = [beijingring!(x) for x in args]
        for i=1:nargs
            NiLang.AD.rings[i][1:i] .= h[1:i,i,j]
            NiLang.AD.rings[i][end:-1:end-i+1] .= h[i,1:i,j]
        end
        @instr (~f)(largs...)
        hes[:,:,j] .= collect_hessian()
    end
    hes
end

function rand_hes(n)
    h = randn(n,n,n)
    h + permutedims(h, (2,1,3))
end

@testset "hessian propagate" begin
    for op in [⊕(*), ⊕(/), ⊕(^), ROT]
        @show op
        h1 = local_hessian(op, (0.3, 0.4, 2.0))
        h2 = local_nhessian(op, (0.3, 0.4, 2.0))
        @test h1 ≈ h2

        h = rand_hes(3)
        h1 = hessian_propagate(copy(h), op, (0.3, 0.4, 2.0))
        h2 = hessian_propagate2(copy(h), op, (0.3, 0.4, 2.0))
        @show h1 - h2
        @test h1 ≈ h2
    end

    for op in [⊕(identity), ⊕(abs), SWAP, ⊕(exp), ⊕(log), ⊕(sin), ⊕(cos)]
        h1 = local_hessian(op, (0.3, 0.4))
        h2 = local_nhessian(op, (0.3, 0.4))
        @test h1 ≈ h2

        h = rand_hes(2)
        h1 = hessian_propagate(copy(h), op, (0.3, 0.4))
        h2 = hessian_propagate2(copy(h), op, (0.3, 0.4))
        @test h1 ≈ h2
    end

    for op in [NEG, CONJ]
        h1 = local_hessian(op, (0.3,))
        h2 = local_nhessian(op, (0.3,))
        @test h1 ≈ h2

        h = rand_hes(1)
        h1 = hessian_propagate(copy(h), op, (0.3,))
        h2 = hessian_propagate2(copy(h), op, (0.3,))
        @test h1 ≈ h2
    end
end


text:
# MIT License
#
# Copyright (c) 2020 Joshua E Gentry

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

@inline function build_rows(
            input::CHAINING,
            counter::UInt64,
            block_len::UInt32,
            flags::UInt32
        )::MVector{4, VECTOR}
    return MVector{4, VECTOR}((
        input[1], input[2],
        VECTOR((IV[1], IV[2], IV[3], IV[4])),
        VECTOR((UInt32(counter & 0xffffffff), UInt32(counter >> 32),  block_len, flags))
    ))
end

@inline function build_m(
            input::AbstractArray{UInt32}
        )::MVector{4, VECTOR}
    @inbounds return MVector{4, VECTOR}((
        vload(VECTOR, input, 1),
        vload(VECTOR, input, 5),
        vload(VECTOR, input, 9),
        vload(VECTOR, input, 13)
    ))
end

@inline function build_m(
            input::BLOCK
        )::MVector{4, VECTOR}
    return MVector{4, VECTOR}((input[1], input[2], input[3], input[4]))
end

@inline function g1(
            rows::MVector{4, VECTOR},
            m::VECTOR
        )::Nothing
    rows[1] = rows[1] + m + rows[2]
    rows[4] = rows[4] ⊻ rows[1]
    rows[4] = rot16(rows[4])
    rows[3] = rows[3] + rows[4]
    rows[2] = rows[2] ⊻ rows[3]
    rows[2] = rot12(rows[2])
    nothing
end

@inline function g2(
            rows::MVector{4, VECTOR},
            m::VECTOR
        )::Nothing
    rows[1] = rows[1] + m + rows[2]
    rows[4] = rows[4] ⊻ rows[1]
    rows[4] = rot8(rows[4])
    rows[3] = rows[3]  + rows[4]
    rows[2] = rows[2] ⊻ rows[3]
    rows[2] = rot7(rows[2])
    nothing
end

@inline function diagonalize(
            rows::MVector{4, VECTOR}
        )::Nothing
    rows[1] = shuffle(rows[1], (2, 1, 0, 3))
    rows[4] = shuffle(rows[4], (1, 0, 3, 2))
    rows[3] = shuffle(rows[3], (0, 3, 2, 1))
    nothing
end

@inline function undiagonalize(
            rows::MVector{4, VECTOR}
        )::Nothing
    rows[1] = shuffle(rows[1], (0, 3, 2, 1))
    rows[4] = shuffle(rows[4], (1, 0, 3, 2))
    rows[3] = shuffle(rows[3], (2, 1, 0, 3))
    nothing
end

@inline function compress_pre(
            rows::MVector{4, VECTOR},
            m::MVector{4, VECTOR},
        )::Nothing
    # Round 1
    t0 = shuffle(m[1], m[2], (2, 0, 2, 0))
    g1(rows, t0)
    t1 = shuffle(m[1], m[2], (3, 1, 3, 1))
    g2(rows, t1)
    diagonalize(rows)
    t2 = shuffle(m[3], m[4], (2, 0, 2, 0))
    t2 = shuffle(t2, (2, 1, 0, 3))
    g1(rows, t2)
    t3 = shuffle(m[3], m[4], (3, 1, 3, 1))
    t3 = shuffle(t3, (2, 1, 0, 3))
    g2(rows, t3)
    undiagonalize(rows)
    m[1], m[2], m[3], m[4] = t0, t1, t2, t3

    # Round 2
    t0 = shuffle(m[1], m[2], (3, 1, 1, 2))
    t0 = shuffle(t0, (0, 3, 2, 1))
    g1(rows, t0)
    t1 = shuffle(m[3], m[4], (3, 3, 2, 2))
    tt = shuffle(m[1], (0, 0, 3, 3))
    t1 = blend(tt, t1, 0xcc)
    g2(rows, t1)
    diagonalize(rows)
    t2 = unpacklo_64(m[4], m[2])
    tt = blend(t2, m[3], 0xc0)
    t2 = shuffle(tt, (1, 3, 2, 0))
    g1(rows, t2)
    t3 = unpackhi_32(m[2], m[4])
    tt = unpacklo_32(m[3], t3)
    t3 = shuffle(tt, (0, 1, 3, 2))
    g2(rows, t3)
    undiagonalize(rows)
    m[1], m[2], m[3], m[4] = t0, t1, t2, t3

    # Round 3
    t0 = shuffle(m[1], m[2], (3, 1, 1, 2))
    t0 = shuffle(t0, (0, 3, 2, 1))
    g1(rows, t0)
    t1 = shuffle(m[3], m[4], (3, 3, 2, 2))
    tt = shuffle(m[1], (0, 0, 3, 3))
    t1 = blend(tt, t1, 0xcc)
    g2(rows, t1)
    diagonalize(rows)
    t2 = unpacklo_64(m[4], m[2])
    tt = blend(t2, m[3], 0xc0)
    t2 = shuffle(tt, (1, 3, 2, 0))
    g1(rows, t2)
    t3 = unpackhi_32(m[2], m[4])
    tt = unpacklo_32(m[3], t3)
    t3 = shuffle(tt, (0, 1, 3, 2))
    g2(rows, t3)
    undiagonalize(rows)
    m[1], m[2], m[3], m[4] = t0, t1, t2, t3

    # Round 4
    t0 = shuffle(m[1], m[2], (3, 1, 1, 2))
    t0 = shuffle(t0, (0, 3, 2, 1))
    g1(rows, t0)
    t1 = shuffle(m[3], m[4], (3, 3, 2, 2))
    tt = shuffle(m[1], (0, 0, 3, 3))
    t1 = blend(tt, t1, 0xcc)
    g2(rows, t1)
    diagonalize(rows)
    t2 = unpacklo_64(m[4], m[2])
    tt = blend(t2, m[3], 0xc0)
    t2 = shuffle(tt, (1, 3, 2, 0))
    g1(rows, t2)
    t3 = unpackhi_32(m[2], m[4])
    tt = unpacklo_32(m[3], t3)
    t3 = shuffle(tt, (0, 1, 3, 2))
    g2(rows, t3)
    undiagonalize(rows)
    m[1], m[2], m[3], m[4] = t0, t1, t2, t3

    # Round 5
    t0 = shuffle(m[1], m[2], (3, 1, 1, 2))
    t0 = shuffle(t0, (0, 3, 2, 1))
    g1(rows, t0)
    t1 = shuffle(m[3], m[4], (3, 3, 2, 2))
    tt = shuffle(m[1], (0, 0, 3, 3))
    t1 = blend(tt, t1, 0xcc)
    g2(rows, t1)
    diagonalize(rows)
    t2 = unpacklo_64(m[4], m[2])
    tt = blend(t2, m[3], 0xc0)
    t2 = shuffle(tt, (1, 3, 2, 0))
    g1(rows, t2)
    t3 = unpackhi_32(m[2], m[4])
    tt = unpacklo_32(m[3], t3)
    t3 = shuffle(tt, (0, 1, 3, 2))
    g2(rows, t3)
    undiagonalize(rows)
    m[1], m[2], m[3], m[4] = t0, t1, t2, t3

    # Round 6
    t0 = shuffle(m[1], m[2], (3, 1, 1, 2))
    t0 = shuffle(t0, (0, 3, 2, 1))
    g1(rows, t0)
    t1 = shuffle(m[3], m[4], (3, 3, 2, 2))
    tt = shuffle(m[1], (0, 0, 3, 3))
    t1 = blend(tt, t1, 0xcc)
    g2(rows, t1)
    diagonalize(rows)
    t2 = unpacklo_64(m[4], m[2])
    tt = blend(t2, m[3], 0xc0)
    t2 = shuffle(tt, (1, 3, 2, 0))
    g1(rows, t2)
    t3 = unpackhi_32(m[2], m[4])
    tt = unpacklo_32(m[3], t3)
    t3 = shuffle(tt, (0, 1, 3, 2))
    g2(rows, t3)
    undiagonalize(rows)
    m[1], m[2], m[3], m[4] = t0, t1, t2, t3

    # Round 7
    t0 = shuffle(m[1], m[2], (3, 1, 1, 2))
    t0 = shuffle(t0, (0, 3, 2, 1))
    g1(rows, t0)
    t1 = shuffle(m[3], m[4], (3, 3, 2, 2))
    tt = shuffle(m[1], (0, 0, 3, 3))
    t1 = blend(tt, t1, 0xcc)
    g2(rows, t1)
    diagonalize(rows)
    t2 = unpacklo_64(m[4], m[2])
    tt = blend(t2, m[3], 0xc0)
    t2 = shuffle(tt, (1, 3, 2, 0))
    g1(rows, t2)
    t3 = unpackhi_32(m[2], m[4])
    tt = unpacklo_32(m[3], t3)
    t3 = shuffle(tt, (0, 1, 3, 2))
    g2(rows, t3)
    undiagonalize(rows)
    nothing
end

function compress(
            chaining::CHAINING,
            block::Union{AbstractVector{UInt32}, BLOCK},
            block_len::UInt32,
            counter::UInt64,
            flags::UInt32
        )::CHAINING
    # @printf("CV1: 0x%08x 0x%08x 0x%08x 0x%08x\n", chaining_value[1][1], chaining_value[1][2], chaining_value[1][3], chaining_value[1][4])
    # @printf("CV2: 0x%08x 0x%08x 0x%08x 0x%08x\n", chaining_value[2][1], chaining_value[2][2], chaining_value[2][3], chaining_value[2][4])
    # @printf("BL1: 0x%08x 0x%08x 0x%08x 0x%08x\n", block[1], block[2], block[3], block[4])
    # @printf("BL2: 0x%08x 0x%08x 0x%08x 0x%08x\n", block[5], block[6], block[7], block[8])
    # @printf("BL3: 0x%08x 0x%08x 0x%08x 0x%08x\n", block[9], block[10], block[11], block[12])
    # @printf("BL4: 0x%08x 0x%08x 0x%08x 0x%08x\n", block[13], block[14], block[15], block[16])
    # @printf("block_len: %d counter: %d flags: %d\n", block_len, counter, flags)
    rows = build_rows(chaining, counter, block_len, flags)
    m    = build_m(block)

    compress_pre(rows, m)

    return CHAINING((rows[1] ⊻ rows[3], rows[2] ⊻ rows[4]))
end

function compress_xof(
            chaining::CHAINING,
            block::BLOCK,
            block_len::UInt32,
            counter::UInt64,
            flags::UInt32,
            out::Vector{UInt8},
            offset::Int = 1
        )::Nothing
    rows = build_rows(chaining, counter, block_len, flags)
    m    = build_m(block)
    compress_pre(rows, m)

    out = reinterpret(UInt32, view(out, offset:offset+64-1))
    vstore(rows[1] ⊻ rows[3], out, 1)
    vstore(rows[2] ⊻ rows[4], out, 5)
    vstore(rows[3] ⊻ chaining[1], out, 9)
    vstore(rows[4] ⊻ chaining[2], out, 13)

    nothing
end


text:
D8fields = [x^8 - 3*x^7 + 4*x^6 - 3*x^5 + 3*x^4 - 3*x^3 + 4*x^2 - 3*x + 1,
x^8 - x^7 + 2*x^6 + 3*x^5 - x^4 + 3*x^3 + 2*x^2 - x + 1,
x^8 - 2*x^7 + 2*x^6 + 2*x^5 - 2*x^4 + 2*x^3 + 2*x^2 - 2*x + 1,
x^8 - 2*x^7 + 2*x^6 - 2*x^5 + 7*x^4 - 10*x^3 + 8*x^2 - 4*x + 1,
x^8 + 6*x^4 + 1,
x^8 - 3*x^7 + 9*x^6 - 13*x^5 + 18*x^4 - 11*x^3 + 11*x^2 - 4*x + 1,
x^8 + 4*x^6 + x^4 - 6*x^2 + 4,
x^8 - 4*x^6 + 8*x^4 - 4*x^2 + 1,
x^8 + 5*x^6 + 4*x^4 + 5*x^2 + 1,
x^8 - 2*x^7 + 2*x^6 + 4*x^5 - x^4 + 4*x^3 + 2*x^2 - 2*x + 1,
x^8 + 3*x^4 + 1,
x^8 + 3*x^4 + 9,
x^8 - 4*x^7 + 8*x^6 - 10*x^5 + 11*x^4 - 10*x^3 + 2*x^2 + 2*x + 1,
x^8 - x^7 + 8*x^6 + 7*x^5 + 19*x^4 + 7*x^3 + 8*x^2 - x + 1,
x^8 - 3*x^7 + x^6 + 6*x^5 + 6*x^4 + 6*x^3 + x^2 - 3*x + 1,
x^8 - 4*x^7 + 8*x^6 - 6*x^5 + 7*x^4 - 10*x^3 + 2*x^2 - 6*x + 9,
x^8 - x^7 + 2*x^6 - 3*x^5 - x^4 + 6*x^3 + 8*x^2 + 8*x + 16,
x^8 - 5*x^6 + 28*x^4 + 15*x^2 + 9,
x^8 + 2*x^6 + 10*x^4 + 2*x^2 + 1,
x^8 - 2*x^7 + 2*x^6 - 6*x^5 + 14*x^4 - 30*x^3 + 50*x^2 - 50*x + 25]

q8fields = [ x^8 + 12*x^6 + 36*x^4 + 36*x^2 + 9,
  x^8 - 12*x^6 + 36*x^4 - 36*x^2 + 9,
  x^8 - x^7 - 34*x^6 + 29*x^5 + 361*x^4 - 305*x^3 - 1090*x^2 + 1345*x - 395,
  x^8 - 3*x^7 + 22*x^6 - 60*x^5 + 201*x^4 - 450*x^3 + 1528*x^2 - 3069*x + 4561,
  x^8 + 60*x^6 + 900*x^4 + 4500*x^2 + 5625,
  x^8 - 60*x^6 + 900*x^4 - 4500*x^2 + 5625,
  x^8 - x^7 - 47*x^6 + 40*x^5 + 581*x^4 - 220*x^3 - 2038*x^2 - 932*x - 109,
  x^8 + 84*x^6 + 1764*x^4 + 12348*x^2 + 21609,
  x^8 - 84*x^6 + 1764*x^4 - 12348*x^2 + 21609,
  x^8 + 44*x^6 + 308*x^4 + 484*x^2 + 121,
  x^8 - 44*x^6 + 308*x^4 - 484*x^2 + 121,
  x^8 + 60*x^6 + 1170*x^4 + 9000*x^2 + 22500,
  x^8 + 60*x^6 + 810*x^4 + 1800*x^2 + 900,
  x^8 - 60*x^6 + 1170*x^4 - 9000*x^2 + 22500,
  x^8 - 60*x^6 + 810*x^4 - 1800*x^2 + 900,
  x^8 - 3*x^7 - 63*x^6 + 90*x^5 + 1311*x^4 - 20*x^3 - 7702*x^2 - 5524*x - 1009,
  x^8 + 68*x^6 + 986*x^4 + 4624*x^2 + 4624,
  x^8 - 68*x^6 + 986*x^4 - 4624*x^2 + 4624,
  x^8 - 3*x^7 - 67*x^6 - 16*x^5 + 863*x^4 + 1276*x^3 + 392*x^2 - 54*x + 1,
  x^8 - x^7 + 50*x^6 + 71*x^5 + 529*x^4 + 2173*x^3 + 842*x^2 + 5545*x + 17623]


text:
module Paraml

import YAML

include("utils.jl")
include("check.jl")
include("parse.jl")

"""
    create_params(
        defn_path::AbstractString,
        param_paths...;
        error_on_unused::Bool = false,
        out_path::AbstractString = "",
    )

Entry function - given the path to the parameter definitions and files, parse and create a params dictionary.

The defn_path and param_paths can be either a single yaml file, or a directory containing yaml files.

# Arguments
- `defn_path`: path to the parameter definitions
- `param_paths...`: paths to parameter files or directories. The files will be merged in the passed order so that item 'a' in the first params will be overwritten by item 'a' in the second params.
- `out_path`: path to directory where computed params will be saved if passed
- `error_on_unused`: throw a hard error if there are unused parameters, otherwise warnings are only printed

# Returns
- dictionary with computed/validated model paramters with defaults filled in where needed
"""
function create_params(
    defn_path::AbstractString,
    param_paths...;
    error_on_unused::Bool = false,
    out_path::AbstractString = "",
)
    defs = build_yaml(defn_path)

    params = Dict{String,Any}()
    for param_path in param_paths
        cur_params = build_yaml(param_path)
        params = merge_rec(params, cur_params)
    end

    classes = parse_classes(defs, params)
    parsed = parse_params(defs, params; classes)

    if !isempty(out_path)
        @info "writing parsed params to $out_path"
        YAML.write_file(out_path, parsed)
    end

    @info "\nChecking for unused parameters..."
    num_unused = warn_unused_params(parsed, params)
    @info "$num_unused unused parameters found"
    if error_on_unused
        @assert num_unused == 0 "There are unused parameters passed to the parser (see print statements)"
    end

    return parsed
end

export create_params

end


text:
# This file contains
# a type for input data

"""
Module containing the type to interact with input data.
"""
module Input

# Export
export InputStruct

# Load extras
using ..Extras

# Arrays with non-standard indexing
using OffsetArrays

# Precisions and formats of numbers
using ..Prec

"""
    InputStruct()

Instantiate this type to interact with input data.

# Data
- `N::`[`IT`](@ref)`=0`: sample size;
- `Δt::`[`RT`](@ref)`=0`: sample step;
- `q::`[`RT`](@ref)`=0`: significance level;
- `t::Vector{`[`RT`](@ref)`}=[]`: time array;
- `x::Vector{`[`RT`](@ref)`}=[]`: values array.

# Methods
- [`read!`](@ref)`(file::AbstractString)`: read input data from a file;
- [`write`](@ref)`(file::AbstractString)`: write input data to a file;
- [`example`](@ref)`(file::AbstractString)`: generate an example of the input/output file;
- [`reset!`](@ref)`()`: reset an instance to default values.

# Note
Data can be also read when calling an instance like so:
```jldoctest; output = false
using Scats
s = Scats.API()
file, _ = mktemp()
s.Input.example(file)
s.Input(file)

# output

```
"""
mutable struct InputStruct

    # Data
    N::IT         # Sample size
    Δt::RT        # Sample step
    q::RT         # significance level
    t::Vector{RT} # Time array
    x::Vector{RT} # Values array

    # Methods
    read!::Function   # Read input data from a file
    write::Function   # Write input data to a file
    example::Function # Generate an example of the input/output file
    reset!::Function  # Reset an instance to default values

    # Construct an object of this type
    function InputStruct()
        this = new(0, 0, 0, [], [])
        this.read! = function (file::AbstractString) read!(this, file) end
        this.write = function (file::AbstractString) write(this, file) end
        this.example = example
        this.reset! = function () reset!(this) end
        return this
    end

    # Read the data when calling an instance
    function (input::InputStruct)(file::AbstractString)
        input.read!(file)
    end

end

# Sources
include("input_exceptions.jl") # Exceptions
include("input_read.jl")       # Read input data from a file
include("input_write.jl")      # Write input data to a file
include("input_example.jl")    # Generate an example of the input/output file
include("input_reset.jl")      # Reset an instance to default values

end


text:
module MeshCore

using StaticArrays

include("utilities.jl")
include("attributes.jl")
include("shapedesc.jl")
include("shapes.jl")
include("increl.jl")

# We can either use/import individual functions from MeshCore like so:
# ```
# using MeshCore: nrelations, nentities
# ```
# or we can bring into our context all exported symbols as
# ```
# using MeshCore.Exports
# ```
include("Exports.jl")

end # module


text:
"""
   SIGNAL TO NOISE MEASUREMENT
"""
function SNRatio( corr::Array{T,N}, method::S, width::I ) where {T<:Real,N}

    max2       = 0.0;
    h, w, d    = size(corr,1), size(corr,2), size(corr,3)
    idx1, max1 = firstPeak( corr )

    if     method == "peak2peak"
		#idx2, max2 = secondPeak( corr, width, idx1[1:N] )
        idx2, max2 = value2ndPeak( corr, max1, width )
    elseif method == "peak2mean"
        max2 = Statistics.mean( corr )
    else
        throw(ErrorException("Error: use 'peak2peak' or 'peak2mean'.\n"))
    end

    # if max2 == 0, sig2noise = inf. Check this at return

	# When using ZNCC, max2 might be negative, then max1/max2 is negative value
	# and we can't take the logarithm. Therefore, shift max1 and max2 up, until
	# max2 = 1 and max1 = max1 - max2 + 1. Then we can divide max1/max2 without problems
	if max2 < 0 && max1 > 0
		max1 = max1 - max2 + 1
		max2 = 1 # max2 - max2 + 1
	end

    sig2noise = max2 !== 0.0 ? max1 / max2 : 0.0

    return sig2noise;
end

function value2ndPeak( cmat::Array{T,N}, max1, rad ) where {T<:AbstractFloat,N}

	h, w, d = size(cmat,1), size(cmat,2), size(cmat,3);
	zrad    = ( N == 2 ) ? 0 : rad;
    max2    = cmat[ 1, 1, 1 ]
    maxid   = ( 1, 1, 1 )

    @inbounds for dep in 1+zrad:d-zrad, col in 1+rad:w-rad, row in 1+rad:h-rad

        ( cmat[row,col,dep] == max1 || cmat[row,col,dep] <= cmat[maxid...] ) && ( continue; );

        ismax = true;
        for z in dep-zrad:dep+zrad, c in col-rad:col+rad, r in row-rad:row+rad
            ( cmat[row,col,dep] < cmat[r,c,z]  ) && ( ismax = false; break; )
        end

        max2  = ( ismax ) ? cmat[row,col,dep] : max2;
		#maxid = ( ismax ) ? ( row, col, dep ) : maxid;
    end

    return maxid, max2
end

max_issues_repo_name: ErichZimmer/quickPIV

max_forks_repo_name: ErichZimmer/quickPIV