text:
export kummer_extension, ray_class_field, hilbert_class_field, prime_decomposition_type, defining_modulus

add_verbose_scope(:ClassField)
add_assert_scope(:ClassField)
#set_assert_level(:ClassField, 1)

###############################################################################
#
#  Number_field interface and reduction to prime power case
#
###############################################################################
@doc Markdown.doc"""
    NumberField(CF::ClassField) -> NfRelNS{nf_elem}

Given a (formal) abelian extension, compute the class field by finding defining
polynomials for all prime power cyclic subfields.

Note, the return type is always a non-simple extension.
"""
function NumberField(CF::ClassField{S, T}; redo::Bool = false, using_norm_relation = false) where {S, T}
  if isdefined(CF, :A) && !redo
    return CF.A
  end

  res = ClassField_pp{S, T}[]
  ord = torsion_units_order(base_field(CF))
  G = codomain(CF.quotientmap)
  @assert issnf(G)
  q = GrpAbFinGenElem[G[i] for i=1:ngens(G)]
  for i=1:ngens(G)
    o = G.snf[i]
    lo = factor(o)
    for (p, e) = lo.fac
      q[i] = p^e*G[i]
      S1, mQ = quo(G, q, false)
      if using_norm_relation && !divides(fmpz(ord), order(S1))[1]
        push!(res, ray_class_field_cyclic_pp_Brauer(CF, mQ))
      else
        push!(res, ray_class_field_cyclic_pp(CF, mQ))
      end
    end
    q[i] = G[i]
  end
  CF.cyc = res
  if isempty(res)
    @assert isone(degree(CF))
    Ky = PolynomialRing(base_field(CF), "y", cached = false)[1]
    CF.A = number_field([gen(Ky)-1], check = false, cached = false)[1]
  else
    CF.A = number_field([x.A.pol for x = CF.cyc], check = false, cached = false)[1]
  end
  return CF.A
end

function ray_class_field_cyclic_pp(CF::ClassField{S, T}, mQ::GrpAbFinGenMap) where {S, T}
  @vprint :ClassField 1 "cyclic prime power class field of degree $(degree(CF))\n"
  CFpp = ClassField_pp{S, T}()
  CFpp.quotientmap = compose(CF.quotientmap, mQ)
  CFpp.rayclassgroupmap = CF.rayclassgroupmap
  @assert domain(CFpp.rayclassgroupmap) == domain(CFpp.quotientmap)

  @vprint :ClassField 1 "computing the S-units...\n"
  @vtime :ClassField 1 _rcf_S_units(CFpp)
  @vprint :ClassField 1 "finding the Kummer extension...\n"
  @vtime :ClassField 1 _rcf_find_kummer(CFpp)
  @vprint :ClassField 1 "reducing the generator...\n"
  @vtime :ClassField 1 _rcf_reduce(CFpp)
  @vprint :ClassField 1 "descending ...\n"
  @vtime :ClassField 1 _rcf_descent(CFpp)
  return CFpp
end

################################################################################
#
#  Using norm relations to get the s-units
#
################################################################################

function ray_class_field_cyclic_pp_Brauer(CF::ClassField{S, T}, mQ::GrpAbFinGenMap) where {S, T}
  @vprint :ClassField 1 "cyclic prime power class field of degree $(order(codomain(mQ)))\n"
  CFpp = ClassField_pp{S, T}()
  CFpp.quotientmap = compose(CF.quotientmap, mQ)
  CFpp.rayclassgroupmap = CF.rayclassgroupmap
  @assert domain(CFpp.rayclassgroupmap) == domain(CFpp.quotientmap)

  e = degree(CFpp)
  v, q = ispower(e)
  k = base_field(CFpp)
  CE = cyclotomic_extension(k, e)

  @vprint :ClassField 1 "computing the S-units...\n"
  @vtime :ClassField 1 _rcf_S_units_using_Brauer(CFpp)
  mr = CF.rayclassgroupmap
  mq = CFpp.quotientmap
  KK = kummer_extension(e, FacElem{nf_elem, AnticNumberField}[CFpp.a])
  ng, mng = norm_group(KK, CE.mp[2], mr)
  attempt = 1
  while !iszero(mng*mq)
    attempt += 1
    @vprint :ClassField 1 "attempt number $(attempt)\n"
    _rcf_S_units_enlarge(CE, CFpp)
    KK = kummer_extension(e, FacElem{nf_elem, AnticNumberField}[CFpp.a])
    ng, mng = norm_group(KK, CE.mp[2], mr)
  end

  @vprint :ClassField 1 "reducing the generator...\n"
  @vtime :ClassField 1 _rcf_reduce(CFpp)
  @vprint :ClassField 1 "descending ...\n"
  @vtime :ClassField 1 _rcf_descent(CFpp)
  return CFpp
end

################################################################################
#
#  S-units computation
#
################################################################################

function _rcf_S_units_enlarge(CE, CF::ClassField_pp)
  lP = CF.sup
  OK = order(lP[1])
  f = maximum(minimum(p) for p in lP)
  for i = 1:5
    f = next_prime(f)
    push!(lP, prime_decomposition(OK, f)[1][1])
  end
  e = degree(CF)
  @vtime :ClassField 3 S, mS = NormRel._sunit_group_fac_elem_quo_via_brauer(nf(OK), lP, e)
  KK = kummer_extension(e, FacElem{nf_elem, AnticNumberField}[mS(S[i]) for i=1:ngens(S)])
  CF.bigK = KK
  lf = factor(minimum(defining_modulus(CF)[1]))
  lfs = Set(collect(keys(lf.fac)))
  CE.kummer_exts[lfs] = (lP, KK)
  _rcf_find_kummer(CF)
  return nothing
end

function _rcf_S_units_using_Brauer(CF::ClassField_pp)

  f = defining_modulus(CF)[1]
  @vprint :ClassField 2 "Kummer extension with modulus $f\n"
  k1 = base_field(CF)

  #@assert Hecke.isprime_power(e)
  @vprint :ClassField 2 "Adjoining the root of unity\n"
  C = cyclotomic_extension(k1, degree(CF))
  @vtime :ClassField 2 automorphisms(C, copy = false)
  G, mG = automorphism_group(C.Ka)
  @vtime :ClassField 3 fl = NormRel.has_useful_brauer_relation(G)
  if fl
    @vtime :ClassField 3 lP, KK = _s_unit_for_kummer_using_Brauer(C, minimum(f))
  else
    lP, KK = _s_unit_for_kummer(C, minimum(f))
  end

  CF.bigK = KK
  CF.sup = lP
  _rcf_find_kummer(CF)
  return nothing
end

#This function finds a set S of primes such that we can find a Kummer generator in it.
function _s_unit_for_kummer_using_Brauer(C::CyclotomicExt, f::fmpz)

  e = C.n
  lf = factor(f)
  lfs = Set(collect(keys(lf.fac)))
  for (k, v) in C.kummer_exts
    if issubset(lfs, k)
      return v
    end
  end
  K = absolute_field(C)
  @vprint :ClassField 2 "Maximal order of cyclotomic extension\n"
  ZK = maximal_order(K)
  if isdefined(ZK, :lllO)
    ZK = ZK.lllO
  end

  lP = Hecke.NfOrdIdl[]

  for p = keys(lf.fac)
    #I remove the primes that can't be in the conductor
    lp = prime_decomposition(ZK, p)
    for (P, s) in lp
      if gcd(norm(P), e) != 1 || gcd(norm(P)-1, e) != 1
        push!(lP, P)
      end
    end
  end

  if length(lP) < 10
    #add some primes
    pp = f
    while length(lP) < 10
      pp = next_prime(pp)
      lpp = prime_decomposition(ZK, pp)
      if !isone(length(lpp))
        push!(lP, lpp[1][1])
      end
    end
  end
  @vprint :ClassField 3 "Computing S-units with $(length(lP)) primes\n"
  @vtime :ClassField 3 S, mS = NormRel._sunit_group_fac_elem_quo_via_brauer(C.Ka, lP, e)
  KK = kummer_extension(e, FacElem{nf_elem, AnticNumberField}[mS(S[i]) for i = 1:ngens(S)])
  C.kummer_exts[lfs] = (lP, KK)
  return lP, KK
end

###############################################################################
#
#  Find small generator of class group
#
###############################################################################

function find_gens(mR::Map, S::PrimesSet, cp::fmpz=fmpz(1))
# mR: SetIdl -> GrpAb (inv of ray_class_group or Frobenius or so)
  ZK = order(domain(mR))
  R = codomain(mR)
  sR = GrpAbFinGenElem[]
  lp = elem_type(domain(mR))[]

  q, mq = quo(R, sR, false)
  s, ms = snf(q)
  if order(s) == 1
    return lp, sR
  end
  for p in S
    if cp % p == 0 || index(ZK) % p == 0
      continue
    end
    @vprint :ClassField 2 "doing $p\n"
    lP = prime_decomposition(ZK, p)

    f=R[1]
    for (P, e) = lP
      try
        f = mR(P)
      catch e
        if !isa(e, BadPrime)
          rethrow(e)
        end
        break # try new prime number
      end
      if iszero(mq(f))
        break
      end
      #At least one of the coefficient of the element
      #must be invertible in the snf form.
      el = ms\f
      to_be = false
      for w = 1:ngens(s)
        if gcd(s.snf[w], el.coeff[w]) == 1
          to_be = true
          break
        end
      end
      if !to_be
        continue
      end
      push!(sR, f)
      push!(lp, P)
      q, mq = quo(R, sR, false)
      s, ms = snf(q)
    end
    if order(q) == 1
      break
    end
  end
  return lp, sR

end

function find_gens_descent(mR::Map, A::ClassField_pp, cp::fmpz)
  ZK = order(domain(mR))
  C = cyclotomic_extension(nf(ZK), degree(A))
  R = codomain(mR)
  Zk = order(codomain(A.rayclassgroupmap))
  sR = GrpAbFinGenElem[]
  lp = elem_type(domain(mR))[]
  q, mq = quo(R, sR, false)
  s, ms = snf(q)

  PPS = A.bigK.frob_gens[1]
  for p in PPS
    P = intersect_prime(C.mp[2], p, Zk)
    local f::GrpAbFinGenElem
    try
      f = mR(P)
    catch e
      if !isa(e, BadPrime)
        rethrow(e)
      end
      break # try new prime number
    end
    if iszero(mq(f))
      continue
    end
    #At least one of the coefficient of the element
    #must be invertible in the snf form.
    el = ms\f
    to_be = false
    for w = 1:ngens(s)
      if gcd(s.snf[w], el.coeff[w]) == 1
        to_be = true
        break
      end
    end
    if !to_be
      continue
    end
    push!(sR, f)
    push!(lp, P)
    q, mq = quo(R, sR, false)
    s, ms = snf(q)
    if order(s) == divexact(order(R), degree(A.K))
      break
    end
  end

  if degree(C.Kr) != 1
    RR = ResidueRing(FlintZZ, degree(A))
    U, mU = unit_group(RR)
    if degree(C.Kr) < order(U)  # there was a common subfield, we
                              # have to pass to a subgroup

      f = C.Kr.pol
      # Can do better. If the group is cyclic (e.g. if p!=2), we already know the subgroup!
      s, ms = sub(U, GrpAbFinGenElem[x for x in U if iszero(f(gen(C.Kr)^Int(lift(mU(x)))))], false)
      ss, mss = snf(s)
      U = ss
      #mg = mg*ms*mss
      mU = mss * ms * mU
    end
    for i = 1:ngens(U)
      l = Int(lift(mU(U[i])))
      S = PrimesSet(100, -1, degree(A), l)
      found = false
      for p in S
        if mod(cp, p) == 0 || mod(index(ZK), p) == 0
          continue
        end
        @vprint :ClassField 2 "Computing Frobenius over $p\n"
        lP = prime_decomposition(ZK, p)

        f = R[1]
        for (P, e) = lP
          lpp = prime_decomposition(C.mp[2], P)
          if divexact(degree(lpp[1][1]), degree(P)) != U.snf[i]
            continue
          end
          try
            f = mR(P)
          catch e
            if !isa(e, BadPrime)
              rethrow(e)
            end
            break
          end
          push!(sR, f)
          push!(lp, P)
          q, mq = quo(R, sR, false)
          s, ms = snf(q)
          found = true
          break
        end
        if found
          break
        end
      end
    end
  end

  if order(q) == 1
    return lp, sR
  end

  @vprint :ClassField 3 "Bad Case in Descent\n"
  S = PrimesSet(300, -1)
  for p in S
    if cp % p == 0 || index(ZK) % p == 0
      continue
    end
    @vprint :ClassField 2 "Computing Frobenius over $p\n"
    lP = prime_decomposition(ZK, p)

    f=R[1]
    for (P, e) = lP
      try
        f = mR(P)
      catch e
        if !isa(e, BadPrime)
          rethrow(e)
        end
        break # try new prime number
      end
      if iszero(mq(f))
        continue
      end
      #At least one of the coefficient of the element
      #must be invertible in the snf form.
      el = ms\f
      to_be = false
      for w = 1:ngens(s)
        if gcd(s.snf[w], el.coeff[w]) == 1
          to_be = true
          break
        end
      end
      if !to_be
        continue
      end
      push!(sR, f)
      push!(lp, P)
      q, mq = quo(R, sR, false)
      s, ms = snf(q)
    end
    if order(q) == 1
      break
    end
  end
  return lp, sR
end

################################################################################
#
#  S-units computation
#
################################################################################

function _rcf_S_units(CF::ClassField_pp)

  f = defining_modulus(CF)[1]
  @vprint :ClassField 2 "Kummer extension with modulus $f\n"
  k1 = base_field(CF)

  #@assert Hecke.isprime_power(e)
  @vprint :ClassField 2 "Adjoining the root of unity\n"
  C = cyclotomic_extension(k1, degree(CF))

  #We could use f, but we would have to factor it.
  @vtime :ClassField 3 lP, KK = _s_unit_for_kummer(C, minimum(f))
  CF.bigK = KK
  CF.sup = lP
  return nothing
end

#This function finds a set S of primes such that we can find a Kummer generator in it.
function _s_unit_for_kummer(C::CyclotomicExt, f::fmpz)

  e = C.n
  lf = factor(f)
  lfs = Set(collect(keys(lf.fac)))
  for (k, v) in C.kummer_exts
    if issubset(lfs, k)
      return v
    end
  end
  K = absolute_field(C)
  @vprint :ClassField 2 "Maximal order of cyclotomic extension\n"
  ZK = maximal_order(K)
  @vprint :ClassField 2 "Class group of cyclotomic extension: $K\n"
  c, mc = class_group(ZK)
  allow_cache!(mc)
  @vprint :ClassField 2 "... $c\n"
  c, mq = quo(c, e, false)
  mc = compose(pseudo_inv(mq), mc)

  lP = Hecke.NfOrdIdl[]

  for p = keys(lf.fac)
     #I remove the primes that can't be in the conductor
     lp = prime_decomposition(ZK, p)
     for (P, s) in lp
       if gcd(norm(P), e) != 1 || gcd(norm(P)-1, e) != 1
         push!(lP, P)
       end
     end
  end

  g = Array{GrpAbFinGenElem, 1}(undef, length(lP))
  for i = 1:length(lP)
    g[i] = preimage(mc, lP[i])
  end

  q, mq = quo(c, g, false)
  mc = compose(pseudo_inv(mq), mc)

  #@vtime :ClassField 3
  lP = vcat(lP, find_gens(pseudo_inv(mc), PrimesSet(100, -1))[1])::Vector{NfOrdIdl}
  @vprint :ClassField 2 "using $lP of length $(length(lP)) for S-units\n"
  if isempty(lP)
    U, mU = unit_group_fac_elem(ZK)
    KK = kummer_extension(e, FacElem{nf_elem, AnticNumberField}[mU(U[i]) for i = 1:ngens(U)])
  else
    #@vtime :ClassField 2
    S, mS = Hecke.sunit_group_fac_elem(lP)
    #@vtime :ClassField 2
    KK = kummer_extension(e, FacElem{nf_elem, AnticNumberField}[mS(S[i]) for i=1:ngens(S)])
  end
  C.kummer_exts[lfs] = (lP, KK)
  return lP, KK
end

###############################################################################
#
#  First step: Find the Kummer extension over the cyclotomic field
#
###############################################################################
#=
  next, to piece things together:
    have a quo of some ray class group in k,
    taking primes in k over primes in Z that are 1 mod n
    then the prime is totally split in Kr, hence I do not need to
      do relative splitting and relative ideal norms. I am lazy
        darn: I still need to match the ideals
    find enough such primes to generate the rcg quotient (via norm)
                       and the full automorphism group of the big Kummer

          Kr(U^(1/n))  the "big" Kummer ext
         /
        X(z) = Kr(x^(1/n)) the "target"
      / /
    X  Kr = k(z)  = Ka
    | /             |
    k               |
    |               |
    Q               Q

    this way we have the map (projection) from "big Kummer" to
    Aut(X/k) = quo(rcg)
    The generator "x" is fixed by the kernel of this map

    Alternatively, "x" could be obtained via Hecke's theorem, ala Cohen

    Finally, X is derived via descent
=#

# mR: GrpAb A -> Ideal in k, only preimage used
# cf: Ideal in K -> GrpAb B, only image
# mp:: k -> K inclusion
# builds a (projection) from B -> A identifying (pre)images of
# prime ideals, the ideals are coprime to cp and ==1 mod n

function build_map(CF::ClassField_pp, K::KummerExt, c::CyclotomicExt)
  #mR should be GrpAbFinGen -> IdlSet
  #          probably be either "the rcg"
  #          or a compositum, where the last component is "the rcg"
  # we need this to get the defining modulus - for coprime testing
  m = defining_modulus(CF)[1]
  ZK = maximal_order(base_ring(K.gen[1]))
  cp = lcm(minimum(m), discriminant(ZK))

  Zk = order(m)
  mp = c.mp[2]
  cp = lcm(cp, index(Zk))
  Sp = Hecke.PrimesSet(100, -1, c.n, 1) #primes = 1 mod n, so totally split in cyclo

  #@vtime :ClassField 3
  lp, sG = find_gens(K, Sp, cp)
  G = K.AutG
  sR = Array{GrpAbFinGenElem, 1}(undef, length(lp))
  #@vtime :ClassField 3
  for i = 1:length(lp)
    p = intersect_nonindex(mp, lp[i], Zk)
    #Since the prime are totally split in the cyclotomic extension by our choice, we can ignore the valuation of the norm
    #sR[i] = valuation(norm(lp[i]), norm(p))*CF.quotientmap(preimage(CF.rayclassgroupmap, p))
    sR[i] = CF.quotientmap(preimage(CF.rayclassgroupmap, p))
  end
  @hassert :ClassField 1 order(quo(G, sG, false)[1]) == 1
       # if think if the quo(G, ..) == 1, then the other is automatic
       # it is not, in general it will never be.
       #example: Q[sqrt(10)], rcf of 16*Zk
  # now the map G -> R sG[i] -> sR[i]
  h = hom(sG, sR, check = false)
  @hassert :ClassField 1 !isone(gcd(fmpz(degree(CF)), minimum(m))) || issurjective(h)
  CF.h = h
  return h
end

function _rcf_find_kummer(CF::ClassField_pp{S, T}) where {S, T}

  #if isdefined(CF, :K)
  #  return nothing
  #end
  f = defining_modulus(CF)[1]
  @vprint :ClassField 2 "Kummer extension with modulus $f\n"
  k1 = base_field(CF)

  #@assert Hecke.isprime_power(e)
  @vprint :ClassField 2 "Adjoining the root of unity\n"
  C = cyclotomic_extension(k1, degree(CF))

  #We could use f, but we would have to factor it.
  #As we only need the support, we save the computation of
  #the valuations of the ideal
  KK = CF.bigK
  lP = CF.sup
  @vprint :ClassField 2 "building Artin map for the large Kummer extension\n"
  @vtime :ClassField 2 h = build_map(CF, KK, C)
  @vprint :ClassField 2 "... done\n"

  #TODO:
  #If the s-units are not large enough, the map might be trivial
  #We could do something better.
  if iszero(h)
    CF.a = FacElem(one(C.Ka))
    return nothing
  end
  k, mk = kernel(h, false)
  G = domain(h)

  # Now, we find the kummer generator by considering the action
  # of the automorphisms on the s-units
  # x needs to be fixed by k
  # that means x needs to be in the kernel:
  # x = prod gen[1]^n[i] -> prod (z^a[i] gen[i])^n[i]
  #                            = z^(sum a[i] n[i]) x
  # thus it works iff sum a[i] n[i] = 0
  # for all a in the kernel
  R = ResidueRing(FlintZZ, C.n, cached=false)
  M = change_base_ring(R, mk.map)
  i, l = right_kernel(M)
  @assert i > 0
  n = lift(l)
  e1 = degree(CF)
  N = GrpAbFinGen(fmpz[fmpz(e1) for j=1:nrows(n)])
  s, ms = sub(N, GrpAbFinGenElem[N(fmpz[n[j, ind] for j=1:nrows(n)]) for ind=1:i], false)
  ms = Hecke.make_domain_snf(ms)
  H = domain(ms)
  @hassert :ClassField 1 iscyclic(H)
  o = Int(order(H))
  c = fmpz(1)
  if o < fmpz(e1)
    c = div(fmpz(e1), o)
  end
  g = ms(H[1])
  @vprint :ClassField 2 "g = $g\n"
  #@vprint :ClassField 2 "final $n of order $o and e=$e\n"
  a = FacElem(Dict{nf_elem, fmpz}(one(C.Ka) => fmpz(1)))
  o2 = div(o, 2)
  for i = 1:ngens(N)
    eeee = div(mod(g[i], fmpz(e1)), c)
    if iszero(eeee)
      continue
    end
    mul!(a, a, KK.gen[i]^eeee)
  end
  #@vprint :ClassField 2 "generator $a\n"
  CF.a = a
  CF.sup_known = true
  CF.o = o
  CF.defect = c
  return nothing
end

###############################################################################
#
#  Descent to K
#
###############################################################################

#This function computes a primitive element for the target extension with the
#roots of unit over the base field and the action of the automorphisms on it.
#The Kummer generator is always primitive! (Carlo and Claus)
function _find_prim_elem(CF::ClassField_pp, AutA)

  AutA_gen = CF.AutG
  A = domain(AutA_gen[1])
  pe = gen(A)
  Auto = Dict{GrpAbFinGenElem, NfRelElem{nf_elem}}(find_orbit(AutA_gen, AutA, pe))
  if degree(CF) != degree(A)
    #In this case, gen(A) might not be primitive...
    while length(Auto) != length(unique(values(Auto)))
      pe += gen(base_field(A))
      Auto = Dict{GrpAbFinGenElem, NfRelElem{nf_elem}}(find_orbit(AutA_gen, AutA, pe))
    end
  end
  @vprint :ClassField 2 "have action on the primitive element!!!\n"
  return pe, Auto
end

function find_orbit(auts, AutG, x)
  @assert issnf(AutG)
  S = gens(AutG)
  t = ngens(AutG)
  order = 1
  elements = Tuple{GrpAbFinGenElem, NfRelElem{nf_elem}}[(id(AutG), x)]
  g = S[1]

  while !iszero(g)
    order = order +1
    push!(elements, (g, auts[1](elements[end][2])))
    g = g + S[1]
  end

  for i in 2:t
    previous_order = order
    for j = 1:previous_order
      order = order + 1
      push!(elements, (S[i]+elements[j][1], auts[i](elements[j][2])))
      for s in 2:Int(AutG.snf[i]-1)
        order = order + 1
        push!(elements, (elements[end][1] + S[i], auts[i](elements[end][2])))
      end
    end
  end
  return elements
end

function _aut_A_over_k(C::CyclotomicExt, CF::ClassField_pp)

  A = CF.K
#=
    now the automorphism group of A OVER k
    A = k(zeta, a^(1/n))
    we have
     tau: A-> A : zeta -> zeta   and a^(1/n) -> zeta a^(1/n)
   sigma: A-> A : zeta -> zeta^l and a^(1/n) -> (sigma(a))^(1/n)
                                                = c a^(s/n)
     for some c in k(zeta) and gcd(s, n) == 1
    Since A is compositum of the class field and k(zeta), A is abelian over
    k, thus sigma * tau = tau * sigma

    sigma * tau : zeta    -> zeta         -> zeta^l
                  a^(1/n) -> zeta a^(1/n) -> zeta^l c a^(s/n)
    tau * sigma : zeta    -> zeta^l       -> zeta^l
                  a^(1/n) -> c a^(s/n)    -> c zeta^s a^(s/n)

    Since A is abelian, the two need to agree, hence l==s and
    c can be computed as root(sigma(a)*a^-s, n)

    This has to be done for enough automorphisms of k(zeta)/k to generate
    the full group. If n=p^k then this is one (for p>2) and n=2, 4 and
    2 for n=2^k, k>2
=#
  e = degree(CF)
  g, mg = unit_group(ResidueRing(FlintZZ, e, cached=false))
  @assert issnf(g)
  @assert (e%8 == 0 && ngens(g)==2) || ngens(g) <= 1

  K = C.Ka
  Kr = C.Kr

  if degree(Kr) < order(g)  # there was a common subfield, we
                              # have to pass to a subgroup
    @assert order(g) % degree(Kr) == 0
    f = Kr.pol
    # Can do better. If the group is cyclic (e.g. if p!=2), we already know the subgroup!
    s, ms = sub(g, GrpAbFinGenElem[x for x in g if iszero(f(gen(Kr)^Int(lift(mg(x)))))], false)
    ss, mss = snf(s)
    g = ss
    #mg = mg*ms*mss
    mg = mss * ms * mg
  end

  @vprint :ClassField 2 "building automorphism group over ground field...\n"
  ng = ngens(g)+1
  AutA_gen = Array{Hecke.NfRelToNfRelMor_nf_elem_nf_elem, 1}(undef, ng)
  AutA_rel = zero_matrix(FlintZZ, ng, ng)
  zeta = C.mp[1]\(gen(Kr))
  n = degree(A)
  @assert e % n == 0

  @vprint :ClassField 2 "... the trivial one (Kummer)\n"
  tau = hom(A, A, zeta^div(e, n)*gen(A), check = false)
  AutA_gen[1] = tau

  AutA_rel[1,1] = n  # the order of tau

  @vprint :ClassField 2 "... need to extend $(ngens(g)) from the cyclo ext\n"
  for i = 1:ngens(g)
    si = hom(Kr, Kr, gen(Kr)^Int(lift(mg(g[i]))), check = false)
    @vprint :ClassField 2 "... extending zeta -> zeta^$(mg(g[i]))\n"
    to_be_ext = hom(K, K, C.mp[1]\(si(image(C.mp[1], gen(K)))), check = false)
    sigma = _extend_auto(A, to_be_ext, Int(lift(mg(g[i]))))
    AutA_gen[i+1] = sigma

    @vprint :ClassField 2 "... finding relation ...\n"
    m = gen(A)
    for j = 1:Int(order(g[i]))
      m = sigma(m)
    end
    # now m SHOULD be tau^?(gen(A)), so sigma^order(g[i]) = tau^?
    # the ? is what I want...
    j = 0
    zeta_i = zeta^mod(div(e, n)*(e-1), e)
    mi = coeff(m, 1)
    @hassert :ClassField 1 m == mi*gen(A)

    while mi != 1
      mul!(mi, mi, zeta_i)
      j += 1
      @assert j <= e
    end
    @vprint :ClassField 2 "... as tau^$(j) == sigma_$i^$(order(g[i]))\n"
    AutA_rel[i+1, 1] = -j
    AutA_rel[i+1, i+1] = order(g[i])
  end
  CF.AutG = AutA_gen
  CF.AutR = AutA_rel
  auts_in_snf!(CF)
  return nothing

end

function auts_in_snf!(CF::ClassField_pp)
  G = abelian_group(CF.AutR)
  S, mS = snf(G)
  auts = CF.AutG
  gens = Vector{NfRelToNfRelMor_nf_elem_nf_elem}(undef, ngens(S))
  for i = 1:ngens(S)
    el = mS(S[i])
    aut = id_hom(domain(CF.AutG[1]))
    for j = 1:length(auts)
      aut *= auts[j]^mod(el[j], exponent(S))
    end
    gens[i] = aut
  end
  CF.AutR = rels(S)
  CF.AutG = gens
  return nothing
end

struct ExtendAutoError <: Exception end

function _extend_auto(K::Hecke.NfRel{nf_elem}, h::Hecke.NfToNfMor, r::Int = -1)
  @hassert :ClassField 1 iskummer_extension(K)
  #@assert iskummer_extension(K)
  k = base_field(K)
  Zk = maximal_order(k)

  if r != -1
    if degree(K) == 2
      r = 1
    else
      #TODO: Do this modularly.
      zeta, ord = Hecke.torsion_units_gen_order(Zk)
      @assert ord % degree(K) == 0
      zeta = k(zeta)^div(ord, degree(K))
      im_zeta = h(zeta)
      r = 1
      z = deepcopy(zeta)
      while im_zeta != z
        r += 1
        mul!(z, z, zeta)
      end
    end
  end

  a = -coeff(K.pol, 0)
  dict = Dict{nf_elem, fmpz}()
  dict[h(a)] = 1
  if r <= div(degree(K), 2)
    add_to_key!(dict, a, -r)
    aa = FacElem(dict)
    @vtime :ClassField 3 fl, b = ispower(aa, degree(K), with_roots_unity = true)
    if !fl
      throw(ExtendAutoError())
    end
    return hom(K, K, h, evaluate(b)*gen(K)^r)
  else
    add_to_key!(dict, a, degree(K)-r)
    aa = FacElem(dict)
    @vtime :ClassField 3 fl, b = ispower(aa, degree(K), with_roots_unity = true)
    if !fl
      throw(ExtendAutoError())
    end
    return hom(K, K, h, evaluate(b)*gen(K)^(r-degree(K)))
  end
end

function _rcf_descent(CF::ClassField_pp)
  if isdefined(CF, :A)
    return nothing
  end

  @vprint :ClassField 2 "Descending ...\n"

  e = degree(CF)
  k = base_field(CF)
  CE = cyclotomic_extension(k, e)
  A = CF.K
  CK = absolute_field(CE)
  if degree(CK) == degree(k)
    #Relies on the fact that, if the cyclotomic extension has degree 1,
    #the absolute field is equal to the base field
    #There is nothing to do! The extension is already on the right field
    CF.A = A
    CF.pe = gen(A)
    return nothing
  end

  Zk = order(codomain(CF.rayclassgroupmap))
  ZK = maximal_order(absolute_field(CE))

  n = degree(A)
  #@vprint :ClassField 2 "Automorphism group (over ground field) $AutA\n"
  _aut_A_over_k(CE, CF)

  AutA_gen = CF.AutG
  AutA = abelian_group(CF.AutR)
  # now we need a primitive element for A/k
  # mostly, gen(A) will do

  @vprint :ClassField 2 "\nnow the fix group...\n"
  if iscyclic(AutA)  # the subgroup is trivial to find!
    @assert issnf(AutA)
    #Notice that this implies that the target field and the cyclotomic extension are disjoint.
    @vprint :ClassField 2 ".. trivial as automorphism group is cyclic\n"
    s, ms = sub(AutA, e, false)
    ss, mss = snf(s)
    ms = mss*ms
    gss = NfRelToNfRelMor_nf_elem_nf_elem[AutA_gen[1]^e]
    @vprint :ClassField 2 "computing orbit of primitive element\n"
    pe = gen(A)
    os = NfRelElem{nf_elem}[x[2] for x in find_orbit(gss, ss, pe)]
  else
    @vprint :ClassField 2 "Computing automorphisms of the extension and orbit of primitive element\n"
    pe, Auto = _find_prim_elem(CF, AutA)

    @vprint :ClassField 2 ".. interesting...\n"
    # want: hom: AutA = Gal(A/k) -> Gal(K/k) = domain(mq)
    # K is the target field.
    # idea: take primes p in k and compare
    # Frob(p, A/k) and preimage(mq, p)
    @assert n == degree(CF.K)

    local canFrob
    let CE = CE, ZK = ZK, n = n, pe = pe, Auto = Auto
      function canFrob(p::NfOrdIdl)
        lP = prime_decomposition(CE.mp[2], p)
        P = lP[1][1]
        F, mF = ResidueFieldSmall(ZK, P)
        Ft = PolynomialRing(F, cached = false)[1]
        mFp = extend_easy(mF, CE.Ka)
        ap = image(mFp, CF.a)
        @vprint :ClassField 1 "projection successful\n"
        polcoeffs = Vector{elem_type(F)}(undef, n+1)
        polcoeffs[1] = -ap
        for i = 2:n
          polcoeffs[i] = zero(F)
        end
        polcoeffs[n+1] = one(F)
        pol = Ft(polcoeffs)
        Ap = ResidueRing(Ft, pol, cached = false)
        xpecoeffs = Vector{elem_type(F)}(undef, n)
        for i = 0:n-1
          xpecoeffs[i+1] = image(mFp, coeff(pe, i))
        end
        xpe = Ft(xpecoeffs)
        imF = Ap(xpe)^norm(p)

        res = GrpAbFinGenElem[]
        for (ky, v) in Auto
          cfs = Vector{fq_nmod}(undef, n)
          for i = 0:n-1
            cfs[i+1] = image(mFp, coeff(v, i))
          end
          xp = Ft(cfs)
          kp = Ap(xp)
          if kp == imF
            push!(res, ky)
            if length(res) >1
              @vprint :ClassField 1 "res has length > 1\n"
              throw(BadPrime(p))
            end
          end
        end
        return res[1]
        error("Frob not found")
      end
    end
    @vprint :ClassField 2 "finding Artin map...\n"
    #TODO can only use non-indx primes, easy primes...
    cp = lcm([minimum(defining_modulus(CF)[1]), index(Zk), index(ZK)])
    #@vtime :ClassField 2 lp, f1 = find_gens(MapFromFunc(canFrob, IdealSet(Zk), AutA),
    #                  PrimesSet(200, -1), cp)
    lp, f1 = find_gens_descent(MapFromFunc(canFrob, IdealSet(Zk), AutA), CF, cp)
    imgs = GrpAbFinGenElem[image(CF.quotientmap, preimage(CF.rayclassgroupmap, p1)) for p1 = lp]
    h = hom(f1, imgs)
    @hassert :ClassField 1 issurjective(h)
    s, ms = kernel(h, false)
    @vprint :ClassField 2 "... done, have subgroup!\n"
    @vprint :ClassField 2 "computing orbit of primitive element\n"
    os = NfRelElem{nf_elem}[Auto[ms(j)] for j in s]
  end

  q, mq = quo(AutA, ms.map, false)
  sq, msq = snf(q)
  @assert iscyclic(sq)
  el_in_q = msq(sq[1])
  gsq = id_hom(A)
  for i = 1:ngens(q)
    if !iszero(el_in_q[i])
      gsq *= AutA_gen[i]^Int(el_in_q[i])
    end
  end
  genssq = NfRelToNfRelMor_nf_elem_nf_elem[gsq]

  @assert Int(order(q)) == degree(CF)

  #now, hopefully either norm or trace will be primitive for the target
  local charpoly
  inc_map = CE.mp[1]
  let inc_map = inc_map, AutA = AutA, e = e

    function charpoly(a::NfRelElem{nf_elem})
      tr_in_K = Vector{nf_elem}(undef, e)
      tr_err = divexact(order(AutA), e)
      el = one(parent(a))
      @vtime :ClassField 2 for i = 2:e+1
        mul!(el, el, a)
        tr_in_K[i-1] = divexact(tr(inc_map(tr(el))), tr_err)
      end
      res = power_sums_to_polynomial(tr_in_K)
      return res
    end
  end

  @vprint :ClassField 2 "trying relative trace\n"
  @assert length(os) > 0
  t = os[1]
  for i = 2:length(os)
    t += os[i]
  end
  CF.pe = t
  #now the minpoly of t - via Galois as this is easiest to implement...
  @vprint :ClassField 2 "char poly...\n"
  @vtime :ClassField 2 f2 = charpoly(t)
  @vprint :ClassField 2 "... done\n"

  if !issquarefree(f2)
    os1 = NfRelElem{nf_elem}[elem for elem in os]
    while !issquarefree(f2)
      @vprint :ClassField 2 "trying relative trace of powers\n"
      for i = 1:length(os)
        os1[i] *= os[i]
      end
      t = os1[1]
      for i = 2:length(os)
        t += os1[i]
      end
      CF.pe = t
      #now the minpoly of t - via Galois as this is easiest to implement...
      @vprint :ClassField 2 "min poly...\n"
      f2 = charpoly(t)
      @vprint :ClassField 2 "... done\n"
    end
  end
  CF.A = number_field(f2, check = false, cached = false)[1]
  return nothing
end

###############################################################################
#
#  Reduction of generators
#
###############################################################################

function _rcf_reduce(CF::ClassField_pp)
  e = CF.o
  if CF.sup_known
    CF.a = reduce_mod_powers(CF.a, e, CF.sup)
    CF.sup_known = false
  else
    CF.a = reduce_mod_powers(CF.a, e)
  end
  CF.K = radical_extension(CF.o, CF.a, check = false, cached = false)[1]
  return nothing
end

###############################################################################
#
#  Auxiliary functions (to be moved)
#
###############################################################################

@doc Markdown.doc"""
    ring_class_group(O::NfAbsOrd)

The ring class group (Picard group) of $O$.
"""
ring_class_group(O::NfAbsOrd) = picard_group(O)

@doc Markdown.doc"""
    ring_class_field(O::NfAbsOrd) -> ClassField

The ring class field of $O$, i.e. the maximal abelian extension ramified
only at primes dividing the conductor with the automorphism group
isomorphic to the Picard group.
"""
function ring_class_field(O::NfAbsOrd)
  M = maximal_order(O)
  f = extend(conductor(O), M)
  R, mR = ray_class_group(f)
  P, mP = picard_group(O)
  g, t = find_gens(mR)
  h = hom(t, [mP \ contract(x, O) for x = g], check = true)
  k = kernel(h, true)[1]
  q, mq = quo(R, k)
  return ray_class_field(mR, mq)
end


text:
export Shamanskii_LS_MA57

"""
A globalized Shamanskii algorithm with Line Search.
This code is an implementation of the idea presented in:

    Global Convergence Technique for the Newton
        Method with Periodic Hessian Evaluation

By F. LAMPARIELLO and M. SCIANDRONE

JOURNAL OF OPTIMIZATION THEORY AND APPLICATIONS:
Vol. 111, No. 2, pp. 341–358, November 2001
"""
function Shamanskii_LS_MA57(nlp          :: AbstractNLPModel,
                 	   		nlp_stop     :: NLPStopping;
                 	   		linesearch   :: Function = shamanskii_line_search,
                 	   		Nwtdirection :: Function = NwtdirectionLDLT_LS,
                 	   		hessian_rep  :: Function = hessian_sparse,
					   		η 		    :: Float64 = 1.5,
					   		ca           :: Float64 = 0.1,
					   		cf           :: Float64 = 0.25,
					   		u            :: Int = 0,
					   		p            :: Int = 1,
                 	   		kwargs...)

    nlp_at_x = nlp_stop.current_state
    n = nlp.meta.nvar
    xt = copy(nlp.meta.x0)
	T = eltype(xt)
    ∇ft = Array{Float64}(undef, n)

    f = obj(nlp, xt)
    ∇f = grad(nlp, xt)
	nlp_stop.meta.optimality0 = norm(copy(∇f))
    OK = update_and_start!(nlp_stop, x = xt, fx = f, gx = ∇f)#, g0 = ∇f)
    ∇fNorm = norm(nlp_at_x.gx)

	global Hhat = hessian_rep(nlp, xt)
    global iter = 0
	global i = 0
	global fm1 = nothing

	global M = Ma57
	global L = SparseMatrixCSC{Float64, Int64}
	global D = SparseMatrixCSC{Float64, Int64}
	global pp = Array{Int64, 1}
	global s = Array{Float64}
	global ρ = nothing
	global ncomp = nothing
	global Q = nothing
	global Γ = nothing
	global ϵ2 = nothing

    @info log_header([:iter, :f, :nrm_g, :α], [Int64, T, T, T])
	@info log_row(Any[iter, nlp_at_x.fx, ∇fNorm, 0.0])

    β = 0.0
    d = zero(nlp_at_x.gx)
    h = LineSearch.LineModel(nlp, xt, d)

    while !OK
		if !(iter == i * p) && (u == 0)
			# println("we do nothing")
		else
			# println("we do something")
			H57 = convert(SparseMatrixCSC{Float64, Int64}, hessian_rep(nlp, xt))
			printstyled("on a H57 \n", color = :yellow)
			try
		        M = Ma57(H57)#, print_level = -1)
		        ma57_factorize(M)
		    catch
		 	println("*******   Problem in MA57_0")
		        M = Ma57(H57,print_level = -1)
		        ma57_factorize(M)
		        res = PDataMA57()
		        res.OK = false
		        return res
		    end
			printstyled("on a M \n", color = :yellow)
			try
		        (L, D, s, pp) = ma57_get_factors(M)
		    catch
		        println("*******   Problem after MA57_1")
		        # println(" Cond(H) = $(cond(full(H)))")
		        res = PDataMA57_0()
		        res.OK = false
		        return res
		    end
			printstyled("on a L, D, s, pp \n", color = :yellow)
			#################  Future object BlockDiag operator?
		    vD1 = diag(D)       # create internal representation for block diagonal D
		    vD2 = diag(D, 1)     #
		    vQ1 = ones(length(vD1))       # vector representation of orthogonal matrix Q
		    vQ2 = zeros(length(vD2))      #
		    vQ2m = zeros(length(vD2))     #
		    veig = copy(vD1)      # vector of eigenvalues of D, initialized to diagonal of D
		                          # if D diagonal, nothing more will be computed
			printstyled("on a vD1, vD2, vQ1, vQ2, vQ2m, veig \n", color = :yellow)
		    i = 1;
		    while i < length(vD1)
		        if vD2[i] == 0.0
		            i += 1
		        else
		            mA = [vD1[i] vD2[i]; vD2[i] vD1[i + 1]] #  2X2 submatrix
		            # DiagmA, Qma = eig(mA)                   #  spectral decomposition of mA
		            X = eigen(mA)
		            DiagmA = X.values
		            Qma = X.vectors
		            veig[i] = DiagmA[1]
		            vQ1[i] = Qma[1, 1]
		            vQ2[i] = Qma[1, 2]
		            vQ2m[i] = Qma[2, 1]
		            vQ1[i + 1] = Qma[2, 2]
		            veig[i + 1] = DiagmA[2]
		            i += 2
		        end
		    end

		    Q = sparse(SparseArrays.spdiagm(0 => vQ1, 1 => vQ2m, -1 => vQ2))
			printstyled("on a Q \n", color = :yellow)
		    Δ = veig
			printstyled("on a  Δ \n", color = :yellow)
			ϵ2 = sqrt(eps(T))
			printstyled("on a  ϵ2 \n", color = :yellow)
			Γ = max.(abs.(Δ), ϵ2)
			printstyled("on a  Γ \n", color = :yellow)
			@show Γ
			@show Γ .== ϵ2
			if (true in (Δ .== ϵ2))
			   u = 1
			else
			   u = 0
			end
			if (iter == i * p)
				i += 1
			end
		end
		# @show L
		# @show D
		# @show Q
		# @show pp
		# @show Γ
		# @show ϵ2
		@show u
        # d, u = Nwtdirection(Hhat, nlp_at_x.gx)
		g = copy(nlp_at_x.gx)
		@show eltype(g)
		@show L
		@show g
		@show g[pp]
		d̃ = L \ g[pp]
		@show eltype(d̃)
		d̂ = L' \ (Q * (Q' * d̃ ./ Γ))
		@show eltype(d̂)
		d = -d̂[invperm(pp)]
		@show eltype(d)

		# @show u
        slope = d' * nlp_at_x.gx

        h = LineSearch.redirect!(h, xt, d)
		printstyled("on a  h \n", color = :yellow)
        ls_at_t = LSAtT(0.0, h₀ = nlp_at_x.fx, g₀ = slope)
		printstyled("on a  ls_at_t \n", color = :yellow)
        stop_ls = LS_Stopping(h, (x, y) -> armijo(x, y, τ₀ = 1e-09), ls_at_t)
		printstyled("on a  stop_ls \n", color = :yellow)
        ls_at_t, good_step_size = linesearch(h, stop_ls; kwargs...)
		printstyled("on a  fini le line search \n", color = :yellow)
        xt  = nlp_at_x.x + ls_at_t.x * d
		printstyled("on a  xt \n", color = :yellow)
		fm1 = copy(nlp_at_x.fx)
        ft  = obj(nlp, xt); ∇ft = grad(nlp, xt)

        ∇ft = grad!(nlp, xt, ∇ft)
		# @show eltype(∇ft)
		# @show eltype(xt)
		# @show eltype(ft)
		OK = update_and_stop!(nlp_stop, x = xt, fx = ft, gx = ∇ft)

		# @show ls_at_t.x
		# @show ca
		# @show fm1
		# @show ft
		# @show cf * abs(ft) + fm1
		# @show cf * abs(ft) + fm1 < ft

		if (fm1 < ft)
			u = 1
		else
			u = 0
		end

        # Move on.
        ∇fNorm = norm(nlp_at_x.gx)
        iter = iter + 1

		@info log_row(Any[iter, nlp_at_x.fx, ∇fNorm, ls_at_t.x])
    end

    return nlp_stop, nlp_stop.meta.optimal
end

ext: jl


text:
#=
- (can't subtype, use fn) Make JsonResponse as a subtype of HTTP response
- Create logger
- macro w/in macro => HTTP.@register (working for now)
- Handle route path variables
-
=#
module QuickAPI

import HTTP
import URIs
import Sockets
import JSON3
import MsgPack
import Dates

const APP = HTTP.Router()

const GET    = :GET
const POST   = :POST
const PUT    = :PUT
const DELETE = :DELETE

const LOCALHOST = Sockets.localhost

const CONTENT_TYPES = Dict{Symbol,String}(
  :text       => "text/plain",
  :json       => "application/json",
  :msgpack    => "application/x-msgpack",
)

# Functions for Response Types

function text(data::Any ; status::Int = 200, header::Array{Pair{String,String},1} = Pair{String,String}[])
    HTTP.Response(
        status,
        ["Content-Type" => CONTENT_TYPES[:text], header...];
        body = string(data) * "\n"
    )
end

function json(data::Any ; status::Int = 200, header::Array{Pair{String,String},1} = Pair{String,String}[])
    HTTP.Response(
        status,
        ["Content-Type" => CONTENT_TYPES[:json], header...];
        body = JSON3.write(data) * "\n"
    )
end

function msgpack(data::Any ; status::Int = 200, header::Array{Pair{String,String},1} = Pair{String,String}[])
    HTTP.Response(
        status,
        ["Content-Type" => CONTENT_TYPES[:msgpack], header...];
        body = MsgPack.pack(data)
    )
end

# Macros for handling GET, POST,
# PUT, DELETE requests

macro get(path::String, response::Expr)
    HTTP.@register(
        APP,
        :GET,
        path,
        response
    )
end

macro post(path::String, response::Expr)
    HTTP.@register(
        APP,
        :POST,
        path,
        response
    )
end

macro put(path::String, response::Expr)
    HTTP.@register(
        APP,
        :PUT,
        path,
        response
    )
end

macro delete(path::String, response::Expr)
    HTTP.@register(
        APP,
        :DELETE,
        path,
        response
    )
end

# Logger function
function log_request()

end

# Functions to Serve the App

function serve(host=LOCALHOST, port=8081; kw...)
    @info "Serving at http://$host:$port/"
    try
        HTTP.serve(
            APP,
            host,
            port;
            kw...
        )
    catch e
        @error e
    end
end

function serve_async(host=LOCALHOST, port=8081; kw...)
    @async HTTP.serve(APP,host,port,kw...)
end

export
    APP,
    LOCALHOST,
    text, json, msgpack,
    @get, @post, @put, @delete,
    serve, serve_async

end # module

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

# using HTTP
# using JSON
using .QuickAPI
using Sockets

@info "Starting module..."

@get "/" (r::HTTP.Request -> begin
    json(Dict(
        "hello" => "world",
        "dogs" => 3
    ))
end)

serve(
    LOCALHOST,
    8081,
    tcpisvalid = sock::Sockets.TCPSocket -> begin
        @show sock
        @show sock.buffer
        true
    end,
)

@info "Done."


text:
using Mocha

data_tr_layer = HDF5DataLayer(name="data-train", source="data/train.txt", batch_size=100, shuffle=@windows ? false : true)
data_tt_layer = HDF5DataLayer(name="data-test", source="data/test.txt", batch_size=100)

conv1_layer = ConvolutionLayer(name="conv1", n_filter=32, kernel=(5,5), pad=(2,2),
    stride=(1,1), filter_init=GaussianInitializer(std=0.0001),
    bottoms=[:data], tops=[:conv1])
pool1_layer = PoolingLayer(name="pool1", kernel=(3,3), stride=(2,2), neuron=Neurons.ReLU(),
    bottoms=[:conv1], tops=[:pool1])
norm1_layer = LRNLayer(name="norm1", kernel=3, scale=5e-5, power=0.75, mode=LRNMode.WithinChannel(),
    bottoms=[:pool1], tops=[:norm1])
conv2_layer = ConvolutionLayer(name="conv2", n_filter=32, kernel=(5,5), pad=(2,2),
    stride=(1,1), filter_init=GaussianInitializer(std=0.01),
    bottoms=[:norm1], tops=[:conv2], neuron=Neurons.ReLU())
pool2_layer = PoolingLayer(name="pool2", kernel=(3,3), stride=(2,2), pooling=Pooling.Mean(),
    bottoms=[:conv2], tops=[:pool2])
norm2_layer = LRNLayer(name="norm2", kernel=3, scale=5e-5, power=0.75, mode=LRNMode.WithinChannel(),
    bottoms=[:pool2], tops=[:norm2])
conv3_layer = ConvolutionLayer(name="conv3", n_filter=64, kernel=(5,5), pad=(2,2),
    stride=(1,1), filter_init=GaussianInitializer(std=0.01),
    bottoms=[:norm2], tops=[:conv3], neuron=Neurons.ReLU())
pool3_layer = PoolingLayer(name="pool3", kernel=(3,3), stride=(2,2), pooling=Pooling.Mean(),
    bottoms=[:conv3], tops=[:pool3])
ip1_layer   = InnerProductLayer(name="ip1", output_dim=10, weight_init=GaussianInitializer(std=0.01),
    weight_regu=L2Regu(250), bottoms=[:pool3], tops=[:ip1])

loss_layer  = SoftmaxLossLayer(name="softmax", bottoms=[:ip1, :label])
acc_layer   = AccuracyLayer(name="accuracy", bottoms=[:ip1, :label])

common_layers = [conv1_layer, pool1_layer, norm1_layer, conv2_layer, pool2_layer, norm2_layer,
                 conv3_layer, pool3_layer, ip1_layer]

backend = DefaultBackend()
init(backend)

net = Net("CIFAR10-train", backend, [data_tr_layer, common_layers..., loss_layer])

lr_policy = LRPolicy.Staged(
  (60000, LRPolicy.Fixed(0.001)),
  (5000, LRPolicy.Fixed(0.0001)),
  (5000, LRPolicy.Fixed(0.00001)),
)
method = SGD()
solver_params = make_solver_parameters(method, max_iter=70000,
    regu_coef=0.004, mom_policy=MomPolicy.Fixed(0.9), lr_policy=lr_policy,
    load_from="snapshots")
solver = Solver(method, solver_params)

# report training progress every 200 iterations
add_coffee_break(solver, TrainingSummary(), every_n_iter=200)

# save snapshots every 5000 iterations
add_coffee_break(solver, Snapshot("snapshots"), every_n_iter=5000)

# show performance on test data every 1000 iterations
test_net = Net("CIFAR10-test", backend, [data_tt_layer, common_layers..., acc_layer])
add_coffee_break(solver, ValidationPerformance(test_net), every_n_iter=1000)

solve(solver, net)

destroy(net)
destroy(test_net)
shutdown(backend)


text:
module Models

export Model, Mixture, HiddenMarkovModel
export ArrivalDistribution, GeometricArrivals, PoissonArrivals
export MixtureParameters, NtlParameters, DpParameters
export DataParameters, GaussianParameters, MultinomialParameters
export GaussianWishartParameters
export ClusterSufficientStatistics, DataSufficientStatistics
export NtlSufficientStatistics, GaussianSufficientStatistics, MultinomialParameters
export HmmSufficientStatistics, MixtureSufficientStatistics
export StationaryHmmSufficientStatistics, NonstationaryHmmSufficientStatistics
export ChangepointSufficientStatistics, Changepoint
export ParametricArrivalsClusterParameters
export PitmanYorArrivals
export SufficientStatistics
export AuxillaryVariables, ArrivalsAuxillaryVariables, DataAuxillaryVariables
export GeometricAuxillaryVariables, PoissonAuxillaryVariables, DpAuxillaryVariables
export GaussianAuxillaryVariables, GaussianWishartAuxillaryVariables, MultinomialAuxillaryVariables
export FiniteTopicModelAuxillaryVariables

using SparseArrays
using DataStructures

abstract type DataParameters end

struct GaussianParameters <: DataParameters
    data_covariance::Matrix{Float64}
    data_precision::Matrix{Float64}
    prior_mean::Vector{Float64}
    prior_covariance::Matrix{Float64}
    dim::Int64
    function GaussianParameters(data_covariance::Matrix{Float64},
                                prior_mean::Vector{Float64},
                                prior_covariance::Matrix{Float64})
        if size(data_covariance) !== size(prior_covariance)
            error("Dimension of data and prior covariance must match")
        end
        dim = length(prior_mean)
        if dim !== size(prior_covariance)[1] || dim !== size(prior_covariance)[2]
            error("Dimension of mean prior must match each dimension of prior covariance")
        end
        data_precision = inv(data_covariance)
        return new(data_covariance, data_precision, prior_mean, prior_covariance, dim)
    end
end

struct GaussianWishartParameters <: DataParameters
    prior_mean::Vector{Float64}
    scale_matrix::Matrix{Float64}
    scale::Float64
    dof::Float64
    dim::Int64
    function GaussianWishartParameters(prior_mean, scale_matrix, scale, dof)
        dim = length(prior_mean)
        return new(prior_mean, scale_matrix, scale, dof, dim)
    end
end

struct MultinomialParameters <: DataParameters
    n::Int64
    dim::Int64
    prior_dirichlet_scale::Vector
    function MultinomialParameters(n::Int64, prior_dirichlet_scale::Vector)
        dim = length(prior_dirichlet_scale)
        if dim < 2
            error("Dirichlet distribution must have dimension at least 2")
        end
        if any(prior_dirichlet_scale .<= 0)
            error("Each value of Dirichlet scale parameter must be greater than 0")
        end
        if n <= 0
            error("n must be greater than 0")
        end
        return new(n, dim, prior_dirichlet_scale)
    end
    function MultinomialParameters(prior_dirichlet_scale::Vector)
        return MultinomialParameters(1, prior_dirichlet_scale)
    end
end

struct FiniteTopicModelParameters <: DataParameters
    num_topics::Int64
    num_words::Int64
    length::Int64
    topic_parameter::Float64
    word_parameter::Float64
    function FiniteTopicModelParameters(;num_topics::Int64=10, num_words::Int64=100, length::Int64=100,
                                        topic_parameter::Float64=0.1, word_parameter::Float64=0.1)
        return new(num_topics, num_words, length, topic_parameter, word_parameter)
    end
end

abstract type ArrivalDistribution end

struct GeometricArrivals <: ArrivalDistribution
    prior::Vector{Float64}
    sample_parameter_posterior::Bool
    function GeometricArrivals(;prior=[1,1], sample_parameter_posterior=false)
        new(prior, sample_parameter_posterior)
    end
end

struct PoissonArrivals <: ArrivalDistribution
    alpha::Float64
    beta::Float64
    sample_parameter_posterior::Bool
    function PoissonArrivals(;alpha=1, beta=10, sample_parameter_posterior=false)
        new(alpha, beta, sample_parameter_posterior)
    end
end

struct PitmanYorArrivals <: ArrivalDistribution
    tau::Float64
    theta::Float64
    function PitmanYorArrivals(;tau::Float64=0., theta::Float64=1.)
        if tau < 0 || tau > 1
            error("tau must be in (0,1)")
        end
        if theta <= -tau
            error("theta must be greater than -tau")
        end
        return new(tau, theta)
    end
end

abstract type ParametricArrivalsClusterParameters end

mutable struct NtlParameters <: ParametricArrivalsClusterParameters
    prior::Vector{Float64}
    arrival_distribution::GeometricArrivals
    sample_parameter_posterior::Bool
    function NtlParameters(prior::Vector{Float64}, arrival_distribution::GeometricArrivals; sample_parameter_posterior::Bool=false)
        return new(prior, arrival_distribution, sample_parameter_posterior)
    end
end

abstract type ClusterParameters end

mutable struct DpParameters <: ClusterParameters
    beta::Float64
    alpha::Float64
    sample_parameter_posterior::Bool
    function DpParameters(beta, alpha; sample_parameter_posterior=false)
        return new(beta, alpha, sample_parameter_posterior)
    end
    function DpParameters(alpha; sample_parameter_posterior=false)
        return new(0, alpha, sample_parameter_posterior)
    end
end

abstract type ClusterSufficientStatistics end

struct MixtureSufficientStatistics <: ClusterSufficientStatistics
    num_observations::Vector{Int64}
    cumulative_num_observations::FenwickTree{Int64}
    clusters::BitArray
end

abstract type HmmSufficientStatistics <: ClusterSufficientStatistics end

# Stationary HMM
mutable struct StationaryHmmSufficientStatistics <: HmmSufficientStatistics
    num_observations::Vector{Int64}
    state_num_observations::Matrix{Int64}
    clusters::BitArray
end

# Nonstationary HMM
mutable struct NonstationaryHmmSufficientStatistics <: HmmSufficientStatistics
    num_observations::Vector{Int64}
    state_num_observations::Matrix{Int64}
    state_assignments::Vector{SparseMatrixCSC}
    clusters::BitArray
end

struct ChangepointSufficientStatistics <: ClusterSufficientStatistics
    num_observations::Vector{Int64}
    clusters::BitArray
end

abstract type DataSufficientStatistics end

struct FiniteTopicModelSufficientStatistics <: DataSufficientStatistics
    cluster_topic_frequencies::Matrix{Int64}
    cluster_topic_posterior::Matrix{Float64}
    topic_token_frequencies::Matrix{Int64}
    topic_token_posterior::Matrix{Float64}
end

struct GaussianSufficientStatistics <: DataSufficientStatistics
    data_means::Matrix{Float64}
    data_precision_quadratic_sums::Vector{Float64}
    posterior_means::Matrix{Float64}
    posterior_covs::Matrix{Float64}
    function GaussianSufficientStatistics(data_means::Matrix{Float64},
                                          data_precision_quadratic_sum::Vector{Float64},
                                          posterior_means::Matrix{Float64},
                                          posterior_covs::Matrix{Float64})
        # Make sure dimensions match
        if size(data_means)[1] !== size(posterior_means)[1]
            error("Dimension of data means and posterior means do not match")
        end
        dim = size(posterior_means)[1]
        if size(posterior_covs)[1] !== dim^2
            error("Dimension of covariances must be square of dimension of means")
        end
        return new(data_means, data_precision_quadratic_sum, posterior_means, posterior_covs)
    end
end

struct GaussianWishartSufficientStatistics <: DataSufficientStatistics
    data_means::Matrix{Float64}
    data_outer_product_sums::Matrix{Float64}
    posterior_means::Matrix{Float64}
    posterior_scale_matrix::Matrix{Float64}
    posterior_scale::Vector{Float64}
    posterior_dof::Vector{Float64}
end

struct MultinomialSufficientStatistics <: DataSufficientStatistics
    total_counts::Matrix{Int64}
    posterior_dirichlet_scale::Matrix{Float64}
    function MultinomialSufficientStatistics(total_counts::Matrix{Int64}, posterior_dirichlet_scale::Matrix{Float64})
        if any(total_counts .< 0)
            error("Each count in vector must be at least 0")
        end
        if any(posterior_dirichlet_scale .<= 0)
            error("Each value in dirichlet scale posterior must be greater than 0")
        end
        return new(total_counts, posterior_dirichlet_scale)
    end
end

mutable struct SufficientStatistics{C<:ClusterSufficientStatistics, D<:DataSufficientStatistics}
    n::Int
    num_assigned::Int64
    cluster::C
    data::D
end

abstract type Model end

mutable struct Mixture{C<:Union{ClusterParameters, ParametricArrivalsClusterParameters, DpParameters}, D<:DataParameters} <: Model
    cluster_parameters::C
    data_parameters::D
end

struct HiddenMarkovModel{C<:Union{ClusterParameters, ParametricArrivalsClusterParameters}, D<:DataParameters} <: Model
    cluster_parameters::C
    data_parameters::D
end

struct Changepoint{C<:Union{ClusterParameters, ParametricArrivalsClusterParameters}, D<:DataParameters} <: Model
    cluster_parameters::C
    data_parameters::D
end

abstract type ArrivalsAuxillaryVariables end

struct GeometricAuxillaryVariables <: ArrivalsAuxillaryVariables end

mutable struct PoissonAuxillaryVariables <: ArrivalsAuxillaryVariables
    last_arrival_time::Int64
end

struct DpAuxillaryVariables <: ArrivalsAuxillaryVariables end

abstract type DataAuxillaryVariables end

struct GaussianAuxillaryVariables <: DataAuxillaryVariables end

struct MultinomialAuxillaryVariables <: DataAuxillaryVariables end

struct GaussianWishartAuxillaryVariables <: DataAuxillaryVariables end

struct FiniteTopicModelAuxillaryVariables <: DataAuxillaryVariables
    token_topic_assignments::Vector{Dict{Int64, Dict{Int64, Int64}}}
    document_topic_frequencies::Matrix
end

struct AuxillaryVariables{A <: ArrivalsAuxillaryVariables, D <: DataAuxillaryVariables}
    arrivals::A
    data::D
end

end

max_stars_repo_name: realseanla/ntl-stick-breaking-julia

max_issues_repo_name: realseanla/ntl-stick-breaking-julia

max_forks_repo_name: realseanla/ntl-stick-breaking-julia


text:
using(DelimitedFiles)

input = DelimitedFiles.readdlm("Day 5/input_day5.txt")

# part 1

function passport(x)

parsed_binary = reduce(replace, ["B"=>1, "F"=>0, "L"=>0, "R"=>1], init=x)
row = parse(Int, parsed_binary[1:7], base = 2)
column = parse(Int, parsed_binary[8:10], base = 2)
seatid = row * 8 + column
return seatid

end

@time begin

seatids = passport.(input) #using . to vectorize bc julia is op

maximum(seatids)

end

# part 2

@time begin

seat_range = minimum(seatids):maximum(seatids)

seat_diff = setdiff(seat_range, seatids)

end

max_stars_repo_path: Day 5/advent2020_day5.jl

max_issues_repo_path: Day 5/advent2020_day5.jl

max_forks_repo_path: Day 5/advent2020_day5.jl


text:
using Test
using Algorithm101
const ALGO = Algorithm101

@testset "Test BinarySearch" begin

    # Test: 1:7:1000 sorted array
    array = collect(1:7:1000)
    target_i = rand(1:length(array))
    target = array[target_i]
    @test ALGO.binary_search(array, target) == target_i

    # Test: If element is not in the array
    target = 6
    @test ALGO.binary_search(array, target) == -1

end


text:
"""

# Method `update_model_file`

Update Stan language model file if necessary

### Method
```julia
CmdStan.update_model_file(
  file::String,
  str::String
)
```
### Required arguments
```julia
* `file::AbstractString`                : File holding existing Stan model
* `str::AbstractString`                 : Stan model string
```

### Related help
```julia
?CmdStan.Stanmodel                 : Create a StanModel
```
"""
function update_model_file(file::AbstractString, str::AbstractString)
  str2 = ""
  if isfile(file)
    resfile = open(file, "r")
    str2 = read(resfile, String)
    str != str2 && rm(file)
  end
  if str != str2
    println("\nFile $(file) will be updated.\n")
    strmout = open(file, "w")
    write(strmout, str)
    close(strmout)
  end
end


text:
using Lint

lintpkg("Lexicon")


text:
function build_cells(::CRISPRi,
                     guides::Vector{Barcode},
                     guide_freq_dist::Categorical,
                     n::Int
                    )
    cells = rand(guide_freq_dist, n)
    phenotypes = Array{Float64}(n)
    @inbounds @fastmath for i in 1:n
        phenotypes[i] = guides[cells[i]].theo_phenotype
    end
    cells, phenotypes
end

function build_cells(behav::CRISPRn,
                     guides::Vector{Barcode},
                     guide_freq_dist::Categorical,
                     n::Int
                    )

    phenotypes = Array{Float64}(n);
    cells = rand(guide_freq_dist, n)
    ko_dist = behav.knockout_dist
    dist = rand(ko_dist, n)
    @inbounds @fastmath for i in 1:n
        phenotypes[i] = guides[cells[i]].theo_phenotype*(dist[i] - 1)/2
    end
    cells, phenotypes
end

"""
$(SIGNATURES)

Simulates the transfection the of the library of guides (`guides`) according to
the parameters in [`Simulation.Library`](@ref) `lib` and the frequencies
represented by the Categorical distribution `guide_freqs_dist`.

This function creates a population of cells post-transfection where each cell
has single guide RNA present according to the frequencies of the guides in the
present in the library. The cells are then grown up so that there are
`setup.bottleneck_representation * num_guides` of them. For FACS screens we
assume a simple linear expansion of the cells, while for Growth screens the
cells grow according as a function of their phenotype and degree of knockdown.
"""
function transfect(setup::FacsScreen,
                   lib::Library,
                   guides::Vector{Barcode},
                   guide_freqs_dist::Categorical)

    moi = setup.moi
    num_guides = length(guides)
    cell_count = num_guides * setup.representation
    expand_to = setup.bottleneck_representation * length(guides)

    cells, cell_phenotypes = build_cells(lib.cas9_behavior, guides, guide_freqs_dist,
                                         round(Int, pdf(Poisson(moi), 1)*cell_count))
    num_cells = length(cells)

    multiples = 1
    if expand_to > num_cells
        multiples = ceil(Int, expand_to/num_cells)
        expansion_c = Array{Int}(num_cells*multiples)
        expansion_p = Array{Float64}(num_cells*multiples)
        for rep in 1:multiples
            rng = (rep-1)*num_cells+1:rep*num_cells
            expansion_c[rng] = cells
            expansion_p[rng] = cell_phenotypes
        end
        cells, cell_phenotypes = expansion_c, expansion_p
    else
        picked = sample(collect(1:num_cells), expand_to, replace=false)
        cells, cell_phenotypes = cells[picked], cell_phenotypes[picked]
    end

    initial_freqs = StatsBase.counts(cells, 1:num_guides) ./ length(cells)

    for i in 1:num_guides
        @inbounds guides[i].initial_freq = initial_freqs[i]
    end
    cells, cell_phenotypes
end

function transfect(setup::GrowthScreen,
                   lib::Library,
                   guides::Vector{Barcode},
                   guide_freqs_dist::Categorical)

    num_guides = length(guides)
    cell_count = num_guides * setup.representation
    initial_cells, cell_phenotypes = build_cells(lib.cas9_behavior, guides, guide_freqs_dist,
                                         round(Int, pdf(Poisson(setup.moi), 1)*cell_count))
    target = num_guides * setup.bottleneck_representation

    if target < length(initial_cells)
        picked = sample(collect(1:length(initial_cells)), target, replace=false)
        cells, cell_phenotypes = initial_cells[picked], cell_phenotypes[picked]
        num_doublings = -1
    else
        cells, cell_phenotypes = copy(initial_cells), copy(cell_phenotypes)
        num_inserted = length(cells)
        num_doublings = 0
        output_c = Array{Int}(target*4)
        output_p = Array{Float64}(target*4)

        while num_inserted < target
            num_inserted = grow!(cells, cell_phenotypes, output_c, output_p, setup)
            cells = copy(view(output_c, 1:num_inserted))
            cell_phenotypes = copy(view(output_p, 1:num_inserted))
            num_doublings += 1
        end
    end

    initial_freqs = counts(cells, 1:length(guides)) ./ length(cells)

    for i in 1:length(guides)
        @inbounds guides[i].initial_freq = initial_freqs[i]
    end
    cells, cell_phenotypes
end


text:
# Copy from https://github.com/JuliaLang/julia/issues/4648#issuecomment-761030051
@generated structEqual(e1, e2) = begin
    if fieldcount(e1) == 0
        return :(true)
    end
    mkEq = fldName -> :(e1.$fldName == e2.$fldName)
    # generate individual equality checks
    eqExprs = map(mkEq, fieldnames(e1))
    # construct &&-expression for chaining all checks
    mkAnd = (expr, acc) -> Expr(:&&, expr, acc)
    # no need in initial accumulator because eqExprs is not empty
    foldr(mkAnd, eqExprs)
end

Base.:(==)(x::m.Node, y::m.Node) = structEqual(x, y)

@testset "Test Modifying AST" begin
    one = () -> m.IntegerLiteral(m.Token(m.INT, "1"), 1)
    two = () -> m.IntegerLiteral(m.Token(m.INT, "2"), 2)
    turn_one_into_two =
        (node::m.Node) -> isa(node, m.IntegerLiteral) && node.value == 1 ? two() : node

    for (input, expected) in [
        (one(), two()),
        (
            m.Program([m.ExpressionStatement(one().token, one())]),
            m.Program([m.ExpressionStatement(one().token, two())]),
        ),
        (
            m.InfixExpression(m.Token(m.PLUS, "+"), one(), "+", two()),
            m.InfixExpression(m.Token(m.PLUS, "+"), two(), "+", two()),
        ),
        (
            m.InfixExpression(m.Token(m.PLUS, "+"), two(), "+", one()),
            m.InfixExpression(m.Token(m.PLUS, "+"), two(), "+", two()),
        ),
        (
            m.PrefixExpression(m.Token(m.MINUS, "-"), "-", one()),
            m.PrefixExpression(m.Token(m.MINUS, "-"), "-", two()),
        ),
        (
            m.IndexExpression(one().token, one(), one()),
            m.IndexExpression(one().token, two(), two()),
        ),
        (
            m.IfExpression(
                m.Token(m.IF, "if"),
                one(),
                m.BlockStatement(one().token, [m.ExpressionStatement(one().token, one())]),
                nothing,
            ),
            m.IfExpression(
                m.Token(m.IF, "if"),
                two(),
                m.BlockStatement(one().token, [m.ExpressionStatement(one().token, two())]),
                nothing,
            ),
        ),
        (
            m.IfExpression(
                m.Token(m.IF, "if"),
                one(),
                m.BlockStatement(one().token, [m.ExpressionStatement(one().token, one())]),
                m.BlockStatement(one().token, [m.ExpressionStatement(one().token, one())]),
            ),
            m.IfExpression(
                m.Token(m.IF, "if"),
                two(),
                m.BlockStatement(one().token, [m.ExpressionStatement(one().token, two())]),
                m.BlockStatement(one().token, [m.ExpressionStatement(one().token, two())]),
            ),
        ),
        (
            m.ReturnStatement(m.Token(m.RETURN, "return"), one()),
            m.ReturnStatement(m.Token(m.RETURN, "return"), two()),
        ),
        (
            m.LetStatement(
                m.Token(m.IDENT, "x"),
                m.Identifier(m.Token(m.IDENT, "x"), "x"),
                one(),
                false,
            ),
            m.LetStatement(
                m.Token(m.IDENT, "x"),
                m.Identifier(m.Token(m.IDENT, "x"), "x"),
                two(),
                false,
            ),
        ),
        (
            m.FunctionLiteral(
                m.Token(m.FUNCTION, "fn"),
                [m.Identifier(m.Token(m.IDENT, "x"), "x")],
                m.BlockStatement(one().token, [m.ExpressionStatement(one().token, one())]),
            ),
            m.FunctionLiteral(
                m.Token(m.FUNCTION, "fn"),
                [m.Identifier(m.Token(m.IDENT, "x"), "x")],
                m.BlockStatement(one().token, [m.ExpressionStatement(one().token, two())]),
            ),
        ),
        (
            m.ArrayLiteral(one().token, [one(), one()]),
            m.ArrayLiteral(one().token, [two(), two()]),
        ),
        (
            m.HashLiteral(one().token, Dict(one() => one())),
            m.HashLiteral(one().token, Dict(two() => two())),
        ),
    ]
        @test begin
            modified = m.modify(input, turn_one_into_two)
            modified == expected
        end
    end
end


text:
using DashBootstrapComponents, DashHtmlComponents

function make_button(placement)
    return dbc_button(
        "Tooltip on $placement",
        id = "tooltip-target-$placement",
        className = "mx-2",
        n_clicks = 0,
    )
end;

function make_tooltip(placement)
    return dbc_tooltip(
        "Tooltip on $placement",
        target = "tooltip-target-$placement",
        placement = placement,
    )
end;

tooltips = html_div(
    [
        [make_button(p) for p in ["top", "left", "right", "bottom"]]
        [make_tooltip(p) for p in ["top", "left", "right", "bottom"]]
    ],
);

max_stars_repo_name: glsdown/dash-bootstrap-components

max_issues_repo_name: glsdown/dash-bootstrap-components

max_forks_repo_name: glsdown/dash-bootstrap-components


text:
using Ion
using Git
using Ion.Internal
using Test

@testset "is_version_number" begin
    @test Internal.is_version_number("0.1.2") == true
    @test Internal.is_version_number("0.1.x") == false
end

# clone some example repos
Ion.Internal.clone("Comonicon", pkgdir(Ion, "test", "repos"))
Ion.Internal.clone("Example", pkgdir(Ion, "test", "repos"))
Ion.Internal.clone("PlutoSliderServer", pkgdir(Ion, "test", "repos"))

@testset "ProjectToRelease" begin
    project = Internal.ProjectToRelease("patch", pkgdir(Ion, "test", "repos", "Comonicon"))
    @test project.name == "Comonicon"
    @test project.path == pkgdir(Ion, "test", "repos", "Comonicon")
    @test project.subdir == "."
    @test project.version.patch+1 == project.release_version.patch
    @test project.version_spec == "patch"
    @test project.branch == "master"

    print(stdin.buffer, 'y')
    project = Internal.ProjectToRelease("patch", pkgdir(Ion, "test", "repos", "Example"))
    readavailable(stdin.buffer)
    @test project.release_version == project.version
    @test Internal.should_update_version(project) == false

    project = Internal.ProjectToRelease("patch", pkgdir(Ion, "test", "repos", "Comonicon", "lib", "ComoniconTestUtils"); note="test note")
    @test project.subdir == "lib/ComoniconTestUtils"
    @test Internal.should_update_version(project) == true

    @test Internal.julia_registrator_comment(project) == """
    Released via [Ion](https://github.com/Roger-luo/Ion.jl)

    @JuliaRegistrator register branch=master subdir=lib/ComoniconTestUtils

    Patch notes:
    test note
    """

    project = Internal.ProjectToRelease("minor", pkgdir(Ion, "test", "repos", "Comonicon", "lib", "ComoniconTestUtils"); note="test note")
    @test Internal.julia_registrator_comment(project) == """
    Released via [Ion](https://github.com/Roger-luo/Ion.jl)

    @JuliaRegistrator register branch=master subdir=lib/ComoniconTestUtils

    Release notes:
    test note
    """

    project = Internal.ProjectToRelease("patch", pkgdir(Ion, "test", "repos", "PlutoSliderServer"))
    @test project.branch == "main"
    @test Internal.should_update_version(project) == true
end

Ion.release("patch", "Comonicon")
Ion.release("current", "../Comonicon")
project = Ion.Internal.ProjectToRelease("current", "../Comonicon")


text:
include("eqOdds.jl")
include("linprog.jl")
include("CalEqOdds.jl")


text:
"""
$(DocStringExtensions.TYPEDSIGNATURES)

A short description of the function

```jldoctest
julia> 1 + 1
2
```
"""
function initialize_transition_probabilities(kmer_graph)

    total_kmers = Graphs.nv(kmer_graph)
    transition_likelihoods = Dict(
        true => SparseArrays.spzeros(total_kmers, total_kmers),
        false => SparseArrays.spzeros(total_kmers, total_kmers)
    )

    for edge in collect(Graphs.edges(kmer_graph))
#         weight = length(kmer_graph.eprops[edge][:evidence])
        weight = kmer_graph.eprops[edge][:weight]
        for o in kmer_graph.eprops[edge][:orientations]
            transition_likelihoods[o.source_orientation][edge.src, edge.dst] = weight
        end
    end

    for source_orientation in (true, false)
        for src in 1:total_kmers
            transition_weights = transition_likelihoods[source_orientation][src, :]
            total_weight = sum(transition_weights)
            dsts, vals = SparseArrays.findnz(transition_weights)
            for (dst, val) in zip(dsts, vals)
                transition_likelihoods[source_orientation][src, dst] = val / total_weight
            end
            normalized_probability = sum(transition_likelihoods[source_orientation][src, :])
            @assert isapprox(normalized_probability, 0) || isapprox(normalized_probability, 1)
        end
    end
    return transition_likelihoods
end

"""
$(DocStringExtensions.TYPEDSIGNATURES)

A short description of the function

```jldoctest
julia> 1 + 1
2
```
"""
function set_initial_state_likelihoods!(
        kmer_graph,
        initial_state,
        kmer_likelihoods,
        error_rate,
        state_likelihoods,
        arrival_paths
    )
    for vertex in collect(Graphs.vertices(kmer_graph))
        hidden_kmer = kmer_graph.vprops[vertex][:kmer]

        fw_alignment =
            BioAlignments.pairalign(
                BioAlignments.LevenshteinDistance(),
                initial_state.fw,
                hidden_kmer)

        fw_probability = kmer_likelihoods[vertex]

        for match in 1:BioAlignments.count_matches(BioAlignments.alignment(fw_alignment))
            fw_probability *= 1 - error_rate
        end

        for edit in 1:fw_alignment.value
            fw_probability *= error_rate
        end

        bw_alignment =
            BioAlignments.pairalign(
                BioAlignments.LevenshteinDistance(),
                initial_state.bw,
                hidden_kmer)

        bw_probability = kmer_likelihoods[vertex]

        for match in 1:BioAlignments.count_matches(BioAlignments.alignment(bw_alignment))
            bw_probability *= 1 - error_rate
        end

        for edit in 1:bw_alignment.value
            bw_probability *= error_rate
        end

        if fw_probability > bw_probability
            state_probability = fw_probability
            state_orientation = true
        elseif fw_probability < bw_probability
            state_probability = bw_probability
            state_orientation = false
        else fw_probability == bw_probability
            state_probability = fw_probability
            state_orientation = missing
        end
        state_likelihoods[vertex, 1] = state_probability
        arrival_paths[vertex, 1] = [vertex => state_orientation]
    end
end

"""
$(DocStringExtensions.TYPEDSIGNATURES)

A short description of the function

```jldoctest
julia> 1 + 1
2
```
"""
function run_viterbi!(
        current_state,
        prior_state,
        observed_nucleotide,
        observed_quality_score,
        observed_error_rate,
        current_vertex,
        prior_vertex,
        state_likelihoods,
        transition_likelihoods,
        shortest_paths,
        arrival_paths,
        kmer_graph,
        kmer_likelihoods
        )
    # if probability of prior state is lower than current probability, skip

#     @show current_state
#     @show prior_state
#     @show current_vertex
#     @show prior_vertex

    current_state_likelihood = state_likelihoods[current_vertex, current_state]
    prior_state_likelihood = state_likelihoods[prior_vertex, prior_state]

    # if we already have a better possible path, skip calculating anything
    if prior_state_likelihood < current_state_likelihood
#         @show prior_state_likelihood < current_state_likelihood
        return
    end

    # take shortest path and assume it's the maximum likelihood path
    # this assumption seems fair because in an ideal situation
    # we're just moving to an adjacent kmer
    # and the shortest path and most likely path should be the same
    shortest_path = shortest_paths[prior_vertex][current_vertex]

#     no path & not considering insertion
    if isempty(shortest_path) && (prior_vertex != current_vertex)
#         @show "no path, skipping"
        return
    end

    # if shortest path isn't viable, exit
    if !isempty(shortest_path)
#         @show "checking if path is viable"

        terminal_orientation_prior_state = last(last(arrival_paths[prior_vertex, prior_state]))
#         @show arrival_paths[prior_vertex, prior_state]
#         @show "we were at vertex $(prior_vertex) in orientation $(terminal_orientation_prior_state)"
        candidate_edge = Graphs.Edge(shortest_path[1], shortest_path[2])

        if !ismissing(terminal_orientation_prior_state) &&
            !any(o -> o.source_orientation == terminal_orientation_prior_state, kmer_graph.eprops[candidate_edge][:orientations])

#             @show "no viable orientation matching edges detected between $(candidate_edge)"
#             @show "full candidate path was $(shortest_path)"
#             @show "orientation options were:"
#             @show kmer_graph.eprops[candidate_edge][:orientations]
            return
        end
    end

    # zero step path - insertion in observed sequence relative to kmer graph
    is_same_vertex = (current_vertex == prior_vertex)
    has_edge = Graphs.has_edge(kmer_graph, Graphs.Edge(prior_vertex, current_vertex))
    if is_same_vertex && has_edge
        shortest_path = [prior_vertex, current_vertex]
    end

    if is_same_vertex
#         @show "same vertex, considering insertion potential"
        emission_likelihood = observed_error_rate
        transition_likelihood = observed_error_rate
        state_likelihood = kmer_likelihoods[current_vertex]
        path_likelihood = prior_state_likelihood * emission_likelihood * transition_likelihood * state_likelihood
        path = [last(arrival_paths[prior_vertex, prior_state])]

        if current_state_likelihood > state_likelihoods[current_vertex, current_state]
#             @show "selecting path"
#             @show path
#             @show path_likelihood
            state_likelihoods[current_vertex, current_state] = path_likelihood
            arrival_paths[current_vertex, current_state] = path
        end
    # one or more step path - match, mismatch, or deletion in observed sequence relative to kmer graph
    elseif !isempty(shortest_path)
#         @show "path is viable!"
#         @show "considering shortest path: $(shortest_path)"

        initial_path_state = last(arrival_paths[prior_vertex, prior_state])

        path = Vector{typeof(initial_path_state)}(undef, length(shortest_path))
        path[1] = initial_path_state

        path_likelihood::Float64 = state_likelihoods[prior_vertex, prior_state]

        for i in 2:length(shortest_path)

            this_vertex = shortest_path[i]
            prior_vertex, prior_orientation = path[i-1]
            edge = Graphs.Edge(prior_vertex, this_vertex)

            possible_edge_orientations::Set{NamedTuple{(:source_orientation, :destination_orientation), Tuple{Bool, Bool}}} = kmer_graph.eprops[edge][:orientations]

#             @show possible_edge_orientations

            if !ismissing(prior_orientation)
                possible_edge_orientations = filter(o -> o.source_orientation == prior_orientation, possible_edge_orientations)
            end

#             @show possible_edge_orientations

            if isempty(possible_edge_orientations)
                path_likelihood *= 0.0
                path = Vector{eltype(path)}()
#                 @show "no possible orientations, bailing early"
                break
            end

#             @show prior_orientation
            if ismissing(prior_orientation)
                if transition_likelihoods[true][prior_vertex, this_vertex] > transition_likelihoods[false][prior_vertex, this_vertex]
                    prior_orientation = true
                    transition_likelihood = transition_likelihoods[true][prior_vertex, this_vertex]::Float64
                elseif transition_likelihoods[true][prior_vertex, this_vertex] < transition_likelihoods[false][prior_vertex, this_vertex]
                    prior_orientation = false
                    transition_likelihood = transition_likelihoods[false][prior_vertex, this_vertex]::Float64
                else transition_likelihoods[true][prior_vertex, this_vertex] == transition_likelihoods[false][prior_vertex, this_vertex]
                    prior_orientation = missing
                    transition_likelihood = transition_likelihoods[true][prior_vertex, this_vertex]::Float64
                end
            else
                transition_likelihood = transition_likelihoods[prior_orientation][prior_vertex, this_vertex]::Float64
            end
            state_likelihood::Float64 = kmer_likelihoods[this_vertex]
            path_likelihood *= transition_likelihood * state_likelihood

            if length(possible_edge_orientations) == 1
                orientation = first(possible_edge_orientations).destination_orientation
                path[i] = this_vertex => orientation
            else
                path[i] = this_vertex => missing
            end
        end

        # see if new nucleotide is a match or mismatch to terminal kmer in path
        if !isempty(path) && path_likelihood > 0
            terminal_kmer_index, terminal_kmer_orientation = last(path)
            terminal_kmer = BioSequences.LongDNASeq(kmer_graph.vprops[terminal_kmer_index][:kmer])::BioSequences.LongDNASeq
            if ismissing(terminal_kmer_orientation)
                fw_is_match = observed_nucleotide == last(terminal_kmer)
                bw_is_match = observed_nucleotide == last(BioSequences.reverse_complement!(terminal_kmer))
                if fw_ismatch && !bw_is_match
                    path[end] = terminal_kmer_index => true
                    path_likelihood *= 1 - observed_error_rate
                elseif !fw_ismatch && bw_is_match
                    path[end] = terminal_kmer_index => false
                    path_likelihood *= 1 - observed_error_rate
                elseif fw_ismatch && bw_is_match
                    path_likelihood *= 1 - observed_error_rate
                elseif !fw_ismatch && !bw_is_match
                    path_likelihood *= observed_error_rate
                end
            elseif terminal_kmer_orientation
                is_match = observed_nucleotide == last(terminal_kmer)
                if is_match
                    path_likelihood *= 1 - observed_error_rate
                else
                    path_likelihood *= observed_error_rate
                end
            else
                terminal_kmer = BioSequences.reverse_complement!(terminal_kmer)
                is_match = observed_nucleotide == last(terminal_kmer)
                if is_match
                    path_likelihood *= 1 - observed_error_rate
                else
                    path_likelihood *= observed_error_rate
                end
            end
        end

        if path_likelihood > state_likelihoods[current_vertex, current_state]
#             @show "selecting path"
#             @show path
#             @show path_likelihood
            state_likelihoods[current_vertex, current_state] = path_likelihood
            arrival_paths[current_vertex, current_state] = path
        end
    end
    return
end

"""
$(DocStringExtensions.TYPEDSIGNATURES)

A short description of the function

```jldoctest
julia> 1 + 1
2
```
"""
function polish_fastq(kmer_graph, fastq_file)

#     @info "Assessing kmer likelihoods"
    kmers = [kmer_graph.vprops[v][:kmer] for v in Graphs.vertices(kmer_graph)]
#     kmer_counts = [length(kmer_graph.vprops[v][:evidence]) for v in Graphs.vertices(kmer_graph)]
    kmer_counts = [kmer_graph.vprops[v][:weight] for v in Graphs.vertices(kmer_graph)]
    kmer_likelihoods = kmer_counts ./ sum(kmer_counts)
    k = kmer_graph.gprops[:k]
    kmer_type = BioSequences.BigDNAMer{k}
    total_kmers = length(kmers)

#     @info "determining shortest paths between kmers"
    shortest_paths = Graphs.enumerate_paths(Graphs.floyd_warshall_shortest_paths(kmer_graph));

    @info "counting the number of records to establish runtime estimate"
    number_of_records = 0
    for fastq_record in FASTX.FASTQ.Reader(open(fastq_file))
        number_of_records += 1
    end
    progress_bar = ProgressMeter.Progress(number_of_records, 1)

    output_fastq_file = replace(fastq_file, ".fastq" => ".k$(kmer_graph.gprops[:k]).fastq")
    fastq_writer = FASTX.FASTQ.Writer(open(output_fastq_file, "w"))
    for fastq_record in FASTX.FASTQ.Reader(open(fastq_file))
        ProgressMeter.next!(progress_bar)

#         @info "Initializing matrices"
        total_states = length(FASTX.sequence(fastq_record))-k+1
        transition_likelihoods = initialize_transition_probabilities(kmer_graph)
        state_likelihoods = zeros(total_kmers, total_states)
        arrival_paths = fill(Pair{Int, Union{Bool, Missing}}[], total_kmers, total_states)

#         @info "Determining Likelihoods of initial states"
        initial_state = first(BioSequences.each(kmer_type, FASTX.sequence(fastq_record)))
        current_state = 1
        # note this is a place for potential improvement, use the q value at each base to guide probability rather than median
        median_q_value = Statistics.median(Int.(FASTX.quality(fastq_record)[1:k]))
        current_error_rate = q_value_to_error_rate(median_q_value)
        # canonical_kmer = BioSequences.canonical(initial_state.fw)
        set_initial_state_likelihoods!(
                kmer_graph,
                initial_state,
                kmer_likelihoods,
                current_error_rate,
                state_likelihoods,
                arrival_paths
            )

#         @info "Determining likelihood of downstream states"

#         non_singleton_states = findall(kmer_counts .> 1)

        for current_state in 2:total_states
            prior_state = current_state - 1

        #     observed_kmer = BioSequences.BigDNAMer{k}(FASTX.sequence(fastq_record)[current_state:current_state+k-1])

        #     @assert observed_kmer == collect(BioSequences.each(kmer_type, FASTX.sequence(fastq_record)))[current_state].fw

        #     canonical_kmer = BioSequences.canonical(observed_kmer)

            observed_nucleotide = FASTX.sequence(fastq_record)[k-1+current_state]
        #     observed_nucleotide = last(observed_kmer)
            observed_quality_score = FASTX.quality(fastq_record)[k-1+current_state]
            observed_error_rate = q_value_to_error_rate(observed_quality_score)

            # we'll assess prior states in order of decreasing likelihood
            # such that we maximize how frequently we are able to utilize the
            # current_state_likelihood > candidate prior state
            # break that won't bother evaluating lower likelihood possibilities
            prior_states_in_decreasing_likelihood = sortperm(state_likelihoods[:, prior_state], rev=true)

            # and skip all prior states with zero probability

            for current_vertex in total_states
                for prior_vertex in prior_states_in_decreasing_likelihood
                    if state_likelihoods[prior_vertex, prior_state] > 0
                        run_viterbi!(
                                current_state,
                                prior_state,
                                observed_nucleotide,
                                observed_quality_score,
                                observed_error_rate,
                                current_vertex,
                                prior_vertex,
                                state_likelihoods,
                                transition_likelihoods,
                                shortest_paths,
                                arrival_paths,
                                kmer_graph,
                                kmer_likelihoods
                                )
                    end
                end
            end
        end

#         try
        maximum_likelihood_path, maximum_likelihood_value =
            determine_maximum_likelihood_path(
                state_likelihoods,
                arrival_paths
                )
#         catch
#             return state_likelihoods, arrival_paths
#         end

        sequence = oriented_path_to_sequence(kmer_graph, maximum_likelihood_path)

#         @info "comparing to original path"
        original_sequence_likelihood = oriented_path_to_likelihood(kmer_graph, kmers, kmer_likelihoods, transition_likelihoods, fastq_record)
        relative_likelihood = maximum_likelihood_value / original_sequence_likelihood
#         relative_likelihood_formatted = NumericIO.formatted(relative_likelihood, ndigits=1, charset=:ASCII)
#         println("relative likelihood of new path to old path is $(relative_likelihood_formatted)")

#         @info "writing updated record"
        identifier = FASTX.identifier(fastq_record) * "_k$(k)"
        description = string(relative_likelihood)
        # because the sequences won't always be the same length, we take an ordered sampling with replacement
        # which introduces some random error but preserves overall patterns and areas of high/low accuracy
        quality_scores = StatsBase.sample(FASTX.quality(fastq_record), length(sequence), ordered=true)

        new_fastq_record = FASTX.FASTQ.Record(
            identifier,
            description,
            sequence,
            quality_scores
        )
        write(fastq_writer, new_fastq_record)
    end
    close(fastq_writer)
    return output_fastq_file
end

"""
$(DocStringExtensions.TYPEDSIGNATURES)

A short description of the function

```jldoctest
julia> 1 + 1
2
```
"""
function determine_maximum_likelihood_path(
    state_likelihoods,
    arrival_paths
    )
    maximum_likelihood_value = maximum(state_likelihoods[:, end])

    maximum_likelihood_path_indices = findall(state_likelihoods[:, end] .== maximum_likelihood_value)

    # if multiple paths are tied, randomly choose one
    maximum_likelihood_path_index = rand(maximum_likelihood_path_indices)

    maximum_likelihood_path = arrival_paths[maximum_likelihood_path_index, end]

    for state_index in size(state_likelihoods, 2)-1:-1:1
        next_kmer, next_orientation = first(maximum_likelihood_path)
        maximum_likelihood_arrival_path = arrival_paths[next_kmer, state_index]

        is_match = last(maximum_likelihood_arrival_path) == (next_kmer => next_orientation)
        if !ismissing(is_match) && !is_match
            error("breaking")
        end
        maximum_likelihood_path = vcat(maximum_likelihood_arrival_path[1:end-1], maximum_likelihood_path)
    end
    return maximum_likelihood_path, maximum_likelihood_value
end

max_stars_repo_name: snystrom/Mycelia

max_issues_repo_name: snystrom/Mycelia

max_forks_repo_name: cjprybol/Eisenia