text:
add_verbose_scope(:AbExt)
add_assert_scope(:AbExt)

add_verbose_scope(:MaxAbExt)

##############################################################################
#
#  Sieves for primes and squarefree numbers
#
##############################################################################

function squarefree_up_to(n::Int; coprime_to::Array{fmpz,1}=fmpz[])

  list = trues(n)
  for x in coprime_to
    t = Int(x)
    while t <= n
      list[t]=false
      t += Int(x)
    end
  end
  i = 2
  b = root(n, 2)
  lp = primes_up_to(b)
  for i = 1:length(lp)
    p2 = lp[i]*lp[i]
    ind = p2
    while ind <= n
      list[ind] = false
      ind += p2
    end
  end
  return findall(list)

end

function SqfSet(n::fmpz; coprime_to::Int=-1)
  return (i for i=-n:n if i!= -1 && i!=0 && issquarefree(i))
end

function primes_up_to(n::Int)

  list= trues(n)
  list[1]=false
  i=2
  s=4
  while s<=n
    list[s]=false
    s+=2
  end
  i=3
  b=sqrt(n)
  while i<=b
    if list[i]
      j=3*i
      s=2*i
      while j<= n
        list[j]=false
        j+=s
      end
    end
    i+=1
  end
  return findall(list)

end

###########################################################################
#
#  Some useful functions
#
###########################################################################

#This function gets a quotient of the ray class group and the action on
# the ray class group
# In output, we get the quotient group as a ZpnGModule

function _action_on_quo(mq::GrpAbFinGenMap, act::Array{GrpAbFinGenMap,1})

  q=mq.header.codomain
  S,mS=snf(q)
  n=Int(S.snf[end])
  R=ResidueField(FlintZZ, n, cached=false)
  quo_action=Array{nmod_mat,1}(undef, length(act))
  for s=1:length(act)
    quo_action[s]= change_base_ring(mS.map*act[i].map*mS.imap, R)
  end
  return ZpnGModule(S, quo_action)

end

# Function that finds an integer in the ideal
function _min_wild(D::Dict{NfOrdIdl, Int})

  res = fmpz(1)
  primes_done=fmpz[]
  for (p,v) in D
    s=minimum(p)
    if s in primes_done
      continue
    end
    push!(primes_done, s)
    res*=s^v
  end
  return res

end

function totally_positive_generators(mr::MapRayClassGrp, wild::Bool=false)

  a = minimum(mr.modulus_fin)
  if isdefined(mr, :tame_mult_grp)
    tmg=mr.tame_mult_grp
    for (p, v) in tmg
      new_p = GrpAbFinGenToNfAbsOrdMap(domain(v), codomain(v), [ make_positive(v.generators[1], a) ], v.discrete_logarithm)
      if isdefined(v, :disc_log)
        new_p.disc_log = v.disc_log
      end
      tmg[p] = new_p
      @hassert :RayFacElem 1 istotally_positive(mr.tame_mult_grp[p].generators[1])
    end
  end
  if wild && isdefined(mr, :wild_mult_grp)
    wld=mr.wild_mult_grp
    for (p,v) in wld
      x=v.generators
      for i=1:length(x)
        x[i]=make_positive(x[i],a)
        @hassert :RayFacElem 1 iscoprime(ideal(parent(x[i]), x[i]), ideal(parent(x[i]), a))
      end
      mr.wild_mult_grp[p] = GrpAbFinGenToNfAbsOrdMap(domain(v), codomain(v), x, v.discrete_logarithm)
    end
  end

end

function make_positive(x::NfOrdElem, a::fmpz)

  els=conjugates_real(elem_in_nf(x))
  m=fmpz(0)
  for i=1:length(els)
    y = BigFloat(midpoint(els[i]/a))
    if y > 0
      continue
    else
      m = max(m,1-ceil(fmpz,y))
    end
  end
  @hassert :RayFacElem 1 iscoprime(ideal(parent(x),x), ideal(parent(x), a))
  @hassert :RayFacElem 1 iscoprime(ideal(parent(x),x+fmpz(m)*a), ideal(parent(x), a))
  @hassert :RayFacElem 1 istotally_positive(x+m*a)
  return x+fmpz(m)*a

end

function _are_there_subs(G::GrpAbFinGen,gtype::Array{Int,1})

  H=DiagonalGroup(gtype)
  H, _=snf(H)
  G1, _=snf(G)
  if length(G1.snf)<length(H.snf)
    return false
  end
  for i=0:length(H.snf)-1
    if !divisible(G1.snf[end-i],H.snf[end-i])
      return false
    end
  end
  return true
end

function issquarefree(n::Union{Int,fmpz})
  if iszero(n)
    throw(error("Argument must be non-zero"))
  end
  return isone(abs(n)) || maximum(values(factor(n).fac)) == 1
end

##########################################################################################################
#
#  Functions for conductors
#
##########################################################################################################

function tame_conductors_degree_2(O::NfOrd, bound::fmpz)

  K = nf(O)
  d = degree(O)
  b1 = Int(root(bound,d))
  ram_primes = ramified_primes(O)
  sort!(ram_primes)
  filter!(x -> x!=2, ram_primes)
  list = squarefree_up_to(b1, coprime_to = vcat(ram_primes,2))

  extra_list = Tuple{Int, fmpz}[(1,fmpz(1))]
  for q in ram_primes
    tr = prime_decomposition_type(O, Int(q))
    e = tr[1][2]
    nq = fmpz(q)^(divexact(d,e))
    if nq > bound
      break
    end
    l=length(extra_list)
    for i = 1:l
      n = extra_list[i][2]*nq
      if n > bound
        continue
      end
      push!(extra_list, (extra_list[i][1]*q, n))
    end
  end

  final_list=Tuple{Int,fmpz}[]
  l = length(list)
  for (el,norm) in extra_list
    for i=1:l
      if fmpz(list[i])^d*norm>bound
        continue
      end
      push!(final_list, (list[i]*el, fmpz(list[i])^d*norm))
    end
  end
  return final_list

end

function squarefree_for_conductors(O::NfOrd, n::Int, deg::Int; coprime_to::Array{fmpz,1}=fmpz[])

  sqf = trues(n)
  primes = trues(n)

  #remove primes that can be wildly ramified or
  #that are ramified in the base field
  for x in coprime_to
    el = Int(x)
    t = el
    while t <= n
      @inbounds sqf[t]=false
      @inbounds primes[t]=false
      t += el
    end
  end

  #sieving procedure

  if !(2 in coprime_to)
    dt = prime_decomposition_type(O,2)
    if isone(gcd(2^dt[1][1]-1, deg))
      j=2
      while j<=n
        @inbounds sqf[j]=false
        @inbounds primes[j]=false
        j += 2
      end
    else
      i=2
      s=4
      while s <= n
        @inbounds primes[s]=false
        s+=2
      end
      s=4
      while s<=n
        @inbounds sqf[s]=false
        s+=4
      end
    end
  end
  i = 3
  b = root(n, 2)
  while i <= b
    if primes[i]
      if gcd(i-1, deg) != 1
        j = i
        while j <= n
          @inbounds primes[j]=false
          j += i
        end
        j = i^2
        t = j
        while j <= n
          @inbounds sqf[j]=false
          j += t
        end
      else
        dt = prime_decomposition_type(O, i)
        if gcd(deg, i^dt[1][1]-1) == 1
          @inbounds primes[i] = false
          @inbounds sqf[i] = false
          j = i
          while j <= n
           @inbounds primes[j] = false
           @inbounds sqf[j] = false
           j+=i
          end
        else
          j=i
          while j <= n
            @inbounds primes[j]=false
            j+=i
          end
          j=i^2
          t=j
          while j<= n
            @inbounds sqf[j]=false
            j+=t
          end
        end
      end
    end
    i+=2
  end
  while i<=n
    if primes[i] && gcd(i-1, deg) == 1
      dt=prime_decomposition_type(O,i)
      if gcd(deg,i^dt[1][1]-1)==1
        @inbounds sqf[i]=false
        j=i
        while j<= n
         @inbounds sqf[j]=false
         j+=i
        end
      end
    end
    i+=2
  end

  if degree(O)==1
    i=2
    while i<=length(sqf)
      @inbounds sqf[i]=false
      i+=4
    end
  end
  return Int[i for i=1:length(sqf) if sqf[i]]

end

function conductors_tame(O::NfOrd, n::Int, bound::fmpz)

  if n == 2
    return tame_conductors_degree_2(O,bound)
  end
  #
  #  First, conductors coprime to the ramified primes and to the
  #  degree of the extension we are searching for.
  #
  d = degree(O)
  K = nf(O)
  wild_ram = collect(keys(factor(fmpz(n)).fac))
  ram_primes = ramified_primes(O)
  filter!(x -> !divisible(fmpz(n),x), ram_primes)
  coprime_to = vcat(ram_primes, wild_ram)
  sort!(ram_primes)
  m = minimum(wild_ram)
  k = divexact(n, m)
  e = Int((m-1)*k)
  b1 = Int(root(bound, degree(O)*e))
  list = squarefree_for_conductors(O, b1, n, coprime_to = coprime_to)

  extra_list = Tuple{Int, fmpz}[(1, fmpz(1))]
  for q in ram_primes
    tr = prime_decomposition_type(O, Int(q))
    f = tr[1][1]
    nq = q^f
    if gcd(nq-1,n) == 1
      continue
    end
    nq = nq^(length(tr)*e)
    if nq > bound
      continue
    end
    l = length(extra_list)
    for i=1:l
      no = extra_list[i][2]*nq
      if no > bound
        continue
      end
      push!(extra_list, (extra_list[i][1]*q, no))
    end
  end

  final_list=Tuple{Int, fmpz}[]
  l=length(list)
  e = Int((m-1)*k)
  for (el,norm) in extra_list
    for i=1:l
      if (list[i]^(e*d)) * norm > bound
        continue
      end
      push!(final_list, (list[i]*el, (fmpz(list[i])^(e*d))*norm))
    end
  end

  return final_list
end

function conductors(O::NfOrd, a::Array{Int, 1}, bound::fmpz, tame::Bool=false)

  #Careful: I am assuming that a is in snf!
  K=nf(O)
  d=degree(O)
  n = prod(a)
  expo = a[end]
  wild_ram=collect(keys(factor(fmpz(n)).fac))

  #
  # First, conductors for tamely ramified extensions
  #
  bound_tame = root(bound, divexact(n, expo))
  list = conductors_tame(O, expo, bound_tame)

  if tame
    return Tuple{Int, Dict{NfOrdIdl, Int}}[(x[1], Dict{NfOrdIdl, Int}()) for x in list]
  end
  #
  # now, we have to multiply the obtained conductors by proper powers of wildly ramified ideals.
  #
  wild_list=Tuple{Int, Dict{NfOrdIdl, Int}, fmpz}[(1, Dict{NfOrdIdl, Int}(), fmpz(1))]
  for q in wild_ram
    lp = prime_decomposition(O, Int(q))
    fq = divexact(d, lp[1][2]*length(lp))
    l = length(wild_list)
    sq = fmpz(q)^(divexact(d,lp[1][2])) #norm of the squarefree part of the integer q
    #=
      we have to use the conductor discriminant formula to understand the maximal possible exponent of q.
      Let ap be the exponent of p in the relative discriminant, let m be the conductor and h_(m,C) the cardinality of the
      quotient of ray class group by its subgroup C.
      Then
          ap= v_p(m)h_(m,C)- sum_{i=1:v_p(m)} h_(m/p^i, C)
      Since m is the conductor, h_(m/p^i, C)<= h_(m,C)/q.
      Consequently, we get
        v_p(m)<= (q*ap)/(h_(m,C)*(q-1))
      To find ap, it is enough to compute a logarithm.
    =#
    v = valuation(expo, q)
    # First, we compute the bound coming from the bound on the discriminant
    boundsubext = root(bound, Int(divexact(n, q^v))) #The bound on the norm of the discriminant on the subextension
                                                     # of order q^v
    #Bound coming from the bound on the discriminant
    obound = flog(boundsubext, sq)

    #Bound coming from the analysis on the different in a local extension
    nbound = q^v + lp[1][2] * v * q^v - 1

    bound_max_ap = min(nbound, obound)  #bound on ap
    bound_max_exp = div(bound_max_ap, (q-1)*q^(v-1)) #bound on the exponent in the conductor

    #Ramification groups bound
    max_nontrivial_ramification_group = div(lp[1][2]*(q^v), q-1)
    if v > 1
      ram_groups_bound = max_nontrivial_ramification_group - sum(q^i for i = 1:v-1) + v
    else
      ram_groups_bound = max_nontrivial_ramification_group + 1
    end
    bound_max_exp = min(ram_groups_bound, bound_max_exp)

    #The prime may be also tamely ramified!
    nisc = gcd(q^(fq)-1, expo)
    if nisc != 1
      fnisc=minimum(keys(factor(nisc).fac))
      nq=sq^((fnisc-1)*(divexact(n, fnisc)))
      for s=1:l
        nn=nq*wild_list[s][3]
        if nn>bound
          continue
        end
        push!(wild_list, (q*wild_list[s][1], wild_list[s][2], nn))
      end
    end
    for i=2:bound_max_exp
      d1=Dict{NfOrdIdl, Int}()
      for j=1:length(lp)
        d1[lp[j][1]]=i
      end
      exp1 = i*(q-1)*divexact(n,q)
      nq= sq^(exp1)
      for s=1:l
        nn=nq*wild_list[s][3]
        if nn>bound
          continue
        end
        d2 = merge(max, d1, wild_list[s][2])
        if nisc!=1
          push!(wild_list, (q*wild_list[s][1], d2, nn))
        else
          push!(wild_list, (wild_list[s][1], d2, nn))
        end
      end
    end
  end

  #the final list
  final_list=Tuple{Int, Dict{NfOrdIdl, Int}}[]
  for (el, nm) in list
    for (q,d,nm2) in wild_list
      if nm*nm2 > bound
        continue
      end
      push!(final_list, (el*q,d))
    end
  end
  return final_list

end

###############################################################################
#
#  Conductors over QQ
#
###############################################################################

function squarefree_for_conductorsQQ(O::NfOrd, n::Int, a::Array{Int, 1}; coprime_to::Array{fmpz,1}=fmpz[])

  G = map(Int, snf(DiagonalGroup(a))[1].snf)
  sqf= trues(n)
  primes= trues(n)
  deg = G[end]
  #remove primes that can be wildly ramified
  for x in coprime_to
    el=Int(x)
    t=el
    while t<= n
      @inbounds sqf[t]=false
      @inbounds primes[t]=false
      t+=el
    end
  end

  single = Array{Int, 1}()
  push!(single, 1)
  multiple = Array{Int, 1}()

  #sieving procedure
  #First, I can remove all the multiples of 2
  if !(2 in coprime_to)
    i=2
    while i<=length(sqf)
      @inbounds sqf[i]=false
      i+=2
    end
  end

  i=3
  b=Base.sqrt(n)
  while i<=b
    if primes[i]
      if gcd(deg,i-1)==1
        @inbounds primes[i]=false
        @inbounds sqf[i]=false
        j=i
        while j<= n
         @inbounds primes[j]=false
         @inbounds sqf[j]=false
         j+=i
        end
      else
        j=i
        while j<= n
          @inbounds primes[j]=false
          j+=i
        end
        j=i^2
        t=2*j
        while j<= n
          @inbounds sqf[j]=false
          j+=t
        end
      end
    end
    i+=2
  end
  while i<=n
    if primes[i]
      if gcd(deg,i-1) == 1
        @inbounds primes[i] = false
        @inbounds sqf[i]=false
        j = i
        while j <= n
         @inbounds sqf[j]=false
         j += i
        end
      end
    end
    i+=2
  end

  if length(G) > 1
    i = 3
    while i < n
      if primes[i]
        push!(single, i)
      elseif sqf[i]
        push!(multiple, i)
      end
      i += 2
    end
  elseif !isprime(deg)
    i = 3
    while i < n
      if primes[i]
        if rem(i-1, deg) == 0
          push!(multiple, i)
        else
          push!(single, i)
        end
      elseif sqf[i]
        push!(multiple, i)
      end
      i += 2
    end
  else
    multiple = Int[i for i = 2:length(sqf) if sqf[i]]
  end

  return single, multiple

end

function conductors_tameQQ(O::NfOrd, a::Array{Int, 1}, bound::fmpz)

  #
  #  First, conductors coprime to the ramified primes and to the
  #  degree of the extension we are searching for.
  #
  n = prod(a)
  wild_ram = collect(keys(factor(fmpz(n)).fac))
  m = minimum(wild_ram)
  k = divexact(n, m)
  b1 = Int(root(fmpz(bound),Int((m-1)*k)))

  return squarefree_for_conductorsQQ(O, b1, a, coprime_to=wild_ram)

end

function conductorsQQ(O::NfOrd, a::Array{Int, 1}, bound::fmpz, tame::Bool=false)

  K = nf(O)
  d = degree(O)
  n = prod(a)
  expo = a[end]
  wild_ram = collect(keys(factor(fmpz(n)).fac))

  #
  # First, conductors for tamely ramified extensions
  #

  single, multiple = conductors_tameQQ(O,a,bound)

  if tame
    return multiple
  end
  #
  # now, we have to multiply the obtained conductors by proper powers of wildly ramified ideals.
  #
  wild_list=Tuple{Int, Int, fmpz}[(1, 1, 1)]
  for q in wild_ram
    l = length(wild_list)
    #=
      we have to use the conductor discriminant formula to understand the maximal possible exponent of q.
      Let ap be the exponent of p in the relative discriminant, let m be the conductor and h_(m,C) the cardinality of the
      quotient of ray class group by its subgroup C.
      Then
          ap= v_p(m)h_(m,C)- sum_{i=1:v_p(m)} h_(m/p^i, C)
      Since m is the conductor, h_(m/p^i, C)<= h_(m,C)/q.
      Consequently, we get
        v_p(m)<= (q*ap)/(h_(m,C)*(q-1))
      To find ap, it is enough to compute a logarithm.
    =#
    v = valuation(expo, q)

    #I don't need to give a bound for a_p on the big extension but only on the maximum extension of q-power order
    #This is the only thing that matters for the exponent of the conductor
    nisc = gcd(q-1,n)
    nbound = q^v + v * q^v - 1
    boundsubext = root(bound, Int(divexact(n, q^v)))
    obound = flog(boundsubext, q)
    nnbound = valuation_bound_discriminant(n, q)
    k = div(expo, Int(q)-1)
    bound_max_ap = min(nbound, obound, nnbound)  #bound on ap
    bound_max_exp = div(bound_max_ap, (q^(v-1))*(q-1)) #bound on the exponent in the conductor
    if nisc != 1
      fnisc=minimum(keys(factor(nisc).fac))
      nq=fmpz(q)^((fnisc-1)*(divexact(n, fnisc)))
      for s=1:l
        nn=nq*wild_list[s][3]
        if nn>bound
          continue
        end
        push!(wild_list, (q*wild_list[s][1], wild_list[s][2], nn))
      end
    end
    t=(q-1)*divexact(n,q)
    for i=2:bound_max_exp
      nq= q^(i*t)
      for s=1:l
        nn=nq*wild_list[s][3]
        if nn>bound
          continue
        end
        push!(wild_list, (wild_list[s][1], wild_list[s][2]*q^i, nn))
      end
    end
  end

  #the final list
  final_list = Int[]
  exps = Int((minimum(wild_ram)-1)*divexact(n, minimum(wild_ram)))
  for el in multiple
    for (q,d,nm2) in wild_list
      if (fmpz(el)^exps)*nm2 > bound
        continue
      end
      push!(final_list, (el*q*d))
    end
  end

  for el in single
    for j = 2:length(wild_list)
      q,d,nm2 = wild_list[j]
      if (fmpz(el)^exps)*nm2 > bound
        continue
      end
      push!(final_list, (el*q*d))
    end
  end
  return final_list

end

###############################################################################
#
#  Abelian extensions
#
###############################################################################

function abelian_extensions(O::Union{FlintIntegerRing, FlintRationalField},
                            gtype::Vector{Int}, discriminant_bound::fmpz;
                            real::Bool = false,
                            tame::Bool = false)

  Qx, x = PolynomialRing(FlintQQ, "x", cached = false)
  K, _ = NumberField(x - 1, "a", cached = false)
  OK = maximal_order(K)
  l = abelian_extensions(OK, gtype, discriminant_bound,
                         real = real,
                         tame = tame)

  newlist = Vector{NfAbsNS}(undef, length(l))
  for j in 1:length(l)
    newlist[j], _ = number_field([Qx([coeff(coeff(f, i), 0) for i in 0:length(f)]) for f in l[j].abs_pol], check = false)
  end
  return newlist
end

function abelian_extensions(O::NfOrd, gtype::Array{Int,1}, absolute_discriminant_bound::fmpz; real::Bool=false, tame::Bool=false, with_autos::Type{Val{T}} = Val{false}) where T

  K=nf(O)
  @assert degree(K)==1
  gtype = map(Int, snf(DiagonalGroup(gtype))[1].snf)
  n = prod(gtype)

  expo = lcm(gtype)

  #Getting conductors
  l_conductors=conductorsQQ(O, gtype, absolute_discriminant_bound, tame)
  @vprint :AbExt 1 "Number of conductors: $(length(l_conductors)) \n"
  if with_autos==Val{true}
    fields = Tuple{NfRel_ns, Array{NfRel_nsToNfRel_nsMor,1}}[]
  else
    fields = NfRel_ns[]
  end

  #Now, the big loop
  for (i, k) in enumerate(l_conductors)
    @vprint :AbExt 1 "Conductor: $k \n"
    @vprint :AbExt 1 "Left: $(length(l_conductors) - i)\n"
    r,mr=Hecke.ray_class_groupQQ(O,k,!real,expo)
    if !Hecke._are_there_subs(r,gtype)
      continue
    end
    ls=subgroups(r,quotype=gtype, fun= (x, y) -> quo(x, y, false)[2])
    for s in ls
      C=ray_class_field(mr, s)
      if Hecke._is_conductor_minQQ(C,n) && Hecke.discriminant_conductorQQ(O,C,k,absolute_discriminant_bound,n)
        @vprint :AbExt 1 "New Field \n"
        L=number_field(C)
        if with_autos==Val{true}
          push!(fields,(L, automorphism_groupQQ(C)))
        else
          push!(fields, L)
        end
      end
    end
  end

  return fields

end

function abelian_normal_extensions(O::NfOrd, gtype::Array{Int,1}, absolute_discriminant_bound::fmpz; ramified_at_infplc::Bool=true, tame::Bool=false, absolute_galois_group::Symbol = :all,  with_autos::Type{Val{T}} = Val{false}, autos::Array{NfToNfMor,1}=NfToNfMor[]) where T

  d=degree(O)
  if d==1
    return abelian_extensions(O, gtype, absolute_discriminant_bound, real=!ramified_at_infplc, tame=tame, with_autos=with_autos)
  end

  K=nf(O)
  n=prod(gtype)
  inf_plc=InfPlc[]

  if with_autos==Val{true}
    fields=Tuple{NfRel_ns,Vector{NfRel_nsToNfRel_nsMor{nf_elem}}}[]
  else
    fields=NfRel_ns[]
  end

  bound = div(absolute_discriminant_bound, abs(discriminant(O))^n)
  if bound == 0
    return fields
  end

  if mod(n,2)==0 && ramified_at_infplc
    inf_plc=real_places(K)
  end

  expo=lcm(gtype)
  C,mC=class_group(O)
  cgrp= gcd(n,order(C))!=1
  if cgrp
    S = prime_ideals_up_to(O, max(100,12*clog(abs(discriminant(O)),10)^2))
  end
  allow_cache!(mC)

  #
  # Getting a small set of generators
  # for the automorphisms group
  #
  if isempty(autos)
    Aut=Hecke.automorphisms(K)
    @assert length(Aut)==degree(O) #The field is normal over Q
    gens = Hecke.small_generating_set(Aut)
  else
    gens = autos
  end

  #Getting conductors
  l_conductors=conductors(O, gtype, bound, tame)
  @vprint :AbExt 1 "Number of conductors: $(length(l_conductors)) \n"

  #Now, the big loop
  for (i, k) in enumerate(l_conductors)
    @vprint :AbExt 1 "Conductor: $k \n"
    @vprint :AbExt 1 "Left: $(length(l_conductors) - i)\n"
    r,mr=ray_class_group_quo(O,expo,k[1], k[2],inf_plc)
    if !_are_there_subs(r,gtype)
      continue
    end
    if cgrp
        mr.prime_ideal_cache = S
    end
    act = induce_action(mr,gens)
    ls = stable_subgroups(r, act, op=(x, y) -> quo(x, y, false)[2], quotype = gtype)
    totally_positive_generators(mr)
    for s in ls
      @hassert :AbExt 1 order(codomain(s))==n
      C = ray_class_field(mr, s)
      if Hecke._is_conductor_min_normal(C) && Hecke.discriminant_conductor(C, bound)
        @vprint :AbExt 1 "New Field \n"
        L = number_field(C)
        if absolute_galois_group != :all
          autabs = absolute_automorphism_group(C, gens)
          G = generic_group(autabs, *)
          if absolute_galois_group == :_16T7
            if isisomorphic_to_16T7(G)
              @vprint :AbExt 1 "I found a field with Galois group 16T7 (C_2 x Q_8)\n"
              @vtime :AbExt 1 push!(fields, L)
            end
          elseif absolute_galois_group == :_8T5
            if isisomorphic_to_8T5(G)
              @vprint :AbExt 1 "I found a field with Galois group 8T5 (Q_8)\n"
              @vtime :AbExt 1 push!(fields, L)
            end
          else
            error("Group not supported (yet)")
          end
        else
          if with_autos==Val{true}
            push!(fields, (L, absolute_automorphism_group(C,gens)))
          else
            push!(fields, L)
          end
        end
      end
    end
  end

  return fields

end

###############################################################################
#
#  Quadratic Extensions of Q
#
###############################################################################

function quadratic_extensions(bound::Int; tame::Bool=false, real::Bool=false, complex::Bool=false, with_autos::Type{Val{T}}=Val{false}) where T

  @assert !(real && complex)
  Qx,x=PolynomialRing(FlintQQ, "x")
  sqf=squarefree_up_to(bound)
  if real
    deleteat!(sqf,1)
  elseif complex
    sqf=Int[-i for i in sqf]
  else
    sqf= vcat(sqf[2:end], Int[-i for i in sqf])
  end
  if tame
    filter!( x -> mod(x,4)==1, sqf)
    return ( number_field(x^2-x+divexact(1-i,4), cached=false, check = false)[1] for i in sqf)
  end
  final_list=Int[]
  for i=1:length(sqf)
    if abs(sqf[i]*4)< bound
      @views push!(final_list,sqf[i])
      continue
    end
    if mod(sqf[i],4)==1
      @views push!(final_list,sqf[i])
    end
  end
  return ( mod(i,4)!=1 ? number_field(x^2-i, cached=false, check = false)[1] : number_field(x^2-x+divexact(1-i,4), cached = false, check = false)[1] for i in final_list)

end

function _quad_ext(bound::Int, only_real::Bool = false)

  Qx, x = PolynomialRing(FlintQQ, "x", cached = false)
  K = NumberField(x-1, cached = false, check = false)[1]
  sqf = squarefree_up_to(bound)
  final_list = Int[]
  for i=2:length(sqf)
    if abs(sqf[i]*4)< bound
      @views push!(final_list, sqf[i])
      continue
    end
    if mod(sqf[i], 4) == 1
      @views push!(final_list, sqf[i])
    end
  end
  if !only_real
    for i = 1:length(sqf)
      if abs(sqf[i]*4)< bound
        @views push!(final_list, -sqf[i])
        continue
      end
      if mod(sqf[i], 4) == 3
        @views push!(final_list, -sqf[i])
      end
    end
  end
  fields_list = Array{Tuple{AnticNumberField, Vector{NfToNfMor}, Vector{NfToNfMor}}, 1}(undef, length(final_list))
  for i = 1:length(final_list)
    if mod(final_list[i],4) != 1
      L, gL = number_field(x^2-final_list[i], cached=false, check = false)
      auts = NfToNfMor[hom(L, L, -gL, check = false)]
      emb = NfToNfMor[hom(K, L, one(L), check = false)]
      fields_list[i] = (L, auts, emb)
    else
      L, gL = number_field(x^2-x+divexact(1-final_list[i], 4), cached=false, check = false)
      auts = NfToNfMor[hom(L, L, 1-gL, check = false)]
      emb = NfToNfMor[hom(K, L, one(L), check = false)]
      fields_list[i] = (L, auts, emb)
    end
  end
  return fields_list

end

###############################################################################
#
#  C2 x C2 extensions of Q
#
###############################################################################

function C22_extensions(bound::Int)

  Qx, x=PolynomialRing(FlintZZ, "x")
  K, _=NumberField(x-1, cached = false)
  Kx,x=PolynomialRing(K,"x", cached=false)
  b1=ceil(Int,Base.sqrt(bound))
  n=2*b1+1
  pairs = _find_pairs(bound)
  return (_ext(Kx,x,i,j) for (i, j) in pairs)

end

function _ext(Ox,x,i,j)

  y=mod(i,4)
  pol1 = x^2
  if y != 1
    pol1 = x^2-i
  else
    pol1 = x^2-x+divexact(1-i,4)
  end
  y=mod(j,4)
  pol2=Ox(1)
  if y!=1
    pol2=x^2-j
  else
    pol2=x^2-x+divexact(1-j,4)
  end
  return number_field([pol1,pol2], cached = false, check = false)

end

function _C22_exts_abexts(bound::Int, only_real::Bool = false)
  Qx, x = PolynomialRing(FlintZZ, "x")
  pairs = _find_pairs(bound, only_real)
  return (_ext_with_autos(Qx, x, i, j) for (i, j) in pairs)
end

function _ext_with_autos(Qx, x, i::Int, j::Int)
  first = i
  second = j
  g = gcd(i, j)
  if g != 1
    third = divexact(i*j, g^2)
    if gcd(first, third) == 1
      second = third
    elseif gcd(second, third) == 1
      first = third
    end
  end
  y1 = mod(first, 4)
  cp1 = Vector{fmpz}(undef, 3)
  cp1[3] = fmpz(1)
  if y1 != 1
    cp1[1] = fmpz(-first)
    cp1[2] = fmpz(0)
  else
    cp1[1] = fmpz(divexact(1-first,4))
    cp1[2] = fmpz(-1)
  end
  y2 = mod(second, 4)
  cp2 = Vector{fmpz}(undef, 3)
  cp2[3] = fmpz(1)
  if y2 != 1
    cp2[1] = fmpz(-second)
    cp2[2] = fmpz(0)
  else
    cp2[1] = fmpz(divexact(1-second, 4))
    cp2[2] = fmpz(-1)
  end
  return Qx(cp1), Qx(cp2)
end

function __get_term(a::fmpq_mpoly, exps::Vector{UInt})
   z = fmpq()
   ccall((:fmpq_mpoly_get_coeff_fmpq_ui, :libflint), Nothing,
         (Ref{fmpq}, Ref{fmpq_mpoly}, Ptr{UInt}, Ref{FmpqMPolyRing}),
         z, a, exps, parent(a))
   return z
end

function _C22_with_max_ord(l)
  list = Vector{Tuple{AnticNumberField, Vector{NfToNfMor}, Vector{NfToNfMor}}}()
  Qx, x = PolynomialRing(FlintQQ, "x", cached = false)
  K = NumberField(x-1, cached = false)[1]
  for (p1, p2) in l
    Kns, g = number_field(fmpz_poly[p1, p2], check = false)
    S, mS = simple_extension(Kns, check = false)
    d1 = discriminant(p1)
    d2 = discriminant(p2)
    cf = gcd(d1, d2)
    B = Vector{nf_elem}(undef, degree(S))
    B[1] = S(1)
    B[2] = mS\(g[1])
    B[3] = mS\(g[2])
    B[4] = B[2] * B[3]
    M = basis_matrix(B, FakeFmpqMat)
    hnf_modular_eldiv!(M.num, M.den, :lowerleft)
    O = NfAbsOrd(S, FakeFmpqMat(M.num, M.den))
    O.disc = d1^2*d2^2
    d3 = numerator(discriminant(S))
    if d3 < 0
      O.disc = -O.disc
    end
    O.index = divexact(d3, O.disc)
    if cf != 1
      fac = factor(cf)
      for (p, v) in fac
        O = pmaximal_overorder(O, p)
      end
    end
    O.ismaximal = 1
    Hecke._set_maximal_order_of_nf(S, O)
    coord1 = __get_term(mS.prim_img.data, UInt[1, 0])
    coord2 = __get_term(mS.prim_img.data, UInt[0, 1])
    auts = Vector{NfToNfMor}(undef, 2)
    if iszero(coeff(p1, 1))
      auts[1] = hom(S, S, (-coord1)*B[2]+coord2*B[3], check = false)
    else
      auts[1] = hom(S, S, 1+(-coord1)*B[2]+coord2*B[3], check = false)
    end
    if iszero(coeff(p2, 1))
      auts[2] = hom(S, S, coord1*B[2]+(-coord2)*B[3], check = false)
    else
      auts[2] = hom(S, S, coord1*B[2]+1+(-coord2)*B[3], check = false)
    end
    cl_auts = Vector{NfToNfMor}(undef, 4)
    cl_auts[1] = id_hom(S)
    cl_auts[2] = auts[1]
    cl_auts[3] = auts[2]
    cl_auts[4] = auts[1] * auts[2]
    Hecke._set_automorphisms_nf(S, cl_auts)
    push!(list, (S, auts, NfToNfMor[hom(K, S, S(1), check = false)]))
  end
  return list
end

function _disc(a::Int, b::Int, c::Int, bound::Int)
  a1 = mod(a, 4)
  b1 = mod(b, 4)
  if a1 == 1 && b1 == 1
    return a*b*c <= bound
  end
  c1 = mod(c, 4)
  if a1 == 1 || b1 == 1 || c1 == 1
    return 16*a*b*c <= bound
  else
    return 64*a*b*c <= bound
  end
end

function _pairs_totally_real(pairs, ls, bound)
  b1=floor(Int, Base.sqrt(bound))
  pairs[1, 1:length(ls)] .= false
  pairs[1:length(ls), 1] .= false
  for j = 2:length(ls)
    for i = j:length(ls)
      pairs[i, j] = false
    end
  end
  for j = 2:length(ls)
    second = ls[j]
    for i = 2:j-1
      if pairs[i, j]
        first = ls[i]
        g = gcd(first, second)
        if isone(g)
          third = first*second
        else
          third = divexact(first*second, g^2)
        end
        if abs(third) > b1
          pairs[i, j] = false
          continue
        end
        k = searchsortedfirst(ls, third)
        if i < k
          pairs[i, k] = false
        else
          pairs[k, i] = false
        end
        if j < k
          pairs[j, k] = false
        else
          pairs[k, j] = false
        end
        if !_disc(first, second, third, bound)
          pairs[i, j] = false
        end
      end
    end
  end
  it = findall(pairs)
  totally_real_exts = Vector{Tuple{Int, Int}}(undef, length(it))
  ind = 1
  for I in it
    totally_real_exts[ind] = (ls[I[1]], ls[I[2]])
    ind += 1
  end
  return totally_real_exts

end

function _find_pairs(bound::Int, only_real::Bool = false)

  #first, we need the squarefree numbers
  b1=ceil(Int, Base.sqrt(bound))
  ls = squarefree_up_to(b1)
  #The first step is to enumerate all the totally real extensions.
  pairs = trues(length(ls), length(ls))
  real_exts = _pairs_totally_real(pairs, ls, bound)
  if only_real
    return real_exts
  end
  ls1 = -ls
  #Now, all the others.
  pairs[1:length(ls), 2:length(ls)] .= true
  for j = 2:length(ls)
    second = ls[j]
    for i = 1:length(ls)
      if pairs[i, j]
        first = ls1[i]
        g = gcd(first, second)
        if isone(g)
          third = first*second
        else
          third = divexact(first*second, g^2)
        end
        abt = -third
        if abt > b1
          pairs[i, j] = false
          continue
        end
        k = searchsortedfirst(ls, abt)
        pairs[k, j] = false
        if !_disc(first, second, third, bound)
          pairs[i, j] = false
        end
      end
    end
  end
  it = findall(pairs)
  ind = 1
  res = Vector{Tuple{Int, Int}}(undef, length(it))
  for I in it
    res[ind] = (ls1[I[1]], ls[I[2]])
    ind += 1
  end
  return vcat(res, real_exts)
end

################################################################################
#
#  Valuation bounds for discriminants
#
################################################################################

function valuation_bound_discriminant(n::Int, p::Union{Integer, fmpz})
  # First compute the p-adic expansion of n
  S = Vector{typeof(p)}()
	q = typeof(p)(n)
  q, r = divrem(q, p)
  push!(S, r)
  while q >= p
    q, r = divrem(q, p)
    push!(S, r)
  end

	if !iszero(q)
		push!(S, q)
	end

	@assert sum(S[i + 1] * p^i for i in 0:length(S)-1) == n

	b = zero(typeof(p))

	for i in 0:length(S) - 1
		b = b + S[i + 1] * (i + 1) * p^i
		if !iszero(S[i + 1])
			b = b - 1
		end
	end

  return b
end

###############################################################################
#
#  From relative to absolute
#
###############################################################################

function _from_relative_to_absQQ(L::NfRel_ns{T}, auts::Array{NfRel_nsToNfRel_nsMor{T}, 1}) where T

  @vprint :AbExt 2 "Computing maximal orders of subfields\n"
  Qx, x = PolynomialRing(FlintQQ, "x")
  polys = Vector{fmpq_poly}(undef, length(L.pol))
  for i = 1:length(L.pol)
    fK = isunivariate(L.pol[i])[2]
    f = Qx([coeff(coeff(fK, j), 0) for j = 0:degree(fK)])
    polys[i] = f
  end
  NS, gNS = number_field(polys, check = false)
  gpols = gens(parent(gNS[1]))
  B, lp = maximal_order_of_components(NS)
  K, mK = simple_extension(NS, check = false)
  BKK = Array{nf_elem, 1}(undef, degree(K))
  ind = degree(polys[1])
  for i = 1:ind
    BKK[i] = mK\(B[1][i])
  end
  for jj = 2:length(polys)
    new_deg = degree(polys[jj])
    for i = 2:new_deg
      el = mK\(B[jj][i])
      for j = 1:ind
        BKK[(i-1)* ind + j] = BKK[j] * el
      end
    end
    ind *= new_deg
  end
  O1 = NfAbsOrd(BKK)
  for p in lp
    if isprime(p)
      O1 = pmaximal_overorder(O1, p)
    else
      fac = factor(p)
      for (k, v) in fac
        O1 = pmaximal_overorder(O1, k)
      end
    end
  end
  #@vtime :AbExt 2 O1 = pmaximal_overorder_at(O1, lp)

  _set_maximal_order(K, O1)
  #Now, we translate the automorphisms.
  imgs = Vector{NfAbsNSElem}(undef, length(auts))
  for i = 1:length(auts)
    fK = isunivariate(auts[i].emb[i].data)[2]
    f = Qx([coeff(coeff(fK, j), 0) for j = 0:degree(fK)])
    imgs[i] = NS(evaluate(f, gpols[i]))
  end
  autsNS = Vector{NfAbsNSToNfAbsNS}(undef, length(auts))
  for t = 1:length(auts)
    imgs = Vector{NfAbsNSElem}(undef, length(polys))
    for s = 1:length(polys)
      fK = isunivariate(auts[t].emb[s].data)[2]
      f = Qx([coeff(coeff(fK, j), 0) for j = 0:degree(fK)])
      imgs[s] = NS(evaluate(f, gpols[s]))
    end
    autsNS[t] = NfAbsNSToNfAbsNS(NS, NS, imgs)
  end
  #Now, to the simple field
  auts_abs = Vector{NfToNfMor}(undef, length(autsNS))
  gK = mK(gen(K))
  for i = 1:length(auts_abs)
    auts_abs[i] = hom(K, K, mK\(autsNS[i](gK)), check = false)
  end

  @vprint :AbExt 2 "Done. Now simplify and translate information\n"
  @vtime :AbExt 2 Ks, mKs = simplify(K)
  #Now, we have to construct the maximal order of this field.
  #I am computing the preimages of mKs by hand, by inverting the matrix.
  arr_prim_img = Array{nf_elem, 1}(undef, degree(Ks))
  arr_prim_img[1] = K(1)
  for i = 2:degree(Ks)
    arr_prim_img[i] = arr_prim_img[i-1]*mKs.prim_img
  end
  M1 = inv(basis_matrix(arr_prim_img, FakeFmpqMat))

  basisO2 = Array{nf_elem, 1}(undef, degree(Ks))
  M = zero_matrix(FlintZZ, 1, degree(Ks))
  for i=1:length(basisO2)
    elem_to_mat_row!(M, 1, denominator(O1.basis_nf[i]), O1.basis_nf[i])
    mul!(M, M, M1.num)
    basisO2[i] = elem_from_mat_row(Ks, M, 1, M1.den*denominator(O1.basis_nf[i]))
  end
  O2 = NfAbsOrd(Ks, basis_matrix(O1, copy = false)*M1)
  O2.ismaximal = 1
  _set_maximal_order_of_nf(Ks, O2)

  #Now, the automorphisms.
  autos = Array{NfToNfMor, 1}(undef, length(auts_abs))
  el = mKs(gen(Ks))
  for i = 1:length(auts)
    x = auts_abs[i](el)
    elem_to_mat_row!(M, 1, denominator(x), x)
    mul!(M, M, M1.num)
    y=Hecke.elem_from_mat_row(Ks, M, 1, M1.den*denominator(x))
    @assert iszero(Ks.pol(y))
    autos[i] = hom(Ks, Ks, y, check = false)
  end
  _set_automorphisms_nf(Ks, closure(autos, degree(Ks)))

  @vprint :AbExt 2 "Finished\n"
  return Ks, autos

end

function _from_relative_to_abs(L::NfRel_ns{T}, auts::Array{NfRel_nsToNfRel_nsMor{T}, 1}) where T

  S, mS = simple_extension(L)
  K, mK, mK2 = absolute_field(S, false)
  #First, we compute the maximal order of the absolute field.
  #Since the computation of the relative maximal order is slow, I bring to the absolute field the elements
  # generating the equation order.
  @vprint :AbExt 2 "Computing the maximal order\n"

  #First, I construct the product basis of the relative extension
  pols = L.pol
  gL = gens(L)
  B = Array{nf_elem, 1}(undef, degree(K))
  B[1] = K(1)
  ind = total_degree(pols[1])
  genjj = mK\(mS\gL[1])
  for i = 2:ind
    B[i] = genjj*B[i-1]
  end
  for jj = 2:length(pols)
    genjj = mK\(mS\gL[jj])
    el = genjj
    new_deg = total_degree(pols[jj])
    for i = 2:new_deg
      for j = 1:ind
        B[(i-1)* ind + j] = B[j]* el
      end
      mul!(el, el, genjj)
    end
    ind *= new_deg
  end

  #Now, I add the elements of the maximal order
  if degree(base_field(L)) > 1
    O = maximal_order(S.base_ring)
    for i = 1:degree(O)
      el = mK\(S(O.basis_nf[i]))
      for j = 1:ind
        B[(i-1)* ind + j] = B[j] * el
      end
    end
  end
  #Now, we compute the maximal order. Then we simplify.
  PO = Order(K, B, check = false, cached = false, isbasis = true)
  @vtime :AbExt 2 O1 = MaximalOrder(PO)
  O1.ismaximal = 1
  _set_maximal_order_of_nf(K, O1)
  @vprint :AbExt 2 "Done. Now simplify and translate information\n"
  @vtime :AbExt 2 Ks, mKs = simplify(K)
  #Now, we have to construct the maximal order of this field.
  #I am computing the preimages of mKs by hand, by inverting the matrix.
  arr_prim_img = Array{nf_elem, 1}(undef, degree(Ks))
  arr_prim_img[1] = K(1)
  for i = 2:degree(Ks)
    arr_prim_img[i] = arr_prim_img[i-1]*mKs.prim_img
  end
  M1 = inv(basis_matrix(arr_prim_img, FakeFmpqMat))
  basisO2 = Array{nf_elem, 1}(undef, degree(Ks))
  M = zero_matrix(FlintZZ, 1, degree(Ks))
  for i=1:length(basisO2)
    elem_to_mat_row!(M, 1, denominator(O1.basis_nf[i]), O1.basis_nf[i])
    mul!(M, M, M1.num)
    basisO2[i] = elem_from_mat_row(Ks, M, 1, M1.den*denominator(O1.basis_nf[i]))
  end
  O2 = Order(Ks, basisO2, check = false, cached = false)
  O2.ismaximal = 1
  _set_maximal_order_of_nf(Ks, O2)

  #Now, the automorphisms.
  autos=Array{NfToNfMor, 1}(undef, length(auts))
  el=mS(mK((mKs(gen(Ks)))))
  for i=1:length(auts)
    x = mK\(mS\(auts[i](el)))
    elem_to_mat_row!(M, 1, denominator(x), x)
    mul!(M, M, M1.num)
    y=Hecke.elem_from_mat_row(Ks, M, 1, M1.den*denominator(x))
    autos[i] = hom(Ks, Ks, y, check = false)
  end
  _set_automorphisms_nf(Ks, closure(autos, degree(Ks)))

  @vprint :AbExt 2 "Finished\n"
  return Ks, autos

end

###############################################################################
#
#  Maximal abelian subfield for fields function
#
###############################################################################

function check_abelian_extensions(class_fields::Vector{Tuple{Hecke.ClassField{Hecke.MapRayClassGrp,GrpAbFinGenMap}, Vector{GrpAbFinGenMap}}}, autos::Array{NfToNfMor, 1}, emb_sub::NfToNfMor)

  @vprint :MaxAbExt 3 "Starting checking abelian extension\n"
  K = base_field(class_fields[1][1])
  d = degree(K)
  G = domain(class_fields[1][2][1])
  expo = G.snf[end]
  com, uncom = ppio(Int(expo), d)
  if com == 1
    return Hecke.ClassField{Hecke.MapRayClassGrp,GrpAbFinGenMap}[x[1] for x in class_fields]
  end
  #I extract the generators that restricted to the domain of emb_sub are the identity.
  #Notice that I can do this only because I know the way I am constructing the generators of the group.
  expG_arr = Int[]
  act_indices = Int[]
  p = 11
  d1 = discriminant(domain(emb_sub))
  d2 = discriminant(K)
  while iszero(mod(d1, p)) || iszero(mod(d2, p))
    p = next_prime(p)
  end
  R = GF(p, cached = false)
  Rx, x = PolynomialRing(R, "x", cached = false)
  fmod = Rx(K.pol)
  mp_pol = Rx(emb_sub.prim_img)
  for i = 1:length(autos)
    pol = Rx(autos[i].prim_img)
    if mp_pol ==  Hecke.compose_mod(mp_pol, pol, fmod)
      push!(act_indices, i)
      #I compute the order of the automorphisms. I need the exponent of the relative extension!
      j = 2
      att = Hecke.compose_mod(pol, pol, fmod)
      while att != x
        att = Hecke.compose_mod(pol, att, fmod)
        j += 1
      end
      push!(expG_arr, j)
    end
  end
  expG = lcm(expG_arr)
  expG1 = ppio(expG, com)[1]
  com1 = ppio(com, expG1)[1]
  @vprint :MaxAbExt 3 "Context for ray class groups\n"

  OK = maximal_order(K)
  rcg_ctx = Hecke.rayclassgrp_ctx(OK, com1*expG1)

  @vprint :MaxAbExt 3 "Ordering the class fields\n"

  mins = Vector{fmpz}(undef, length(class_fields))
  for i = 1:length(mins)
    mins[i] = minimum(defining_modulus(class_fields[i][1])[1])
  end
  ismax = trues(length(mins))
  for i = 1:length(ismax)
    for j = i+1:length(ismax)
      if ismax[j]
        i2 = ppio(mins[i], mins[j])[2]
        if isone(i2)
          ismax[i] = false
          break
        end
        i3 = ppio(mins[j], mins[i])[2]
        if isone(i3)
          ismax[j] = false
        end
      end
    end
  end
  ord_class_fields = Vector{Int}(undef, length(ismax))
  j1 = 1
  j2 = length(ismax)
  for i = 1:length(ismax)
    if ismax[i]
      ord_class_fields[j1] = i
      j1 += 1
    else
      ord_class_fields[j2] = i
      j2 -= 1
    end
  end

  cfields = Hecke.ClassField{Hecke.MapRayClassGrp, GrpAbFinGenMap}[]
  for i = 1:length(class_fields)
    @vprint :MaxAbExt 3 "Class Field $i\n"
    C, res_act = class_fields[ord_class_fields[i]]
    res_act_new = Vector{GrpAbFinGenMap}(undef, length(act_indices))
    for i = 1:length(act_indices)
      res_act_new[i] = res_act[act_indices[i]]
    end
    fl = check_abelian_extension(C, res_act_new, emb_sub, rcg_ctx)
    if fl
      push!(cfields, C)
    end
  end
  return cfields

end

function check_abelian_extension(C::Hecke.ClassField, res_act::Vector{GrpAbFinGenMap}, emb_sub::NfToNfMor, rcg_ctx::Hecke.ctx_rayclassgrp)

  #I consider the action on every P-sylow and see if it is trivial
  G = codomain(C.quotientmap)
  expG = Int(exponent(G))
  fac = factor(rcg_ctx.n)
  prime_to_test = Int[]
  new_prime = false
  for (P, v) in fac
    # I check that the action on the P-sylow is the identity.
    for i = 1:ngens(G)
      if !divisible(G.snf[i], P)
        continue
      end
      pp, q = ppio(G.snf[i], P)
      new_prime = false
      for j = 1:length(res_act)
        if res_act[j](q*G[i]) != q*G[i]
          new_prime = true
          break
        end
      end
      if new_prime
        break
      end
    end
    if !new_prime
      push!(prime_to_test, P)
    end
  end
  if isempty(prime_to_test)
    return true
  end

  n = prod(prime_to_test)
  n1, m = ppio(Int(G.snf[end]), n)
  if m != 1
    # Now, I compute the maximal abelian subextension of the n-part of C
    Q, mQ = quo(G, n1, false)
    C1 = ray_class_field(C.rayclassgroupmap, Hecke.GrpAbFinGenMap(Hecke.compose(C.quotientmap, mQ)))
    #@vtime :MaxAbExt 1
    fl = _maximal_abelian_subfield(C1, emb_sub, rcg_ctx, expG)
  else
    #@vtime :MaxAbExt 1
    fl = _maximal_abelian_subfield(C, emb_sub, rcg_ctx, expG)
  end
  return fl

end

function _bound_exp_conductor_wild(O::NfOrd, n::Int, q::Int, bound::fmpz)
  d = degree(O)
  lp = prime_decomposition_type(O, q)
  f_times_r = divexact(d, lp[1][2])
  sq = fmpz(q)^f_times_r
  nbound = n+n*lp[1][2]*valuation(n,q)-div(fmpz(n), q^(valuation(n,q)))
  obound = flog(bound, sq)
  bound_max_ap = min(nbound, obound)  #bound on ap
  return div(q*bound_max_ap, n*(q-1)) #bound on the exponent in the conductor
end

function minimumd(D::Dict{NfOrdIdl, Int}, deg_ext::Int)
  primes_done = Int[]
  res = 1
  for (P, e) in D
    p = Int(minimum(P))
    if p in primes_done
      continue
    else
      push!(primes_done, p)
    end
    ram_index = P.splitting_type[1]
    s, t = divrem(e, ram_index)
    if iszero(t)
      d = min(s, valuation(deg_ext, p)+2)
      res *= p^d
    else
      d = min(s+1, valuation(deg_ext, p)+2)
      res *= p^d
    end
  end
  return res
end

function _maximal_abelian_subfield(A::Hecke.ClassField, mp::Hecke.NfToNfMor, ctx::Hecke.ctx_rayclassgrp, expG::Int)

  K = base_field(A)
  k = domain(mp)
  ZK = maximal_order(K)
  zk = maximal_order(k)
  expected_order = div(degree(K), degree(k))
  if gcd(expected_order, degree(A)) == 1
    return false
  end
  # disc(ZK/Q) = N(disc(ZK/zk)) * disc(zk)^deg
  # we need the disc ZK/k, well a conductor.
  d = abs(div(discriminant(ZK), discriminant(zk)^expected_order))
  mR1 = A.rayclassgroupmap
  expo = Int(exponent(codomain(A.quotientmap)))

  #First, a suitable modulus for A over k
  #I take the discriminant K/k times the norm of the conductor A/K

  fac1 = factor(d)
  fm0 = Dict{NfOrdIdl, Int}()
  for (p, v) in fac1
    lPp = prime_decomposition(zk, p)
    if divisible(fmpz(expected_order), p)
      theoretical_bound = _bound_exp_conductor_wild(zk, expG, Int(p), d)
      for i = 1:length(lPp)
        fm0[lPp[i][1]] = min(theoretical_bound, Int(v))
      end
    else
      for i = 1:length(lPp)
        fm0[lPp[i][1]] = 1
      end
    end
  end
  #Now, I want to compute f_m0 = merge(max, norm(mR1.fact_mod), fac2)
  primes_done = fmpz[]
  for (P, e) in mR1.fact_mod
    p = minimum(P)
    if p in primes_done
      continue
    else
      push!(primes_done, p)
    end
    lp = prime_decomposition(zk, p)
    if !divisible(fmpz(expected_order * expo), p)
      for i = 1:length(lp)
        fm0[lp[i][1]] = 1
      end
    else
      #I need the relative norm of P expressed as a prime power.
      pm = Hecke.intersect_prime(mp, P)
      fpm = divexact(P.splitting_type[2], pm.splitting_type[2])
      theoretical_bound1 = _bound_exp_conductor_wild(zk, lcm(expo, expG), Int(p), d*norm(defining_modulus(A)[1]))
      for i = 1:length(lp)
        if haskey(fm0, lp[i][1])
          fm0[lp[i][1]] =  min(fm0[lp[i][1]] * fpm * e, theoretical_bound1)
        else
          fm0[lp[i][1]] = min(fpm * e, theoretical_bound1)
        end
      end
    end
  end
  # Now, I extend the modulus to K
  fM0 = Dict{NfOrdIdl, Int}()
  primes_done = fmpz[]
  for (P, e) in fm0
    p = Hecke.minimum(P)
    if p in primes_done
      continue
    else
      push!(primes_done, p)
    end

    lp = prime_decomposition(ZK, p)
    multip = divexact(lp[1][2], P.splitting_type[1])
    if !divisible(fmpz(expected_order * expo), p)
      for i = 1:length(lp)
        fM0[lp[i][1]] = 1
      end
    else
      if isdefined(A, :abs_disc)
        d = A.abs_disc
        ev = prod(w^z for (w,z) in d)
        ev = divexact(ev, discriminant(ZK))
        theoretical_bound2 = _bound_exp_conductor_wild(ZK, expo, Int(p), ppio(ev, p)[1])
        for i = 1:length(lp)
          fM0[lp[i][1]] = min(multip * e, theoretical_bound2)
        end
      else
        for i = 1:length(lp)
          fM0[lp[i][1]] = multip * e
        end
      end
    end
  end
  ind = 0
  #@vtime :MaxAbExt 1
  if isdefined(ctx, :computed)
    flinf = isempty(mR1.modulus_inf)
    for i = 1:length(ctx.computed)
      idmr, ifmr, mRRR = ctx.computed[i]
      if flinf != ifmr
        continue
      end
      contained = true
      for (P, v) in fM0
        if !haskey(idmr, P) || idmr[P] < v
          contained = false
        end
      end
      if contained
        ind = i
        break
      end
    end
  end
  if iszero(ind)
    R, mR = Hecke.ray_class_group_quo(ZK, fM0, mR1.modulus_inf, ctx, check = false)
    if isdefined(ctx, :computed)
      push!(ctx.computed, (fM0, isempty(mR1.modulus_inf), mR))
    else
      ctx.computed = [(fM0, isempty(mR1.modulus_inf), mR)]
    end
  else
    mR = ctx.computed[ind][3]
    R = domain(mR)
  end
  if degree(zk) != 1
    if istotally_real(K) && isempty(mR1.modulus_inf)
      inf_plc = InfPlc[]
    else
      inf_plc = real_places(k)
    end
    #@vtime :MaxAbExt 1
    r, mr = Hecke.ray_class_group_quo(zk, ctx.n, fm0, inf_plc)
  else
    rel_plc = true
    if istotally_real(K) && isempty(mR1.modulus_inf)
      rel_plc = false
    end
    modulo = minimumd(fm0, expo * expected_order)
    #@vtime :MaxAbExt 1
    r, mr = Hecke.ray_class_groupQQ(zk, modulo, rel_plc, ctx.n)
  end
  #@vtime :MaxAbExt 1
  lP, gS = Hecke.find_gens(mR, coprime_to = minimum(mR1.modulus_fin))
  listn = NfOrdIdl[norm(mp, x) for x in lP]
  # Create the map between R and r by taking norms
  preimgs = Vector{GrpAbFinGenElem}(undef, length(listn))
  for i = 1:length(preimgs)
    preimgs[i] = mr\listn[i]
  end
  proj = hom(gS, preimgs)
  #compute the norm group of A in R
  prms = Vector{GrpAbFinGenElem}(undef, length(lP))
  for i = 1:length(lP)
    if haskey(mR1.prime_ideal_preimage_cache, lP[i])
      f = mR1.prime_ideal_preimage_cache[lP[i]]
    else
      f = mR1\lP[i]
      mR1.prime_ideal_preimage_cache[lP[i]] = f
    end
    prms[i] = A.quotientmap(f)
  end
  proj1 = hom(gS, prms)
  S, mS = kernel(proj1, false)
  mS1 = mS*proj
  G, mG = Hecke.cokernel(mS1, false)
  return (order(G) == expected_order)::Bool

end


text:
using HDF5
using Base
using ProgressMeter

function Loadh5Param(dir)
    h5file = dir * "ParamLLEJulia.h5"
    res = Dict()
    sim = Dict()
    sim_name = ["res", "sim"]
    par = [["Qc", "R", "ng", "Qi", "gamma","dispfile"], ["dphi","Pin", "Tscan", "domega_init", "domega_end", "f_pmp", "mu_sim", "debug", "ind_aux"]]
    cnt = 1
    for sim_par = [res, sim]
        for it = par[cnt]
            sim_par[it] = h5read(h5file, sim_name[cnt] *"/" * it)
        end
        cnt += 1
    end

    return res, sim
end

function SaveResults(dir, S)
    h5file = dir * "ResultsJulia.h5"
    # print(h5file)
    # print("\n")
    h5open(h5file, "w") do file
        g = g_create(file, "Results") # create a group
        for ii in S
            g[ii[1]*"Real"] = real(ii[2])              # create a scalar dataset inside the group
            g[ii[1]*"Imag"] = imag(ii[2])
        end
        attrs(g)["Description"] = "This group contains only a single dataset" # an attribute
    end
        # attrs(g)["Description"] = "This group contains only a single dataset" # an attribute
end

# ----------------------------------------------------
# -- STARTING MAIN SCRIPT --
# ----------------------------------------------------
# FFTW.set_num_threads(Sys.CPU_CORES)

tmp_dir = ARGS[1]
tol= float(ARGS[2])
maxiter = Int(float(ARGS[3]))
stepFactor = float(ARGS[4])
res, sim = Loadh5Param(tmp_dir)
logfile =  open(tmp_dir * "Log.log" ,"w")

# -- Collect Data ---
# ----------------------------------------------------
c0 = 299792458
ħ = 6.634e-34/2/pi

# -- res parameters --
ng = res["ng"][1]
R = res["R"][1]
γ = res["gamma"][1]
L = 2*pi*R

Q0 = res["Qi"][1]
Qc = res["Qc"][1]

# -- sim parameters --
debug = Bool(sim["debug"][1])
δω_init = sim["domega_init"][1]
δω_end = sim["domega_end"][1]
t_end = sim["Tscan"][1]
fpmp = sim["f_pmp"][1]
Pin = sim["Pin"][1]
if length(sim["f_pmp"])>1
    faux = sim["f_pmp"][2:length(sim["f_pmp"])]
    Paux = sim["Pin"][2:length(sim["f_pmp"])]
else
    faux = []
    Paux = []
end
ω0 = 2*pi*fpmp
dφ = sim["dphi"]
μ_sim = sim["mu_sim"]
μ = collect(μ_sim[1]:μ_sim[2])
pmp_sim = find(μ .== 0)[1]
ind_aux = sim["ind_aux"]
write(logfile,string(ind_aux))
print(string(ind_aux))
# -- Setup the different Parameters --
# ----------------------------------------------------
tR = L*ng/c0
f0 = c0/ng/L
T = 1*tR
θ = ω0/Qc*tR
α=1/2* (ω0/Q0 + ω0/Qc) * tR

# -- losses --
dω = collect(μ_sim[1]:μ_sim[end])*2*pi/T
Qc_disp=Qc*ones(size(dω))
Qi_disp=Q0*ones(size(dω))

# -- Dispersion --
dβ=dφ/L
dβ=dβ-dβ[pmp_sim]
dφ_Kerr=θ*Pin*γ*L/α^2

# -- Input Energy --
Ein= zeros(size(μ))
Ein[pmp_sim]=sqrt.(Pin)*length(μ)

for ii=1:length(faux)
    Ein[pmp_sim - ind_aux[ii]] = sqrt.(Paux[ii])*length(μ)
end
Ein_couple=sqrt(θ).*Ein

# -- Noise Background --
Ephoton=ħ*(ω0+dω)
phase=2*pi*(rand(1,length(μ)))
array=rand(1,length(μ))
Enoise=array'.*sqrt.(Ephoton/2/tR).*exp.(1im*phase') .*length(μ)

dt=0.1/(sqrt(Pin))
δω_init = δω_init * tR
δω_end = δω_end * tR

t_end = t_end*tR
t_ramp=t_end

# -- Seting up Simulation --
# ----------------------------------------------------

# -- Allocate FFT to go faster --
ifft_plan = plan_ifft(zeros(size(Enoise)))
fft_plan = plan_fft(zeros(size(Enoise)))

# -- Sim length --
tol=1e-3
# tol=1e-20
maxiter=6
Nt= round(t_ramp/tR/dt)
t1=0

# -- Detuning ramp --
xx = collect(1:Nt)
Δω_pmp = δω_init+ xx/Nt * (δω_end- δω_init)

# -- Initiial State --
u0=ifft_plan*Enoise

# -- Output Dict --
S = Dict()
num_probe=1000
S["u_probe"] = 1im*zeros(num_probe, length(u0))
S["Em_probe"] = 1im*zeros(num_probe, length(u0))
S["comb_power"] = zeros(num_probe,)
S["detuning"] = zeros(num_probe,)

# -- Progress Bar --
probe=0
# pb = Progress(Int(Nt),0.5, "Calculating temporal loop:", Int(25))

space = "."^100
logfile =  open(tmp_dir * "log.log" ,"a")
write(logfile,string(Int(round(0))) * "\n")
close(logfile)
probe_pbar = 1/100
cnt_pbar = 0
# -- Main Solver --
# ----------------------------------------------------
for it = 1:1:Nt
    # print("coucouc\n")
    if it/Nt >=probe_pbar
        cnt_pbar = cnt_pbar + 1
        # print("\b"^103)
        # print("coucou")
        # print("="^cnt_pbar)
        # print("."^(space-cnt_pbar))
        # print(" ")
        logfile =  open(tmp_dir * "log.log" ,"a")
        write(logfile,string(Int(round(probe_pbar*100))) * "\n")
        close(logfile)
        probe_pbar = probe_pbar + 0.01
    end
    u0=ifft_plan*(fft_plan*(u0) + Ein_couple*dt)
    u1=u0
    cbeta = -α  + 1im*Δω_pmp[Int(it)] +  1im*L*dβ
    halfprop = exp.(cbeta*dt/2)
    cnt = 0
    uhalf = ifft_plan*(halfprop.*(fft_plan*(u0)))
    half1 = ifft_plan*(γ*(fft_plan*( abs.(u0).^2.*u0) ) )./u0
    for ii = 1:maxiter
        half2 = ifft_plan*(γ*(fft_plan*( abs.(u1).^2.*u1)))./u1
        uv = uhalf .* exp.(1im*L.*(half1 + half2)*dt/2)
        uv2 = ifft_plan*(halfprop.*(fft_plan*(uv)))
        if (norm(uv2-u1,2)/norm(u1,2) < tol)
            u1 = uv2
            break
        else
            u1 = uv2
        end
        cnt += 1
    end

    if (cnt == maxiter)
        logfile =  open(tmp_dir * "log.log" ,"a")
        write(logfile,"Failed to converge to N=" * string(it) * " over Nt=" * string(Nt))
        close(logfile)
    end
    u0=u1

   if (it*num_probe/Nt > probe)
        probe += 1
        S["u_probe"][probe,:]=u0[:,1]
        Em_probe = (fft_plan*(u0))/length(μ)
        Em_probe = Em_probe[:,1]
        S["Em_probe"][probe,:]= Em_probe
        deleteat!(Em_probe,pmp_sim)
        S["comb_power"][probe]=(sum(abs.(Em_probe).^2))/Pin
        S["detuning"][probe] = Δω_pmp[Int(it)]/tR
   end

end
S["Ewg"] = 1im*zeros(size(S["Em_probe"]))
for ii=1:size(S["Em_probe"],1)
    S["Ewg"][ii,:] = Ein/length(μ)-S["Em_probe"][ii,:].*sqrt(θ)
end
S["ω"] = (ω0 + dω)

# print("\n")
# print("Congrats, everything has been solved.\n")
# print("Now we are going to save the results in a .h5 file\n")
# print("Everything will be in: ")
logfile =  open(tmp_dir * "log.log" ,"a")
write(logfile,string(Int(round(100))) * "\n")
close(logfile)

SaveResults(tmp_dir, S)
# write(logfile,"Great we have saved everything... Moving back to Python\n")

ext: jl

max_stars_repo_name: jiatianxu/pyLLE

max_issues_repo_name: jiatianxu/pyLLE

max_forks_repo_name: jiatianxu/pyLLE


text:
#=
1. Create a right cylinder C:
	(a) the base of C is a regular polygon with N sides.
	(b) the base is parallel to the xy-plane.
2. Create a copy C′ of C.
3. Rotate C′ around its vertical centerline by α degrees.
4. Measure time for C ∪ C′.
=#

using ViewerGL; GL = ViewerGL
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation;


text:
using Printf

"""
Iterator for the sifting algorithm.
"""
mutable struct SiftIterable{solT <: AbstractVector, vecT <: AbstractVector, numT <: Real}
    "Measurements/IMF"
    y::solT
    "Time series values"
    x::vecT
    "Positions of the local maxima"
    maxs::Vector{Int}
    "Positions of the local minima"
    mins::Vector{Int}
    "Indices of the zero crossings"
    crosses::Vector{Int}
    "Tolerance"
    tol::numT
    "Residuals used to check stopping criteria"
    residual::numT
    "Maximum number of iteration"
    maxiter::Int
    "Number of steps of a stable sift after which the sifting is aborted"
    stop_steps::Int
    "Number of steps on which the number of zero crossings was fix"
    fix_steps::Int
end

@inline converged(it::SiftIterable) = it.residual ≤ var(it.y) * it.tol

@inline start(it::SiftIterable) = 0

@inline done(it::SiftIterable, iteration::Int) = iteration ≥ it.maxiter || converged(it)

# kernel sifting algorithm
function iterate(it::SiftIterable, iteration::Int=start(it))
    # Check for termination first
    if done(it, iteration)
        return nothing
    end

    localmaxmin!(it.y, it.maxs, it.mins)
    maxlen = length(it.maxs)
    minlen = length(it.mins)

    it.fix_steps == it.stop_steps && return nothing

    zerocrossing!(it.y, it.crosses)

    abs(length(it.crosses) - maxlen - minlen) <= 1 && (it.fix_steps += 1)
    if maxlen < 4 || minlen < 4
        return nothing
    end
    smin = get_edgepoint(it.y, it.x, it.mins, first, isless)
    smax = get_edgepoint(it.y, it.x, it.maxs, first, !isless)
    emin = get_edgepoint(it.y, it.x, it.mins, last, isless)
    emax = get_edgepoint(it.y, it.x, it.maxs, last, !isless)

    # TODO: Format for clarity?
    maxTS = interpolate([first(it.x); it.x[it.maxs]; last(it.x)], [smax; it.y[it.maxs]; emax], it.x, DierckXInterp())
    minTS = interpolate([first(it.x); it.x[it.mins]; last(it.x)], [smin; it.y[it.mins]; emin], it.x, DierckXInterp())
    subs = 0.5*(maxTS + minTS)
    it.residual = sum(abs, subs)

    it.y .-= subs

    # Return the residual at item and iteration number as state
    it.residual, iteration + 1
end

# Utility functions

function sift_iterator(y, x;
                      reltol::Real = 0.1,
                      maxiter::Int = 6,
                      stop_steps::Int = 4)
    maxs, mins = Int[], Int[]
    crosses = Int[]
    residual = sum(abs, y)
    tolerance = 0.1
    fix_steps = 0
    # Return the iterable
    return SiftIterable(y, x, maxs, mins, crosses,
        tolerance, residual, maxiter, stop_steps, fix_steps)
end

"""
    sift(y, x=1:length(y); kwargs...) -> imf

Same as [`sift!`](@ref), but allocates a solution vector `imf` initialized with zeros.
"""
sift(y, x=1:length(y); kwargs...) = sift!(zeros(eltype(y), size(y)), y, x; kwargs...)

"""
    sift!(imf, y, x; kwargs...) -> imf

# Arguments

- `imf`: intrinsic mode function, will be updated in-place;
- `y`: data measurements;
- `x`: time stamp series.

## Keywords

- `reltol::Real = 0.1`: relative tolerance for the stopping condition;
- `maxiter::Int = 10`: maximum number of iterations;
- `verbose::Bool = false`: print method information;
- `log::Bool = false`: keep track of the residual in each iteration.

# Output

**if `log` is `false`**

- `imf`: intrinsic mode function.
"""
function sift!(imf, y, x;
             reltol::Real = 0.1,
             maxiter::Int = 6,
             stop_steps::Int = 4,
             log::Bool = false,
             verbose::Bool = false,
             kwargs...)

    imf .= y
    # Actually perform sifting
    iterable = sift_iterator(imf, x; reltol, maxiter, stop_steps)

    for (iteration, item) = enumerate(iterable)
        verbose && @printf("%3d\t%1.2e\n", iteration, iterable.residual)
    end

    verbose && println()

    iterable.y
end


text:
function x = cell_set(x,i,v)

    """
     cell_set - assign value

        x = cell_set(x,i,v);
            is equivalent to
        x{i}=v;

     Copyright (c) 2008 Gabriel Peyre
     """
     x[i] = v
 end


text:
# ------------------------------------------------------------------------------------------
# ## Functions
#
#
# In the last notebook, we talked about modeling data with functions. A **function** is one
# of the most fundamental concepts in computing (and also in mathematics).
#
# A function is a piece of a program that receives **input arguments**, processes them by
# doing certain calculations on them, and returns **outputs**.
#
# For example, we might have a function `g` that takes a number as an input and returns the
# square of that number as an output. How can we define this function `g` on a computer?
# Julia gives us a few different ways to do this.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ### Defining functions
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# Firstly, we could write `g` as follows:
# ------------------------------------------------------------------------------------------

g(x) = x^2

a = "Machine is learning is fun "
g(a)

# ------------------------------------------------------------------------------------------
# Alternatively, we could declare this function using the `function` and `end` keywords:
# ------------------------------------------------------------------------------------------

function g1(x)
    return x^2
end

# ------------------------------------------------------------------------------------------
# The third way we could have declared this function is as an "anonymous" or "lambda"
# function. "Anonymous" functions are functions that truly don't need names! For example, we
# could have declared a function that squares its input as
# ------------------------------------------------------------------------------------------

(x -> x^2)("I ♡ Julia. ") # \heartsuit + <tab>

# ------------------------------------------------------------------------------------------
# Now that we've done that, we can't access the function `x -> x^2` again because we have no
# name to call! That seems a little silly, doesn't it?
#
# Actually, there are times where functions without names are useful to us. We'll see that
# later in this notebook. For now, note that you have the *option* to access an "anonymous"
# function later by binding a variable to it when you declare it. For example,
# ------------------------------------------------------------------------------------------

g = x -> x^2

# ------------------------------------------------------------------------------------------
# This syntax says, "I want to use the variable g to access a function that takes some input
# called `x` and maps that input to the square of `x`.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## An important sigmoidal function
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# A particular function that is used a lot in machine learning is a so-called "sigmoidal"
# function (meaning a function that is S-shaped, i.e. the graph of the function looks like
# an `S`).
#
# The sigmoid function that we will use is given the name $\sigma$, and is defined by the
# following mathematical expression:
#
# $$\sigma(x) := \frac{1}{1 + \exp(-x)}.$$
# ------------------------------------------------------------------------------------------

σ(w*x) = 1/ (1+e^(-w*x))

# ------------------------------------------------------------------------------------------
# #### Exercise 1
#
# Use the first syntax given above to define the function `σ` in Julia. Note that Julia
# actually allows us to use the symbol σ as a variable name! To do so, type `\sigma<TAB>` in
# the code cell.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## Plotting functions
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# Let's draw the function σ to see what it looks like. Throughout this course, we'll use the
# Julia package `Plots.jl` for all of the graphics. This package provides a flexible syntax
# for plotting, in which options to change attributes like the width of the lines used in
# the figure are given as named keyword arguments.
#
# In addition, it allows us to use different "backends", which are the other libraries that
# actually carry out the plotting following the instructions from `Plots.jl`.
# ------------------------------------------------------------------------------------------

using Plots
gr()   # use the PlotlyJS "backend" (plotting library)

plot(σ, -5, 5)

hline!([0, 1], ls=:dash, lw=3)  # add horizontal lines at 0 and 1, with dashed style and linewidth 3
vline!([0], ls=:dash, lw=3)     # add a vertical line at 0

# ------------------------------------------------------------------------------------------
# We can think of $\sigma$ as a smooth version of a step or threshold function (often called
# a "Heaviside" function). To see this, let's modify the steepness of the jump in $\sigma$
# and compare it to the Heaviside function; we'll see how all this works in more detail
# later:
# ------------------------------------------------------------------------------------------

using Interact

heaviside(x) = x < 0 ? 0.0 : 1.0

@manipulate for w in 0.1:0.1:20
    plot(x -> σ(w*x), -5, 5, label="sigma", lw=2)
    plot!(heaviside, ls=:dash, label="step")
end

# ------------------------------------------------------------------------------------------
# This particular function takes any real number as input, and gives an output between $0$
# and $1$. It is continuous and smooth.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# #### Exercise 2
#
# Declare the sigmoid function above as an anonymous function with a different name.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ### Mutating functions: `...!`
#
# To generate our plot of σ above, we used some functions that end with `!`. What does a `!`
# at the end of a function name mean in Julia?
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# Functions that change or modify their inputs are called **mutating functions**. But wait,
# don't all functions do that?
#
# No, actually. Functions typically take *inputs* and use those *inputs* to generate
# *outputs*, but the inputs themselves usually don't actually get changed by a function. For
# example, copy and execute the following code:
#
# ```julia
# v1 = [9, 4, 7, 11]
# v2 = sort(v1)
# ```
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# `v2` is a sorted version of `v1`, but after calling `sort`, `v1` is still unsorted.
#
# However, now trying adding an exclamation point after `sort` and executing the following
# code:
#
# ```julia
# sort!(v1)
# ```
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# Look at the values in `v1` now!
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# This time, the original vector itself was changed (mutated), and is now sorted. Unlike
# `sort`, `sort!` is a mutating function. Did the `!` make `sort!` mutating? Well, no, not
# really. In Julia, `!` indicates mutating functions by convention. When the author of
# `sort!` wrote `sort!`, they added a `!` to let you to know that `sort!` is mutating, but
# the `!` isn't what makes a function mutating or non-mutating in the first place.
#
# #### Exercise
#
# Some of our plotting commands end with `!`. Copy and execute the following code:
#
# ```julia
# r = -5:0.1:5
# g(x) = x^2
# h(x) = x^3
# plot(r, g, label="g")
# plot!(r, h, label="h")
# ```
#
# Then change the code slightly to remove the `!` after `plot!(r, h)`. How does this change
# your output? What do you think it means to add `!` after plotting commands?
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## Pointwise application of functions, `.(...)` (known as "broadcasting")
#
# We saw in a previous notebook that we needed to add `.` after the names of some functions,
# as in
#
# ```julia
# green_amount = mean(Float64.(green.(apple)))
# ```
#
# What are those extra `.`s really doing?
#
# When we add a `.` after a function's name, we are telling Julia that we want to
# "**broadcast**" that function over the inputs passed to the function. This means that we
# want to apply that function *element-wise* over the inputs; in other words, it will apply
# the function to each element of the input, and return an array with the newly-calculated
# values.
#
# For example, copy and execute the following code:
# ```julia
# g.(r)
# ```
# Since the function `g` squares it's input, this squares all the elements of the range `r`.
#
# What happens if instead we just call `g` on `r` via
#
# ```julia
# g(r)
# ```
# ? Try this out and see what happens.
# ------------------------------------------------------------------------------------------

x = [1 2 3;4 5 6]
f(t) = √t + 5
f.(x)

# ------------------------------------------------------------------------------------------
# You should see an error message after calling `g(r)`, which says that Julia cannot
# multiply two vectors. When we call `g(r)`, we ask Julia to multiply `r` by `r`. When we
# call `g.(r)`, we ask Julia to multiply *each element* in `r` by itself.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# #### Exercise 3
#
# Copy and execute the following code to get the type of the object `numbers = [1, 2,
# "three", 4.0]`:
#
# ```julia
# numbers = [1, 2, "three", 4.0]
# typeof(numbers)
# ```
#
# What is the type of `numbers`?
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# #### Exercise 4
#
# Broadcast `typeof` over `numbers` to see what the types of the elements stored inside
# `numbers` are.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# #### Exercise 5
#
# Write a `for` loop that applies `g` to each of the elements of `r` and prints the results.
# Verify that the numbers printed by this `for` loop are equal to the entries of the output
# of `g.(r)`.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# #### Exercise 6
#
# Define a range `xs` between -5 and 5 with steps of 0.5.
# Apply the $\sigma$ function pointwise to this range and define `ys` as the result.
# What does the result look like? Plot these as points and join them with lines.
#
# Make the plot interactive where you can vary the step size. Fix the range of the plot in
# the `x` and `y` directions using the functions `xlims!` and `ylims!`.
# ------------------------------------------------------------------------------------------

max_stars_repo_path: .nbexports/introductory-tutorials/broader-topics-and-ecosystem/intro-to-ml/04. Tools - Functions.jl

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
println('\n', " "^4, "> Loading the packages...")

using LinearAlgebra # Norm
using Plots # Plotting

# Use the GR backend for plots
gr()

# Change the default font for plots
default(fontfamily="Computer Modern", dpi=300, legend=:topright)

# Define the value of ϰ
const ϰ = -1

# Integrate equations of motion using the
# Euler's method, return the values of position
# and velocity on each step
function motion_euler(
    r₀::Vector{F},
    v₀::Vector{F},
    h::F,
    n::I,
)::Tuple{Matrix{F},Matrix{F}} where
{F <: AbstractFloat,I <: Unsigned}
    # Determine the length of the input vectors
    N = length(r₀)
    # Prepare output matrixes
    r = Matrix{F}(undef, n + 1, N)
    v = Matrix{F}(undef, n + 1, N)
    r[1, :] = copy(r₀)
    v[1, :] = copy(v₀)
    # Compute the solutions
    for i in 2:(n + 1)
        ρ = norm(r[i - 1, :])^3
        for k in 1:N
            r[i, k] = r[i - 1, k] + h * v[i - 1, k]
            v[i, k] = v[i - 1, k] + h * ϰ * r[i - 1, k] / ρ
        end
    end
    return r, v
end

# Integrate equations of motion using the
# two-step Adams–Bashforth's method, return
# the values of position and velocity on each step
function motion_ab2(
    r₀::Vector{F},
    v₀::Vector{F},
    h::F,
    n::I,
)::Tuple{Matrix{F},Matrix{F}} where
{F <: AbstractFloat,I <: Unsigned}
    # Determine the length of the input vectors
    N = length(r₀)
    # Prepare output matrixes
    r = Matrix{F}(undef, n + 1, N)
    v = Matrix{F}(undef, n + 1, N)
    r[1, :] = copy(r₀)
    v[1, :] = copy(v₀)
    # Compute the second value of the solution
    # by using the one-step Euler's method
    ρ = norm(r[1, :])^3
    for k in 1:N
        r[2, k] = r[1, k] + h * v[1, k]
        v[2, k] = v[1, k] + h * ϰ * r[1, k] / ρ
    end
    # Define a couple of independent coefficients
    k₁ = 3 / 2 * h
    k₂ = -1 / 2 * h
    # Compute the rest in two steps
    for i in 3:(n + 1)
        ρₖ = norm(r[i - 1, :])^3
        ρₖ₋₁ = norm(r[i - 2, :])^3
        # Define a couple of dependent coefficients
        k₃ = k₁ * ϰ / ρₖ
        k₄ = k₂ * ϰ / ρₖ₋₁
        for k in 1:N
            r[i, k] = r[i - 1, k] + k₁ * v[i - 1, k] + k₂ * v[i - 2, k]
            v[i, k] = v[i - 1, k] + k₃ * r[i - 1, k] + k₄ * r[i - 2, k]
        end
    end
    return r, v
end

# Integrate equations of motion using all available
# methods and create a comparison plot between the solutions
function motion_integrate(rₐ::Matrix{F}, r₀::Vector{F}, v₀::Vector{F}, h::F, n::I) where {F <: AbstractFloat,I <: Unsigned}
    rₑ, _ = motion_euler(r₀, v₀, h, n)
    rₐ₂, _ = motion_ab2(r₀, v₀, h, n)

    # Plot the exact solution
    p = plot(
        rₐ[1:1000:end, 1],
        rₐ[1:1000:end, 2],
        label="Т",
        xlabel="x",
        ylabel="y",
        xlims=extrema(rₐ[:, 1]) .+ (-0.1, 0.1),
        ylims=extrema(rₐ[:, 2]) .+ (-0.1, 0.1),
    );

    # Add the solution achieved by
    # the Euler's method to the plot
    plot!(p, rₑ[:, 1], rₑ[:, 2]; label="Э");

    # Add the solution achieved by the two-step
    # Adams–Bashforth's method to the plot
    plot!(p, rₐ₂[:, 1], rₐ₂[:, 2]; label="АБ-2");

    # Save the figure
    savefig(p, "$(@__DIR__)/../plots/static/orbit_$(h).png")
    savefig(p, "$(@__DIR__)/../plots/static/orbit_$(h).pdf")

    println(" "^6, "* The figures `orbit_$(h).png` and `orbit_$(h).pdf` are saved. *")

    # Create an animation of the same thing but with moving circles
    anim = @animate for i in 1:10^(length("$(n)") - 3):n
        # Recreate the static plot
        plot(
            rₐ[1:1000:end, 1],
            rₐ[1:1000:end, 2],
            label="Т",
            xlabel="x",
            ylabel="y",
            xlims=extrema(rₐ[:, 1]) .+ (-0.1, 0.1),
            ylims=extrema(rₐ[:, 2]) .+ (-0.1, 0.1),
        );
        plot!(rₑ[:, 1], rₑ[:, 2]; label="Э")
        plot!(rₐ₂[:, 1], rₐ₂[:, 2]; label="АБ-2")
        # Add moving points
        scatter!([rₑ[i, 1],], [rₑ[i, 2],]; label="")
        scatter!([rₐ₂[i, 1],], [rₐ₂[i, 2],]; label="")
    end

    # Save the animation
    gif(anim, "$(@__DIR__)/../plots/animated/orbit_$(h).gif"; fps=30, show_msg=false)

    println(" "^6, "* The animation `orbit_$(h).gif` is saved. *")

end

println(" "^4, "> Integrating the equations of motion...")

# Define the initial values of the position and velocity
r₀ = [1.0, 0.0]
v₀ = [0.0, 0.5]

# Calculate the (close-to) analytical solution
rₐ, _ = motion_ab2(r₀, v₀, 0.000001, UInt(2714081))

println()
motion_integrate(rₐ, r₀, v₀, 0.01, UInt(271))
motion_integrate(rₐ, r₀, v₀, 0.001, UInt(2714))
motion_integrate(rₐ, r₀, v₀, 0.0001, UInt(27141))
motion_integrate(rₐ, r₀, v₀, 0.00001, UInt(271408))
println()


text:
struct KFolds{T <: Persa.AbstractDataset} <: ResamplingStrategy{T}
    dataset::T
    index::Array{Int}
    k::Int
end

function KFolds(dataset::Persa.AbstractDataset; k::Int = 10, is_shuffle::Bool = true)
  return KFolds(dataset, splitKFold(length(dataset), k, is_shuffle), k)
end

function getTrainData(kfold::KFolds{T}, fold::Int) :: T where T
    index = findall(r -> r != fold, kfold.index)
    return Persa.sample(kfold.dataset, index)
end

function getTestData(kfold::KFolds{T}, fold::Int) :: T where T
    index = findall(r -> r == fold, kfold.index)
    return Persa.sample(kfold.dataset, index)
end

Base.length(kf::KFolds) = kf.k

function splitKFold(y, num_folds, is_shuffle)
  if is_shuffle
    i = shuffle(collect(1:y));
  else
     i = collect(1:y);
  end

  fold_size = round(Int, floor(y/num_folds));
  remainder = y-num_folds*fold_size;
  groups = zeros(Int, y);
  cursor = 1;
  group = 1;

  while cursor<=y
    this_fold_size = group <= remainder ? fold_size+1 : fold_size;
    groups[i[cursor:cursor+this_fold_size-1]] .= group;
    group += 1;
    cursor += this_fold_size;
  end

  return groups;
end


text:
"""
    LatLonAccel
Acceleration in the frenet frame
"""
immutable LatLonAccel
    a_lat::Float64
    a_lon::Float64
end
Base.show(io::IO, a::LatLonAccel) = @printf(io, "LatLonAccel(%6.3f, %6.3f)", a.a_lat, a.a_lon)
Base.length(::Type{LatLonAccel}) = 2
Base.convert(::Type{LatLonAccel}, v::Vector{Float64}) = LatLonAccel(v[1], v[2])
function Base.copy!(v::Vector{Float64}, a::LatLonAccel)
    v[1] = a.a_lat
    v[2] = a.a_lon
    v
end
function propagate{D<:Union{VehicleDef, BicycleModel}}(veh::Entity{VehicleState, D, Int}, action::LatLonAccel, roadway::Roadway, ΔT::Float64)

    a_lat = action.a_lat
    a_lon = action.a_lon

     v = veh.state.v
     ϕ = veh.state.posF.ϕ
    ds = v*cos(ϕ)
     t = veh.state.posF.t
    dt = v*sin(ϕ)

    ΔT² = ΔT*ΔT
    Δs = ds*ΔT + 0.5*a_lon*ΔT²
    Δt = dt*ΔT + 0.5*a_lat*ΔT²

    ds₂ = ds + a_lon*ΔT
    dt₂ = dt + a_lat*ΔT
    speed₂ = sqrt(dt₂*dt₂ + ds₂*ds₂)
    v₂ = sqrt(dt₂*dt₂ + ds₂*ds₂) # v is the magnitude of the velocity vector
    ϕ₂ = atan2(dt₂, ds₂)

    roadind = move_along(veh.state.posF.roadind, roadway, Δs)
    footpoint = roadway[roadind]

    posG = convert(VecE2, footpoint.pos) + polar(t + Δt, footpoint.pos.θ + π/2)

    posG = VecSE2(posG.x, posG.y, footpoint.pos.θ + ϕ₂)

    return VehicleState(posG, roadway, v₂)
end
function Base.get(::Type{LatLonAccel}, rec::SceneRecord, roadway::Roadway, vehicle_index::Int, pastframe::Int=0)
    accel_lat = get(ACCFT, rec, roadway, vehicle_index, pastframe)
    accel_lon = get(ACCFS, rec, roadway, vehicle_index, pastframe)
    LatLonAccel(accel_lat, accel_lon)
end


text:
## TRAITS FOR MEASURES

is_measure_type(::Any) = false

const MEASURE_TRAITS =
    [:name, :target_scitype, :supports_weights, :prediction_type, :orientation,
     :reports_each_observation, :aggregation, :is_feature_dependent, :docstring,
     :distribution_type]

# already defined in model_traits.jl:
# name              - fallback for non-MLJType is string(M) where M is arg
# target_scitype    - fallback value = Unknown
# supports_weights  - fallback value = false
# prediction_type   - fallback value = :unknown (also: :deterministic,
#                                           :probabilistic, :interval)
# docstring         - fallback value is value of `name` trait.

# specfic to measures:
orientation(::Type) = :loss  # other options are :score, :other
reports_each_observation(::Type) = false
aggregation(::Type) = Mean()  # other option is Sum() or callable object
is_feature_dependent(::Type) = false

# extend to instances:
orientation(m) = orientation(typeof(m))
reports_each_observation(m) = reports_each_observation(typeof(m))
aggregation(m) = aggregation(typeof(m))
is_feature_dependent(m) = is_feature_dependent(typeof(m))

# specific to probabilistic measures:
distribution_type(::Type) = missing

## AGGREGATION

abstract type AggregationMode end

struct Sum <: AggregationMode end
(::Sum)(v) = sum(v)

struct Mean <: AggregationMode end
(::Mean)(v) = mean(v)

# for rms and it's cousins:
struct RootMeanSquare <: AggregationMode end
(::RootMeanSquare)(v) = sqrt(mean(v.^2))

aggregate(v, measure) = aggregation(measure)(v)

# aggregation is no-op on scalars:
const MeasureValue = Union{Real,Tuple{<:Real,<:Real}} # number or interval
aggregate(x::MeasureValue, measure) = x

## DISPATCH FOR EVALUATION

# yhat - predictions (point or probabilisitic)
# X - features
# y - target observations
# w - per-observation weights

value(measure, yhat, X, y, w) = value(measure, yhat, X, y, w,
                                      Val(is_feature_dependent(measure)),
                                      Val(supports_weights(measure)))

## DEFAULT EVALUATION INTERFACE

#  is feature independent, weights not supported:
value(measure, yhat, X, y, w, ::Val{false}, ::Val{false}) = measure(yhat, y)

#  is feature dependent:, weights not supported:
value(measure, yhat, X, y, w, ::Val{true}, ::Val{false}) = measure(yhat, X, y)

#  is feature independent, weights supported:
value(measure, yhat, X, y, w, ::Val{false}, ::Val{true}) = measure(yhat, y, w)
value(measure, yhat, X, y, ::Nothing, ::Val{false}, ::Val{true}) = measure(yhat, y)

#  is feature dependent, weights supported:
value(measure, yhat, X, y, w, ::Val{true}, ::Val{true}) = measure(yhat, X, y, w)
value(measure, yhat, X, y, ::Nothing, ::Val{true}, ::Val{true}) = measure(yhat, X, y)

## HELPERS

"""
    check_dimension(ŷ, y)

Check that two vectors have compatible dimensions
"""
function check_dimensions(ŷ::AbstractVector, y::AbstractVector)
    length(y) == length(ŷ) ||
        throw(DimensionMismatch("Differing numbers of observations and "*
                                "predictions. "))
    return nothing
end

function check_pools(ŷ, y)
    levels(y) == levels(ŷ[1]) ||
        error("Conflicting categorical pools found "*
              "in observations and predictions. ")
    return nothing
end

## FOR BUILT-IN MEASURES

abstract type Measure <: MLJType end
is_measure_type(::Type{<:Measure}) = true
is_measure(m) = is_measure_type(typeof(m))

## DISPLAY AND INFO

Base.show(stream::IO, ::MIME"text/plain", m::Measure) =
    print(stream, "$(name(m)) (callable Measure)")
Base.show(stream::IO, m::Measure) = print(stream, name(m))

function MLJBase.info(M, ::Val{:measure_type})
    values = Tuple(@eval($trait($M)) for trait in MEASURE_TRAITS)
    return NamedTuple{Tuple(MEASURE_TRAITS)}(values)
end

MLJBase.info(m, ::Val{:measure}) = info(typeof(m))

## INCLUDE SPECIFIC MEASURES AND TOOLS

include("continuous.jl")
include("confusion_matrix.jl")
include("finite.jl")
include("loss_functions_interface.jl")

## DEFAULT MEASURES

default_measure(model::M) where M<:Supervised =
    default_measure(model, target_scitype(M))
default_measure(model, ::Any) = nothing
default_measure(model::Deterministic,
                ::Type{<:Union{AbstractVector{Continuous},
                               AbstractVector{Count}}}) = rms
# default_measure(model::Probabilistic,
#                 ::Type{<:Union{AbstractVector{Continuous},
#                                AbstractVector{Count}}}) = rms
default_measure(model::Deterministic,
                ::Type{<:AbstractVector{<:Finite}}) =
                    misclassification_rate
default_measure(model::Probabilistic,
                ::Type{<:AbstractVector{<:Finite}}) =
                    cross_entropy


text:
#!/usr/bin/env julia

home = homedir()
namespace = "pandola"

find_file(x) = begin
    (root, dirs, files) = x
    map(a -> "$root/$a", files)
end

ff = map(find_file, walkdir("$home/$namespace/bin/testdir"))

#strArr = reduce(append!, [], ff) |> println

strArr = map(x -> open(readall, x), reduce(append!, [], ff))

re = r"class=\"?\'?([\w-\s]+)\"?\'?"

pick_class(re, str) = begin
    next_class(rm, acc = []) = begin
        if rm == nothing
            return acc;
        else
            next_class(match(re, str, rm.offset + 1), push!(acc, rm.match))
        end
    end

    next_class(match(re, str, 1))
end
match_class(x) = match(re, x).captures[1]
split_class(x) = split(x, r"\s")
find_class(x) = x |> match_class |> split_class
reduce_class(acc, x) = begin
    fill_acc(acc, x) = begin
        if x == "" || x ∈ acc
            print("OK, find -> \"$x\". But existed. ")
            print("\e[0;33;49mSkip\e[0m\n")
            acc
        else
            print("OK, find -> \"$x\". ")
            print("\e[0;32;49mThunk\e[0m\n")
            push!(acc, x)
        end
    end
    foreach(y -> fill_acc(acc, y), x)
    acc
end

res = pick_class(re, reduce(string, strArr))
res1 = reduce(reduce_class, [], map(find_class, res)) |> sort #|> println
println("\nFound $(length(res1)) roles.")

map(x -> println(x), res1)

css = open(readall, "index.css")
css_content = "a-z0-9\"\'"
css_symbol = "#\\(\\)\\!π\\.%-:,;"
css_comment = "\\*\\/"
css_space = "\\n\\t\\s"
cons_re(classname) = begin
    role = "\\.$classname[:a-z,\\s]*.*\\s*"
    content = "{[$css_content$css_symbol$css_comment$css_space]*}"
    Regex("$role$content", "i")
end
#res2 = map(x -> cons_re(x), res1)
res2 = map(x -> match(cons_re(x), css, 1), res1)
println("\n")
map(x -> println(x), res2)

concat_style(acc, curr) = begin
    if curr != nothing
        acc * curr.match * "\n\n"
    else
        acc
    end
end
reduce(concat_style, "", res2) |> println

max_stars_repo_name: yuffiy/pandola

max_issues_repo_name: HairyRabbit/pandola

max_forks_repo_name: HairyRabbit/pandola


text:
using LinearAlgebra, StatsBase

# Input data
data = map(
    x -> map(
        rawpoint -> (raw = split(rawpoint, ","); Tuple(map(v -> parse(Int64, v), raw))),
        split(x, " -> "),
    ),
    readlines("Day05.input"),
);

# Part 1
allpoints =
    (line) -> (filter(
        p -> (norm(line[1] .- p) + norm(line[2] .- p) ≈ norm(line[1] .- line[2])),
        vec(
            collect(
                Iterators.product(
                    range(extrema(first.(line))...),
                    range(extrema(last.(line))...),
                ),
            ),
        ),
    ));

getans = (data) -> count(>(1), values(countmap(Iterators.flatten(map(allpoints, data)))));

ans = getans(filter((line) -> any(line[1] .== line[2]), data));
println("Part 1: $ans")

# Part 2
ans = getans(data);
println("Part 2: $ans")


text:
import Base: atan2, cross, dot, normalize
export norm2, perp, project, projectN, reflect, reflectN, rotate
export projectNextrema, voronoi_region, overlapping, ininterval, minmaxV, wrap1

norm2(v::SVector{2}) = dot(v,v)
perp(v::SVector{2}) = SVector(v[2],-v[1])
cross(v1::SVector{2}, v2::SVector{2}) = v1[1]*v2[2] - v1[2]*v2[1]
project(v1::SVector{2}, v2::SVector{2}) = (dot(v1,v2) / norm2(v2)) * v2
projectN(v1::SVector{2}, v2::SVector{2}) = dot(v1,v2) * v2
reflect(v1::SVector{2}, ax::SVector{2}) = v1 - 2*project(v1,ax)
reflectN(v1::SVector{2}, ax::SVector{2}) = v1 - 2*projectN(v1,ax)
function rotate(v::SVector{2}, θ)
    c = cos(θ)
    s = sin(θ)
    SVector(v[1]*c - v[2]*s, v[1]*s + v[2]*c)
end
atan2(v::SVector{2}) = atan2(v[2], v[1])
function projectNextrema{T}(pts::AbstractVector{SVector{2,T}}, n::SVector{2,T})
    dmin = T(Inf)
    dmax = -T(Inf)
    for p in pts
        d = dot(p,n)
        d < dmin && (dmin = d)
        d > dmax && (dmax = d)
    end
    SVector(dmin, dmax)
end
function voronoi_region(p::SVector{2}, v::SVector{2})
    d = dot(p,v)
    d < 0 && return -1
    d > norm2(v) && return 1
    return 0
end
overlapping(I1::SVector{2}, I2::SVector{2}) = I1[1] <= I2[2] && I2[1] <= I1[2]
ininterval(x, I::SVector{2}) = I[1] <= x <= I[2]
minmaxV(x, y) = x < y ? SVector(x,y) : SVector(y,x)  # Base.minmax returns tuple
wrap1(i, N) = i < 1 ? N : i > N ? 1 : i      # marginally faster than mod1


text:
####
#### PCA
####

"""
    PCA(; kwargs...)

$PCA_DESCR

# Keyword Parameters

- `maxoutdim::Int=0`: maximum number of output dimensions, uses the smallest dimension of
    training feature matrix if 0 (default).
- `method::Symbol=:auto`: method to use to solve the problem, one of `:auto`,`:cov`
    or `:svd`
- `pratio::Float64=0.99`: ratio of variance preserved
- `mean::Union{Nothing, Real, Vector{Float64}}=nothing`: if set to nothing(default)
    centering will be computed and applied, if set to `0` no
    centering(assumed pre-centered), if a vector is passed, the centering is done with
    that vector.
"""
@mlj_model mutable struct PCA <: MMI.Unsupervised
    maxoutdim::Int = 0::(_ ≥ 0)
    method::Symbol = :auto::(_ in (:auto, :cov, :svd))
    pratio::Float64 = 0.99::(0.0 < _ ≤ 1.0)
    mean::Union{Nothing, Real, Vector{Float64}} = nothing::(_check_typeof_mean(_))
end

function _check_typeof_mean(x)
    return x isa Vector{<:Real} || x === nothing || (x isa Real && iszero(x))
end

function MMI.fit(model::PCA, verbosity::Int, X)
    Xarray = MMI.matrix(X)
    mindim = minimum(size(Xarray))
    maxoutdim = model.maxoutdim == 0 ? mindim : model.maxoutdim
    # NOTE: copy/transpose
    fitresult = MS.fit(
        MS.PCA, transpose(Xarray);
        method=model.method,
        pratio=model.pratio,
        maxoutdim=maxoutdim,
        mean=model.mean
    )
    cache = nothing
    report = (
        indim=MS.indim(fitresult),
        outdim=MS.outdim(fitresult),
        tprincipalvar=MS.tprincipalvar(fitresult),
        tresidualvar=MS.tresidualvar(fitresult),
        tvar=MS.tvar(fitresult),
        mean=copy(MS.mean(fitresult)),
        principalvars=copy(MS.principalvars(fitresult))
    )
    return fitresult, cache, report
end

metadata_model(PCA,
    input=Table(Continuous),
    output=Table(Continuous),
    weights=false,
    descr=PCA_DESCR,
    path="$(PKG).PCA"
)

####
#### KernelPCA
####

"""
    KernelPCA(; kwargs...)

$KPCA_DESCR

# Keyword Parameters

- `maxoutdim::Int = 0`: maximum number of output dimensions, uses the smallest
    dimension of training feature matrix if 0 (default).
- `kernel::Function=(x,y)->x'y`: kernel function of 2 vector arguments x and y, returns a
    scalar value
- `solver::Symbol=:auto`: solver to use for the eigenvalues, one of `:eig`(default),
    `:eigs`
- `inverse::Bool=false`: perform calculation for inverse transform
- `beta::Real=1.0`: strength of the ridge regression that learns the inverse transform
    when inverse is true
- `tol::Real=0.0`: Convergence tolerance for eigs solver
- `maxiter::Int=300`: maximum number of iterations for eigs solver
"""
@mlj_model mutable struct KernelPCA <: MMI.Unsupervised
    maxoutdim::Int = 0::(_ ≥ 0)
    kernel::Union{Nothing, Function} = default_kernel
    solver::Symbol = :eig::(_ in (:eig, :eigs))
    inverse::Bool = false
    beta::Real = 1.0::(_ ≥ 0.0)
    tol::Real = 1e-6::(_ ≥ 0.0)
    maxiter::Int = 300::(_ ≥ 1)
end

function MMI.fit(model::KernelPCA, verbosity::Int, X)
    Xarray = MMI.matrix(X)
    mindim = minimum(size(Xarray))
    # default max out dim if not given
    maxoutdim = model.maxoutdim == 0 ? mindim : model.maxoutdim
    fitresult = MS.fit(
        MS.KernelPCA,
        permutedims(Xarray);
        kernel=model.kernel,
        maxoutdim=maxoutdim,
        solver=model.solver,
        inverse=model.inverse,
        β=model.beta,tol=model.tol,
        maxiter=model.maxiter
    )
    cache  = nothing
    report = (
        indim=MS.indim(fitresult),
        outdim=MS.outdim(fitresult),
        principalvars=copy(MS.principalvars(fitresult))
    )
    return fitresult, cache, report
end

metadata_model(
    KernelPCA,
    input=Table(Continuous),
    output=Table(Continuous),
    weights=false,
    descr=KPCA_DESCR,
    path="$(PKG).KernelPCA"
)

####
#### ICA
####

"""
    ICA(; kwargs...)

$ICA_DESCR

# Keyword Parameters

- `k::Int=0`: number of independent components to recover, set automatically if `0`
- `alg::Symbol=:fastica`: algorithm to use (only `:fastica` is supported at the moment)
- `fun::Symbol=:tanh`: approximate neg-entropy functor, via the function
    `MultivariateStats.icagfun`, one of `:tanh` and `:gaus`
- `do_whiten::Bool=true`: whether to perform pre-whitening
- `maxiter::Int=100`: maximum number of iterations
- `tol::Real=1e-6`: convergence tolerance for change in matrix W
- `mean::Union{Nothing, Real, Vector{Float64}}=nothing`: mean to use, if nothing (default)
    centering is computed andapplied, if zero, no centering, a vector of means can
    be passed
- `winit::Union{Nothing,Matrix{<:Real}}=nothing`: initial guess for matrix `W` either
    an empty matrix (random initilization of `W`), a matrix of size `k × k` (if `do_whiten`
    is true), a matrix of size `m × k` otherwise. If unspecified i.e `nothing` an empty
    `Matrix{<:Real}` is used.
"""
@mlj_model mutable struct ICA <: MMI.Unsupervised
    k::Int = 0::(_ ≥ 0)
    alg::Symbol = :fastica::(_ in (:fastica,))
    fun::Symbol = :tanh::(_ in (:tanh, :gaus))
    do_whiten::Bool = true
    maxiter::Int=100::(_ ≥ 1)
    tol::Real = 1e-6::(_ ≥ 0.0)
    winit::Union{Nothing, Matrix{<:Real}}= nothing
    mean::Union{Nothing, Real, Vector{Float64}} = nothing::(_check_typeof_mean(_))
end

function MMI.fit(model::ICA, verbosity::Int, X)
    Xarray = MMI.matrix(X)
    n, p = size(Xarray)
    m = min(n, p)
    k = ifelse(model.k ≤ m, model.k, m)
    fitresult = MS.fit(
        MS.ICA, transpose(Xarray), k;
        alg=model.alg,
        fun=MS.icagfun(model.fun, eltype(Xarray)),
        do_whiten=model.do_whiten,
        maxiter=model.maxiter,
        tol=model.tol,
        mean=model.mean,
        winit=model.winit === nothing ? zeros(eltype(Xarray), 0, 0) : model.winit
    )
    cache = nothing
    report = (
        indim=MS.indim(fitresult),
        outdim=MS.outdim(fitresult),
        mean=copy(MS.mean(fitresult))
    )
    return fitresult, cache, report
end

metadata_model(
    ICA,
    input=Table(Continuous),
    output=Table(Continuous),
    weights=false,
    descr=ICA_DESCR,
    path="$(PKG).ICA"
)

####
#### PPCA
####

"""
    PPCA(; kwargs...)

$PPCA_DESCR

# Keyword Parameters

- `maxoutdim::Int=0`: maximum number of output dimensions, uses max(no_of_features - 1, 1)
    if 0 (default).
- `method::Symbol=:ml`: method to use to solve the problem, one of `:ml`, `:em`, `:bayes`.
- `maxiter::Int=1000`: maximum number of iterations.
- `tol::Real=1e-6`: convergence tolerance.
- `mean::Union{Nothing, Real, Vector{Float64}}=nothing`: if set to nothing(default)
    centering will be computed and applied, if set to `0` no
    centering(assumed pre-centered), if a vector is passed, the centering is done with
    that vector.
"""
@mlj_model mutable struct PPCA <: MMI.Unsupervised
    maxoutdim::Int = 0::(_ ≥ 0)
    method::Symbol = :ml::(_ in (:ml, :em, :bayes))
    maxiter::Int = 1000
    tol::Real = 1e-6::(_ ≥ 0.0)
    mean::Union{Nothing, Real, Vector{Float64}} = nothing::(_check_typeof_mean(_))
end

function MMI.fit(model::PPCA, verbosity::Int, X)
    Xarray = MMI.matrix(X)
    def_dim = max(1, size(Xarray, 2) - 1)
    maxoutdim = model.maxoutdim == 0 ? def_dim : model.maxoutdim
    # NOTE: copy/transpose
    fitresult = MS.fit(
        MS.PPCA, transpose(Xarray);
        method=model.method,
        tol=model.tol,
        maxiter=model.maxiter,
        maxoutdim=maxoutdim,
        mean=model.mean
    )
    cache = nothing
    report = (
        indim=MS.indim(fitresult),
        outdim=MS.outdim(fitresult),
        tvar=MS.var(fitresult),
        mean=copy(MS.mean(fitresult)),
        loadings=MS.loadings(fitresult)
    )
    return fitresult, cache, report
end

metadata_model(PPCA,
    input=Table(Continuous),
    output=Table(Continuous),
    weights=false,
    descr=PPCA_DESCR,
    path="$(PKG).PPCA"
)

####
#### FactorAnalysis
####

"""
    FactorAnalysis(; kwargs...)

$PPCA_DESCR

# Keyword Parameters

- `method::Symbol=:cm`: Method to use to solve the problem, one of `:ml`, `:em`, `:bayes`.
- `maxoutdim::Int=0`: Maximum number of output dimensions, uses max(no_of_features - 1, 1)
    if 0 (default).
- `maxiter::Int=1000`: Maximum number of iterations.
- `tol::Real=1e-6`: Convergence tolerance.
- `eta::Real=tol`: Variance lower bound
- `mean::Union{Nothing, Real, Vector{Float64}}=nothing`: If set to nothing(default)
    centering will be computed and applied, if set to `0` no
    centering(assumed pre-centered), if a vector is passed, the centering is done with
    that vector.
"""
@mlj_model mutable struct FactorAnalysis <: MMI.Unsupervised
    method::Symbol=:cm::(_ in (:em, :cm))
    maxoutdim::Int=0::(_ ≥ 0)
    maxiter::Int=1000::(_ ≥ 1)
    tol::Real=1e-6::(_ ≥ 0.0)
    eta::Real=tol::(_ ≥ 0.0)
    mean::Union{Nothing, Real, Vector{Float64}} = nothing::(_check_typeof_mean(_))
end

function MMI.fit(model::FactorAnalysis, verbosity::Int, X)
    Xarray = MMI.matrix(X)
    def_dim = max(1, size(Xarray, 2) - 1)
    maxoutdim = model.maxoutdim == 0 ? def_dim : model.maxoutdim
    # NOTE: copy/transpose
    fitresult = MS.fit(
        MS.FactorAnalysis, transpose(Xarray);
        method=model.method,
        maxiter=model.maxiter,
        tol=model.tol,
        η=model.eta,
        maxoutdim=maxoutdim,
        mean=model.mean
    )
    cache = nothing
    report = (
        indim=MS.indim(fitresult),
        outdim=MS.outdim(fitresult),
        variance=MS.var(fitresult),
        covariance_matrix=MS.cov(fitresult),
        mean=MS.mean(fitresult),
        loadings=MS.loadings(fitresult)
    )
    return fitresult, cache, report
end

metadata_model(FactorAnalysis,
    input=Table(Continuous),
    output=Table(Continuous),
    weights=false,
    descr=FactorAnalysis_DESCR,
    path="$(PKG).FactorAnalysis"
)

####
#### Common interface
####

model_types = [
    (PCA, PCAFitResultType),
    (KernelPCA, KernelPCAFitResultType),
    (ICA, ICAFitResultType),
    (PPCA, PPCAFitResultType),
    (FactorAnalysis, FactorAnalysisResultType)
]

for (M, MFitResultType) in model_types
    @eval function MMI.fitted_params(::$M, fr)
        return (projection=copy(MS.projection(fr)),)
    end

    @eval function MMI.transform(::$M, fr::$MFitResultType, X)
        # X is n x d, need to transpose twice
        Xarray = MMI.matrix(X)
        Xnew = transpose(MS.transform(fr, transpose(Xarray)))
        return MMI.table(Xnew, prototype=X)
    end

    if hasmethod(MS.reconstruct, Tuple{MFitResultType{Float64}, Matrix{Float64}})
        @eval function MMI.inverse_transform(::$M, fr::$MFitResultType, Y)
            # X is n x p, need to transpose twice
            Yarray = MMI.matrix(Y)
            Ynew = transpose(MS.reconstruct(fr, transpose(Yarray)))
            return MMI.table(Ynew, prototype=Y)
        end
    end
end