text:
parr = [
    ["Q", "u[1]", 12, 12, 12, 12],
    ["Q", "u[2]", 12, 12, 12, 12],
    ["Q", :η, 12, 12, 12, 12],
    ["Q", :θ, 12, 12, 12, 0],
    ["s_aux", :y, 12, 12, 12, 12],
    ["s_aux", :w, 12, 12, 12, 12],
    ["s_aux", :pkin, 12, 12, 12, 12],
    ["s_aux", :wz0, 12, 12, 10, 12],
    ["s_aux", "uᵈ[1]", 12, 12, 12, 12],
    ["s_aux", "uᵈ[2]", 12, 12, 12, 12],
    ["s_aux", "ΔGᵘ[1]", 12, 12, 12, 12],
    ["s_aux", "ΔGᵘ[2]", 12, 12, 12, 12],
]

#! format: off
imex = [
 [     "Q",   "u[1]", -3.79026877309739870792e-02,  3.77520487392760603629e-02, -1.40391503423503398694e-06,  5.33572816364154094149e-03 ],
 [     "Q",   "u[2]", -7.13128520899887746037e-03,  6.54348134460300615012e-03, -5.44654209452686030146e-06,  9.83670484102603564017e-04 ],
 [     "Q",       :η, -5.75156897888061836160e-03,  5.06823909948793567548e-03, -1.61345029895201705567e-05,  1.58599716799624735307e-03 ],
 [     "Q",       :θ,  2.00000000000000000000e+01,  2.00000000000003126388e+01,  2.00000000000000532907e+01,  6.62241227636115094309e-14 ],
 [ "s_aux",       :y,  0.00000000000000000000e+00,  1.00000000000000011642e+06,  5.00000000000000000000e+05,  2.92779390974978974555e+05 ],
 [ "s_aux",       :w, -1.93125254703897220022e-05,  1.89413928193982494038e-05,  6.96287201343549477191e-08,  1.82962329428399789732e-06 ],
 [ "s_aux",    :pkin, -1.60000000000008335554e+00,  0.00000000000000000000e+00, -8.00000000000004152234e-01,  4.68447025559968477371e-01 ],
 [ "s_aux",     :wz0, -1.78460606500926580929e-06,  1.42096532424030984111e-06, -3.15575167805227919492e-09,  7.00511144160315581441e-07 ],
 [ "s_aux",  "uᵈ[1]",  0.00000000000000000000e+00,  0.00000000000000000000e+00,  0.00000000000000000000e+00,  0.00000000000000000000e+00 ],
 [ "s_aux",  "uᵈ[2]",  0.00000000000000000000e+00,  0.00000000000000000000e+00,  0.00000000000000000000e+00,  0.00000000000000000000e+00 ],
 [ "s_aux", "ΔGᵘ[1]",  0.00000000000000000000e+00,  0.00000000000000000000e+00,  0.00000000000000000000e+00,  0.00000000000000000000e+00 ],
 [ "s_aux", "ΔGᵘ[2]",  0.00000000000000000000e+00,  0.00000000000000000000e+00,  0.00000000000000000000e+00,  0.00000000000000000000e+00 ],
]

explicit = [
 [     "Q",   "u[1]", -3.74270752639266485184e-02,  3.72763215301358599718e-02, -1.40392694289570846012e-06,  5.29521849931501031872e-03 ],
 [     "Q",   "u[2]", -7.13776376792429369245e-03,  6.54949226335153813355e-03, -5.45004642316166964817e-06,  9.84283148803937564184e-04 ],
 [     "Q",       :η, -5.75146380759241834801e-03,  5.06819905743131991077e-03, -1.61399463682884620498e-05,  1.58594255118522063641e-03 ],
 [     "Q",       :θ,  1.99999999999998934186e+01,  2.00000000000014495072e+01,  2.00000000000005542233e+01,  2.51921479392847498778e-13 ],
 [ "s_aux",       :y,  0.00000000000000000000e+00,  1.00000000000000011642e+06,  5.00000000000000000000e+05,  2.92779390974978974555e+05 ],
 [ "s_aux",       :w, -1.91903846873016302693e-05,  1.88201368002799860018e-05,  6.88331165233151285677e-08,  1.81657738735039775970e-06 ],
 [ "s_aux",    :pkin, -1.60000000000012065904e+00,  0.00000000000000000000e+00, -8.00000000000023470115e-01,  4.68447025559979413067e-01 ],
 [ "s_aux",     :wz0, -1.79359066321377860068e-06,  1.42978124247214314059e-06, -3.01455602003261245701e-09,  6.97246485534473377005e-07 ],
 [ "s_aux",  "uᵈ[1]",  0.00000000000000000000e+00,  0.00000000000000000000e+00,  0.00000000000000000000e+00,  0.00000000000000000000e+00 ],
 [ "s_aux",  "uᵈ[2]",  0.00000000000000000000e+00,  0.00000000000000000000e+00,  0.00000000000000000000e+00,  0.00000000000000000000e+00 ],
 [ "s_aux", "ΔGᵘ[1]",  0.00000000000000000000e+00,  0.00000000000000000000e+00,  0.00000000000000000000e+00,  0.00000000000000000000e+00 ],
 [ "s_aux", "ΔGᵘ[2]",  0.00000000000000000000e+00,  0.00000000000000000000e+00,  0.00000000000000000000e+00,  0.00000000000000000000e+00 ],
]

long = [
 [     "Q",   "u[1]", -1.26047572708769956007e-01,  8.73949752192766937675e-02,  4.77293205332110539161e-05,  7.01620954203984511371e-03 ],
 [     "Q",   "u[2]", -7.74039642997905313893e-02,  1.17079049312642666236e-01, -1.81789027476037060851e-04,  1.05959299008931383257e-02 ],
 [     "Q",       :η, -7.88344390091725993885e-02,  3.98295576514583818395e-02, -9.41579559913946891004e-05,  2.80408708312551574471e-02 ],
 [     "Q",       :θ,  1.99999999999542765750e+01,  2.00000000000509103870e+01,  2.00000000000029913849e+01,  2.16653930768717257459e-12 ],
 [ "s_aux",       :y,  0.00000000000000000000e+00,  4.00000000000000046566e+06,  2.00000000000000000000e+06,  1.15573129229947482236e+06 ],
 [ "s_aux",       :w, -2.39010107629734968942e-05,  6.21800814931180984403e-05,  2.57411032687584406218e-07,  4.26276984125372136047e-06 ],
 [ "s_aux",    :pkin, -1.60000000000374065223e+00,  0.00000000000000000000e+00, -8.00000000000172128978e-01,  4.61946285916578103237e-01 ],
 [ "s_aux",     :wz0, -2.51713525691602801182e-06,  9.22709404870064716511e-07, -8.84750036597309359132e-10,  5.52781925121154085940e-07 ],
 [ "s_aux",  "uᵈ[1]",  0.00000000000000000000e+00,  0.00000000000000000000e+00,  0.00000000000000000000e+00,  0.00000000000000000000e+00 ],
 [ "s_aux",  "uᵈ[2]",  0.00000000000000000000e+00,  0.00000000000000000000e+00,  0.00000000000000000000e+00,  0.00000000000000000000e+00 ],
 [ "s_aux", "ΔGᵘ[1]",  0.00000000000000000000e+00,  0.00000000000000000000e+00,  0.00000000000000000000e+00,  0.00000000000000000000e+00 ],
 [ "s_aux", "ΔGᵘ[2]",  0.00000000000000000000e+00,  0.00000000000000000000e+00,  0.00000000000000000000e+00,  0.00000000000000000000e+00 ],
]

#! format: on

refVals =
    (explicit = (explicit, parr), imex = (imex, parr), long = (long, parr))

max_stars_repo_name: mwarusz/CLIMA

max_issues_repo_name: mwarusz/CLIMA

max_forks_repo_name: mwarusz/CLIMA


text:
t(letter) = Int(letter) - 96

function counting_sort_letters(A, position, max_length=30)
  B = Array{String}(undef,length(A))
  C = zeros(Int64,max_length)
  for letters in A
    letter_position = t(letters[position])
    C[letter_position]+=1
  end
  for i=2:length(C)
    C[i] += C[i - 1]
  end
  for j=1:length(A)
    i = length(A) - j + 1
    a_index = t(A[i][position])
    B[C[a_index]] = A[i]
    C[a_index]-=1
  end
 return B
end

accepted_response = ["ca", "cba", "ab", "abc", "ccc"]

# @time counting_sort_letters(["ccc", "cba", "ca", "ab", "abc"], 2)==accepted_response

function counting_sort_length(A)
  # Returnerer direkte om det bare er ett (eller ingen) ord
  if length(A) <= 1
      return A
  end

  # Finner det lengste ordet
  lw = 0
  for i in A
      if length(i) >= lw
          lw = length(i)
      end
  end

  a_lenght = length(A)
  B = Array{String}(undef, length(A))
  C = zeros(Int64, lw + 1)

  # Setter inn i liste C
  for i in A
      p = length(i) + 1
      C[p] += 1
  end

  # Klargjør liste C for innsettingen i liste B
  for i = 2:length(C)
      C[i] += C[i-1]
  end

  # Setter inn i liste B i sortert rekkefølge

  for i = length(A) : -1 : 1
      p = length(A[i])+1
      B[C[p]] = A[i]
      C[p] -= 1
  end
  return B
end

# @info "COUNTING SORT!" counting_sort_length(["bbbb","aaaaaa", "aa", "aaaa", "ccc", "aaaaa"]) == ["aa", "ccc", "bbbb", "aaaa", "aaaaa","aaaaaa"]
# @info "Counting sort" counting_sort_length(["bbbb", "aa", "aaaa", "ccc", "aaaaa"])
# @info "Counting 2" counting_sort_length(["a"])

te(letter) = Int(letter) - 96

function counting_sort_letters_two(A, position, max_length=3)
  B = Array{String}(undef,length(A))
  C = zeros(Int64,max_length)
  for letters in A
    letter_position = te(letters[position])
    C[letter_position]+=1
  end
  for i=2:length(C)
    C[i] += C[i - 1]
  end
  for j=1:length(A)
    i = length(A) - j + 1
    a_index = te(A[i][position])
    B[C[a_index]] = A[i]
    C[a_index]-=1
  end
 return B
end

function flexradix(A, max_length)
for i in 1:length(A)
    i_length = length(A[i])
    if i_length < max_length
        for j=1:(max_length - i_length)
            A[i]*= " "
        end
    end

end
    println(A)
  for j = 1:max_length
    i = max_length - j + 1
    A = counting_sort_letters(A,i)
  end
  for i in 1:length(A)
  A[i] = rstrip(A[i])
  end
  return A
end

#flexradix(["kobra", "aggie", "agg", "kort", "hyblen"],6)

ext: jl


text:
export kummer_extension, ray_class_field

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

function kummer_extension(n::Int, gen::Array{nf_elem, 1})
  g = [FacElem(x) for x=gen]
  return kummer_extension(n, g)
end

function _extend_auto(K::Hecke.NfRel{nf_elem}, h::Hecke.NfToNfMor)
  @assert iskummer_extension(K)
  k = base_ring(K)
  Zk = maximal_order(k)
  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 = zeta
  while im_zeta != z
    r += 1
    z *= zeta
  end

  a = -coeff(K.pol, 0)
  a = a^r//h(a) # this assumes K/k to be abelian
  #global last_rt = (a, degree(K))
  fl, b = hasroot(a, degree(K))
  @assert fl

  return NfRelToNfRelMor(K, K, h, 1//b*gen(K)^r)
end

###############################################################################
#
#  Computation of Frobenius automorphisms
#
###############################################################################

# the Frobenius at p in K:
#K is an extension of k, p a prime in k,
#returns a vector in (Z/nZ)^r representing the Frob
function can_frobenius(p::NfOrdIdl, K::KummerExt)
  if haskey(K.frob_cache, p)
    return K.frob_cache[p]
  end
  Zk = order(p)
  if index(Zk) % minimum(p) == 0
    #index divisors and residue class fields don't agree
    # ex: x^2-10, rcf of 29*Zk, 7. 239 is tricky...
    throw(BadPrime(p))
  end

  if nbits(minimum(p)) > 64
    error("Oops")
  end

  F, mF = ResidueFieldSmall(Zk, p)
  mF = extend_easy(mF, number_field(Zk))

  # K = sqrt[n](gen), an automorphism will be
  # K[i] -> zeta^? K[i]
  # Frob(sqrt[n](a), p) = sqrt[n](a)^N(p) (mod p) = zeta^r sqrt[n](a)
  # sqrt[n](a)^N(p) = a^(N(p)-1 / n) = zeta^r mod p

  z_p = inv(mF(Zk(K.zeta)))
  @assert norm(p) % K.n == 1
  ex = div(norm(p)-1, K.n)
  aut = Array{fmpz,1}(length(K.gen))
  for j=1:length(K.gen)
    mu = mF(K.gen[j])^ex  # can throw bad prime!
    i = 0
    while !isone(mu)
      i += 1
      @assert i <= K.n
      mul!(mu, mu, z_p)
    end
    aut[j]= fmpz(i)
  end
  z = K.AutG(aut)
  K.frob_cache[p] = z
  return z
end

function can_frobenius(p::NfOrdIdl, K::KummerExt, g::FacElem{nf_elem})
  Zk = order(p)
  if index(Zk) % minimum(p) == 0
    #index divisors and residue class fields don't agree
    # ex: x^2-10, rcf of 29*Zk, 7. 239 is tricky...
    throw(BadPrime(p))
  end

  if nbits(minimum(p)) > 64
    error("Oops")
  end

  F, mF = ResidueFieldSmall(Zk, p)
  mF = extend_easy(mF, number_field(Zk))

  #K = sqrt[n](gen), an automorphism will be
  # K[i] -> zeta^? K[i]
  # Frob(sqrt[n](a), p) = sqrt[n](a)^N(p) (mod p) = zeta^r sqrt[n](a)
  # sqrt[n](a)^N(p) = a^(N(p)-1 / n) = zeta^r mod p

  z_p = inv(mF(Zk(K.zeta)))
  @assert norm(p) % K.n == 1
  ex = div(norm(p)-1, K.n)
  aut = fmpz[]

  mu = mF(g)^ex  # can throw bad prime!
  i = 0
  while !isone(mu)
    i += 1
    @assert i <= K.n
    mul!(mu, mu, z_p)
  end
  return i
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))[]

  st = start(S)
  np = 0
  extra = 1

  q, mq = quo(R, sR, false)
  while true
    p, st = next(S, st)
    if cp % p == 0
      continue
    end
    @vprint :ClassField 2 "doin` $p\n"
    np += 1

    if np % ngens(R) == 0
      @vprint :ClassField 3 "need many generators... $np for group with $(ngens(R))"
    end
    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
        continue
      end
      if iszero(mq(f))
        continue
      end
      push!(sR, f)
      push!(lp, P)
      q, mq = quo(R, sR, false)
    end
    if order(q) == 1
      break
    end
  end
  return lp, sR

end

###############################################################################
#
#  Ray Class Field, number_field interface and reduction to prime power case
#
###############################################################################

#doc"""
#    ray_class_field(m::Map, p::Int=0) -> ClassField
#> Creates the (formal) abelian extension defined by the map $m: A \to I$
#> where $I$ is the set of ideals coprime to the modulus defining $m$ and $A$
#> is a quotient of the ray class group (or class group).
#> If $p$ is given and non-zero, only the quotient modulo $p$-th powers is
#> created.
#"""
function ray_class_field(m::Union{MapClassGrp, MapRayClassGrp})
  return ray_class_field(m, GrpAbFinGenMap(domain(m)))
end

function ray_class_field(m::Union{MapClassGrp, MapRayClassGrp}, quomap::GrpAbFinGenMap)
  CF = ClassField()
  CF.rayclassgroupmap = m
  S, mS = snf(codomain(quomap))
  CF.quotientmap = Hecke._compose(inv(mS), quomap)
  #CF.mq = Hecke.make_snf(Hecke._compose(m, inv(quomap)))
  return CF
end

doc"""
    number_field(CF::ClassField) -> Array{GenPol{nf_elem}, 1}
> Given a (formal) abelian extension, compute defining polynomials
> for all prime power cyclic subfields.
"""
function number_field(CF::ClassField)
  if isdefined(CF, :A)
    return CF.A
  end

  res = []
  G = codomain(CF.quotientmap)
  @assert issnf(G)
  q = [G[i] for i=1:ngens(G)]
  for i=1:ngens(G)
    o = order(G[i])
    lo = factor(o)
    for (p, e) = lo.fac
      q[i] = p^e*G[i]
      S, mQ = quo(G, q, false)
      push!(res, ray_class_field_cyclic_pp(CF, mQ))
    end
    q[i] = G[i]
  end
  CF.cyc = res
  CF.A = number_field([x.A.pol for x = CF.cyc])[1]
  return CF.A
end

function ray_class_field_cyclic_pp(CF::ClassField, mQ::GrpAbFinGenMap)
  @vprint :ClassField 1 "cyclic prime power class field of degree $(degree(CF))\n"
  CFpp = ClassField_pp()
  #CFpp.mq = _compose(CF.mq, inv(mQ))
  CFpp.quotientmap = _compose(mQ, CF.quotientmap)
  CFpp.rayclassgroupmap = CF.rayclassgroupmap
  @assert domain(CFpp.rayclassgroupmap) == domain(CFpp.quotientmap)
  @vprint :ClassField 1 "finding the Kummer extension...\n"
  _rcf_find_kummer(CFpp)
  @vprint :ClassField 1 "reducing the generator...\n"
  _rcf_reduce(CFpp)
  @vprint :ClassField 1 "descending ...\n"
  _rcf_descent(CFpp)
  return CFpp
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(mR::Map, K::KummerExt, c::CyclotomicExt)
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)
  ZK = maximal_order(base_ring(K.gen[1]))
  cp = lcm(minimum(m), discriminant(ZK))
  cf = Hecke.MapFromFunc(x->can_frobenius(x, K), IdealSet(ZK), K.AutG)

  mp = c.mp[2]
  ZK = maximal_order(c.Ka)
  @hassert :ClassField 1 order(domain(cf)) == ZK
  Zk = order(m)
  Id_Zk = Hecke.NfOrdIdlSet(Zk)
  k = nf(Zk)
  @hassert :ClassField 1 k == domain(mp)

  Sp = Hecke.PrimesSet(200, -1, c.n, 1) #primes = 1 mod n, so totally split in cyclo

  lp, sG = find_gens(cf, Sp, cp)

  G = codomain(cf)
  sR = Array{GrpAbFinGenElem, 1}(length(lp))

  for i=1:length(lp)
    p = Id_Zk(intersect_nonindex(mp, lp[i]))
    sR[i]= valuation(norm(lp[i]), norm(p))*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)
  return h
end

function _rcf_find_kummer(CF::ClassField_pp)
  #mq = CF.mq
  if isdefined(CF, :K)
    return CF.K
  end
  #f = _modulus(mq)
  f = defining_modulus(CF)
  @vprint :ClassField 2 "Kummer extension ... with conductor(?) $f\n"
  k1 = nf(order(f))
  e = degree(CF)
  @assert Hecke.isprime_power(e)

  @vprint :ClassField 2 "Adjoining the root of unity\n"
  C = cyclotomic_extension(k1, e)
  K = C.Ka

  @vprint :ClassField 2 "Maximal order of cyclotomic extension\n"
  ZK = maximal_order(K)

  @vprint :ClassField 2 "Class group of cyclotomic extension\n"
  c, mc = class_group(ZK)
  allow_cache!(mc)
  @vprint :ClassField 2 "... $c\n"
  c, mq = quo(c, e, false)
  mc = _compose(mc, inv(mq))

  lf = factor(minimum(f)*e)
  lP = Hecke.NfOrdIdl[]
  #Why am I considering e? It does not appear in the theory.
  #Furthermore, some of the factors of the modulus can be ignored, since
  #I am only considering the prime power part.
  for p = keys(lf.fac)
    lp = prime_decomposition(ZK, p)  #TODO: make it work for fmpz
    lP = vcat(lP, [x[1] for x = lp])
  end
  g = GrpAbFinGenElem[preimage(mc, x) for x = lP]

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

  lP = vcat(lP, Hecke.find_gens(inv(mc), PrimesSet(100, -1))[1])
  @vprint :ClassField 2 "using $lP of length $(length(lP)) for S-units\n"
#if false
#  println("enlarging to be Galois closed - just in case...", length(lP))
#  lP = Set(lP)
#  _lp = Set(minimum(x) for x = lP)
#  for p = _lp
#    fp = prime_decomposition(ZK, Int(p))
#    for I = fp
#      push!(lP, I[1])
#    end
#  end
#  lP = collect(lP)
#  println("finally:", length(lP))
#end
  @vprint :ClassField 2 "using $lP of length $(length(lP)) for s-units\n"
  @vtime :ClassField 2 S, mS = Hecke.sunit_group_fac_elem(lP)
  @vprint :ClassField 2 "... done\n"
  @vtime :ClassField 2 KK = kummer_extension(e, [mS(S[i]) for i=1:ngens(S)])
  CF.bigK = KK

  @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"

  k, mk = kernel(h)
  G = domain(h)

  # 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
  # Carlo: The ray class group should already have the rigth exponent.
  #        Why do I need to repeat the kernel mod n?
  R = ResidueRing(FlintZZ, C.n, cached=false)
  M = MatrixSpace(R, ngens(k), ngens(G), false)(mk.map)
  #=
  M = zero_matrix(R, ngens(k), ngens(G))
  for i=1:ngens(k)
    ki = mk(k[i])
    for j=1:ngens(G)
      M[i, j] = ki[j]
    end
  end
  =#
  l, i = nullspace(M)
  @assert i>0
  n = lift(l)
  N = GrpAbFinGen([e for j=1:rows(n)])
  s, ms = sub(N, GrpAbFinGenElem[sum([n[j, k]*N[j] for j=1:rows(n)]) for k=1:i], false)
  ms = Hecke.make_snf(ms)
  @hassert :ClassField 1 iscyclic(domain(ms))
  o = Int(order(domain(ms)[1]))
  c = 1
  if o < e
    c = div(e, o)
  end
  g = ms(domain(ms)[1])
  @vprint :ClassField 2 "final $n of order $o and e=$e\n"
  a = prod([KK.gen[i]^div(mod(n[i], e), c) for i=1:ngens(N)])
  @vprint :ClassField 2 "generator $a\n"
  CF.a = a
  CF.sup = lP
  CF.sup_known = true
  CF.o = o
#  CF.K = pure_extension(Int(o), a)[1] #needs to evaluate a - too expensive!
  return nothing
end

function _rcf_reduce(CF::ClassField_pp)
  #e = order(codomain(CF.quotientmap))
  e = degree(CF)
  if CF.sup_known
    CF.a = reduce_mod_powers(CF.a, CF.o, CF.sup)
    CF.sup_known = false
  else
    CF.a = reduce_mod_powers(CF.a, CF.o)
  end
  CF.K = pure_extension(CF.o, CF.a)[1]
  return nothing
end

function _find_prim_elem(A::NfRel, AutA_gen::Array{NfRelToNfRelMor{nf_elem,  nf_elem},1}, AutA::GrpAbFinGen, oA::fmpz, C::CyclotomicExt)
  pe = gen(A)# + 0*gen(C.Ka)
  Auto=Dict{Hecke.GrpAbFinGenElem, Any}()
  cnt = 0
  while true
    @vprint :ClassField 3 "candidate: $pe\n"
    Im = Set{Hecke.NfRelElem{nf_elem}}()
    for j = AutA
      im = grp_elem_to_map(AutA_gen, j, pe)
      if im in Im
        pe += gen(C.Ka)
        cnt += 1
        if cnt > 100
          error("Too many attempts to find primitive elements")
        end
        break
      else
        push!(Im, im)
      end
      Auto[j]=im
    end
    if length(Im) == oA
      break
    end
  end
  @vprint :ClassField 2 "have primitive element!!!  "
  @vprint :ClassField 3 " $pe"

  return pe, Auto
end

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

#=
    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 autpmorphisms 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)
  #mg = Hecke.make_snf(mg)
  @assert (e%8 == 0 && ngens(domain(mg))==2) ||
           ngens(domain(mg)) <= 1

  if degree(C.Kr) < order(g)  # there was a common subfield, we
                              # have to pass to a subgroup
    @assert order(g) % degree(C.Kr) == 0
    f = C.Kr.pol
    s, ms = sub(g, [x for x = g if iszero(f(gen(C.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"

  AutA_gen = Hecke.NfRelToNfRelMor{nf_elem,  nf_elem}[]
  AutA_rel = zero_matrix(FlintZZ, ngens(g)+1, ngens(g)+1)
  zeta = C.mp[1](gen(C.Kr))
  n = degree(A)
  @assert e % n == 0

  @vprint :ClassField 2 "... the trivial one (Kummer)\n"
  tau = Hecke.NfRelToNfRelMor{nf_elem,  nf_elem}(A, A, zeta^div(e, n)*gen(A))
  push!(AutA_gen, tau)

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

  K = C.Ka
  @vprint :ClassField 2 "... need to extend $(ngens(g)) from the cyclo ext\n"
  for i=1:ngens(g)
    si = Hecke.NfRelToNfRelMor{nf_elem, nf_elem}(C.Kr, C.Kr, gen(C.Kr)^Int(lift(mg(g[i]))))
    @vprint :ClassField 2 "... extending zeta -> zeta^$(mg(g[i]))\n"
    sigma = _extend_auto(A, Hecke.NfToNfMor(K, K, C.mp[1](si(preimage(C.mp[1], gen(K))))))
    push!(AutA_gen, sigma)

#    pe = 17*gen(K) + gen(A)
#    @assert sigma(tau(pe)) - tau(sigma(pe)) == 0
#    @assert sigma(tau(pe)) == tau(sigma(pe))

    @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
    local zeta_i::nf_elem = (inv(zeta)^div(e, n))::nf_elem
    local mi::nf_elem = coeff(m, 1)
    @hassert :ClassField 1 m == mi*gen(A)

    while mi != 1
      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
  return nothing

end

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

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

  e = degree(CF)
  k = nf(order(codomain(CF.rayclassgroupmap)))
  C = cyclotomic_extension(k, e)
  A = CF.K
  n = degree(A)
  @vprint :ClassField 2 "Automorphism group (over ground field) $AutA\n"
  _aut_A_over_k(A, C, CF)

  AutA_gen = CF.AutG
  AutA = GrpAbFinGen(CF.AutR)
  AutA_snf, mp = snf(AutA)
  # now we need a primitive element for A/k
  # mostly, gen(A) will do
  @vprint :ClassField 2 "Searching for primitive element...\n"
  pe, Auto = _find_prim_elem(A::NfRel, AutA_gen, AutA, order(AutA_snf), C)
  @vprint :ClassField 2 "\nnow the fix group...\n"

  if iscyclic(AutA_snf)  # the subgroup is trivial to find!
    @vprint :ClassField 2 ".. trivial as automorphism group is cyclic\n"
    s, ms = sub(AutA_snf, [e*AutA_snf[1]], false)
  else
    @vprint :ClassField 2 ".. interesting...\n"
    #want: hom: AutA = Gal(A/k) -> Gal(K/k) = domain(mq)
    # idea: take primes p in k and compare
    #  Frob(p, A/k) and preimage(mq, p)
    @assert n == degree(CF.K)
    Zk = order(defining_modulus(CF))
    function canFrob(p::NfOrdIdl)
      lP = Hecke.prime_decomposition_nonindex(C.mp[2], p)
      P = lP[1][1]
      Q = lP[end][1]
      F, mF = ResidueFieldSmall(order(P), P)
      Ft, t = PolynomialRing(F, cached=false)
      mFp = extend_easy(mF, C.Ka)
      ap = mFp(CF.a)
      Ap = ResidueRing(Ft, t^n-ap)

      function toAp(x) # x in CF.K
        xp = zero(Ft)
        @assert coeff(x, n) == 0
        for i=0:n-1
          setcoeff!(xp, i, mFp(coeff(x, i)))
        end
        return Ap(xp)
      end

#      @assert length(Set([toAp(v) for (k,v) = Auto])) == length(Auto)
#TODO: bad prime possible if set is too small (pe might interfere with p)
      imF = toAp(pe)^norm(p)
      for (ky,v) = Auto
        kp = toAp(v)
        if kp == imF
          return ky
        end
      end
      error("Frob not found")
    end
    @vprint :ClassField 2 "finding Artin map...\n"
#TODO can only use non-indx primes, easy primes...
    @vtime :ClassField 2 lp, f = find_gens(MapFromFunc(canFrob, IdealSet(Zk), AutA),
                      PrimesSet(200, -1), minimum(defining_modulus(CF)))
    h = hom(f, [image(CF.quotientmap, preimage(CF.rayclassgroupmap, p)) for p = lp])
    @hassert :ClassField 1 issurjective(h)
    h = _compose(mp, h)
    h = hom(AutA_snf, [h(AutA_snf[i]) for i=1:ngens(AutA_snf)])
    s, ms = kernel(h)
    @vprint :ClassField 2 "... done, have subgroup!\n"
  end

  #now, hopefully either norm or trace will be primitive for the target
  #norm, trace relative to s, the subgroup

  @vprint :ClassField 2 "computing orbit of primitive element\n"
  os = [Auto[mp(ms(j))] for j=s]

  function coerce_down(a)
    @assert a.data.length <= 1
    b = coeff(a, 0)
    c = preimage(C.mp[1], b)
    @assert c.data.length <= 1
    return coeff(c, 0)
  end

  function minpoly(a)
    @vtime :ClassField 2 o = [grp_elem_to_map(AutA_gen, mp(mq(j)), t) for j = q]
    @vtime :ClassField 2 f = prod(T-x for x=o)
    @assert degree(f) == length(o)
    @assert length(o) == e
    @vtime :ClassField 2 g = [coerce_down(coeff(f, i)) for i=0:Int(e)]
    return PolynomialRing(parent(g[1]), cached = false)[1](g)
  end

#  n = prod(os) # maybe primitive??
  @vprint :ClassField 2 "trying relative trace\n"
  t = sum(os)
  CF.pe = t
  #now the minpoly of t - via Galois as this is easiest to implement...
  q, mq = quo(AutA_snf, [ms(s[i]) for i=1:ngens(s)], false)
  @assert order(q) == degree(CF)
  AT, T = PolynomialRing(A, "T", cached = false)
  @vprint :ClassField 2 "char poly...\n"
  f = minpoly(t)
  @vprint :ClassField 2 "... done\n"
  if !issquarefree(f)
    @vprint :ClassField 2 "trying relative trace of squares\n"
    for i=1:length(os)
      os[i] *= os[i]
    end
    t = sum(os)
    CF.pe = t
    #now the minpoly of t - via Galois as this is easiest to implement...
    q, mq = quo(AutA_snf, [ms(s[i]) for i=1:ngens(s)], false)
    @assert order(q) == degree(CF)
    AT, T = PolynomialRing(A, "T", cached = false)
    @vprint :ClassField 2 "char poly...\n"
    f = minpoly(t)
    @vprint :ClassField 2 "... done\n"
    @hassert :ClassField 1 issquarefree(f)
  end

  CF.A = number_field(f)[1]
  return nothing
end

function grp_elem_to_map(A::Array, b::Hecke.GrpAbFinGenElem, pe)
  res = pe
  for i=1:length(A)
    if b[i] == 0
      continue
    end
    for j=1:b[i]
      res = A[i](res)
    end
  end
  return res
end

function _reduce(a::fq_nmod)
  A = parent(a)
  if a.length < 2*degree(A)
    ccall((:fq_nmod_reduce, :libflint), Void, (Ptr{fq_nmod}, Ptr{FqNmodFiniteField}), &a, &A)
  else
    ccall((:nmod_poly_rem, :libflint), Void, (Ptr{fq_nmod}, Ptr{fq_nmod}, Ptr{Void}, Ptr{Void}), &a, &a, pointer_from_objref(A)+6*sizeof(Int) + 2*sizeof(Ptr{Void}), pointer_from_objref(A)+sizeof(fmpz))
  end
end

#TODO: move elsewhere - and use. There are more calls to nmod_set/reduce
function (A::FqNmodFiniteField)(x::nmod_poly)
  u = A()
  ccall((:fq_nmod_set, :libflint), Void,
                     (Ptr{fq_nmod}, Ptr{nmod_poly}, Ptr{FqNmodFiniteField}),
                                     &u, &x, &A)
  _reduce(u)
#  ccall((:fq_nmod_reduce, :libflint), Void, (Ptr{fq_nmod}, Ptr{FqNmodFiniteField}), &u, &A)
  return u
end

function _nf_to_fq!(a::fq_nmod, b::nf_elem, K::FqNmodFiniteField, a_tmp::nmod_poly)
  nf_elem_to_nmod_poly!(a_tmp, b)
  ccall((:fq_nmod_set, :libflint), Void,
                     (Ptr{fq_nmod}, Ptr{nmod_poly}, Ptr{FqNmodFiniteField}),
                                     &a, &a_tmp, &K)
  _reduce(a)
#  ccall((:fq_nmod_reduce, :libflint), Void, (Ptr{fq_nmod}, Ptr{FqNmodFiniteField}), &a, &K)
end

doc"""
   extend_easy(f::Hecke.NfOrdToFqNmodMor, K::AnticNumberField) -> NfToFqNmodMor
> For a residue field map from a prime ideal, extend the domain of the map
> to the entire field.
> Requires the prime ideal to be coprime to the index, unramified and
> over a small integer.
> Will throw a {{{BadPrime}}} exception if applied to an element in the
> field with a $p$ in the denominator or with valuation $>0$.
> Also applies to {{{FacElem{nf_elem}}}}.
"""
function extend_easy(f::Hecke.NfOrdToFqNmodMor, K::AnticNumberField)
  nf(domain(f)) != K && error("Number field is not the number field of the order")

  O = domain(f) #used for the hassert and thus the testing
  z = Hecke.NfToFqNmodMor()
  z.header.domain = K
  z.header.codomain = f.header.codomain

  p = characteristic(z.header.codomain)
  y = f(NfOrdElem(domain(f), gen(K)))
  Ft, t = PolynomialRing(ResidueRing(FlintZZ, UInt(p), cached=false), cached=false)
  K = number_field(domain(f))
#  g = gcd(Ft(K.pol), Ft(K(f.P.gen_two)))
#it would be cool to assert that g is the poly defining the codomain
  qm1 = size(codomain(f))-1

  function _image(x::NfOrdElem)
    return f(x)
  end

  Fq = codomain(f)
  s = Fq()
  t = Ft()

  function _image(x::nf_elem)
    m = denominator(x)
    if m %p == 0
      throw(BadPrime(p))
    end
    return Fq(Ft(x))
  end

  function _image(x::FacElem{nf_elem, AnticNumberField})
    r = one(Fq)
    for (k, v) = x.fac
      if v == 0 || v%qm1 == 0
        continue
      end
      if denominator(k) % p == 0
        throw(BadPrime(p))
      end
      _nf_to_fq!(s, k, Fq, t)
      if iszero(s)
        throw(BadPrime(p))
      end
      vr = v % qm1
      if vr < 0
        vr = (-vr) %qm1
        ccall((:fq_nmod_inv, :libflint), Void, (Ptr{fq_nmod}, Ptr{fq_nmod}, Ptr{FqNmodFiniteField}), &s, &s, &Fq)
      end
      ccall((:fq_nmod_pow, :libflint), Void, (Ptr{fq_nmod}, Ptr{fq_nmod}, Ptr{fmpz}, Ptr{FqNmodFiniteField}), &s, &s, &vr, &Fq)
      mul!(r, r, s)
    end
#too expensive - and automatically used in test
#    @hassert :ClassField 3 r == f(O(evaluate(x)))
    return r
  end

  function _preimage(x::fq_nmod)
    return elem_in_nf(preimage(f, x))
  end

  z.header.image = _image
  z.header.preimage = _preimage

  return z
end

doc"""
    reduce_mod_powers(a::nf_elem, n::Int) -> nf_elem
    reduce_mod_powers(a::nf_elem, n::Int, primes::Array{NfOrdIdl, 1}) -> nf_elem
> Given some non-zero algebraic integeri $\alpha$, try to find  $\beta$ s.th.
> $\beta$ is "small" and $\alpha/\beta$ is an $n$-th power.
> If the factorisation of $a$ into prime ideals is known, the ideals
> should be passed in.
"""
function reduce_mod_powers(a::nf_elem, n::Int)
  M = maximal_order(parent(a))
  p = factor(M(a)*M)
  return reduce_mod_powers(a, n, collect(keys(p)))
end

function reduce_mod_powers(a::nf_elem, n::Int, primes::Array{NfOrdIdl, 1})
  # works quite well if a is not too large. There has to be an error
  # somewhere in the precision stuff...
  @vprint :ClassField 2 "reducing modulo $(n)-th powers\n"
  @vprint :ClassField 3 "starting with $a\n"
  Zk = maximal_order(parent(a))
  val = [ div(valuation(a, x), n) for x = primes if !isone(x)]
  if all(x->x==0, val)
    I = ideal(maximal_order(parent(a)), 1)
  else
    I = prod([primes[x]^val[x] for x=1:length(primes)])
  end
  Iinv = inv(I)

  p = 100
  old = precision(BigFloat)
  r1, r2 = signature(parent(a))
  no = abs(norm(a))
  l = Int[]
  while true
   if p > 40000
      error("Something wrong in reduce_mod_powers")
    end
    setprecision(BigFloat, p)
    l = Int[]
    try
      fn = log2(BigFloat(no))/n/degree(parent(a))
      m = Hecke.minkowski(a, p)
      for i=1:r1
        push!(l, Int(BigInt(round(fn - 1/n*log2(abs(m[i]))))))
      end
      for i=1:r2
        v = Int(BigFloat(round(fn - 1/2*1/n*log2((m[r1+2*i-1]^2 + m[r1+2*i]^2)))))
        push!(l, v)
        push!(l, v)
      end
    catch e
      p *= 2
      continue
    end
    if abs(sum(l)) <= length(l)
      try
        b = Hecke.short_elem(Iinv, -matrix(FlintZZ, 1, length(l), l), prec=p)
      catch e
        if e != Hecke.LowPrecisionLLL()
          rethrow(e)
        end
        p *= 2
        if p> 40000
          println("\n\nELT\n\n", a, "\n\nn: $n $primes\n\n")
          error("too much prec")
        end
        continue
      end
#=
  N(x) = prod x_i
  N(x)^(2/n) <= 1/n sum x_i^2 <= 1/n 2^((n-1)/2) disc^(1/n) (LLL)
 so
  N(x) <= (2^((n-1)/4)/n)^(n/2) disc^(1/2)
=#
      if abs(norm(b)//norm(Iinv)) <= abs(discriminant(Zk))
        c = a*b^n
        @vprint :ClassField 3 "leaving with $c\n"
        return c
      end
      println("bad norm")
      println(abs(norm(b)//norm(Iinv)))
      println("should be")
      println(abs(discriminant(Zk)))
    end
    p *= 2
    if p> 40000
      println("abs sum ", l)
      error("too much prec")
    end
  end
  b = Hecke.short_elem(Iinv, -matrix(FlintZZ, 1, length(l), l), prec=p)
  c = a*b^n
  @vprint :ClassField 3 "leaving with $c\n"
  return c
end

function reduce_mod_powers(a::FacElem{nf_elem, AnticNumberField}, n::Int, decom::Dict{NfOrdIdl, Int})
  b = compact_presentation(a, n, decom = decom)
  b = prod([k^(v % n) for (k,v) = b.fac if !iszero(v % n)])
  b *= denominator(b, maximal_order(parent(b)))^n  #non-optimal, but integral...
  return FacElem(b)
end

function reduce_mod_powers(a::FacElem{nf_elem, AnticNumberField}, n::Int, primes::Array{NfOrdIdl, 1})
  lp = Dict((p, Int(valuation(a, p))) for p = primes)
  return reduce_mod_powers(a, n, lp)
  return b
end

function reduce_mod_powers(a::FacElem{nf_elem, AnticNumberField}, n::Int)
  Zk = maximal_order(base_ring(a))
  lp = factor_coprime(a, IdealSet(Zk))
  return reduce_mod_powers(a, n, Dict((p, Int(v)) for (p, v) = lp))
end

function rel_auto(A::ClassField_pp)
  # sqrt[n](a) -> zeta sqrt[n](a) on A.A
  @assert isdefined(A, :A)
  #on A.K, the Kummer: sqrt[n](a) = gen(A.K) -> zeta gen(A.K)
  #we have the embedding A.A -> A.K : gen(A.A) -> A.pe
  M = SMat(base_ring(A.K))
  b = A.K(1)
  push!(M, SRow(b))
  for i=2:degree(A)
    b *= A.pe
    push!(M, SRow(b))
  end
  tau = NfRelToNfRelMor(A.K, A.K, A.bigK.zeta*gen(A.K))
  N = SRow(tau(A.pe))
  C = cyclotomic_extension(base_field(A), degree(A))
  Mk = _expand(M, C.mp[1])
  Nk = _expand(N, C.mp[1])
  s = solve(Mk, Nk) # will not work, matrix non-square...
  im = A.A()
  r = degree(C.Kr)
  for (i, c) = s
    setcoeff!(im, i-1, c)
  end

  return NfRelToNfRelMor(A.A, A.A, im)
end

function rel_auto(A::ClassField)
  aut = [rel_auto(x) for x = A.cyc]
  K = number_field(A)
  g = gens(K)
  Aut = NfRel_nsToNfRel_nsMor[]
  for i=1:length(aut)
    push!(Aut, NfRel_nsToNfRel_nsMor(K, K, [j==i ? aut[i].prim_img.data(g[j]) : g[j] for j=1:length(aut)]))
  end
  return Aut
end

function extend_aut(A::ClassField, tau::T) where T <: Map
  # tau: k       -> k
  #global last_extend = (A, tau)
  k = domain(tau)
  @assert k == codomain(tau)
  @assert k == base_field(A)
  lp = factor(fmpz(degree(A)))
  all_h = [A.A() for x=A.cyc]
  for (p, v) = lp.fac
#    println("doin' $p^$v")
    Cp = [Ap for Ap = A.cyc if degree(Ap) % Int(p) == 0]
    i = 1
    om = 0
    im = 0
    while i <= length(Cp)
      if degree(Cp[i]) > om
        om = degree(Cp[i])
        im = i
      end
      i += 1
    end
    # now Cp[im] is of maximal exponent - hence, it should have the maximal
    # big Kummer extension. By construction (above), the set of s-units
    # SHOULD guarantee this....
    # Now I want all generators in terms of this large Kummer field.
    #
    # Idea: similar to pSelmer in Magma:
    #  look at frob(p, k(zeta_n)(sqrt[n](a)))(sqrt[n](a)):
    #  sqrt[n](a) -> zeta_n^i sqrt[n](a) = sqrt[n](a)^N(p) mod p
    # so N(p) = en+f => f = 1 (otherwise the congruence cannot work)
    # zeta_n^i = a^e mod p for any prime p in k(zeta_n)
    # since a = S-Unit (mod n-th powers) I can compare frob(p, a) to
    # the frob(sqrt[n](U_S)) and find the presentation:
    # next, we want the "same" for sqrt[n](tau(a)).
    # Given, that can_frob is kind of non-cheap, we want to be clever
    # zeta_n^i = tau(a)^e mod p =>
    # tau^-1(zeta_n)^i = a^e mod tau^-1(p) or
    # zeta_n^(ij) = a^e mod tau^-1(p)
    # So: need j = j(tau) and the permutation on lp...
    # in particular, tau need to be extened to k(zeta_n)
    # Cheat: g = Kr.pol => g(zeta) = 0
    #  tau(g)(zeta^r) = 0 for some suitable r
    # then zeta -> zeta^r should be a valid extension....
    #
    # TODO: Don't use the bigK.gen (the full s-units), just use the
    # group generated by the .a, the various Kummer generators, lifted
    # to the large cyclotomic extension.
    C = cyclotomic_extension(k, Int(om))
    g = C.Kr.pol
    tau_g = parent(g)([tau(coeff(g, i)) for i=0:degree(g)])
#    println("g: $g")
#    println("tau(g): $tau_g")
    i = 1
    z = gen(C.Kr)
    while gcd(i, om) != 1 || !iszero(tau_g(z))
      i *= 1
      z *= gen(C.Kr)
    end
    z_i = i

    z_i_inv = invmod(z_i, om)

    Tau = NfRelToNfRelMor(C.Kr, C.Kr, tau, z)
    tau_Ka = NfToNfMor(C.Ka, C.Ka, C.mp[1](Tau(C.mp[1]\gen(C.Ka))))
    #TODO: need the inverse of this or similar...
    # currently, this is not used as it did not work.

    lp = collect(keys(Cp[im].bigK.frob_cache))
    pp = maximum(minimum(x) for x = lp)
    S = Base.Iterators.flatten((lp, PrimeIdealsSet(order(lp[1]), pp, fmpz(-1), indexdivisors=false, ramified=false, degreebound = 1)))

    @assert C.Ka == base_ring(Cp[im].K)

    all_s = []
    all_tau_s = []
    all_emb = []
    for c in Cp
#      println("om: $om -> ", degree(c))
      Cs = cyclotomic_extension(k, Int(degree(c)))
      Emb = NfRelToNfRelMor(Cs.Kr, C.Kr, gen(C.Kr)^div(om, degree(c)))
      emb = inv(Cs.mp[1]) * Emb * C.mp[1]
      a = FacElem(Dict(emb(k) => v for (k,v) = c.a.fac))
      tau_a = FacElem(Dict(tau_Ka(k) => v for (k,v) = a.fac))
      push!(all_emb, (a, tau_a, emb))

      z = matrix(FlintZZ, 0, length(Cp[im].bigK.gen), fmpz[])
      za = Int[]
      zta = Int[]
      G = DiagonalGroup([om for i=1:length(Cp[im].bigK.gen)])
      el_z = elem_type(G)[]
      q, mq = quo(G, el_z, false)

      for p = S
        local f
        local fa
        local tfa
        try
          f = can_frobenius(p, Cp[im].bigK).coeff
          fa = can_frobenius(p, Cp[im].bigK, a)
          tfa = can_frobenius(p, Cp[im].bigK, tau_a)
        catch e
          if typeof(e) != BadPrime
            rethrow(e)
          end
          continue
        end
        el = G(f)
        if iszero(mq(el))
          continue
        end
        push!(el_z, el)
        q, mq = quo(G, el_z, false)
        z = vcat(z, f)
        push!(za, fa)
        push!(zta, tfa)
        if order(q) == 1
          break
        end
      end

      z = hcat(z, om*identity_matrix(FlintZZ, rows(z)))
      fl, s = cansolve(z, matrix(FlintZZ, rows(z), 1, za))
      @assert fl
      s = [mod(s[x, 1], om) for x=1:length(Cp[im].bigK.gen)]
      #println("s: $s")

      fl, tau_s = cansolve(z, matrix(FlintZZ, rows(z), 1, zta))
      @assert fl
      tau_s = [(z_i_inv*tau_s[x, 1]) % om for x=1:length(Cp[im].bigK.gen)]
#      println("tau(s): $tau_s")
      # so a = s -> z_i^-1 * tau_s = tau(a) (mod n) :
      push!(all_s, s)
      push!(all_tau_s, tau_s)
    end
    sG, msG = sub(Cp[im].bigK.AutG, [Cp[im].bigK.AutG(x) for x=all_s], false)
    all_b = []
    # if the above is correct, then tau_s in <s>
    for j = 1:length(Cp)
      sG, msG = sub(Cp[im].bigK.AutG, [Cp[im].bigK.AutG(all_s[x]) for x=1:length(all_s)], false)
      sG, msG = sub(Cp[im].bigK.AutG,
          vcat([(degree(Cp[j]) > om ? div(om, degree(Cp[i])):1) * Cp[im].bigK.AutG(all_s[x])
                       for x=1:length(all_s)],
               [degree(Cp[j])*Cp[im].bigK.AutG[x] for x=1:ngens(Cp[im].bigK.AutG)]), false)
      ts = all_tau_s[j]
      x = Cp[im].bigK.AutG(ts)
      fl, y = haspreimage(msG, x)
#      println(fl, " -> $(x.coeff) -> $(y.coeff)")
      @assert fl
      #need to establish the embeddings (which are also needed below)
#      println([Int(y[i]) for i=1:length(Cp)])
#      println(div(om, degree(Cp[j])))
      mu = prod(all_emb[i][1]^Int(y[i]) for i=1:length(Cp)) * inv(all_emb[j][2])
      mmu = evaluate(mu)
      #global last_rt = (mmu, degree(Cp[j]))
      rt = root(mmu, degree(Cp[j]))
      push!(all_b, (rt, y))
    end
    Ka = C.Ka
    KaT, X = PolynomialRing(Ka, "T", cached = false)
    KK, gKK = number_field([X^degree(Cp[j]) - evaluate(all_emb[j][1]) for j=1:length(Cp)])
    h = NfRel_nsToNfRel_nsMor(KK, KK, tau_Ka, [inv(all_b[i][1])*prod(gKK[j]^Int(divexact(all_b[i][2][j], div(om, degree(Cp[j])))) for j=1:length(Cp)) for i=1:length(Cp)])

    # now "all" that remains is to restrict h to the subfield, using lin. alg..
    # .. and of course move away form the Cp stuff.

    #TODO: better (more efficient) maps from NfRel -> NfRel_ns in case
    # we're using only one generator
    #Similar: NfRel_ns -> NfRel_ns small gens set -> large gens set
    all_pe =[]
    for j=1:length(Cp)
      emb = NfRelToNfRel_nsMor(Cp[j].K, KK, all_emb[j][3], gens(KK)[j])
#      println("start:")
#      println(gen(Cp[j].K), " -> ", emb(gen(Cp[j].K)))
#      println(Cp[j].K.pol, " -> ", minpoly(emb(gen(Cp[j].K))))
      pe = emb(Cp[j].pe)
      tau_pe = h(pe)
#      println("$(Cp[j].pe) pe: $pe -> $tau_pe")
#      println(norm(minpoly(Cp[j].pe)))
#      println(norm(minpoly(pe)))
#      println(norm(minpoly(tau_pe)))
#      println("=======")
      push!(all_pe, (pe, tau_pe))
    end

    B = [KK(1), all_pe[1][1]]
    d = degree(Cp[1])
    while length(B) < degree(Cp[1])
      push!(B, B[end]*all_pe[1][1])
    end

    for j=2:length(Cp)
      d *= degree(Cp[j])
      D = copy(B)
      while length(B) < d
        D = [x*all_pe[j][1] for x = D]
        append!(B, D)
      end
    end
    M = SMat(Ka)
    for i=1:d
      push!(M, SRow(B[i]))
    end
    AA, gAA = number_field([c.A.pol for c = Cp])
    @assert d == degree(AA)
    @assert d == length(B)
    b_AA = basis(AA)
    Mk = _expand(M, C.mp[1])
    #@hassert :ClassField 2 nullspace(Mk')[1] == 0
    all_im = NfRel_nsElem{nf_elem}[]
    for j=1:length(Cp)
      N = SRow(all_pe[j][2])
      Nk = _expand(N, C.mp[1])
       n = solve(Mk, Nk)
      im = sum(v*b_AA[i] for (i, v) = n)
      push!(all_im, im)
    end
    im = NfRel_nsElem{nf_elem}[]
    i = 1
    j = 1
    while j<=length(A.cyc)
      if i<= length(Cp) && degree(A.cyc[j]) == degree(Cp[i])
        push!(im, gens(A.A)[j])
        i += 1
        j += 1
      else
        j += 1
      end
    end
    emb = NfRel_nsToNfRel_nsMor(KK, A.A, im)
    i = 1
    j = 1
    while j<=length(A.cyc)
      if i<= length(Cp) && degree(A.cyc[j]) == degree(Cp[i])
        all_h[j] = emb(all_im[i])
        i += 1
        j += 1
      else
        j += 1
      end
    end
  end
  return NfRel_nsToNfRel_nsMor(A.A, A.A, tau, all_h)
end

#M is over K, mp: K -> K/k, expand M into a matrix over k
function _expand(M::Generic.Mat{nf_elem}, mp::Map)
  Kr = domain(mp)
  Ka = codomain(mp)
  k = base_ring(Kr)
  d = degree(Kr)
  N = zero_matrix(k, rows(M), cols(M) * d)
  for i=1:rows(M)
    for j = 1:cols(M)
      a = mp\M[i,j]
      for k=0:d-1
        N[i, (j-1)*d+k+1] = coeff(a, k)
      end
    end
  end
  return N
end

function _expand(M::SRow{nf_elem}, mp::Map)
  Kr = domain(mp)
  k = base_ring(Kr)
  d = degree(Kr)
  sr = SRow(k)
  for (j, v) = M
    a = mp\v
    for l=0:d-1
      c = coeff(a, l)
      if !iszero(c)
        push!(sr.pos, (j-1)*d+1+l+1)
        push!(sr.values, c)
      end
    end
  end
  return sr
end

function _expand(M::SMat{nf_elem}, mp::Map)
  Kr = domain(mp)
  k = base_ring(Kr)
  N = SMat(k)
  for i=1:rows(M)
    sr = _expand(M[i], mp)
    push!(N, sr)
  end
  return N
end

###############################################################################
#
#  Modulus function
#
###############################################################################

function defining_modulus(CF::ClassField)
  return _modulus(CF.rayclassgroupmap)
end

function defining_modulus(CF::ClassField_pp)
  return _modulus(CF.rayclassgroupmap)
end

function _modulus(mq::MapRayClassGrp)
  return mq.modulus_fin
end

function _modulus(mq::MapClassGrp)
  return ideal(order(codomain(mq)), 1)
end

###############################################################################
#
#  Auxiliary functions (to be moved)
#
###############################################################################
#=
function (I_Zk::NfOrdIdlSet)(a::NfOrdIdl)
  if parent(a) == I_Zk
    return a
  end
  Zk = order(I_Zk)
  Zl = order(a)
  @assert has_2_elem(a)
  b = ideal(Zk, a.gen_one, Zk(Zk.nf(Zl.nf(a.gen_two))))
  for i in [:gens_normal, :gens_weakly_normal, :iszero, :minimum]
    if isdefined(a, i)
      setfield!(b, i, getfield(a, i))
    end
  end
  n = divexact(degree(Zk.nf), degree(Zl.nf))
  if isdefined(a, :norm)
    b.norm = a.norm^n
  end
  if isdefined(a, :princ_gen)
    b.princ_gen = Zk(Zk.nf(Zl.nf(a.princ_gen)))
  end
  if isdefined(a, :isprime) && Zk.nf == Zl.nf && Zk.ismaximal == 1 &&
    Zl.ismaximal == 1
    b.isprime = a.isprime
    if isdefined(a, :splitting_type)
      b.splitting_type = a.splitting_type
    end
  end
  return b
end
=#

doc"""
  base_ring(A::ClassField)
> The maximal order of the field that $A$ is defined over.
"""
function base_ring(A::ClassField)
  return order(defining_modulus(A))
end

doc"""
  base_field(A::ClassField)
> The number field that $A$ is defined over.
"""
function base_field(A::ClassField)
  return number_field(base_ring(A))
end

function base_ring(A::ClassField_pp)
  return order(defining_modulus(A))
end

function base_field(A::ClassField_pp)
  return number_field(base_ring(A))
end

doc"""
  degree(A::ClassField)
> The degree of $A$ over its base field, ie. the size of the defining ideal group.
"""
function degree(A::ClassField)
  if A.degree==-1
    A.degree=Int(order(codomain(A.quotientmap)))
  end
  return A.degree
end

function degree(A::ClassField_pp)
  if A.degree==-1
    A.degree=Int(order(codomain(A.quotientmap)))
  end
  return A.degree
end

doc"""
   factor_coprime(a::FacElem{nf_elem, AnticNumberField}, I::NfOrdIdlSet) -> Dict{NfOrdIdl, fmpz}
> Factors the rincipal ideal generated by $a$ into coprimes by computing a coprime
> basis from the principal ideals in the factorisation of $a$.
"""
function factor_coprime(a::FacElem{nf_elem, AnticNumberField}, I::NfOrdIdlSet)
  Zk = order(I)
  A = Dict{NfOrdIdl, fmpz}()
  for (e,v) = a.fac
    N, D = integral_split(ideal(Zk, e))
    if !isone(N)
      if haskey(A, N)
        A[N] += v
      else
        A[N] = v
      end
    end
    if !isone(D)
      if haskey(A, D)
        A[D] -= v
      else
        A[D] = -v
      end
    end
  end
  if length(A) == 0
    A[ideal(Zk, 1)] = 1
  end
  return factor_coprime(FacElem(A))
end

doc"""
  factor(a::nf_elem, I::NfOrdIdlSet) -> Dict{NfOrdIdl, fmpz}
> Factors the principal ideal generated by $a$.
"""
function factor(a::nf_elem, I::NfOrdIdlSet)
  return factor(ideal(order(I),  a))
end

doc"""
  factor(a::FacElem{nf_elem, AnticNumberField}, I::NfOrdIdlSet) -> Dict{NfOrdIdl, fmpz}
> Factors the principal ideal generated by $a$ by refinind a coprime factorisation.
"""
function factor(a::FacElem{nf_elem, AnticNumberField}, I::NfOrdIdlSet)
  cp = factor_coprime(a, I)
  f = Dict{NfOrdIdl, fmpz}()
  for (I, v) = cp
    lp = factor(I)
    for (p, e) = lp
      f[p] = e*v
    end
  end
  return f
end

doc"""
   factor(a::fmpq, ::FlintIntegerRing) -> Fac{fmpz}
> Factor the rational number $a$ into prime numbers
"""
function factor(a::fmpq, ::FlintIntegerRing)
  fn = factor(numerator(a))
  fd = factor(denominator(a))
  for (p,e) = fd.fac
    fn.fac[p] = -e
  end
  return fn
end

doc"""
  number_field(I::NfOrd)
> Return the number fields containing $I$.
"""
@inline function number_field(O::NfOrd)
  return O.nf
end

function order(A::FacElemMon{IdealSet})
  return order(A.base_ring)
end

function order(A::FacElemMon{NfOrdIdlSet})
  return order(A.base_ring)
end


text:
"""
    handle_packets(link::Link, T)

Handle incoming packets on `link`, returning a type `T`.
"""
function handle_packets(link::Link, T)
    while true
        pkt = nextpacket(link)
        if pkt == PKT_RETURN
            return get(link, T)
        elseif pkt == PKT_TEXT
            print(get(link, String))
        elseif pkt == PKT_MESSAGE
            newpacket(link)
            @assert nextpacket(link) == PKT_TEXT
            @warn get(link, String)
        else
            error("Unsupported packet type $pkt")
        end
    end
end

function weval(link::Link, T, expr; vars...)
    if !isempty(vars)
        expr = W"With"([W"Set"(WSymbol(k), v) for (k,v) in vars], expr)
    end
    put(link, W"EvaluatePacket"(expr))
    endpacket(link)
    handle_packets(link, T)
end

"""
    weval([T,] expr; vars...)

Evaluate expression `expr`, returning a value of type `T` (default = `Any`).

If keyword arguments `vars` are provided, then `expr` is wrapped by a
[`With`](https://reference.wolfram.com/language/ref/With.html) block with `vars` assigned
as local constants.

```
julia> weval(W`Sin[x+2]`)
W"Sin"(W"Plus"(2, W"x"))

julia> weval(W`Sin[x+2]`; x=3)
W"Sin"(5)

julia> weval(W`Sin[x+2]`; x=3.0)
-0.9589242746631385
```
"""
weval(T, expr; vars...) = weval(_defaultlink(), T, expr; vars...)
weval(expr; vars...) = weval(Any, expr; vars...)

function wevalstr(link, T, str::AbstractString)
    weval(link, W"ToExpression"(str), T)
end
wevalstr(T, expr) = wevalstr(_defaultlink(), T, expr)
wevalstr(expr) = wevalstr(Any, expr)

function parseexpr(str::AbstractString)
    r = weval(W"ToExpression"(str, W"StandardForm", W"Hold"))
    r.args[1]
end

macro W_cmd(str)
    quote parseexpr($(esc(Meta.parse("\"$(escape_string(str))\"")))) end
end


text:
include_dpo()
include(joinpath(@__DIR__, "pendulum_minimum_time.jl"))

function fd(model::Pendulum{Midpoint, FreeTime}, x⁺, x, u, w, h, t)
	h = u[end]
    x⁺ - (x + h * f(model, 0.5 * (x + x⁺), u, w) + w)
end

# Nominal solution
x̄, ū = unpack(z̄, prob)
prob_nom = prob.prob

# DPO
β = 1.0
δ = 5.0e-3

# Initial samples
x1_sample = resample(x1, Diagonal([1.0; 1.0]), 1.0e-3)

# Mean problem
prob_mean = trajectory_optimization(
				model,
				EmptyObjective(),
				T,
				ul = control_bounds(model, T, [Inf; 0.0], [Inf; 0.0])[1],
				uu = control_bounds(model, T, [Inf; 0.0], [Inf; 0.0])[2],
				dynamics = false)

# Sample problems
prob_sample = [trajectory_optimization(
				model,
				EmptyObjective(),
				T,
				xl = state_bounds(model, T, x1 = x1_sample[i])[1],
				xu = state_bounds(model, T, x1 = x1_sample[i])[2],
				ul = ul,
				uu = uu,
				dynamics = false,
				con = con_free_time) for i = 1:2 * model.n]

# Sample objective
Q = [(t < T ? Diagonal(10.0 * ones(model.n))
	: Diagonal(100.0 * ones(model.n))) for t = 1:T]
R = [Diagonal(1.0e-1 * [1.0; 10.0]) for t = 1:T-1]

obj_sample = sample_objective(Q, R)
policy = linear_feedback(model.n, model.m - 1)
dist = disturbances([Diagonal(δ * ones(model.d)) for t = 1:T-1])
sample = sample_params(β, T)

prob_dpo = dpo_problem(
	prob_nom, prob_mean, prob_sample,
	obj_sample,
	policy,
	dist,
	sample)

# TVLQR policy
K, P = tvlqr(model, x̄, ū, 0.0, Q, R)

z0 = pack(z̄, K, prob_dpo)

# Solve
if true
	include_snopt() # set path
	z, info = solve(prob_dpo, copy(z0),
		nlp = :SNOPT7,
		tol = 1.0e-3, c_tol = 1.0e-3,
		time_limit = 60 * 2)
	@save joinpath(@__DIR__, "sol_dpo.jld2") z
else
	println("Loading solution...")
    @load joinpath(@__DIR__, "sol_dpo.jld2") z
end


text:
abstract QuPropagatorMethod

@doc """
QuPropagator and QuStateEvolution types are the same.

The central piece which dispatches to various solvers. The various outer constructions, allow to
give in various inputs. For example if the Hamiltonian is `QuBase.AbstractQuMatrix` and the  initial state
is `QuBase.AbstractQuVector`, we construct a `QuSchrodingerEq` with these parameters and then construct
`QuPropagator`.

### Fields :

Inputs :
* eq <: QuEquation

  Type of the equation to be solved, for example : `QuSchrodingerEq`, `QuLindbladMasterEq`
* init_state <: Union(QuBase.AbstractQuVector, QuBase.AbstractQuMatrix)

  Initial state of the system.
* tlist

  Time step array/range
* method <: QuPropagatorMethod

  Method to be used to solve the equation.

Output :
* QuPropagator construct depending on the input.
""" ->
immutable QuStateEvolution{QPM<:QuPropagatorMethod, QVM<:@compat(Union{QuBase.AbstractQuVector,QuBase.AbstractQuMatrix}), QE<:QuEquation}
    eq::QE
    init_state::QVM
    tlist
    method::QPM
    QuStateEvolution(eq, init_state, tlist, method) = new(eq, init_state, tlist, method)
end

QuStateEvolution{QPM<:QuPropagatorMethod, QV<:QuBase.AbstractQuVector, QE<: QuEquation}(eq::QE, init_state::QV, tlist, method::QPM) = QuStateEvolution{QPM,QV,QE}(eq, init_state, tlist, method)

QuStateEvolution{QPM<:QuPropagatorMethod, QM<:QuBase.AbstractQuMatrix, QE<: QuEquation}(eq::QE, init_state::QM, tlist, method::QPM) = QuStateEvolution{QPM,QM,QE}(eq, init_state, tlist, method)

QuStateEvolution{QPM<:QuPropagatorMethod, QV<:QuBase.AbstractQuVector}(eq::QuSchrodingerEq, init_state::QV, tlist, method::QPM) = QuStateEvolution{QPM,QV,QuSchrodingerEq}(eq, init_state, tlist, method)

QuStateEvolution{QPM<:QuPropagatorMethod, QV<:QuBase.AbstractQuVector}(hamiltonian::QuBase.AbstractQuMatrix, init_state::QV,  tlist, method::QPM) = QuStateEvolution{QPM,QV,QuSchrodingerEq}(QuSchrodingerEq(hamiltonian),init_state, tlist, method)

QuStateEvolution{QPM<:QuPropagatorMethod, QM<:QuBase.AbstractQuMatrix}(eq::QuLiouvillevonNeumannEq, init_state::QM, tlist, method::QPM) = QuStateEvolution{QPM,QM,QuLiouvillevonNeumannEq}(eq, init_state, tlist, method)

QuStateEvolution{QPM<:QuPropagatorMethod, QM<:QuBase.AbstractQuMatrix}(hamiltonian::QuBase.AbstractQuMatrix, init_state::QM,  tlist, method::QPM) = QuStateEvolution{QPM,QM,QuLiouvillevonNeumannEq}(QuLiouvillevonNeumannEq(liouvillian_op(hamiltonian)),init_state, tlist, method)

QuStateEvolution{QPM<:QuPropagatorMethod, QM<:QuBase.AbstractQuMatrix}(eq::QuLindbladMasterEq, init_state::QM, tlist, method::QPM) = QuStateEvolution{QPM,QM,QuLindbladMasterEq}(eq, init_state, tlist, method)

QuStateEvolution{QPM<:QuPropagatorMethod, QM<:QuBase.AbstractQuMatrix, COT<:QuBase.AbstractQuMatrix}(hamiltonian::QuBase.AbstractQuMatrix, collapse_ops::Vector{COT}, init_state::QM, tlist, method::QPM) = QuStateEvolution{QPM,QM,QuLindbladMasterEq}(QuLindbladMasterEq(hamiltonian,collapse_ops), init_state, tlist, method)

# QuLindbladMasterEqUncached is used as the construction of Lindblad operator is not required for every tracjectory.
QuStateEvolution{QPM<:QuPropagatorMethod, QV<:QuBase.AbstractQuVector, COT<:QuBase.AbstractQuMatrix}(hamiltonian::QuBase.AbstractQuMatrix, collapse_ops::Vector{COT}, init_state::QV, tlist, method::QPM) = QuStateEvolution{QPM,QV,QuLindbladMasterEq}(QuLindbladMasterEqUncached(hamiltonian,collapse_ops), init_state, tlist, method)

# QuPropagator and QuStateEvolution types are the same.
typealias QuPropagator QuStateEvolution

@doc """
QuPropagator State

Iterator version for QuPropagator used to propagate for a single time-step.

### Fields :

Inputs :
* psi

  Current state to be evolved
* t

  Current time
* t_state

  Index of next time from time step array

Output :
* QuPropagatorState construct.
""" ->
immutable QuPropagatorState
    state
    t
    t_state
end

@doc """
Iterator for QuPropagator

Iterator `start` method for QuPropagator

### Arguments

Inputs :
* prob :: QuPropagator

  The  `QuPropagator` is used to get the initial state of the system.

Output :
* QuPropagatorState construct with parameters as initial state, start time and time state.
""" ->
function Base.start(prob::QuStateEvolution)
    init_state = prob.init_state
    t_state = start(prob.tlist)
    t,t_state = next(prob.tlist,t_state)
    return QuPropagatorState(init_state,t,t_state)
end

@doc """
Iterator for QuPropagator

Iterator `next` method for QuPropagator

### Arguments

Inputs :
* prob :: QuPropagator{QuPropagatorMethod}

  The  `QuPropagator` contruct which involves the parameters of the system.
* qustate :: QuPropagatorState

  Current QuPropagatorState

Output :
* Next state, Time corresponding to the current time and related QuPropagatorState
""" ->
function Base.next{QPM<:QuPropagatorMethod}(prob::QuStateEvolution{QPM}, qustate::QuPropagatorState)
    current_qustate = qustate.state
    current_t = qustate.t
    t,t_state = next(prob.tlist, qustate.t_state)
    next_qustate = propagate(prob.method, prob.eq, t, current_t, current_qustate)
    return (t, next_qustate), QuPropagatorState(next_qustate, t, t_state)
end

@doc """
Iterator for QuPropagator

Iterator `done` method for QuPropagator

### Arguments

Inputs :
* prob :: QuPropagator

* qustate :: QuPropagatorState

Output :
* True if the current state is final state, else false.
""" ->
Base.done(prob::QuStateEvolution, qustate::QuPropagatorState) = done(prob.tlist, qustate.t_state)

@doc """
Propagation of the time evolution operator.

### Arguments

Inputs :
* op <: QuBase.AbstractQuMatrix

  Operator to be evolved.
* dt :: Float64

  Time step.

* tf, ti :: Float64

  Final Time, Initial Time.

Output :
* Evolved operator.
""" ->
QuEvolutionOp{QM<:QuBase.AbstractQuMatrix}(op::QM, dt::Float64) = expm(-im*op*dt)

QuEvolutionOp{QM<:QuBase.AbstractQuMatrix}(op::QM, tf::Float64,  ti::Float64) = QuEvolutionOp(op, tf-ti)

QuEvolutionOp{QE<:QuEquation}(eq::QE, dt::Float64) = QuEvolutionOp(operator(eq), dt)

QuEvolutionOp{QE<:QuEquation}(eq::QE, tf::Float64, ti::Float64) = QuEvolutionOp(eq, tf-ti)

function Base.show(io::IO, qprop::QuPropagator)
    field_params = fieldnames(qprop.eq)
    println(io, "Summarizing the system :")
    if :lindblad in field_params && :hamiltonian in field_params
        println(io, "Equation type : $(typeof(qprop.eq))")
        println(io, "Size of the Lindblad operator of the system : $(size(coeffs(qprop.eq.lindblad)))")
        println(io, "Size of the Hamiltonian of the system : $(size(coeffs(qprop.eq.hamiltonian)))")
        println(io, "Number of collapse operators : $(length(qprop.eq.collapse_ops))")
        println(io, "Size of the Density matrix : $(size(coeffs(qprop.init_state)))")
    elseif :hamiltonian in field_params
        println(io, "Equation type : $(typeof(qprop.eq))")
        println(io, "Size of the Hamiltonian of the system : $(size(coeffs(qprop.eq.hamiltonian)))")
        println(io, "Size of the Initial state : $(size(coeffs(qprop.init_state)))")
    elseif :liouvillian in field_params
        println(io, "Equation type : $(typeof(qprop.eq))")
        println(io, "Size of the Liouvillian of the system : $(size(coeffs(qprop.eq.liouvillian)))")
        println(io, "Size of the Density matrix : $(size(coeffs(qprop.init_state)))")
    end
    println(io, "Time steps used : $(qprop.tlist)")
    println(io, "Solver used : $(qprop.method)")
end

export  QuStateEvolution,
      QuPropagator,
      QuEvolutionOp,
      QuPropagatorState


text:
function ALpreserve(; verbose = true)
    setpreserver!(activehistory)
    verbose && println("Preserving State")
end

function ALrestore(; verbose = true)
    (preserver === missing) && (println("Nothing to restore"); return)

    ((length(replset.commands)==0) || (replset.commands[end]=="#Session Started"))
        popsessionhistory!()

    setreplset!(activelogicset(preserver))

    setactivehistory!(preserver)

    pushsessionhistory!(replset)

    verbose && println("Restoring State - " * reportfeasible())
end


text:
function venkat_limiter(qtilde, loc_globaldata, idx, configData, phi)
    # smallest_dist(globaldata, idx)
    VL_CONST = configData["core"]["vl_const"]::Float64
    ds = loc_globaldata[idx].short_distance
    epsi = VL_CONST * ds
    epsi = epsi ^ 3
    # phi = zeros(Float64, 4)
    del_pos = zero(Float64)
    del_neg = zero(Float64)
    for i in 1:4
        q = copy(loc_globaldata[idx].q[i])
        del_neg = qtilde[i] - q
        if abs(del_neg) <= 1e-5
            phi[i] = 1.0
        elseif abs(del_neg) > 1e-5
            if del_neg > 0.0
                # maximum(globaldata, idx, i)
                del_pos = loc_globaldata[idx].max_q[i] - q
            elseif del_neg < 0.0
                # minimum(globaldata, idx, i)
                del_pos = loc_globaldata[idx].min_q[i] - q
            end
            num = (del_pos*del_pos) + (epsi*epsi)
            num = (num*del_neg) + 2 * (del_neg*del_neg*del_pos)

            den = (del_pos*del_pos) + (2 *del_neg*del_neg)
            den = den + (del_neg*del_pos) + (epsi*epsi)
            den = den*del_neg

            temp = num/den
            if temp < 1.0
                phi[i] = temp
            else
                phi[i] = 1.0
            end
        end
    end
    return nothing
end

# @inline function maximum(globaldata, idx, i)
#     globaldata[idx].max_q[i] = copy(globaldata[idx].q[i])
#     for itm in globaldata[idx].conn
#         if globaldata[idx].max_q[i] < globaldata[itm].q[i]
#             globaldata[idx].max_q[i] = globaldata[itm].q[i]
#         end
#     end
#     return  nothing
# end

# @inline function minimum(globaldata, idx, i)
#     globaldata[idx].min_q[i] = copy(globaldata[idx].q[i])
#     for itm in globaldata[idx].conn
#         if globaldata[idx].min_q[i] > globaldata[itm].q[i]
#             globaldata[idx].min_q[i] = globaldata[itm].q[i]
#         end
#     end
#     return nothing
# end

@inline function smallest_dist(globaldata, idx)
    min_dist = 1000.0
    for itm in globaldata[idx].conn
        ds = hypot(globaldata[idx].x - globaldata[itm].x, globaldata[idx].y - globaldata[itm].y)
        if ds < min_dist
            min_dist = ds
        end
    end
    globaldata[idx].short_distance = min_dist
    return nothing
end

@inline function max_q_values(globaldata, idx)
    maxq = copy(globaldata[idx].q)
    for itm in globaldata[idx].conn
        currq = globaldata[itm].q
        for i in 1:4
            if maxq[i] < currq[i]
                maxq[i] = currq[i]
            end
        end
    end
    return maxq
end

@inline function min_q_values(globaldata, idx)
    minq = copy(globaldata[idx].q)
    for itm in globaldata[idx].conn
        currq = globaldata[itm].q
        for i in 1:4
            if minq[i] > currq[i]
                minq[i] = currq[i]
            end
        end
    end
    return minq
end

@inline function qtilde_to_primitive(result::Array{Float64,1}, qtilde::Array{Float64,1}, configData)
    gamma::Float64 = configData["core"]["gamma"]
    beta = -qtilde[4]*0.5
    temp = 0.5/beta
    u1 = qtilde[2]*temp
    u2 = qtilde[3]*temp

    temp1 = qtilde[1] + beta*(u1*u1 + u2*u2)
    temp2 = temp1 - (log(beta)/(gamma-1))
    rho = exp(temp2)
    pr = rho*temp
    result[1] = u1
    result[2] = u2
    result[3] = rho
    result[4] = pr
    return  nothing
end

max_stars_repo_name: TestSubjector/HybridMeshfreeSolver

max_issues_repo_name: TestSubjector/HybridMeshfreeSolver

max_forks_repo_name: TestSubjector/HybridMeshfreeSolver


text:
using BinaryBuilder, Pkg

# The version of this JLL is decoupled from the upstream version.
# Whenever we package a new upstream release, we initially map its
# version X.Y.Z to X00.Y00.Z00 (i.e., multiply each component by 100).
# So for example version 2.6.3 would become 200.600.300.
#
# Together, this allows to increment the patch level of the JLL for minor tweaks.
# If a rebuild of the JLL is needed which keeps the upstream version identical
# but breaks ABI compatibility for any reason, one can increment the minor or major
# version (depending on whether package using this JLL use `~` or `^` compat entries)
# e.g. go from 200.600.300 to 200.601.300 or 201.600.300
# Similar tricks can also be used to package prerelease versions; e.g. one might
# map a prerelease of 2.7.0 to 200.690.000.

name = "METIS4"
upstream_version = v"4.0.3"
version_offset = v"0.0.0" # reset to 0.0.0 once the upstream version changes
version = VersionNumber(upstream_version.major * 100 + version_offset.major,
                        upstream_version.minor * 100 + version_offset.minor,
                        upstream_version.patch * 100 + version_offset.patch)

# Collection of sources required to build METIS
sources = [
    ArchiveSource("http://glaros.dtc.umn.edu/gkhome/fetch/sw/metis/OLD/metis-$(upstream_version).tar.gz",
                  "5efa35de80703c1b2c4d0de080fafbcf4e0d363a21149a1ad2f96e0144841a55"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
if [[ "${target}" == *-mingw* ]]; then
  COPTIONS="${COPTIONS} -D__VC__"  # to resolve missing srand48/drand48 symbols
fi

# build libmetis.a
cd $WORKSPACE/srcdir/metis-*
for f in ${WORKSPACE}/srcdir/patches/*.patch; do
  atomic_patch -p1 ${f}
done
cd Lib
make -j${nproc} COPTIONS="${COPTIONS} -fPIC"
cd ..

# We copy the .a files into ${prefix}/lib since the main purpose is to link them in other builds.
# Specifically this is in a separate location than the typical location for libraries on Windows.
mkdir -p ${prefix}/lib
mv libmetis.a ${prefix}/lib
mkdir -p ${prefix}/include
cp Lib/metis.h ${prefix}/include
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    FileProduct("lib/libmetis.a", :libmetis_a)
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[]

# Build the tarballs
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies;
               julia_compat="1.6", preferred_gcc_version=v"6")

max_stars_repo_name: MichelJuillard/Yggdrasil

max_issues_repo_name: MichelJuillard/Yggdrasil

max_forks_repo_name: MichelJuillard/Yggdrasil


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Find crossing of a function.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

"""
    find_crossing(f::Function, t₀::Number, t₁::Number, s₀, s₁; Δ = 1e-3, max = 100)

Return the crossing time `tc` in which the function `f(t)` goes from the state
`s₀` to the state `s₁`. It is assumed that `f(t₀) = s₀` and `f(t₁) = s₁`.

If the computed interval is smalled than `Δ`, or if the number of iterations is
higher than `max`, then the algorithm stops.

# Examples

```julia-repl
julia> SatelliteToolbox.find_crossing(
    t -> (sin(t) > 0),
    -0.3,
    0.3,
    false,
    true;
    Δ = 1e-10
)
6.984919309616089e-11
```
"""
function find_crossing(
    f::Function,
    t₀::Number,
    t₁::Number,
    s₀,
    s₁,
    vargs...;
    Δ = 1e-3,
    max = 100
)
    it = 0

    T = typeof((t₁ + t₀) / 2)

    while it <= max
        # Call the function at the middle of the interval.
        ti = (t₁ + t₀) / 2
        si = f(ti, vargs...)

        # Compute the new interval.
        if si == s₀
            t₀ = ti
        elseif si == s₁
            t₁ = ti
        else
            error("The function `f` returned an unexpected state.")
        end

        # If the interval is small enough, then just return.
        (t₁ - t₀) < Δ && break

        it += 1
    end

    return T(t₁)
end


text:
function erdos_renyi(n::Int64, p::Float64; seed::Integer=-1)
    if seed > 0
        Random.seed!(seed)
    end
    g = Graph();
    add_vertex!(g, n);
    m = round(Int, n * (n - 1) / 2);
    ne = rand(Binomial(m, p));
    g = Graph(n, ne);
    return g
end


text:
using GeoMakie, MakieGallery, AbstractPlotting.MakieLayout, GLMakie, Test

const OPENGL = GLMakie.WORKING_OPENGL

@show OPENGL

OPENGL && begin @info "OpenGL detected"; using GLMakie; end
OPENGL || begin @warn "No OpenGL detected!  Using CairoMakie instead."; using CairoMakie; CairoMakie.activate!(); end

# Download reference images from master
MakieGallery.current_ref_version[] = "master"

# Load only the GeoMakie examples
MakieGallery.load_database(["geomakie.jl"]);

filter!(MakieGallery.database) do entry
    !(entry.unique_name ∈ (:air_particulates,))
end

tested_diff_path = joinpath(@__DIR__, "tested_different")
test_record_path = joinpath(@__DIR__, "test_recordings")

isdir(tested_diff_path) && rm(tested_diff_path; force = true, recursive = true)
mkpath(tested_diff_path)

isdir(test_record_path) && rm(test_record_path; force = true, recursive = true)
mkpath(test_record_path)

examples = MakieGallery.record_examples(test_record_path)

@test length(examples) == length(MakieGallery.database)

# MakieGallery.generate_preview(test_record_path, joinpath(homedir(), "Desktop", "index.html"))
# MakieGallery.generate_thumbnails(test_record_path)
# MakieGallery.gallery_from_recordings(test_record_path, joinpath(test_record_path, "index.html"))

printstyled("Running ", color = :green, bold = true)
println("visual regression tests")

MakieGallery.run_comparison(test_record_path, tested_diff_path)


text:
cd(@__DIR__)

using Pkg

using Test, TestSetExtensions, SafeTestsets
using SearchLight
using SearchLightMySQL
using Dates

# @testset ExtendedTestSet "SearchLight PostgreSQL adapter tests" begin
#   @includetests ARGS
# end

# run a simple connect test for the first time

@testset ExtendedTestSet "SearchLight tests" begin
    @includetests ARGS
end

max_stars_repo_name: GenieFramework/SearchLightMySQL

max_issues_repo_name: GenieFramework/SearchLightMySQL

max_forks_repo_name: GenieFramework/SearchLightMySQL


text:
"""
    distance(x, y)

Computes euclidian distance between `x` and `y`.
"""
distance(x, y) = norm(x .- y, 2)

"""
    direction(x, y; normed = true)

Computes direction from `x` to `y`.
"""
function direction(x, y; normed = true)
    d = y .- x
    return normed ? d./norm(d) : d
end

"""
    direction_angle(d)

Computes angle in radians of the direction `d`.
"""
direction_angle(d) = atan(d[2], d[1])

"""
    direction_angle(x, y)

Computes angle in radians of the direction from `x` to `y`.
"""
direction_angle(x, y) = atan(direction(x, y))

"""
    euclidian(r::Real, ϕ::Real, origin = (0, 0))

Computes euclidian coordinates from the given polar coordinates coordinates `(r, ϕ)`.
"""
function euclidian(r::Real, ϕ::Real, origin = (0, 0))
    r >= 0 || throw(ArgumentError("r must be greater or equal to 0."))
    ϕ = mod1(ϕ + π, 2π) - π
    return origin .+ r .* (cos(ϕ), sin(ϕ))
end

"""
    polar(x::Real, y::Real)

Computes polar coordinates from the given 2D euclidian coordinates `(x, y)`.
"""
polar(x::Real, y::Real) = (hypot(x, y), atan(y, x))


text:
module Neuron

using CUDA

include("lif.jl")
include("srm0.jl")

end


text:
# Visual demonstration of the tracing scheme in different analytical field.
using PyPlot, FieldTracer

include("../../test/utility/dipole.jl")

"""
    test_trace_asymptote(issingle=false)

Test streamline tracing by plotting vectors and associated streamlines through a simple
velocity field where Vx=x, Vy=-y. Support for single and double precision data.
"""
function test_trace_asymptote(issingle=false)

   # Start by creating a velocity vector field.
   if issingle
      xmax, ymax = 200.0f0, 20.0f0
      x = -10.0f0:0.25f0:xmax
      y = -10.0f0:0.25f0:ymax
   else
      xmax, ymax = 200.0, 20.0
      x = -10.0:0.25:xmax
      y = -10.0:0.25:ymax
   end

   xgrid = [i for j in y, i in x]
   ygrid = [j for j in y, i in x]

   if issingle
      vx = xgrid * 1.0f0
      vy = ygrid * -1.0f0

      xstart = 1.0f0
      ystart = 10.0f0
   else
      vx = xgrid * 1.0
      vy = ygrid * -1.0

      xstart = 1.0
      ystart = 10.0
   end

   gridtype = "meshgrid"

   x1, y1 = trace(vx, vy, xstart, ystart, x, y; alg=RK4(), ds=0.1, gridtype)
   x2, y2 = trace(vx, vy, xstart, ystart, x, y; alg=RK4(), ds=0.5, gridtype)
   x3, y3 = trace(vx, vy, xstart, ystart, x, y; alg=RK4(), ds=1.0, gridtype)
   x4, y4 = trace(vx, vy, xstart, ystart, x, y; alg=Euler(), ds=0.1, gridtype)
   x5, y5 = trace(vx, vy, xstart, ystart, x, y; alg=Euler(), ds=0.5, gridtype)
   x6, y6 = trace(vx, vy, xstart, ystart, x, y; alg=Euler(), ds=1.0, gridtype)

   # analytical solution const = x*y
   c = xstart * ystart
   if issingle
      x_anly = 1.0f0:0.001f0:xmax
   else
      x_anly = 1.0:0.001:xmax
   end
   y_anly = c ./ x_anly

   fig = plt.figure(figsize=(8,6))
   ax1 = plt.subplot(111)
   ax1.plot(x_anly, y_anly, "r--", label="Analytic",linewidth=3.0)
   ax1.plot(x1, y1, "b",   label="RK4 ds=0.1", linewidth=1.5)
   ax1.plot(x2, y2, "b--", label="RK4 ds=0.5", linewidth=1.5)
   ax1.plot(x3, y3, "b:",  label="RK4 ds=1.0", linewidth=1.5)
   ax1.plot(x4, y4, "g", label="Euler ds=0.1", linewidth=.75)
   ax1.plot(x5, y5, "g--", label="Euler ds=0.5", linewidth=.75)
   ax1.plot(x6, y6, "g:",  label="Euler ds=1.0", linewidth=.75)
   ax1.legend(loc="upper right")
   if issingle
      ax1.set_title("Runge Kutta 4 vs Euler's: Asymptotic Field, Single Precision")
   else
      ax1.set_title("Runge Kutta 4 vs Euler's: Asymptotic Field, Double Precision")
   end
   ax1.set_xlabel("Normalized X Coordinate")
   ax1.set_ylabel("Normalized Y Coordinate")
   ax1.set_xlim([3, 30])
   ax1.set_ylim([0.25, 2.5 ])

   ax1.annotate("Euler's method diverges strongly\nalong curves except when"*
      " taking \nvery small steps.  RK4 is\nfar more accurate for all dS.",
      xy=(10.5,0.83), xycoords="data", xytext=(9,.3),
      arrowprops=Dict(:fc=>"black",:shrink=>0.05),
      horizontalalignment="center")
   ax1.annotate("Tracing begins at x=1, y=10.",
      xy=(4.8,2.45), xycoords="data", xytext=(6,1.9),
      arrowprops=Dict(:fc=>"black",:shrink=>0.05))
   ax1.annotate("The trace continues until \n"*
      "x=200.  At that point, Euler \n"*
      "dS=0.1 and RK4 dS=1.0 converge\n"*
      " at the same point, despite \n"*
      "the 10X difference in step size.",
      xy=(29,.5), xycoords="data", xytext=(20,1),
      arrowprops=Dict(:fc=>"black",:shrink=>0.005),
      horizontalalignment="center")
end

"Trace field lines through a Earth like magnetic dipole field."
function test_trace_dipole()
   # Start by creating a field of unit vectors.
   x = -100.0:5.0:101.0
   y = -100.0:5.0:101.0

   bx, by = b_hat(x, y)

   fig2 = plt.figure(figsize=(12,6))
   fig2.subplots_adjust(wspace=0.15, left=0.08, right=0.94)
   ax2 = plt.subplot(121)
   ax3 = plt.subplot(322)
   ax4a= plt.subplot(347)
   ax4b= plt.subplot(348)
   ax5 = plt.subplot(326)
   ax2.quiver(x,y, bx, by, units="x", pivot="middle")
   ax3.quiver(x,y, bx, by, units="x", pivot="tip")
   ax4a.quiver(x,y,bx, by, units="x", pivot="tip")
   ax4b.quiver(x,y,bx, by, units="x", pivot="tip")
   ax5.quiver(x,y, bx, by, units="x", pivot="tip")

   # Trace through this field.
   xstart = 10.0 # ystart = 25.0
   ds = 0.1
   gridtype = "meshgrid"
   for ystart in 0.:5.:31.
      x1, y1 = trace(bx, by, xstart, ystart, x, y; alg=RK4(), ds, gridtype)
      l1 = ax2.plot(x1,y1,"b")[1]
      ax3.plot(x1,y1,"b"); ax4b.plot(x1,y1,"b")
      ax5.plot(x1,y1,"b"); ax4a.plot(x1,y1,"b")
      x2, y2 = trace(bx, by, xstart, ystart, x, y; alg=Euler(), ds, gridtype)
      l2 = ax2.plot(x2,y2,"r")[1]
      ax3.plot(x2,y2,"r"); ax4b.plot(x2,y2,"r")
      ax5.plot(x2,y2,"r"); ax4a.plot(x2,y2,"r")
      x3, y3 = b_line(xstart, ystart, npoints=300)
      l3 = ax2.plot(x3,y3,"k--")[1]
      ax3.plot(x3,y3,"k--"); ax4b.plot(x3,y3,"k--")
      ax5.plot(x3,y3,"k--"); ax4a.plot(x3,y3,"k--")

      ax2.set_xlim([-2,  100])
      ax2.set_ylim([-30, 100])
      ax2.set_title("Full View")
      ax2.set_xlabel("Normalized X Coordinate")
      ax2.set_ylabel("Normalized Y Coordinate")
      ax2.legend((l1, l2, l3),("RK4", "Euler", "Analytical"), loc="upper left")

      ax3.set_title("Zoomed Views")
      ax3.set_xlim([8.5, 17.5])
      ax3.set_ylim([3, 33])
      goodpos = ax3.get_position()

      ax4a.set_xlim([20,30])
      ax4a.set_ylim([-12,12])
      pos = ax4a.get_position()
      pos.x0 = goodpos.x0
      pos.x1 = pos.x0 + (goodpos.x1-goodpos.x0)/2.0 -0.01
      ax4a.set_position(pos)

      ax4b.set_xlim([50,60])
      ax4b.set_ylim([-12,12])
      pos = ax4b.get_position()
      pos.x0 = goodpos.x0 + (goodpos.x1-goodpos.x0)/2.0 +0.01
      ax4b.set_position(pos)
      ax4b.yaxis.set_ticklabels([])

      ax5.set_xlim([1,7])
      ax5.set_ylim([-7,-3])
      ax5.set_xlabel("Normalized X Coordinate")

      fig2.suptitle("RK4 vs Euler's Method: Dipole Field for dS="*
         string(round(ds,digits=3)))
   end
end


text:
export sNES, snes_init, snes_populate, snes_generation, populate, evaluate, generation

# adapted from
# https://github.com/francescoalemanno/NaturalES.jl/blob/master/src/separable_nes.jl

import Cambrian: populate, evaluate, generation

mutable struct sNESState <: ESState
    μ::Array{Float64}
    σ::Array{Float64}
    u::Array{Float64}
    s::Array{Float64}
end

function sNESState(d::Int, n::Int;
                   μ::Array{Float64}=randn(d),
                   σ::Array{Float64}=ones(d))
    u = zeros(n)
    u = max.(0, log(n/2+1) .- log.(1:n))
    u .= u ./ sum(u) .- 1/n
    s = randn(d, n)
    sNESState(μ, σ, u, s)
end

function snes_config(d::Int; n_population = 4+3*ceil(Int, log(d)),
                     ημ = 1.0, ησ = (3 + log(d)) / (5 * sqrt(d)))
    (n_population = n_population, ημ = ημ, ησ = ησ)
end

mutable struct sNES <: Cambrian.AbstractEvolution
    config::NamedTuple
    logger::CambrianLogger
    population::Array{AbstractESIndividual}
    elites::Array{AbstractESIndividual}
    state::sNESState
    fitness::Function
    gen::Int
end

function snes_init(cfg::NamedTuple, state::ESState; T::Type=ESIndividual)
    population = Array{AbstractESIndividual}(undef, cfg.n_population)
    for i in 1:cfg.n_population
        genes = state.μ .+ state.σ .* view(state.s, :, i)
        population[i] = T(genes, -Inf*ones(cfg.d_fitness))
    end
    population
end

function sNES(cfg::NamedTuple, fitness::Function, state::sNESState; T::Type=ESIndividual, logfile=string("logs/", cfg.id, ".csv"))
    logger = CambrianLogger(logfile)
    population = snes_init(cfg, state, T=T)
    elites = deepcopy([population[i] for i in 1:cfg.n_elite])
    sNES(cfg, logger, population, elites, state, fitness, 0)
end

"""
Separable Natural Evolutionary Strategies default constructor
Will use a random starting point using N(0, 1)
If cfg contains keys from snes_config, these will be used
To provide initial state, see sNES(cfg, fitness, state)
"""
function sNES(cfg::NamedTuple, fitness::Function; T::Type=ESIndividual, logfile=string("logs/", cfg.id, ".csv"))
    logger = CambrianLogger(logfile)
    cfg = merge(snes_config(cfg.n_genes), cfg)
    state = sNESState(cfg.n_genes, cfg.n_population)
    sNES(cfg, fitness, state; T=T, logfile=logfile)
end

"generate next population"
function snes_populate(e::sNES)
    for i in eachindex(e.population)
        e.population[i].genes .= e.state.μ .+ e.state.σ .* view(e.state.s, :, i)
        e.population[i].fitness .= -Inf
    end
end

"update NES state, called after populate and evaluate"
function snes_generation(e::sNES)
    d = e.config.n_genes
    n = e.config.n_population

    # copy population information
    F = zeros(n)
    for i in eachindex(e.population)
        F[i] = -e.population[i].fitness[1]
    end
    idx = sortperm(F)

    # compute gradients
    ∇μ = zeros(d)
    ∇σ = zeros(d)
    for i in 1:n
        j = idx[i]
        ∇μ .+= e.state.u[i] .* e.state.s[:, j]
        ∇σ .+= e.state.u[i] .* (e.state.s[:, j].^2 .- 1.0)
    end

    # update state variables
    e.state.μ .+= e.config.ημ .* e.state.σ .* ∇μ
    e.state.σ .*= exp.(e.config.ησ/2 .* ∇σ)
    randn!(e.state.s)
    Cambrian.elites_generation(e)
end

populate(e::sNES) = snes_populate(e)
evaluate(e::sNES) = fitness_evaluate(e, e.fitness)
generation(e::sNES) = snes_generation(e)


text:
#region import and models
using Pkg
Pkg.activate(".")
Pkg.instantiate()
using DrWatson
using JLD2
using Distributed
using SharedArrays
using ScoreDrivenERGM
using Logging

begin
nWorkers = 8
addprocs(nWorkers - nprocs())
@sync @everywhere begin
    using Pkg
    Pkg.activate(".")
    Pkg.instantiate()
    using ScoreDrivenERGM
    import ScoreDrivenERGM:StaticNets, DynNets
    import ScoreDrivenERGM.DynNets:SdErgm,SdErgmDirBin0Rec0, simulate_and_estimate_parallel
    using ScoreDrivenERGM.Utilities

    model_mle = DynNets.SdErgmDirBin0Rec0_mle(scoreScalingType="FISH_D")

    end
end

#endregion

model_mle = DynNets.SdErgmPml(staticModel = StaticNets.NetModeErgmPml("edges +  gwesp(decay = 0.25,fixed = TRUE)", true), indTvPar = [true, true], scoreScalingType="FISH_D")

#region simulations
dgpSetARlowlow, dgpSetARlow, dgpSetARmed, dgpSetARhigh, dgpSetSIN, dgpSetSDlow, dgpSetSD, dgpSetSDhigh = ScoreDrivenERGM.DynNets.list_example_dgp_settings(model_mle)

sigmaVals = round.((10).^(range(log10(0.005), log10(0.1), length=6)), sigdigits=1)
dgpList = [deepcopy(dgpSetARlowlow) for i in 1:length(sigmaVals)]
setindex!.(getfield.(getfield.(dgpList, :opt),:sigma), sigmaVals)

c= Dict{String, Any}()
c["model"] =[DynNets.SdErgmDirBin0Rec0_mle(scoreScalingType="FISH_D")]
c["T"] = [500]
c["N"] = [500]
c["dgpSettings"] = dgpList
c["nSample"] = 50

list = sort(dict_list(c), by=x->(string(x["model"])))

for d in list

    timeSim = @elapsed allObsT, allvEstSdResPar, allfVecT_filt, allParDgpT, allConvFlag, allftot_0,  allfVecT_filt_SS =  ScoreDrivenERGM.DynNets.simulate_and_estimate_parallel(d["model"], d["dgpSettings"], d["T"], d["N"],  d["nSample"];)

    modelTag = string(d["model"])

    res1 =  (;modelTag, allObsT, allvEstSdResPar, allfVecT_filt, allParDgpT, allConvFlag, allftot_0, allfVecT_filt_SS) |>DrWatson.ntuple2dict |> DrWatson.tostringdict

    estDict = merge(res1, d)

    saveName = replace.( savename(d, "jld2";allowedtypes = (Real, String, Symbol, NamedTuple, Tuple, ScoreDrivenERGM.DynNets.SdErgm) ), r"[\"]" => "")

    timeSave = @elapsed save( datadir("sims", "dgpAR_Fil_SS", saveName), estDict)

    Logging.@info("Time sim = $timeSim ,  time save = $timeSave ")

    # res2 =  (;allAT ) |>DrWatson.ntuple2dict |> DrWatson.tostringdict
    # matsDict = merge(res2, d)
    # @tagsave( datadir("sims", "sampleDgpFilterSD_sampled_mats", saveName), matsDict)

end

#endregion

# #region load and plot estimates
begin
using Statistics
using DataFrames
using PyPlot
import ScoreDrivenERGM.Utilities.AReg

rolling_mean(x::Array{Float64}, size = 3) = vcat(fill(mean(x[i:i+size]), size), [mean(x[i:i+size]) for i =1:length(x)-size])

@elapsed df = collect_results( datadir("sims", "dgpAR_Fil_SS"))
df["modelTag"] = string.(df["model"])

N = 100
modTag = string(model_mle)

df = df[df.N .== N, :]
df = transform(df, :allfVecT_filt => ByRow(x-> dropdims(std(x, dims=2), dims=2)) => :std_SD )

df = transform(df, :allfVecT_filt_SS => ByRow(y-> dropdims(mapslices( (x-> AReg.fitARp(Float64.(x), 1)[2] ), y, dims=(2) ), dims=2)) => :std_SS )

df = transform(df, :allfVecT_filt => ByRow(y-> dropdims(mapslices( (x-> AReg.fitARp(Float64.(x), 1)[2] ), y, dims=(2) ), dims=2)) => :std_SD )

df = transform(df, :allParDgpT => ByRow(y-> dropdims(mapslices( (x-> AReg.fitARp(Float64.(x), 1)[2] ), y, dims=(2) ), dims=2)) => :std_DGP_est )

df = transform(df, :dgpSettings => ByRow(x-> getfield(getfield(x,:opt), :sigma)[1]) => :sigma )

end

begin
df = sort(df, :sigma)
qVals = unique(getfield.(getfield.(df.dgpSettings,:opt), :sigma))
fig, ax = plt.subplots(nrows=2, ncols=length(qVals), figsize=(12, 8), sharey = "row")

labNames = ["theta", "eta"]
for i=1:nrow(df)
    for p=1:2
        all_data = ([df.std_SS[i][p,:] df.std_SD[i][p,:] df.std_DGP_est[i][p,:]] .-df.sigma[i])#./df.sigma[i]
        # plot violin plot
        ax[p, i].boxplot(all_data,
                        showmeans=true, showfliers=false)
        xlims = ax[p, i].get_xlim()
        ax[p, i].hlines(0, xlims[1], xlims[2], linestyle=":" , colors = "r")

        ax[p,i].set_xticklabels(["SS", "SD", "DGP"])
        ax[p,i].grid()
        ax[p,1].set_ylabel(labNames[p])
        ax[1,i].set_title("var dgp = $(qVals[i][1]) ")
    end
end

fig.suptitle("Estimate Variance of Latent AR1 in reciprocity p* model N = $N")
plt.tight_layout()

plt.savefig("C:\\Users\\digan\\Dropbox\\Dynamic_Networks\\repos\\ScoreDrivenExponentialRandomGraphs\\_research\\analysis_for_paper_revision\\plots_for_revision\\estimate_var_latent_AR1_N_$(N)_$modTag")
end
#endregion

max_stars_repo_name: domenicodigangi/SDERGM

max_issues_repo_name: domenicodigangi/SDERGM

max_forks_repo_name: domenicodigangi/SDERGM


text:
# This file is a part of SimilaritySearch.jl

export CosineDistance, AngleDistance, NormalizedCosineDistance, NormalizedAngleDistance
using LinearAlgebra
import Distances: evaluate

"""
   CosineDistance()

The cosine is defined as:
```math
\\cos(u, v) = \\frac{\\sum_i u_i v_i}{\\sqrt{\\sum_i u_i^2} \\sqrt{\\sum_i v_i^2}}
```

The cosine distance is defined as ``1 - \\cos(u,v)``
"""
struct CosineDistance <: PreMetric end

"""
   AngleDistance()

The angle distance is defined as:
```math
∠(u, v)= \\arccos(\\cos(u, v))
```

"""
struct AngleDistance <: PreMetric end
"""
    NormalizedCosineDistance()

Similar to [`CosineDistance`](@ref) but suppose that input vectors are already normalized

```math
1 - \\sum_i {u_i v_i}
```

"""
struct NormalizedCosineDistance <: PreMetric end

"""
    NormalizedAngleDistance()

Similar to [`AngleDistance`](@ref) but suppose that input vectors are already normalized

```math
\\arccos \\sum_i {u_i v_i}
```

"""
struct NormalizedAngleDistance <: PreMetric end

const π_2 = π / 2

function fastacos(d)
    if d <= -1.0
        π
    elseif d >= 1.0
        0.0
    elseif d == 0  # turn around for zero vectors, in particular for denominator=0
        π_2
    else
        acos(d)
    end
end

"""
    evaluate(::NormalizedCosineDistance, a, b)

Computes the cosine distance between two vectors, it expects normalized vectors (see [normalize!](@ref) method).
Please use NormalizedAngleDistance if you are expecting a metric function (cosine_distance is a faster
alternative whenever the triangle inequality is not needed)
"""
evaluate(::NormalizedCosineDistance, a, b) = one(eltype(a)) - dot(a, b)

"""
    evaluate(::AngleDistance, a, b)

Computes the angle  between twovectors. It supposes that all vectors are normalized (see `normalize!` function)

"""
evaluate(::NormalizedAngleDistance, a, b) = fastacos(dot(a, b))

"""
    evaluate(::CosineDistance, a, b)

Computes the cosine distance between two vectors.
Please use AngleDistance if you are expecting a metric function (cosine_distance is a faster
alternative whenever the triangle inequality is not needed)
"""
evaluate(::CosineDistance, a, b) = one(eltype(a)) - dot(a, b) / (norm(a) * norm(b))

"""
    evaluate(::AngleDistance, a, b)

Computes the angle  between twovectors.

"""
function evaluate(::AngleDistance, a, b)
    d = dot(a, b) / (norm(a) * norm(b))
    fastacos(d)
end


text:
"""
Type for storing data acquired with a transient response (TR) experimental paradigm.

#### Example

Put an example here

```julia
s = read_TR("filename")
s = rereference(s, "Cz")
```

"""
mutable struct TR <: EEG
    data::Array
    sensors::Array{Sensor}
    triggers::Dict
    system_codes::Dict
    samplingrate::typeof(1.0u"Hz")
    reference_channel::Array{AbstractString,1}
    file_path::AbstractString
    file_name::AbstractString
    processing::Dict
    header::Dict
end

#######################################
#
# Show
#
#######################################

import Base.show
function Base.show(io::IO, a::TR)
    time_length = round.(size(a.data, 1) / samplingrate(Float64, a) / 60)
    println(
        io,
        "Transient response measurement of $time_length mins with $(size(a.data,2)) channels sampled at $(samplingrate(a))",
    )

end

#######################################
#
# Read data
#
#######################################

"""
    read_TR(fname::AbstractString)
    read_TR(args...)

Read a file or IO stream and store the data in an `TR` type.

# Arguments

- `fname`: Name of the file to be read
- `min_epoch_length`: Minimum epoch length in samples. Shorter epochs will be removed (0)
- `max_epoch_length`: Maximum epoch length in samples. Longer epochs will be removed (0 = all)
- `valid_triggers`: Triggers that are considered valid, others are removed ([1,2])
- `stimulation_amplitude`: Amplitude of stimulation (NaN)
- `remove_first`: Number of epochs to be removed from start of recording (0)
- `max_epochs`: Maximum number of epochs to retain (0 = all)

# Supported file formats

- BIOSEMI (.bdf)
"""
function read_TR(
    fname::AbstractString;
    valid_triggers::Array{Int} = [1, 2],
    min_epoch_length::Int = 0,
    max_epoch_length::Int = 0,
    remove_first::Int = 0,
    max_epochs::Int = 0,
    kwargs...,
)

    @info("Importing EEG from file: $fname")
    file_path, file_name, ext = fileparts(fname)

    #
    # Read file data
    #

    # Import raw data
    if ext == "bdf"
        data, triggers, system_codes, samplingrate, reference_channel, header =
            import_biosemi(fname; kwargs...)
    else
        warn("File type $ext is unknown")
    end

    # Create electrodes
    elecs = Electrode[]
    for e in header["chanLabels"]
        push!(elecs, Electrode(e, Talairach(NaN * u"m", NaN * u"m", NaN * u"m"), Dict()))
    end

    # Create EEG type
    if unit(samplingrate) == unit(1.0u"Hz")
        #nothing
    else
        samplingrate = samplingrate * 1.0u"Hz"
    end

    a = TR(
        data,
        elecs,
        triggers,
        system_codes,
        samplingrate,
        [reference_channel],
        file_path,
        file_name,
        Dict(),
        header,
    )

    #
    # Clean up
    #

    # Remove status channel information
    remove_channel!(a, "Status")

    # Clean epoch index
    a.triggers = clean_triggers(
        a.triggers,
        valid_triggers,
        min_epoch_length,
        max_epoch_length,
        remove_first,
        max_epochs,
    )

    # Try and match sensor names to known locations
    locs = read_elp(joinpath(datadep"BioSemi64Locations", "biosemi64.elp"))
    new_sens, idx = match_sensors(locs, labels(sensors(a)))
    if length(new_sens) == length(sensors(a))
        @debug("Sucsessfully matches location of all sensors, using new locations")
        a.sensors = new_sens
    end

    return a
end


text:
# This file was generated, do not modify it. # hide
using Statistics #hide
@show trace.σ #hide
@show trace.β #hide
@show mean(y) #hide
@show std(y)  #hide

max_stars_repo_name: cscherrer/InformativePrior

max_issues_repo_name: cscherrer/InformativePrior

max_forks_repo_name: cscherrer/InformativePrior


text:
include("bilog.jl")
include("colorize_field.jl")
using ImageView
using Parameters

function imshowall(images, beam_energy_map, transmit_time_map, sim_params)
    @unpack dbrange, orientation = sim_params

    vertical_enabled = orientation == :vertical
    if vertical_enabled
        images = mapslices(rotr90, images; dims=[1, 2])
        beam_energy_map = rotr90(beam_energy_map)
        transmit_time_map = rotr90(transmit_time_map)
    end

    ImageView.imshow(colorize_field(bilog(images, dbrange)))
    ImageView.imshow(bilog(beam_energy_map, dbrange))
    ImageView.imshow(transmit_time_map)

    return  # nothing
end


text:
using CSV, CategoricalArrays, StatsPlots; pyplot()

df = CSV.read("../data/companyData.csv")
years = levels(df.Year)
data  = reshape(df.MarketCap, 5, 3)

p1 = groupedbar(years, data, bar_position=:stack)
p2 = groupedbar(years, data, bar_position=:dodge)
plot(p1, p2, bar_width=0.7, fill=[:blue :red :green], label=["A" "B" "C"],
	ylims=(0,6), xlabel="Year", ylabel="Market Cap (MM)",
	legend=:topleft, size=(800,400))

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
# using Revise, FLPDevelopment, BrowseTables

# Directory = "/home/beatriz/mainen.flipping.5ht@gmail.com/Flipping/run_task_2/Test"
# Directory = "/home/beatriz/mainen.flipping.5ht@gmail.com/Flipping/run_task_2/Dev_raw_data"
pathMac = "/Volumes/GoogleDrive/My Drive/Flipping/Datasets/Pharmacology/BeatrizDevelopment"
pathLin = "/home/beatriz/mainen.flipping.5ht@gmail.com/Flipping/Datasets/Pharmacology/BeatrizDevelopment"
if ispath(pathMac)
    path = pathMac
elseif ispath(pathLin)
    path = pathLin
else
    error("unknown path")
end
##
Age_d, Age_p, Age_b, Age_s = process_dataset(path)
######### Filtering Datas #############
for df in (Age_p, Age_b, Age_s)
    df[!,:Exp_type] = [x in age_exp ? "Development" : "Projections" for x in df.MouseID]
    filter!(r->r.Exp_type == "Development",df)
    df[!,:Age] = [x in youngs ? "Juveniles" : "Adults" for x in df.MouseID]
    transform!(df, :Age => categorical, renamecols=false)
    gd = groupby(df,:Session)
    transform!(gd, :Streak => maximum => :Performance)
end

levels!(Age_p.Age,["Adults", "Juveniles"])
levels!(Age_b.Age,["Adults", "Juveniles"])
levels!(Age_s.Age,["Adults", "Juveniles"])

Afreq = countmap(Age_s.AfterLast)
Aprob = Dict()
for (a,f) in Afreq
    Aprob[a] = round(f/nrow(Age_s),digits = 5)
end
Age_s[!,:P_AfterLast] = [Aprob[a] for a in Age_s.AfterLast]
gd = groupby(Age_s,:Session)
transform!(gd, :AfterLast => frequency)

# for (df,name) in zip([Age_p, Age_b, Age_s],["Pokes", "Bouts", "Streaks"])
#     filename = joinpath(path, "Results_" * string(today()),"FullInfo" * name * ".csv")
#     CSV.write(filename,df)
# end


text:
function rediscretizet!(system::LegSystem,n::Int64)
    # time remaining
    trem = system.solverparams.th*(system.solverparams.numbeats+1) - system.t[n+1];
    # number of steps remaining
    nrem = ceil(trem/system.solverparams.h);
    # update numsteps
    system.solverparams.numsteps = n+nrem;
    # update discrete times
    ttoadd = [0:1:nrem;]*system.solverparams.h + system.t[n+1];
    resize!(system.t,n+1)
    append!(system.t,ttoadd[2:end]);
    # resize allocated space for solution variables
    for i = 1:length(system.branches.ID)
        system.branches.A[i] = system.branches.A[i][:,1:n+1];
        system.branches.A[i] = [system.branches.A[i] zeros(system.solverparams.JL,nrem)];
        system.branches.Q[i] = system.branches.Q[i][:,1:n+1];
        system.branches.Q[i] = [system.branches.Q[i] zeros(system.solverparams.JL,nrem)];
        system.branches.P[i] = system.branches.P[i][:,1:n+1];
        system.branches.P[i] = [system.branches.P[i] zeros(system.solverparams.JL,nrem)];
        if isempty(system.branches.children[i])
            system.branches.term[i].P = system.branches.term[i].P[1:n+1];
            system.branches.term[i].P = [system.branches.term[i].P;zeros(length(ttoadd[2:end]),1)];
            system.branches.term[i].Q = system.branches.term[i].Q[1:n+1];
            system.branches.term[i].Q = [system.branches.term[i].Q;zeros(length(ttoadd[2:end]),1)];
        end
    end
    return system
end

max_stars_repo_name: dcanuto/leg-hemodynamics

max_issues_repo_name: dcanuto/leg-hemodynamics

max_forks_repo_name: dcanuto/leg-hemodynamics


text:
using Documenter, Word2Vec

makedocs()