text:
################################################################################
#
#  Ramification at infinite places
#
################################################################################

function _real_level(L::GAP.GapObj)
  #I compute the index of the first derived subgroup containing all the elements of order 2
  #First, I compute the subgroup generated by them.
  lElem = GAP.GapObj[]
  G = L[1]
  ElemsG = GAP.Globals.Elements(G)
  for i = 1:length(ElemsG)
    g = ElemsG[i]
    if GAP.Globals.Order(g) == 2
      push!(lElem, g)
    end
  end
  S = GAP.Globals.Subgroup(G, GAP.julia_to_gap(lElem))
  #Now, I check containment.
  k = 0
  for i = 2:length(L)
    k += 1
    if !GAP.Globals.IsSubgroup(L[i], S)
      break
    end
  end
  return k
end

###############################################################################
#
#  Abelian extensions of QQ
#
###############################################################################

function abelian_extensionsQQ(gtype::Vector{Int}, bound::fmpz, only_real::Bool = false; unramified_outside::Vector{fmpz} = fmpz[])

  gtype = map(Int, snf(abelian_group(gtype))[1].snf)
  #Easy case: quadratic and biquadratic extensions
  if gtype == [2]
    lq = Hecke._quad_ext(Int(bound), only_real, unramified_outside = unramified_outside)
    res = Vector{FieldsTower}(undef, length(lq))
    for i = 1:length(lq)
      x = lq[i]
      E = EquationOrder(x[1])
      E.ismaximal = 1
      E.index = fmpz(1)
      E.gen_index = fmpq(1)
      E.disc = discriminant(E)
      set_attribute!(x[1], :maximal_order => E)
      auts = Vector{NfToNfMor}(undef, 2)
      auts[1] = hom(x[1], x[1], gen(x[1]), check = false)
      auts[2] = x[2][1]
      set_automorphisms(x[1], auts)
      res[i] = FieldsTower(x[1], x[2], x[3])
      res[i].isabelian = true
    end
    return res
  end
  if gtype == Int[2,2]
    l = Hecke._C22_exts_abexts(Int(bound), only_real, unramified_outside = unramified_outside)
    @vprint :Fields 1 "Computing maximal orders\n"
    @vprint :FieldsNonFancy 1 "Computing maximal orders\n"
    lf = Hecke._C22_with_max_ord(l)
    res1 = Vector{FieldsTower}(undef, length(lf))
    for i = 1:length(lf)
      res1[i] = FieldsTower(lf[i][1], lf[i][2], lf[i][3])
      res1[i].isabelian = true
    end
    return res1
  end
  l1 = _abelian_extensionsQQ(gtype, bound, only_real, unramified_outside = unramified_outside)
  list1 = Vector{FieldsTower}(undef, length(l1))
  @vprint :Fields 1 "Computing maximal orders\n\n"
  @vprint :FieldsNonFancy 1 "Computing maximal orders\n"
  for i = 1:length(l1)
    @vprint :Fields 1 "\e[1FComputing maximal order $(i) /$(length(l1)) \n"
    x = l1[i]
    K, auts = _relative_to_absoluteQQ(x[1], x[2])
    if length(gtype) == 1
      #If the group is cyclic, I prefer to have a generator!
      new_auts = Vector{NfToNfMor}(undef, 1)
      new_auts[1] = auts[1]
      for j = 2:length(auts)
        new_auts[1] *= auts[j]
      end
      list1[i] = FieldsTower(K, new_auts, [hom(base_field(x[1]), K, K(1), check = false)])
    else
      list1[i] = FieldsTower(K, auts, [hom(base_field(x[1]), K, K(1), check = false)])
    end
    list1[i].isabelian = true
  end
  @vprint :Fields 1 "\e[1F$(Hecke.set_cursor_col())$(Hecke.clear_to_eol())"
  return list1

end

function _abelian_extensionsQQ(gtype::Vector{Int}, absolute_discriminant_bound::fmpz, only_real::Bool = false; unramified_outside::Vector{fmpz} = fmpz[])

  Qx, x = PolynomialRing(FlintQQ, "x", cached = false)
  K, _ = NumberField(x-1, "a", cached = false)
  O = maximal_order(K)
  n = prod(gtype)
  expo = lcm(gtype)

  #Getting conductors
  l_conductors = Hecke.conductorsQQ(O, gtype, absolute_discriminant_bound; unramified_outside = unramified_outside)
  sort!(l_conductors, rev = true)
  len = length(l_conductors)
  @vprint :Fields 1 "Number of conductors: $(len) \n\n"
  @vprint :FieldsNonFancy 1 "Number of conductors: $(len) \n"

  complex = iseven(expo) && !only_real

  #Now, the big loop
  class_fields = Vector{Hecke.ClassField{MapRayClassGrp, GrpAbFinGenMap}}()
  for (i, k) in enumerate(l_conductors)
    if iszero(mod(i, 1000))
      pt = len - i
      @vprint :Fields 1 "\e[1F$(Hecke.set_cursor_col())$(Hecke.clear_to_eol()) Conductors to test: $(pt)\n"
    end
    r, mr = Hecke.ray_class_groupQQ(O, k, complex, expo)
    if !has_quotient(r, gtype)
      continue
    end
    ls = subgroups(r, quotype = gtype, fun = (x, y) -> quo(x, y, false)[2])
    for s in ls
      s::GrpAbFinGenMap
      C = ray_class_field(mr, s)::ClassField{MapRayClassGrp, GrpAbFinGenMap}
      if Hecke._is_conductor_minQQ(C, n) && Hecke.discriminant_conductorQQ(O, C, k, absolute_discriminant_bound)
        push!(class_fields, C)
      end
    end
  end
  fields = Vector{Tuple{Hecke.NfRelNS{nf_elem}, Vector{Hecke.NfRelNSToNfRelNSMor_nf_elem}}}(undef, length(class_fields))
  for i = 1:length(class_fields)
    @vprint :Fields 1 "\e[1FComputing class field $(i) /$(length(class_fields)) \n"
    C = class_fields[i]
    fields[i] = (number_field(C), Hecke.automorphism_groupQQ(C))
  end
  @vprint :Fields 1 "\e[1F$(Hecke.set_cursor_col())$(Hecke.clear_to_eol())Number of fields found: $(length(fields)) \n"
  @vprint :FieldsNonFancy 1 "Number of fields found: $(length(fields)) \n"
  return fields

end

###############################################################################
#
#  General abelian layer
#
###############################################################################

function check_extension(C::Hecke.ClassField, bound::fmpz, Dcond::Dict, Ddisc::Dict)
  @vtime :Fields 3 fl2 = Hecke._is_conductor_min_normal(C, lwp = Dcond)
  if !fl2
    return false
  end
  @vtime :Fields 3 fl3 = Hecke.discriminant_conductor(C, bound, lwp = Ddisc)
  return fl3

end

function _construct_grp(IdH::GAP.GapObj, uncom::Int)
  G = GAP.Globals.SmallGroup(IdH)
  L = GAP.Globals.DerivedSeries(G)
  mH1 = GAP.Globals.NaturalHomomorphismByNormalSubgroup(G, L[end-1])
  K = GAP.Globals.Kernel(mH1)
  gens = GAP.Globals.MinimalGeneratingSet(K)
  new_gens = []
  for i = 1:length(gens)
    o = GAP.Globals.Order(gens[i])
    ex = ppio(o, uncom)[1]
    push!(new_gens, gens[i]^ex)
  end
  S = GAP.Globals.Subgroup(G, GAP.julia_to_gap(new_gens))
  Q = GAP.Globals.FactorGroup(G, S)
  IdCheck = GAP.Globals.IdGroup(Q)
  return IdCheck
end

function max_ramified_prime(O::NfOrd, gtype::Vector{Int}, bound::fmpz)
  n = prod(gtype)
  fac = factor(n)
  m = Int(minimum(keys(fac.fac)))
  k = divexact(n, m)
  b1 = Int(root(bound, degree(O)*(m-1)*k))
  return b1
end

function _abelian_normal_extensions(F::FieldsTower, gtype::Vector{Int}, absbound::fmpz, IdCheck::GAP.GapObj, only_real::Bool, IdG::GAP.GapObj; unramified_outside::Vector{fmpz} = fmpz[])
  K = F.field
  O = maximal_order(K)
  n = prod(gtype)
  inf_plc = Vector{InfPlc}()
  if abs(discriminant(O))^n > absbound
    return Vector{Hecke.ClassField{Hecke.MapRayClassGrp, GrpAbFinGenMap}}[]
  end
  bound = div(absbound, abs(discriminant(O))^n)
  @vprint :Fields 2 "\n"
  @vtime :Fields 2 l_conductors = conductors_with_restrictions(F, gtype, IdG, bound, unramified_outside = unramified_outside)
  @vprint :Fields 1 "   Number of conductors: $(length(l_conductors))\n"
  @vprint :FieldsNonFancy 1 "Number of conductors: $(length(l_conductors))\n"
  if length(l_conductors) == 0
    @vprint :Fields 1 "\n"
    return Vector{Hecke.ClassField{Hecke.MapRayClassGrp, GrpAbFinGenMap}}[]
  end
  @vprint :Fields 2 "\n"
  @vprint :Fields 1 "Computing class group of $(K.pol)"
  @vprint :Fields 2 "\n"
  @vtime :Fields 2 Cl, mCl = class_group(O, use_aut = true)
  @vprint :Fields 1 "$(Hecke.set_cursor_col())$(Hecke.clear_to_eol())"
  if mod(n, 2) == 0 && !only_real
    inf_plc = real_places(K)
  end
  expo = lcm(gtype)
  if length(l_conductors) == 1 && isone(l_conductors[1][1]) && isempty(l_conductors[1][2]) && !divisible(order(Cl)* (2^length(inf_plc)), fmpz(n))
    @vprint :Fields 1 "\n"
    return Vector{Hecke.ClassField{Hecke.MapRayClassGrp, GrpAbFinGenMap}}[]
  end
  Hecke.allow_cache!(mCl)
  @vtime :Fields 3 rcg_ctx = Hecke.rayclassgrp_ctx(O, expo)
  @vprint :Fields 1 "\n"
  j = -1
  #first_group = true
  autos = F.generators_of_automorphisms
  class_fields_with_act = Tuple{ClassField{MapRayClassGrp, GrpAbFinGenMap}, Vector{GrpAbFinGenMap}}[]
  fun_sub = (x, y) -> quo(x, y, false)[2]
  for k in l_conductors
    j += 1
    @vprint :Fields 1 "\e[1FConductors to test: $(length(l_conductors)-j) \n"
    @vprint :Fields 3 "\n\n"
    @vtime :Fields 3 r, mr = Hecke.ray_class_group_quo(O, k[1], k[2], inf_plc, rcg_ctx)
    S, mS = snf(r)
    if !Hecke.has_quotient(S, gtype)
      continue
    end
    #=
    if first_group
      B1 = maximum(fmpz[next_prime(max_ramified_prime(O, gtype, bound)), 211])
      rfF = ramified_primes(F)
      B1 = max(B1, maximum(rfF))+1
      lP = Hecke.find_gens(pseudo_inv(rcg_ctx.class_group_map), PrimesSet(B1, fmpz(-1)))[1]
      rcg_ctx.class_group_map.small_gens = lP
      first_group = false
    end
    =#
    #mr.clgrpmap.small_gens = rcg_ctx.class_group_map.small_gens
    @vtime :Fields 3 act = Hecke.induce_action(mr, autos)
    act_S = induce_action_on_subgroup(mS, act)
    imS = inv(mS)
    @vtime :Fields 3 ls = stable_subgroups(S, act_S, op = fun_sub, quotype = gtype)
    Dcond = Dict{Int, Vector{GrpAbFinGenElem}}()
    Ddisc = Dict{Tuple{Int, Int}, Vector{GrpAbFinGenElem}}()
    for s in ls
      s::GrpAbFinGenMap
      @hassert :Fields 1 order(codomain(s)) == n
      codomain(s).exponent = expo
      @vtime :Fields 4 C = ray_class_field(mr, imS*s)::ClassField{MapRayClassGrp, GrpAbFinGenMap}
      fl = check_extension(C, bound, Dcond, Ddisc)
      if fl
        res_act = _action(s, act_S)
        @vtime :Fields 3 fl1 = check_group_extension(IdCheck, autos, res_act)
        if fl1
          push!(class_fields_with_act, (C, res_act))
        end
      end
    end
    @vprint :Fields 3 "\n\n"
  end
  if isempty(class_fields_with_act)
    return Vector{Hecke.ClassField{Hecke.MapRayClassGrp, GrpAbFinGenMap}}[]
  end
  @vprint :Fields 1 "\e[1F$(Hecke.clear_to_eol())Sieving $(length(class_fields_with_act)) abelian extensions"
  candidates = check_abelian_extensions(class_fields_with_act, F, IdG)
  @vprint :Fields 1 "$(Hecke.set_cursor_col())$(Hecke.clear_to_eol())\n"
  t_candidates = findall(candidates)
  res_cfields = ClassField{MapRayClassGrp, GrpAbFinGenMap}[class_fields_with_act[t_candidates[i]][1] for i = 1:length(t_candidates)]
  return res_cfields
end

################################################################################
#
#  From the class fields to the computation of the defining polynomials
#
################################################################################

function from_class_fields_to_fields(class_fields::Vector{ClassField{MapRayClassGrp, GrpAbFinGenMap}}, autos::Vector{NfToNfMor}, grp_to_be_checked::Dict{Int, GAP.GapObj}, target_group::GAP.GapObj)

  if isempty(class_fields)
    @vprint :Fields 1 "\e[1F$(Hecke.set_cursor_col())$(Hecke.clear_to_eol())"
    return Tuple{Hecke.NfRelNS{nf_elem}, Vector{Hecke.NfRelNSToNfRelNSMor_nf_elem}}[]
  end
  K = base_ring(class_fields[1])
  divisors_of_n = collect(keys(grp_to_be_checked))
  sort!(divisors_of_n)
  pclassfields = Vector{Vector{ClassField{MapRayClassGrp, GrpAbFinGenMap}}}(undef, length(divisors_of_n))
  right_grp = trues(length(class_fields))
  ind = 1
  for p in divisors_of_n
    it = findall(right_grp)
    if iszero(length(it))
      break
    end
    cfieldsp = Vector{ClassField{MapRayClassGrp, GrpAbFinGenMap}}(undef, length(class_fields))
    for i in it
      cfieldsp[i] = Hecke.maximal_p_subfield(class_fields[i], Int(p))
    end
    idE = grp_to_be_checked[p]
    if iszero(mod(order(torsion_unit_group(maximal_order(K))[1]), p^(valuation(exponent(class_fields[1]), p))))
      #In this case, since we already know the class group of K and the RCF requires K(zeta_p^v) = K,
      #we compute the class field directly over K
      compute_fields(cfieldsp, autos, idE, right_grp)
      pclassfields[ind] = cfieldsp
      ind += 1
      continue
    end
    #Now, we check whether the abelian extensions can be realized over a subfield from a group theoretical point of view.
    E = GAP.Globals.SmallGroup(idE)
    S, H, ab_inv = max_ab_norm_sub_containing(E)
    if S == H
      #No, it can't.
      compute_fields(cfieldsp, autos, idE, right_grp)
      pclassfields[ind] = cfieldsp
      ind += 1
      continue
    end
    @vprint :Fields 1 "$(Hecke.set_cursor_col())$(Hecke.clear_to_eol())"
    # I can compute the fields over a subfield.
    #First, I need the subfields.
    #ab_inv gives the abelian invariants of the maximal abelian subextension of K over the target subfield.
    K = base_field(class_fields[1])
    assure_automorphisms(K, autos)
    subfields = compute_subfields(K, E, H, S)
    computing_over_subfields(cfieldsp, subfields, idE, autos, right_grp, ab_inv, GAP.Globals.IsNormal(E, S))
    pclassfields[ind] = cfieldsp
    ind += 1
  end
  it = findall(right_grp)
  if isempty(it)
    return Vector{Tuple{NfRelNS{nf_elem}, Vector{Hecke.NfRelNSToNfRelNSMor_nf_elem}}}()
  end
  if length(divisors_of_n) == 1 || iscoprime(degree(class_fields[it[1]]), degree(K))
    fields = Vector{Tuple{NfRelNS{nf_elem}, Vector{Hecke.NfRelNSToNfRelNSMor_nf_elem}}}(undef, length(it))
    ind = 1
    for i in it
      res = Vector{typeof(class_fields[1])}(undef, length(divisors_of_n))
      for j = 1:length(divisors_of_n)
        res[j] = pclassfields[j][i]
      end
      fields[ind] = _ext_and_autos(res, autos)
      ind += 1
    end
  else
    #I need to check the isomorphism class of the Galois group
    fields = Vector{Tuple{NfRelNS{nf_elem}, Vector{Hecke.NfRelNSToNfRelNSMor_nf_elem}}}()
    for i in it
      res = Vector{typeof(class_fields[1])}(undef, length(divisors_of_n))
      for j = 1:length(divisors_of_n)
        res[j] = pclassfields[j][i]
      end
      fres = _ext_and_autos(res, autos)
      Cperm = permutation_group(fres[2])
      if GAP.Globals.IdGroup(Cperm) == target_group
        push!(fields, fres)
      end
    end
  end
  return fields

end

function compute_fields(class_fields::Vector{Hecke.ClassField{Hecke.MapRayClassGrp, GrpAbFinGenMap}}, autos::Vector{NfToNfMor}, grp_to_be_checked::GAP.GapObj, right_grp)

  use_brauer = true
  it = findall(right_grp)
  K = base_field(class_fields[it[1]])
  OK = maximal_order(K)
  if divisible(torsion_units_order(K), exponent(class_fields[it[1]]))
    use_brauer = false
  end

  fields = Tuple{Hecke.NfRelNS{nf_elem}, Vector{Hecke.NfRelNSToNfRelNSMor_nf_elem}}[]
  expo = Int(exponent(codomain(class_fields[it[1]].quotientmap)))

  if !use_brauer
    set_up_cycl_ext(K, expo, autos)
  end
  @vprint :Fields 3 "Computing the fields directly\n"
  for i in it
    C = class_fields[i]
    L = number_field(C, using_norm_relation = use_brauer)
    autL = Hecke.absolute_automorphism_group(C, autos)
    if !isone(gcd(degree(K), expo))
      Cpperm = permutation_group(autL)
      if GAP.Globals.IdGroup(Cpperm) != grp_to_be_checked
        right_grp[i] = false
      end
    end
  end
  return right_grp
end

function _ext_and_autos(resul::Vector{Hecke.ClassField{S, T}}, autos::Vector{NfToNfMor}) where S where T
  if length(resul) == 1
    return resul[1].A, resul[1].AbsAutGrpA
  end
  K = domain(autos[1])
  Kx, x = PolynomialRing(K, "x", cached = false)
  pols = typeof(x)[]
  for i = 1:length(resul)
    append!(pols, [to_univariate(Kx, resul[i].A.pol[w]) for w = 1:length(resul[i].A.pol)])
  end
  L, gL = number_field(pols, cached = false, check = false)
  autL = Vector{Hecke.NfRelNSToNfRelNSMor_nf_elem}()
  imgs_auts_base_field = Vector{Vector{Hecke.NfRelNSElem{nf_elem}}}(undef, length(autos))
  for i = 1:length(autos)
    imgs_auts_base_field[i] = gens(L)
  end
  #Now, the automorphisms...
  for i = 1:length(resul)
    Cp = resul[i]
    p = ispower(degree(Cp))[2]
    w = 1
    while mod(degree(pols[w]), p) != 0
      w += 1
    end
    auts = Cp.AbsAutGrpA
    for phi in auts
      if (isdefined(phi.image_data.base_field_map_data, :prim_image) ? phi.image_data.base_field_map_data.prim_image : codomain(phi)(gen(base_field(codomain(phi))))) == codomain(phi)(gen(K))
        imgsphi = gens(L)
        for ind = w:(w+length(Cp.A.pol)-1)
          imgsphi[ind] = evaluate(to_univariate(Kx, image_generators(phi)[ind-w+1].data), imgsphi[ind])
        end
        push!(autL, hom(L, L, imgsphi))
      else
        ind_aut = 1
        #while autos[ind_aut] != phi.coeff_aut
        while !_isequal(base_field(domain(phi)), codomain(phi), __convert_map_data(autos[ind_aut].image_data, codomain(phi)), phi.image_data.base_field_map_data)
          ind_aut += 1
        end
        for ind = w:(w+length(Cp.A.pol)-1)
          imgs_auts_base_field[ind_aut][ind] = Hecke.msubst(image_generators(phi)[ind-w+1].data, gL[w:(w+length(Cp.A.pol)-1)])
        end
      end
    end
  end
  for i = 1:length(imgs_auts_base_field)
    push!(autL, hom(L, L, autos[i], imgs_auts_base_field[i]))
  end
  return L, autL

end

function set_up_cycl_ext(K::AnticNumberField, n::Int, autK::Vector{NfToNfMor})
  fac = factor(n)
  for (p, v) in fac
    e = Int(p)^v
    if e == 2
      continue
    end
    C = cyclotomic_extension(K, e)
    if degree(C.Kr) == 1
      continue
    end
    auts = automorphisms(C, gens = autK, copy = false)
    @vprint :Fields 1 ": computing class group of cyclotomic extension of order $e"
    @vprint :FieldsNonFancy 1 "computing class group of cyclotomic extension of order $e\n"
    Cl, mCl = class_group(maximal_order(C.Ka), use_aut = true)
    Hecke.allow_cache!(mCl)
    @vprint :Fields 1 "$(Hecke.set_cursor_col())$(Hecke.clear_to_eol())"
  end
  return nothing

end

function _action(t::Hecke.GrpAbFinGenMap, act::Vector{Hecke.GrpAbFinGenMap})

  T = codomain(t)
  S, mS = snf(T)
  new_act = Vector{Hecke.GrpAbFinGenMap}(undef, length(act))
  for i = 1:length(act)
    res = mS.map*act[i].map*mS.imap
    new_act[i] = Hecke.GrpAbFinGenMap(S, S, res)
  end
  return new_act

end

################################################################################
#
#  Maximal subgroups which is in direct product with the derived subgroup
#
################################################################################

function max_ab_norm_sub_containing(G::GAP.GapObj)

  D = GAP.Globals.DerivedSeries(G)
  H = D[end-1]
  G1 = GAP.Globals.Centralizer(G, H)
  if G1 == H
    return H, H, Int[]
  end
  #First, I check if the centralizer split directly as a direct product.
  sc = GAP.Globals.ComplementClassesRepresentatives(G1, H)
  if !GAP.Globals.IsEmpty(sc)
    return G1, H, GAP.gap_to_julia(Vector{Int}, GAP.Globals.AbelianInvariants(sc[1]))
  end
  lS = GAP.Globals.ConjugacyClassesSubgroups(G1)
  #TODO: Subgroups in the quotient by H and not in the full group
  candidate = H
  sizecandidate = GAP.Globals.Size(H)
  ab_invariants = GAP.gap_to_julia(Vector{Int}, GAP.Globals.AbelianInvariants(H))
  for i = 1:length(lS)
    S = lS[i]
    s = GAP.Globals.Representative(S)
    if !GAP.Globals.IsSubset(s, H)
      continue
    end
    sc = GAP.Globals.ComplementClassesRepresentatives(s, H)
    if GAP.Globals.IsEmpty(sc)
      continue
    end
    sncandidate = GAP.Globals.Size(s)
    if sncandidate > sizecandidate
      candidate = s
      sizecandidate = sncandidate
      ab_invariants = GAP.gap_to_julia(Vector{Int}, GAP.Globals.AbelianInvariants(sc[1]))
    end
  end
  return candidate, H, ab_invariants
end

################################################################################
#
#  Computing the extension over subfields
#
################################################################################

function computing_over_subfields(class_fields, subfields, idE, autos, right_grp, ab_invariants::Vector{Int}, is_normal_subfield::Bool)

  it = findall(right_grp)
  new_class_fields, subs, to_be_done = translate_class_field_down(subfields, class_fields, it, ab_invariants)
  for x in to_be_done
    C = class_fields[x]
    L = number_field(C, using_norm_relation = true)
    auts = absolute_automorphism_group(C, autos)
    Cpperm = permutation_group(auts)
    if GAP.Globals.IdGroup(Cpperm) == idE
      error("Something went wrong")
    end
    right_grp[x] = false
  end
  it = findall(right_grp)
  if isempty(it)
    return Vector{Tuple{Hecke.NfRelNS{nf_elem}, Vector{Hecke.NfRelNSToNfRelNSMor_nf_elem}}}()
  end
  use_brauer = true
  if !is_normal_subfield || !iszero(mod(order(torsion_unit_group(base_ring(new_class_fields[it[1]]))[1]), exponent(new_class_fields[it[1]])))
    use_brauer = false
  end
  for i in it
    number_field(new_class_fields[i], using_norm_relation = use_brauer)
  end
  translate_fields_up(class_fields, new_class_fields, subs, it)
  #Now, finally, the automorphisms computation and the isomorphism check
  for i in it
    C = class_fields[i]
    C1 = new_class_fields[i]
    mL = subs[1]
    indsubf = 1
    while domain(mL) != base_field(C1)
      indsubf += 1
      mL = subs[indsubf]
    end
    maprel = hom(C1.A, C.A, mL, gens(C.A))
    autsrelC1 = Hecke.rel_auto(C1)
    autsrelC = Vector{Hecke.NfRelNSToNfRelNSMor_nf_elem}(undef, length(autsrelC1))
    for s = 1:length(autsrelC1)
      el = autsrelC1[s]
      autsrelC[s] = hom(C.A, C.A, [maprel(x) for x in image_generators(el)])
      #@hassert :Fields 1 isconsistent(autsrelC[s])
    end
    rel_extend = Hecke.new_extend_aut(C, autos)
    autsA = vcat(rel_extend, autsrelC)
    C.AbsAutGrpA = autsA
    if !iscoprime(degree(C), degree(base_field(C)))
      Cpperm = permutation_group(autsA)
      if GAP.Globals.IdGroup(Cpperm) != idE
        right_grp[i] = false
      end
    end
  end
  it = findall(right_grp)
  fields = Vector{Tuple{Hecke.NfRelNS{nf_elem}, Vector{Hecke.NfRelNSToNfRelNSMor_nf_elem}}}(undef, length(it))
  ind = 1
  for i in it
    C = class_fields[i]
    fields[ind] = (C.A, C.AbsAutGrpA)
    ind += 1
  end
  return fields
end

function translate_extensions(mL::NfToNfMor, class_fields, new_class_fields, ctxK, it, ab_invariants::Vector{Int})

  to_be_done = Int[]
  L = domain(mL)
  OL = maximal_order(L)
  K = codomain(mL)
  OK = maximal_order(K)
  n = Int(exponent(class_fields[it[1]]))
  ordext = Int(degree(class_fields[it[1]]))
  ctx = Hecke.rayclassgrp_ctx(OL, n)
  d = divexact(discriminant(maximal_order(K)), discriminant(OL)^(divexact(degree(K), degree(L))))
  f = factor(ideal(OL, d))
  F = factor(ideal(OK, d))
  ab_invariants_mod = Int[]
  for i = 1:length(ab_invariants)
    r = gcd(ab_invariants[i], n)
    if isone(r)
      continue
    end
    push!(ab_invariants_mod, r)
  end
  for indclf in it
    if isassigned(new_class_fields, indclf)
      continue
    end
    C = class_fields[indclf]
    #First, I need a modulus.
    #I take the intersection of the modulus of C with L
    mR = C.rayclassgroupmap
    fM0 = copy(mR.fact_mod)
    fm0 = Dict{NfOrdIdl, Int}()
    for (p, v) in fM0
      p1 = Hecke.intersect_prime(mL, p)
      if !haskey(fm0, p1)
        if iscoprime(minimum(p1, copy = false), n)
          fm0[p1] = 1
        else
          fm0[p1] = v
        end
      end
      ep = divexact(ramification_index(p), ramification_index(p1))
      fM0[p] = ep*v
    end
    #Now, I have problems, so I need to add the ramification of the other extension.
    for (p, v) in f
      if !haskey(fm0, p)
        if isone(gcd(minimum(p), n))
          fm0[p] = 1
        else
          fm0[p] = v
        end
      else
        if !isone(gcd(minimum(p), n))
          fm0[p] = max(v, fm0[p])
        end
      end
      lPP = prime_decomposition(mL, p)
      for (P, vP) in lPP
        if haskey(fM0, P)
          fM0[P] = max(fM0[P], 2*vP*fm0[p])
        else
          fM0[P] = vP*fm0[p]*2
        end
      end
    end
    infplc = InfPlc[]
    if iszero(mod(n, 2))
      infplc = real_places(L)
    end
    @vprint :Fields 3 "Checking if I can compute $(indclf) over a subfield\n\n "
    @vtime :Fields 3 r, mr = Hecke.ray_class_group_quo(OL, fm0, infplc, ctx, check = false)
    if exponent(r) < n || order(r) < degree(C)
      push!(to_be_done, indclf)
      continue
    end
    #Now, the norm group of K over L
    @vtime :Fields 3 ngL, mngL = Hecke.norm_group(mL, mr, prod(ab_invariants_mod))
    @hassert :Fields 1 divisible(divexact(fmpz(degree(codomain(mL))), degree(domain(mL))), divexact(order(r), order(ngL)))
    if !divisible(order(ngL), degree(C)) || !divisible(exponent(C), n)
      push!(to_be_done, indclf)
      continue
    end
    #Finally, the norm group of C over L
    #I use the usual strategy, as in check_abelian_extension
    for (p, v) in F
      if iscoprime(minimum(p, copy = false), degree(C))
        fM0[p] = 1
      else
        if haskey(fM0, p)
          fM0[p] = max(v, fM0[p])
        else
          fM0[p] = v
        end
      end
    end
    inf_plc2 = InfPlc[]
    if !isempty(infplc)
      inf_plc2 = real_places(K)
    end
    @vtime :Fields 3 RM, mRM = Hecke.ray_class_group_quo(OK, fM0, inf_plc2, ctxK, check = false)
    if !isdefined(ctxK.class_group_map, :small_gens)
      ctxK.class_group_map.small_gens = find_gens(pseudo_inv(ctxK.class_group_map), PrimesSet(101, -1))[1]
    end
    @vtime :Fields 3 lP, gS = Hecke.find_gens(mRM, coprime_to = minimum(defining_modulus(mR)[1]))
    listn = NfOrdIdl[norm(mL, 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 C in RM
    prms = Vector{GrpAbFinGenElem}(undef, length(lP))
    for i = 1:length(lP)
      prms[i] = C.quotientmap(mR\lP[i])
    end
    RMtoR = hom(gS, prms)
    k, mk = kernel(RMtoR, false)
    @hassert :Fields 1 isisomorphic(cokernel(mk, false)[1], codomain(C.quotientmap))
    mp = mk*proj
    ck, mck = cokernel(mp, false)
    #If everything could work, then ck should be the direct product of the abelian extension I am searching for and
    #the maximal abelian subextension of K/L
    G1 = snf(cokernel(mngL, false)[1])[1]
    G2 = snf(codomain(C.quotientmap))[1]
    if !has_quotient(ck, map(Int, vcat(G1.snf, G2.snf)))
      push!(to_be_done, indclf)
      continue
    end
    @vprint :Fields 1 "$(Hecke.set_cursor_col())$(Hecke.clear_to_eol())Doing field $(indclf) / $(length(class_fields))"
    s, ms = image(mngL*mck, false)
    fl, ms1 = has_complement(ms)
    @assert fl
    mq1 = cokernel(ms1, false)[2]
    mqq = mck * mq1
    @hassert :Fields 1 domain(mqq) == r
    C1 = ray_class_field(mr, mqq)
    number_field(C1)
    @hassert :Fields 1 degree(C1) == degree(C)
    new_class_fields[indclf] = C1
  end
  @vprint :Fields 1 "$(Hecke.set_cursor_col())$(Hecke.clear_to_eol())"
  return to_be_done

end

function create_sub(ss, iso, PermGAP, auts, K)
  lS = GAP.Globals.MinimalGeneratingSet(ss)
  lS1 = GAP.GapObj[GAP.Globals.Image(iso, lS[i]) for i in 1:length(lS)]
  inds = Vector{Int}(undef, length(lS1))
  for j = 1:length(inds)
    el = lS1[j]
    inds[j] = findfirst(x -> x == el, PermGAP)
  end
  mL = Hecke.fixed_field(K, NfToNfMor[auts[i] for i in inds])[2]
  return mL
end

function compute_subfields(K::AnticNumberField, E, H, S)

  proj = GAP.Globals.NaturalHomomorphismByNormalSubgroup(E, H)
  Hn = GAP.Globals.ImagesSource(proj)
  imgsS = GAP.Globals.Image(proj, S)
  autsHn = GAP.Globals.Elements(GAP.Globals.AutomorphismGroup(Hn))
  orbitS = [imgsS]
  for i = 1:length(autsHn)
    varphi = autsHn[i]
    Sn = GAP.Globals.Image(varphi, imgsS)
    if !(Sn in orbitS)
      push!(orbitS, Sn)
    end
  end
  auts = automorphisms(K, copy = false)
  Hperm = _from_autos_to_perm(auts)
  Hauts = _perm_to_gap_grp(Hperm)
  PermGAP = Vector{GAP.GapObj}(undef, length(Hperm))
  for w = 1:length(Hperm)
    PermGAP[w] = _perm_to_gap_perm(Hperm[w])
  end
  iso = GAP.Globals.IsomorphismGroups(Hn, Hauts)
  return (create_sub(ss, iso, PermGAP, auts, K) for ss in orbitS)

end

function translate_class_field_down(subfields, class_fields, it, ab_invariants)
  new_class_fields = similar(class_fields)
  #Now, I translate the fields over the subfields.
  to_be_done = Int[i for i in it]
  created_subfields = NfToNfMor[]
  K = base_field(class_fields[to_be_done[1]])
  OK = maximal_order(K)
  ctxK = Hecke.rayclassgrp_ctx(OK, exponent(class_fields[to_be_done[1]]))
  for mL in subfields
    push!(created_subfields, mL)
    to_be_done_new = translate_extensions(mL, class_fields, new_class_fields, ctxK, it, ab_invariants)
    if length(to_be_done_new) == 0
      return new_class_fields, created_subfields, to_be_done_new
    end
    to_be_done = to_be_done_new
  end
  return new_class_fields, created_subfields, to_be_done
end

function translate_fields_up(class_fields, new_class_fields, subfields, it)
  K = base_field(class_fields[it[1]])
  Ky = PolynomialRing(K, "y", cached = false)[1]
  for i in it
    C = class_fields[i]
    C1 = new_class_fields[i]
    mL = subfields[1]
    indsubf = 1
    while domain(mL) != base_field(C1)
      indsubf += 1
      mL = subfields[indsubf]
    end
    L = domain(mL)
    D = Dict{Int, NfToNfMor}()
    for j = 1:length(new_class_fields[i].cyc)
      d = degree(new_class_fields[i].cyc[j])
      if !haskey(D, d)
        CEK = cyclotomic_extension(K, d)
        CEL = cyclotomic_extension(L, d)
        img = gen(CEK.Kr)
        if degree(CEK.Kr) != euler_phi(d)
          pp = map_coefficients(mL, CEL.Kr.pol, cached = false)
          while !iszero(pp(img))
            mul!(img, img, gen(CEK.Kr))
          end
        end
        mrel = hom(CEL.Kr, CEK.Kr, mL, img)
        #@hassert :Fields 1 isconsistent(mrel)
        g = mrel(CEL.mp[1](gen(CEL.Ka)))
        mp = hom(CEL.Ka, CEK.Ka, CEK.mp[1]\(g), check = false)
        D[d] = mp
      end
    end
    cyc = Vector{Hecke.ClassField_pp{Hecke.MapRayClassGrp, GrpAbFinGenMap}}(undef, length(C1.cyc))
    for j = 1:length(cyc)
      Ccyc = C1.cyc[j]
      d = degree(Ccyc)
      #First, easy: the degree
      Cpp = Hecke.ClassField_pp{Hecke.MapRayClassGrp, GrpAbFinGenMap}()
      Cpp.rayclassgroupmap = C.rayclassgroupmap
      Cpp.degree = d
      #Then, the fac elem corresponding to the generator of the Kummer Extension
      Cpp.a = FacElem(Dict{nf_elem, fmpz}(D[d](x) => v for (x, v) in Ccyc.a))
      #Now, the Kummer extension
      Lzeta = codomain(D[d])
      Lt = PolynomialRing(Lzeta, "t", cached = false)[1]
      d1 = degree(Ccyc.K)
      coeffs = Vector{nf_elem}(undef, d1 + 1)
      coeffs[1] = D[d](coeff(Ccyc.K.pol, 0))
      for s = 2:length(coeffs)-1
        coeffs[s] = zero(Lzeta)
      end
      coeffs[end] = one(Lzeta)
      Cpp.K = number_field(Lt(coeffs), cached = false, check = false)[1]
      #The target extension
      fdef = map_coefficients(mL, Ccyc.A.pol, parent = Ky, cached = false)
      Cpp.A = number_field(fdef, cached = false, check = false)[1]
      #Now, the primitive element of the target extension seen in Cpp.K
      mrel2 = hom(Ccyc.K, Cpp.K, D[d], gen(Cpp.K))
      #@hassert :Fields 1 isconsistent(mrel2)
      @hassert :Fields 1 parent(Ccyc.pe) == domain(mrel2)
      Cpp.pe = mrel2(Ccyc.pe)
      CEKK = cyclotomic_extension(K, d)
      @hassert :Fields 1 iszero(map_coefficients(CEKK.mp[2], fdef, cached = false)(Cpp.pe))
      Cpp.o = d1
      cyc[j] = Cpp
    end
    C.cyc = cyc
    C.A = number_field([cyc[j].A.pol for j = 1:length(cyc)], cached = false, check = false)[1]
  end
  return nothing

end


text:
using StatisticalRethinking, Optim
#gr(size=(600,600));

p_grid = range(0, step=0.001, stop=1)
prior = ones(length(p_grid))
likelihood = [pdf(Binomial(9, p), 6) for p in p_grid]
posterior = likelihood .* prior
posterior = posterior / sum(posterior)
samples = sample(p_grid, Weights(posterior), length(p_grid));
samples[1:5]

N = 10000
samples = sample(p_grid, Weights(posterior), N);

chn = MCMCChains.Chains(reshape(samples, N, 1, 1), ["toss"]);

MCMCChains.describe(chn)

plot(chn)

p = Vector{Plots.Plot{Plots.GRBackend}}(undef, 2)
p[1] = scatter(1:N, samples, markersize = 2, ylim=(0.0, 1.3), lab="Draws")

w = 6
n = 9
x = 0:0.01:1
p[2] = density(samples, ylim=(0.0, 5.0), lab="Sample density")
p[2] = plot!( x, pdf.(Beta( w+1 , n-w+1 ) , x ), lab="Conjugate solution")

plot(p..., layout=(1, 2))

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl

ext: jl


text:
using OrdinaryDiffEq, PoissonRandom

function rate_to_proportion(r::Float64,t::Float64)
    1-exp(-r*t)
end;

function sir_abm!(du,u,p,t)
    (β,c,γ,δt) = p
    N = length(u)
    # Initialize du to u
    for i in 1:N
        du[i] = u[i]
    end
    for i in 1:N # loop through agents
        # If recovered
        if u[i]==Recovered
            continue
        # If susceptible
        elseif u[i]==Susceptible
            ncontacts = pois_rand(c*δt)
            while ncontacts > 0
                j = rand(1:N)
                if j==i
                    continue
                end
                a = u[j]
                if a==Infected && rand() < β
                    du[i] = Infected
                    break
                end
                ncontacts -= 1
            end
        # If infected
        else u[i]==Infected
            if rand() < γ
                du[i] = Recovered
            end
        end
    end
    nothing
end

δt = 0.1
nsteps = 400
tf = nsteps*δt
tspan = (0.0,nsteps)
t = 0:δt:tf;
β = 0.05
c = 10.0
γ = rate_to_proportion(0.25,δt)
p = [β,c,γ,δt]

@enum InfectionStatus Susceptible Infected Recovered
Base.zero(::Type{<:InfectionStatus}) = Susceptible

N = 1000
I0 = 10
u0 = Array{InfectionStatus}(undef,N);
for i in 1:N
    if i <= I0
        s = Infected
    else
        s = Susceptible
    end
    u0[i] = s
end

prob_map = DiscreteProblem(sir_abm!,u0,tspan,p)
sol_map = solve(prob_map,FunctionMap())


text:
#!/usr/bin/env julia

const Tlist = [0.9, 0.8, 0.7]
const tw1list = map(Int64,[1e4,1e6,1e8])

const reps=8

function rectify_times!(array)
    L = length(array)
    for i in 2:L
        if array[i] == array[i-1]
            for j in i:L
                array[j] += 1
            end
        end
    end
end

function generate_times()
    tws = logspace(log10(1),log10(3e8+10),100)
    tws = map(x->floor(Int64,x), tws)
    tws = vcat(tws, tw1list) |> sort
    rectify_times!(tws)

    tws = unique(tws) |> sort
    ts = logspace(4,8,10)
    ts = map(x->floor(Int64,x), ts) |> sort

    allt = [tw + t for tw in tws, t in [0; ts]] |> vec |> unique |> sort

    writedlm("tiempos_log.dat", allt)
    writedlm("tiempos_lineal.dat",[])
end

if isdir("jobs")
    error("Cowardly refusing to use existing directory 'jobs'")
else
    mkdir("jobs")
end

function generate_jobs(;T=0.9)
    suffix = "quench_T_$(T)"
    mkdir("jobs/$suffix")
    for i in 1:reps
        jobid  = join(rand('a':'z',5))
        open("jobs/$(suffix)/$(jobid).job", "w") do f
            println(f, "/home/sergio/output/$(suffix)_replica_$(i)    Nombre del directorio")
            println(f, "400000001                                Timewindow (maximo numero pasos MC simulacion) ")
            println(f, "0                                        FlagJ                                          ")
            println(f, "0                                        FlagS                                          ")
            println(f, "666                                      Flag_h (now, flag_dummy: not used)             ")
            println(f, "0                                        Flag_continua                                  ")
            println(f, "$T                                      T1                                             ")
            println(f, "$T                                      T2                                             ")
            println(f, "1                                     tw1                                            ")
            println(f, "2                                    tw2                                            ")
            println(f, "0                                        Hfield                                         ")
            println(f, "1                                        TimeHon                                        ")
        end
    end
end

generate_times()

info("Generating jobs")
f = 0
for T in Tlist
    generate_jobs(T=T)
    f+=1
    print("$f*$reps files generated  \r")
end
println()

max_stars_repo_name: RedPointyJackson/tfg

max_issues_repo_name: RedPointyJackson/tfg

max_forks_repo_name: RedPointyJackson/tfg


text:
using SymEngine

"""
Saves info given by the user for a given stencil.
There are 2 "modes":
    + Expression with optional coefficients.
    + A cubic matrix with coefficient and uses_vsq flag.
"""
mutable struct StencilDefinition
    expr::Union{Expr, Nothing}
    coefs::Union{Array{Float64, N} where N ,Array{Basic, N} where N ,Nothing}
    max_radius::Integer
    uses_vsq::Bool
end

"""
Saves the generated function and other revelant info.
"""
mutable struct StencilInstance
    stencil::Union{StencilDefinition, NTuple{N,StencilDefinition} where N}
    stencil_sym::Basic
    max_radius::Int
    front_z_max::Int
    behind_z_max::Int
    bdim::Int
    #syms::Tuple{Symbol,3}
    #data::Array{Float32,3}
    #coefs::Array{Float32,N} where N
    kernel::Function
    kernel_expr::Expr
    uses_vsq::Bool
    combined_time_step::Union{Bool, Int}
    m_kernel::Union{Nothing, Function}
    m_kernel_exp::Union{Nothing, Expr}
    m_kernel_sym::Union{Nothing, Basic}
    m_max_radius::Union{Bool, Int}
    m_front_z_max::Union{Bool, Int}
    m_behind_z_max::Union{Bool, Int}
end

mutable struct stencil_info
    exists::Bool
    sym_math::Basic
    front_max::Integer
    behind_max::Integer
    max_radius::Integer
    stencil_info() = new(false, symbols(:a), 0,0,0)
end
nothing


text:
import ApproxFun: weight, dimension

export JacobiQWeight

"""
`JacobiQWeight`
weights a basis on `ℂ\\𝕀` weighted by `(z+1)^α*(z-1)^β`.
Note the inconsistency of the parameters with `JacobiQ`.
When the domain is `[a,b]` the weight is inferred by mapping to `[-1,1]`
"""
struct JacobiQWeight{S,DD} <: WeightSpace{S,DD,Float64}
    α::Float64
    β::Float64
    space::S
    function JacobiQWeight{S,DD}(α::Float64,β::Float64,space::S) where {S,DD}
        if isa(space,JacobiQWeight)
            JacobiQWeight(α+space.α,β+space.β,space.space)
        else
            new{S,DD}(α,β,space)
        end
    end
end

JacobiQWeight(a::Number,b::Number,d::RealUnivariateSpace) =
    JacobiQWeight{typeof(d),typeof(domain(d))}(Float64(a),Float64(b),d)
JacobiQWeight(a::Number,b::Number,d::IntervalOrSegmentDomain) =
    JacobiQWeight(Float64(a),Float64(b),Space(d))
JacobiQWeight(a::Number,b::Number,d::Vector) =
    JacobiQWeight(Float64(a),Float64(b),Space(d))
JacobiQWeight(a::Number,b::Number) = JacobiQWeight(a,b,Chebyshev())

JacobiQWeight(a::Number,b::Number,s::PiecewiseSpace) = PiecewiseSpace(JacobiQWeight(a,b,vec(s)))

spacescompatible(A::JacobiQWeight,B::JacobiQWeight) =
    A.α==B.α && A.β == B.β && spacescompatible(A.space,B.space)

transformtimes(f::Fun{JW1},g::Fun{JW2}) where {JW1<:JacobiQWeight,JW2<:JacobiQWeight}=
            Fun(JacobiQWeight(f.space.α+g.space.α,f.space.β+g.space.β,f.space.space),
                coefficients(transformtimes(Fun(f.space.space,f.coefficients),
                                            Fun(g.space.space,g.coefficients))))
transformtimes(f::Fun{JW},g::Fun) where {JW<:JacobiQWeight} = Fun(f.space,coefficients(transformtimes(Fun(f.space.space,f.coefficients),g)))
transformtimes(f::Fun,g::Fun{JW}) where {JW<:JacobiQWeight} = Fun(g.space,coefficients(transformtimes(Fun(g.space.space,g.coefficients),f)))

##  α and β are opposite the convention for JacobiQ polynomials
# Here, α is the left algebraic singularity and β is the right algebraic singularity.

jacobiQweight(α,β,x) = (x+1)^α*(x-1)^β
jacobiQweight(α,β,d::Domain) = Fun(JacobiQWeight(α,β,ConstantSpace(d)),[1.])
jacobiQweight(α,β) = jacobiQweight(α,β,ChebyshevInterval())

weight(sp::JacobiQWeight,x) = jacobiQweight(sp.α,sp.β,tocanonical(sp.space.domain,x))
dimension(sp::JacobiQWeight) = dimension(sp.space)

Base.first(f::Fun{JW}) where {JW<:JacobiQWeight} =
    space(f).α > 0 ? zero(cfstype(f)) : f(first(domain(f)))
Base.last(f::Fun{JW}) where {JW<:JacobiQWeight} =
    space(f).β > 0 ? zero(cfstype(f)) : f(last(domain(f)))

setdomain(sp::JacobiQWeight,d::Domain) = JacobiQWeight(sp.α,sp.β,setdomain(sp.space,d))


text:
using IOMCovidData

include("user\\paths.jl")

IOMCovidData.updatedatacsv(coviddatafolder)

max_stars_repo_name: dannywinrow/iom-govt-data

max_issues_repo_name: dannywinrow/iom-govt-data

max_forks_repo_name: dannywinrow/iom-govt-data


text:
using SafeTestsets

@safetestset "Kalman Filter tests" begin include("kalman_filter_tests.jl") end
@safetestset "Likelihood tests" begin include("likelihood_tests.jl") end
@safetestset "ODE tests" begin include("ode_tests.jl") end


text:
# calculate crc checksums for data.  this process is closely related
# to the routines in IntModN (particularly GF2Poly), but the
# traditional implementations are optimized to the point where
# implementing them from "first principles" makes little sense
# (although we can use IntModN to check the results here).

# the basic idea is to do a polynomial division along the byte stream.
# it's all explained very nicely at
# http://en.wikipedia.org/wiki/Computation_of_cyclic_redundancy_checks
# and http://www.zlib.net/crc_v3.txt

module CRC

include("crc_calc.jl")
include("crc_cmdline.jl")

end


text:
using SparseArrays
using Profile

"""
    assemble!(dofmap, a_kernel!, L_kernel!, bc_dofs, bc_vals)

Assemble bilinear and linear kernels, `a_kernel`, `L_kernel` and
applies boundary contitions in symmetric way.

Allocates and returns sparse CSC matrix and dense vector
of rhs.

"""
function assemble!(dofmap::DofMap, a_kernel!::Function, L_kernel!::Function,
    bc_dofs::Array{Bool, 1}, bc_vals::Array{Float64, 1})

    mesh = dofmap.mesh
    dofsize = dofmap.dim
    tdim = mesh.reference_cell.get_dimension()::Int64

    dofsize_local = size(dofmap.cell_dofs)[2]
    num_values = dofsize_local ^ 2
    num_cells = num_entities(mesh, tdim)
    cell_vert_conn = mesh.topology[(tdim, 0)]

    vert_per_cell = size(cell_vert_conn)[2]
    gdim = size(mesh.vertices)[2]

    # Vector of non-zero values in the global matrix
    # These values will be scattered according to "sparsity pattern"
    values = Array{Float64, 1}()

    # Storage for coordinates of cell vertices
    cell_coords = zeros(Float64, vert_per_cell, gdim)

    # Storage for insertion index sets
    I = Array{Int64, 1}()
    J = Array{Int64, 1}()

    b = zeros(Float64, dofsize)

    A_local = zeros(Float64, dofsize_local, dofsize_local)
    b_local = zeros(Float64, dofsize_local)

    for cell_id in 1:num_cells

        fill!(A_local, 0.0)
        fill!(b_local, 0.0)

        # Pack vertex coordinates of the cell
        cell_coords[1:vert_per_cell, 1:gdim] = mesh.vertices[cell_vert_conn[cell_id, 1:vert_per_cell], 1:gdim]

        a_kernel!(A_local, cell_coords)
        L_kernel!(b_local, cell_coords)

        # Append values to COO format
        for (i, iglobal) in enumerate(dofmap.cell_dofs[cell_id, 1:dofsize_local])

            if bc_dofs[iglobal]
                A_local[i, 1:dofsize_local] .= 0.0
                b_local[1:dofsize_local] -= A_local[1:dofsize_local, i] * bc_vals[iglobal]
                A_local[1:dofsize_local, i] .= 0.0
                A_local[i, i] = 1.0
                b_local[i] = bc_vals[iglobal]
            end

            b[iglobal] += b_local[i]
            for j in 1:dofsize_local
                # Push only nonzero values to avoid later "eliminate zeros" call
                if A_local[i, j] != 0.0
                    push!(values, A_local[i, j])
                    push!(I, iglobal)
                    push!(J, dofmap.cell_dofs[cell_id, j])
                end
            end
        end
    end

    # Return CSC sparse array
    A = sparse(I, J, values)

    return A, b
end


text:
# Copyright (c) 2021, Joseph M. Rutherford
# All rights reserved.
#
# Code provided under the license contained in the LICENSE file.
#
# This file includes the Sudoku module source code.
# This file imports all sudoku_*_test.jl scripts into a separate namespace.
# Tests are executed using runtests.jl.

module SudokuTester

include("../src/sudoku.jl")

include("sudoku_common_test.jl")
include("sudoku_permute_test.jl")
include("sudoku_compose_test.jl")
include("sudoku_valid_test.jl")
include("sudoku_solve_test.jl")

end # end module

max_stars_repo_name: joseph-m-rutherford/sudoku-julia

max_issues_repo_name: joseph-m-rutherford/sudoku-julia

max_forks_repo_name: joseph-m-rutherford/sudoku-julia


text:
#
# VUMPS code
#

struct Hᶜ
  ∑h::InfiniteITensorSum
  Hᴸ::InfiniteMPS
  Hᴿ::InfiniteMPS
  ψ::InfiniteCanonicalMPS
  n::Int
end

struct Hᴬᶜ
  ∑h::InfiniteITensorSum
  Hᴸ::InfiniteMPS
  Hᴿ::InfiniteMPS
  ψ::InfiniteCanonicalMPS
  n::Int
end

function Base.:*(H::Hᶜ, v::ITensor)
  ∑h = H.∑h
  Hᴸ = H.Hᴸ
  Hᴿ = H.Hᴿ
  ψ = H.ψ
  ψ′ = dag(ψ)'
  n = H.n
  l = linkinds(only, ψ.AL)
  l′ = linkinds(only, ψ′.AL)
  r = linkinds(only, ψ.AR)
  r′ = linkinds(only, ψ′.AR)
  δˡ = δ(l[n], l′[n])
  δˡ⁻¹ = δ(l[n - 1], l′[n - 1])
  δʳ = δ(r[n], r′[n])
  δʳ⁺¹ = δ(r[n + 1], r′[n + 1])
  Hᶜᴸv = v * Hᴸ[n] * dag(δʳ)
  Hᶜᴿv = v * δˡ * Hᴿ[n]
  Hᶜʰv =
    v * ψ.AL[n] * δˡ⁻¹ * ψ′.AL[n] * ∑h[(n, n + 1)] * ψ.AR[n + 1] * dag(δʳ⁺¹) * ψ′.AR[n + 1]
  Hᶜv = Hᶜᴸv + Hᶜʰv + Hᶜᴿv
  return Hᶜv * dag(δˡ) * δʳ
end

function Base.:*(H::Hᴬᶜ, v::ITensor)
  ∑h = H.∑h
  Hᴸ = H.Hᴸ
  Hᴿ = H.Hᴿ
  ψ = H.ψ
  ψ′ = dag(ψ)'
  n = H.n
  l = linkinds(only, ψ.AL)
  l′ = linkinds(only, ψ′.AL)
  r = linkinds(only, ψ.AR)
  r′ = linkinds(only, ψ′.AR)
  s = siteinds(only, ψ)
  s′ = siteinds(only, ψ′)

  δˢ(n) = δ(s[n], s′[n])
  δˡ(n) = δ(l[n], l′[n])
  δʳ(n) = δ(r[n], r′[n])

  Hᴬᶜᴸv = v * Hᴸ[n - 1] * dag(δˢ(n)) * dag(δʳ(n))
  Hᴬᶜᴿv = v * δˡ(n - 1) * dag(δˢ(n)) * Hᴿ[n]
  Hᴬᶜʰ¹v = v * ψ.AL[n - 1] * δˡ(n - 2) * ψ′.AL[n - 1] * ∑h[(n - 1, n)] * dag(δʳ(n))
  Hᴬᶜʰ²v = v * δˡ(n - 1) * ψ.AR[n + 1] * dag(δʳ(n + 1)) * ψ′.AR[n + 1] * ∑h[(n, n + 1)]
  Hᴬᶜv = Hᴬᶜᴸv + Hᴬᶜʰ¹v + Hᴬᶜʰ²v + Hᴬᶜᴿv
  return Hᴬᶜv * dag(δˡ(n - 1)) * δˢ(n) * δʳ(n)
end

function (H::Hᶜ)(v)
  Hᶜv = H * v
  ## return Hᶜv * δˡ * δʳ
  return noprime(Hᶜv)
end

function (H::Hᴬᶜ)(v)
  Hᴬᶜv = H * v
  ## return Hᶜv * δˡ⁻¹ * δˢ * δʳ
  return noprime(Hᴬᶜv)
end

# Struct for use in linear system solver
struct Aᴸ
  ψ::InfiniteCanonicalMPS
  n::Int
end

function (A::Aᴸ)(x)
  ψ = A.ψ
  ψᴴ = dag(ψ)
  ψ′ = ψᴴ'
  ψ̃ = prime(linkinds, ψᴴ)
  n = A.n

  N = length(ψ)
  #@assert n == N

  l = linkinds(only, ψ.AL)
  l′ = linkinds(only, ψ′.AL)
  r = linkinds(only, ψ.AR)
  r′ = linkinds(only, ψ′.AR)

  xT = translatecell(x, -1)
  for k in (n - N + 1):n
    xT = xT * ψ.AL[k] * ψ̃.AL[k]
  end
  δˡ = δ(l[n], l′[n])
  δʳ = δ(r[n], r′[n])
  xR = x * ψ.C[n] * ψ′.C[n] * dag(δʳ) * denseblocks(δˡ)
  return xT - xR
end

function left_environment(hᴸ, 𝕙ᴸ, ψ; tol=1e-15)
  ψ̃ = prime(linkinds, dag(ψ))
  N = nsites(ψ)

  Aᴺ = Aᴸ(ψ, N)
  Hᴸᴺ¹, info = linsolve(Aᴺ, 𝕙ᴸ[N], 1, -1; tol=tol)
  # Get the rest of the environments in the unit cell
  Hᴸ = InfiniteMPS(Vector{ITensor}(undef, N))
  Hᴸ[N] = Hᴸᴺ¹
  Hᴸᴺ¹ = translatecell(Hᴸᴺ¹, -1)
  for n in 1:(N - 1)
    Hᴸ[n] = Hᴸ[n - 1] * ψ.AL[n] * ψ̃.AL[n] + hᴸ[n]
  end
  # Compute more accurate environments
  # Not currently working
  #for n in 1:(N - 1)
  #  Aⁿ = Aᴸ(ψ, n)
  #  Hᴸ[n], info = linsolve(Aⁿ, 𝕙ᴸ[n], Hᴸ[n], 1, -1; tol=tol)
  #end
  return Hᴸ
end

# Struct for use in linear system solver
struct Aᴿ
  hᴿ::InfiniteMPS
  ψ::InfiniteCanonicalMPS
  n::Int
end

function (A::Aᴿ)(x)
  hᴿ = A.hᴿ
  ψ = A.ψ
  ψᴴ = dag(ψ)
  ψ′ = ψᴴ'
  ψ̃ = prime(linkinds, ψᴴ)
  n = A.n

  N = length(ψ)
  @assert n == N

  l = linkinds(only, ψ.AL)
  l′ = linkinds(only, ψ′.AL)
  r = linkinds(only, ψ.AR)
  r′ = linkinds(only, ψ′.AR)

  xT = x
  for k in reverse(1:N)
    xT = xT * ψ.AR[k] * ψ̃.AR[k]
  end
  xT = translatecell(xT, 1)
  δˡ = δ(l[n], l′[n])
  δʳ = δ(r[n], r′[n])
  xR = x * ψ.C[n] * ψ′.C[n] * δˡ * denseblocks(dag(δʳ))
  return xT - xR
end

function right_environment(hᴿ, ψ; tol=1e-15)
  ψ̃ = prime(linkinds, dag(ψ))
  # XXX: replace with `nsites`
  #N = nsites(ψ)
  N = length(ψ)

  A = Aᴿ(hᴿ, ψ, N)
  Hᴿᴺ¹, info = linsolve(A, hᴿ[N], 1, -1; tol=tol)
  # Get the rest of the environments in the unit cell
  Hᴿ = InfiniteMPS(Vector{ITensor}(undef, N))
  Hᴿ[N] = Hᴿᴺ¹
  for n in reverse(1:(N - 1))
    Hᴿ[n] = Hᴿ[n + 1] * ψ.AR[n + 1] * ψ̃.AR[n + 1] + hᴿ[n]
  end
  return Hᴿ
end

function vumps_iteration(args...; multisite_update_alg="sequential", kwargs...)
  if multisite_update_alg == "sequential"
    return vumps_iteration_sequential(args...; kwargs...)
  elseif multisite_update_alg == "parallel"
    return vumps_iteration_parallel(args...; kwargs...)
  else
    error(
      "Multisite update algorithm multisite_update_alg = $multisite_update_alg not supported, use \"parallel\" or \"sequential\"",
    )
  end
end

function vumps_iteration_sequential(
  ∑h::InfiniteITensorSum,
  ψ::InfiniteCanonicalMPS;
  (ϵᴸ!)=fill(1e-15, nsites(ψ)),
  (ϵᴿ!)=fill(1e-15, nsites(ψ)),
  eigsolve_tol=(x -> x / 100),
)
  Nsites = nsites(ψ)
  ϵᵖʳᵉˢ = max(maximum(ϵᴸ!), maximum(ϵᴿ!))
  krylov_tol = eigsolve_tol(ϵᵖʳᵉˢ)
  ψᴴ = dag(ψ)
  ψ′ = ψᴴ'
  # XXX: make this prime the center sites
  ψ̃ = prime(linkinds, ψᴴ)

  # TODO: replace with linkinds(ψ)
  l = CelledVector([commoninds(ψ.AL[n], ψ.AL[n + 1]) for n in 1:Nsites])
  l′ = CelledVector([commoninds(ψ′.AL[n], ψ′.AL[n + 1]) for n in 1:Nsites])
  r = CelledVector([commoninds(ψ.AR[n], ψ.AR[n + 1]) for n in 1:Nsites])
  r′ = CelledVector([commoninds(ψ′.AR[n], ψ′.AR[n + 1]) for n in 1:Nsites])

  ψ = copy(ψ)
  C̃ = InfiniteMPS(Vector{ITensor}(undef, Nsites))
  Ãᶜ = InfiniteMPS(Vector{ITensor}(undef, Nsites))
  Ãᴸ = InfiniteMPS(Vector{ITensor}(undef, Nsites))
  Ãᴿ = InfiniteMPS(Vector{ITensor}(undef, Nsites))
  eᴸ = Vector{Float64}(undef, Nsites)
  eᴿ = Vector{Float64}(undef, Nsites)
  for n in 1:Nsites
    hᴸ = InfiniteMPS([
      δ(only(l[k - 2]), only(l′[k - 2])) *
      ψ.AL[k - 1] *
      ψ.AL[k] *
      ∑h[(k - 1, k)] *
      ψ′.AL[k - 1] *
      ψ′.AL[k] for k in 1:Nsites
    ])
    hᴿ = InfiniteMPS([
      δ(only(dag(r[k + 2])), only(dag(r′[k + 2]))) *
      ψ.AR[k + 2] *
      ψ.AR[k + 1] *
      ∑h[(k + 1, k + 2)] *
      ψ′.AR[k + 2] *
      ψ′.AR[k + 1] for k in 1:Nsites
    ])
    eᴸ = [
      (hᴸ[k] * ψ.C[k] * δ(only(dag(r[k])), only(dag(r′[k]))) * ψ′.C[k])[] for k in 1:Nsites
    ]
    eᴿ = [(hᴿ[k] * ψ.C[k] * δ(only(l[k]), only(l′[k])) * ψ′.C[k])[] for k in 1:Nsites]
    for k in 1:Nsites
      # TODO: remove `denseblocks` once BlockSparse + DiagBlockSparse is supported
      hᴸ[k] -= eᴸ[k] * denseblocks(δ(inds(hᴸ[k])))
      hᴿ[k] -= eᴿ[k] * denseblocks(δ(inds(hᴿ[k])))
    end

    function left_environment_cell(ψ, ψ̃, hᴸ, n)
      Nsites = nsites(ψ)
      𝕙ᴸ = copy(hᴸ)
      for k in reverse((n - Nsites + 2):n)
        𝕙ᴸ[k] = 𝕙ᴸ[k - 1] * ψ.AL[k] * ψ̃.AL[k] + 𝕙ᴸ[k]
      end
      return 𝕙ᴸ[n]
    end

    #for k in 2:Nsites
    #  hᴸ[k] = hᴸ[k - 1] * ψ.AL[k] * ψ̃.AL[k] + hᴸ[k]
    #end
    𝕙ᴸ = copy(hᴸ)
    for k in 1:Nsites
      𝕙ᴸ[k] = left_environment_cell(ψ, ψ̃, hᴸ, k)
    end
    Hᴸ = left_environment(hᴸ, 𝕙ᴸ, ψ; tol=krylov_tol)
    for k in 2:Nsites
      hᴿ[k] = hᴿ[k + 1] * ψ.AR[k + 1] * ψ̃.AR[k + 1] + hᴿ[k]
    end
    Hᴿ = right_environment(hᴿ, ψ; tol=krylov_tol)

    Cvalsₙ₋₁, Cvecsₙ₋₁, Cinfoₙ₋₁ = eigsolve(
      Hᶜ(∑h, Hᴸ, Hᴿ, ψ, n - 1), ψ.C[n - 1], 1, :SR; ishermitian=true, tol=krylov_tol
    )
    Cvalsₙ, Cvecsₙ, Cinfoₙ = eigsolve(
      Hᶜ(∑h, Hᴸ, Hᴿ, ψ, n), ψ.C[n], 1, :SR; ishermitian=true, tol=krylov_tol
    )
    Avalsₙ, Avecsₙ, Ainfoₙ = eigsolve(
      Hᴬᶜ(∑h, Hᴸ, Hᴿ, ψ, n), ψ.AL[n] * ψ.C[n], 1, :SR; ishermitian=true, tol=krylov_tol
    )
    C̃[n - 1] = Cvecsₙ₋₁[1]
    C̃[n] = Cvecsₙ[1]
    Ãᶜ[n] = Avecsₙ[1]

    function ortho_overlap(AC, C)
      AL, _ = polar(AC * dag(C), uniqueinds(AC, C))
      return noprime(AL)
    end

    function ortho_polar(AC, C)
      UAC, _ = polar(AC, uniqueinds(AC, C))
      UC, _ = polar(C, commoninds(C, AC))
      return noprime(UAC) * noprime(dag(UC))
    end

    Ãᴸ[n] = ortho_polar(Ãᶜ[n], C̃[n])
    Ãᴿ[n] = ortho_polar(Ãᶜ[n], C̃[n - 1])

    # Update state for next iteration
    #ψ = InfiniteCanonicalMPS(Ãᴸ, C̃, Ãᴿ)
    ψ.AL[n] = Ãᴸ[n]
    ψ.AR[n] = Ãᴿ[n]
    ψ.C[n - 1] = C̃[n - 1]
    ψ.C[n] = C̃[n]
    ψᴴ = dag(ψ)
    ψ′ = ψᴴ'
    # XXX: make this prime the center sites
    ψ̃ = prime(linkinds, ψᴴ)

    # TODO: replace with linkinds(ψ)
    l = CelledVector([commoninds(ψ.AL[n], ψ.AL[n + 1]) for n in 1:Nsites])
    l′ = CelledVector([commoninds(ψ′.AL[n], ψ′.AL[n + 1]) for n in 1:Nsites])
    r = CelledVector([commoninds(ψ.AR[n], ψ.AR[n + 1]) for n in 1:Nsites])
    r′ = CelledVector([commoninds(ψ′.AR[n], ψ′.AR[n + 1]) for n in 1:Nsites])
  end
  for n in 1:Nsites
    ϵᴸ![n] = norm(Ãᶜ[n] - Ãᴸ[n] * C̃[n])
    ϵᴿ![n] = norm(Ãᶜ[n] - C̃[n - 1] * Ãᴿ[n])
  end
  return ψ, (eᴸ, eᴿ)
end

function vumps_iteration_parallel(
  ∑h::InfiniteITensorSum,
  ψ::InfiniteCanonicalMPS;
  (ϵᴸ!)=fill(1e-15, nsites(ψ)),
  (ϵᴿ!)=fill(1e-15, nsites(ψ)),
  eigsolve_tol=(x -> x / 100),
)
  Nsites = nsites(ψ)
  ϵᵖʳᵉˢ = max(maximum(ϵᴸ!), maximum(ϵᴿ!))
  krylov_tol = ϵᵖʳᵉˢ / 100
  ψᴴ = dag(ψ)
  ψ′ = ψᴴ'
  # XXX: make this prime the center sites
  ψ̃ = prime(linkinds, ψᴴ)

  # TODO: replace with linkinds(ψ)
  l = CelledVector([commoninds(ψ.AL[n], ψ.AL[n + 1]) for n in 1:Nsites])
  l′ = CelledVector([commoninds(ψ′.AL[n], ψ′.AL[n + 1]) for n in 1:Nsites])
  r = CelledVector([commoninds(ψ.AR[n], ψ.AR[n + 1]) for n in 1:Nsites])
  r′ = CelledVector([commoninds(ψ′.AR[n], ψ′.AR[n + 1]) for n in 1:Nsites])

  hᴸ = InfiniteMPS([
    δ(only(l[n - 2]), only(l′[n - 2])) *
    ψ.AL[n - 1] *
    ψ.AL[n] *
    ∑h[(n - 1, n)] *
    ψ′.AL[n - 1] *
    ψ′.AL[n] for n in 1:Nsites
  ])

  hᴿ = InfiniteMPS([
    δ(only(dag(r[n + 2])), only(dag(r′[n + 2]))) *
    ψ.AR[n + 2] *
    ψ.AR[n + 1] *
    ∑h[(n + 1, n + 2)] *
    ψ′.AR[n + 2] *
    ψ′.AR[n + 1] for n in 1:Nsites
  ])

  eᴸ = [
    (hᴸ[n] * ψ.C[n] * δ(only(dag(r[n])), only(dag(r′[n]))) * ψ′.C[n])[] for n in 1:Nsites
  ]
  eᴿ = [(hᴿ[n] * ψ.C[n] * δ(only(l[n]), only(l′[n])) * ψ′.C[n])[] for n in 1:Nsites]

  for n in 1:Nsites
    # TODO: use these instead, for now can't subtract
    # BlockSparse and DiagBlockSparse tensors
    #hᴸ[n] -= eᴸ[n] * δ(inds(hᴸ[n]))
    #hᴿ[n] -= eᴿ[n] * δ(inds(hᴿ[n]))
    hᴸ[n] -= eᴸ[n] * denseblocks(δ(inds(hᴸ[n])))
    hᴿ[n] -= eᴿ[n] * denseblocks(δ(inds(hᴿ[n])))
  end

  # Sum the Hamiltonian terms in the unit cell
  for n in 2:Nsites
    hᴸ[n] = hᴸ[n - 1] * ψ.AL[n] * ψ̃.AL[n] + hᴸ[n]
  end
  Hᴸ = left_environment(hᴸ, ψ; tol=krylov_tol)

  for n in 2:Nsites
    hᴿ[n] = hᴿ[n + 1] * ψ.AR[n + 1] * ψ̃.AR[n + 1] + hᴿ[n]
  end
  Hᴿ = right_environment(hᴿ, ψ; tol=krylov_tol)

  C̃ = InfiniteMPS(Vector{ITensor}(undef, Nsites))
  Ãᶜ = InfiniteMPS(Vector{ITensor}(undef, Nsites))
  for n in 1:Nsites
    Cvalsₙ, Cvecsₙ, Cinfoₙ = eigsolve(
      Hᶜ(∑h, Hᴸ, Hᴿ, ψ, n), ψ.C[n], 1, :SR; ishermitian=true, tol=krylov_tol
    )
    Avalsₙ, Avecsₙ, Ainfoₙ = eigsolve(
      Hᴬᶜ(∑h, Hᴸ, Hᴿ, ψ, n), ψ.AL[n] * ψ.C[n], 1, :SR; ishermitian=true, tol=krylov_tol
    )
    C̃[n] = Cvecsₙ[1]
    Ãᶜ[n] = Avecsₙ[1]
  end

  # TODO: based on minimum singular values of C̃, use more accurate
  # method for finding Ãᴸ, Ãᴿ
  Ãᴸ = InfiniteMPS(Vector{ITensor}(undef, Nsites))
  Ãᴿ = InfiniteMPS(Vector{ITensor}(undef, Nsites))
  for n in 1:Nsites
    Ãᴸⁿ, X = polar(Ãᶜ[n] * dag(C̃[n]), uniqueinds(Ãᶜ[n], C̃[n]))
    Ãᴿⁿ, _ = polar(Ãᶜ[n] * dag(C̃[n - 1]), uniqueinds(Ãᶜ[n], C̃[n - 1]))
    Ãᴸⁿ = noprime(Ãᴸⁿ)
    Ãᴿⁿ = noprime(Ãᴿⁿ)
    Ãᴸ[n] = Ãᴸⁿ
    Ãᴿ[n] = Ãᴿⁿ
  end

  for n in 1:Nsites
    ϵᴸ![n] = norm(Ãᶜ[n] - Ãᴸ[n] * C̃[n])
    ϵᴿ![n] = norm(Ãᶜ[n] - C̃[n - 1] * Ãᴿ[n])
  end
  return InfiniteCanonicalMPS(Ãᴸ, C̃, Ãᴿ), (eᴸ, eᴿ)
end

function vumps(
  ∑h,
  ψ;
  maxiter=10,
  tol=1e-8,
  outputlevel=1,
  multisite_update_alg="sequential",
  eigsolve_tol=(x -> x / 100),
)
  N = nsites(ψ)
  (ϵᴸ!) = fill(tol, nsites(ψ))
  (ϵᴿ!) = fill(tol, nsites(ψ))
  outputlevel > 0 &&
    println("Running VUMPS with multisite_update_alg = $multisite_update_alg")
  for iter in 1:maxiter
    ψ, (eᴸ, eᴿ) = vumps_iteration(
      ∑h,
      ψ;
      (ϵᴸ!)=(ϵᴸ!),
      (ϵᴿ!)=(ϵᴿ!),
      multisite_update_alg=multisite_update_alg,
      eigsolve_tol=eigsolve_tol,
    )
    ϵᵖʳᵉˢ = max(maximum(ϵᴸ!), maximum(ϵᴿ!))
    maxdimψ = maxlinkdim(ψ[0:(N + 1)])
    outputlevel > 0 && println(
      "VUMPS iteration $iter (out of maximum $maxiter). Bond dimension = $maxdimψ, energy = $((eᴸ, eᴿ)), ϵᵖʳᵉˢ = $ϵᵖʳᵉˢ, tol = $tol",
    )
    if ϵᵖʳᵉˢ < tol
      println(
        "Precision error $ϵᵖʳᵉˢ reached tolerance $tol, stopping VUMPS after $iter iterations (of a maximum $maxiter).",
      )
      break
    end
  end
  return ψ
end

##################################################################
# Old functionality, only used for testing

## function left_environment_recursive(hᴸ, ψ; niter=10)
##   ψ̃ = prime(linkinds, dag(ψ))
##   # XXX: replace with `nsites`
##   #N = nsites(ψ)
##   N = length(ψ)
##   Hᴸᴺ¹ = hᴸ[N]
##   for _ in 1:niter
##     Hᴸᴺ¹ = translatecell(Hᴸᴺ¹, -1)
##     for n in 1:N
##       Hᴸᴺ¹ = Hᴸᴺ¹ * ψ.AL[n] * ψ̃.AL[n]
##     end
##     # Loop over the Hamiltonian terms in the unit cell
##     for n in 1:N
##       hᴸⁿ = hᴸ[n]
##       for k in (n + 1):N
##         hᴸⁿ = hᴸⁿ * ψ.AL[k] * ψ̃.AL[k]
##       end
##       Hᴸᴺ¹ += hᴸⁿ
##     end
##   end
##   # Get the rest of the environments in the unit cell
##   Hᴸ = InfiniteMPS(Vector{ITensor}(undef, N))
##   Hᴸ[N] = Hᴸᴺ¹
##   Hᴸᴺ¹ = translatecell(Hᴸᴺ¹, -1)
##   for n in 1:(N - 1)
##     Hᴸ[n] = Hᴸ[n - 1] * ψ.AL[n] * ψ̃.AL[n] + hᴸ[n]
##   end
##   return Hᴸ
## end
##
## function right_environment_recursive(hᴿ, ψ; niter=10)
##   ψ̃ = prime(linkinds, dag(ψ))
##   # XXX: replace with `nsites`
##   #N = nsites(ψ)
##   N = length(ψ)
##   Hᴿᴺ¹ = hᴿ[0]
##   for _ in 1:niter
##     Hᴿᴺ¹ = translatecell(Hᴿᴺ¹, 1)
##     for n in reverse(1:N)
##       Hᴿᴺ¹ = Hᴿᴺ¹ * ψ.AR[n] * ψ̃.AR[n]
##     end
##     # Loop over the Hamiltonian terms in the unit cell
##     for n in reverse(0:(N - 1))
##       hᴿⁿ = hᴿ[n]
##       for k in reverse(1:n)
##         hᴿⁿ = hᴿⁿ * ψ.AR[k] * ψ̃.AR[k]
##       end
##       Hᴿᴺ¹ += hᴿⁿ
##     end
##   end
##   Hᴿᴺ¹ = translatecell(Hᴿᴺ¹, 1)
##   # Get the rest of the environments in the unit cell
##   Hᴿ = InfiniteMPS(Vector{ITensor}(undef, N))
##   Hᴿ[N] = Hᴿᴺ¹
##   for n in reverse(1:(N - 1))
##     Hᴿ[n] = Hᴿ[n + 1] * ψ.AR[n + 1] * ψ̃.AR[n + 1] + hᴿ[n]
##   end
##   return Hᴿ
## end


text:
export param_first, param_double, param_increment, param_inbetween, hasparam_inbetween
"""
    param_first(platform::Platform)
"""
param_first(platform::Platform) = 1
"""
    param_double(platform::Platform, param)
"""
param_double(platform::Platform, param) =
    extensionsize(dictionary(platform, param))
"""
    param_increment(platform::Platform, param)
"""
param_increment(platform::Platform, n) = addone(n)
"""
    param_inbetween(platform::Platform, param1, param2)
"""
param_inbetween(platform::Platform, param1, param2) =
    inbetween(param1, param2)
"""
    hasparam_inbetween(platform::Platform, param1, param2, stoptolerance)
"""
hasparam_inbetween(platform::Platform, param1, param2, stoptolerance::Real) =
    hasparam_inbetween(param1, param2, stoptolerance)

addone(x::Number) = x+1
addone(x) = map(addone, x)
inbetween(n1::Int, n2::Int) = (n1+n2) >> 1
hasparam_inbetween(n1::Int, n2::Int, stoptolerance) = inbetween(n1, n2) != n1
inbetween(n1::T, n2::T) where {T<:Tuple} = map(inbetween, n1, n2)
hasparam_inbetween(n1::T, n2::T, stoptolerance) where {T<:Tuple} =
    any(map(hasparam_inbetween, n1, n2, Ref(stoptolerance))...)


text:
### Proximal gradient method
export ProxGradParams, fit!

type ProxGradParams<:AbstractParams
    stepsize::Float64 # initial stepsize
    max_iter::Int # maximum number of outer iterations
    inner_iter_X::Int # how many prox grad steps to take on X before moving on to Y (and vice versa)
    inner_iter_Y::Int # how many prox grad steps to take on Y before moving on to X (and vice versa)
    abs_tol::Float64 # stop if objective decrease upon one outer iteration is less than this * number of observations
    rel_tol::Float64 # stop if objective decrease upon one outer iteration is less than this * objective value
    min_stepsize::Float64 # use a decreasing stepsize, stop when reaches min_stepsize
end
function ProxGradParams(stepsize::Number=1.0; # initial stepsize
				        max_iter::Int=100, # maximum number of outer iterations
                inner_iter_X::Int=1, # how many prox grad steps to take on X before moving on to Y (and vice versa)
                inner_iter_Y::Int=1, # how many prox grad steps to take on Y before moving on to X (and vice versa)
                inner_iter::Int=1,
                abs_tol::Number=0.00001, # stop if objective decrease upon one outer iteration is less than this * number of observations
                rel_tol::Number=0.0001, # stop if objective decrease upon one outer iteration is less than this * objective value
				        min_stepsize::Number=0.01*stepsize) # stop if stepsize gets this small
    stepsize = convert(Float64, stepsize)
    inner_iter_X = max(inner_iter_X, inner_iter)
    inner_iter_Y = max(inner_iter_Y, inner_iter)
    return ProxGradParams(convert(Float64, stepsize),
                          max_iter,
                          inner_iter_X,
                          inner_iter_Y,
                          convert(Float64, abs_tol),
                          convert(Float64, rel_tol),
                          convert(Float64, min_stepsize))
end

### FITTING
function fit!(glrm::GLRM, params::ProxGradParams;
			  ch::ConvergenceHistory=ConvergenceHistory("ProxGradGLRM"),
			  verbose=true,
			  kwargs...)
	### initialization
	A = glrm.A # rename these for easier local access
	losses = glrm.losses
	rx = glrm.rx
	ry = glrm.ry
	X = glrm.X; Y = glrm.Y
    # check that we didn't initialize to zero (otherwise we will never move)
    if vecnorm(Y) == 0
    	Y = .1*randn(k,d)
    end
	k = glrm.k
    m,n = size(A)

    # find spans of loss functions (for multidimensional losses)
    yidxs = get_yidxs(losses)
    d = maximum(yidxs[end])
    # check Y is the right size
    if d != size(Y,2)
        warn("The width of Y should match the embedding dimension of the losses.
            Instead, embedding_dim(glrm.losses) = $(embedding_dim(glrm.losses))
            and size(glrm.Y, 2) = $(size(glrm.Y, 2)).
            Reinitializing Y as randn(glrm.k, embedding_dim(glrm.losses).")
            # Please modify Y or the embedding dimension of the losses to match,
            # eg, by setting `glrm.Y = randn(glrm.k, embedding_dim(glrm.losses))`")
        glrm.Y = randn(glrm.k, d)
    end

    XY = @compat Array{Float64}((m, d))
    gemm!('T','N',1.0,X,Y,0.0,XY) # XY = X' * Y initial calculation

    # step size (will be scaled below to ensure it never exceeds 1/\|g\|_2 or so for any subproblem)
    alpharow = params.stepsize*ones(m)
    alphacol = params.stepsize*ones(n)
    # stopping criterion: stop when decrease in objective < tol, scaled by the number of observations
    scaled_abs_tol = params.abs_tol * mapreduce(length,+,glrm.observed_features)

    # alternating updates of X and Y
    if verbose println("Fitting GLRM") end
    update_ch!(ch, 0, objective(glrm, X, Y, XY, yidxs=yidxs))
    t = time()
    steps_in_a_row = 0
    # gradient wrt columns of X
    g = zeros(k)
    # gradient wrt column-chunks of Y
    G = zeros(k, d)
    # rowwise objective value
    obj_by_row = zeros(m)
    # columnwise objective value
    obj_by_col = zeros(n)

    # cache views for better memory management
    # make sure we don't try to access memory not allocated to us
    @assert(size(Y) == (k,d))
    @assert(size(X) == (k,m))
    # views of the columns of X corresponding to each example
    ve = [view(X,:,e) for e=1:m]
    # views of the column-chunks of Y corresponding to each feature y_j
    # vf[f] == Y[:,f]
    vf = [view(Y,:,yidxs[f]) for f=1:n]
    # views of the column-chunks of G corresponding to the gradient wrt each feature y_j
    # these have the same shape as y_j
    gf = [view(G,:,yidxs[f]) for f=1:n]

    # working variables
    newX = copy(X)
    newY = copy(Y)
    newve = [view(newX,:,e) for e=1:m]
    newvf = [view(newY,:,yidxs[f]) for f=1:n]

    for i=1:params.max_iter
# STEP 1: X update
        # XY = X' * Y was computed above

        # reset step size if we're doing something more like alternating minimization
        if params.inner_iter_X > 1 || params.inner_iter_Y > 1
            for ii=1:m alpharow[ii] = params.stepsize end
            for jj=1:n alphacol[jj] = params.stepsize end
        end

        for inneri=1:params.inner_iter_X
        for e=1:m # for every example x_e == ve[e]
            scale!(g, 0) # reset gradient to 0
            # compute gradient of L with respect to Xᵢ as follows:
            # ∇{Xᵢ}L = Σⱼ dLⱼ(XᵢYⱼ)/dXᵢ
            for f in glrm.observed_features[e]
                # but we have no function dLⱼ/dXᵢ, only dLⱼ/d(XᵢYⱼ) aka dLⱼ/du
                # by chain rule, the result is: Σⱼ (dLⱼ(XᵢYⱼ)/du * Yⱼ), where dLⱼ/du is our grad() function
                curgrad = grad(losses[f],XY[e,yidxs[f]],A[e,f])
                if isa(curgrad, Number)
                    axpy!(curgrad, vf[f], g)
                else
                    # on v0.4: gemm!('N', 'T', 1.0, vf[f], curgrad, 1.0, g)
                    gemm!('N', 'N', 1.0, vf[f], curgrad, 1.0, g)
                end
            end
            # take a proximal gradient step to update ve[e]
            l = length(glrm.observed_features[e]) + 1 # if each loss function has lipshitz constant 1 this bounds the lipshitz constant of this example's objective
            obj_by_row[e] = row_objective(glrm, e, ve[e]) # previous row objective value
            while alpharow[e] > params.min_stepsize
                stepsize = alpharow[e]/l
                # newx = prox(rx[e], ve[e] - stepsize*g, stepsize) # this will use much more memory than the inplace version with linesearch below
                ## gradient step: Xᵢ += -(α/l) * ∇{Xᵢ}L
                axpy!(-stepsize,g,newve[e])
                ## prox step: Xᵢ = prox_rx(Xᵢ, α/l)
                prox!(rx[e],newve[e],stepsize)
                if row_objective(glrm, e, newve[e]) < obj_by_row[e]
                    copy!(ve[e], newve[e])
                    alpharow[e] *= 1.05
                    break
                else # the stepsize was too big; undo and try again only smaller
                    copy!(newve[e], ve[e])
                    alpharow[e] *= .7
                    if alpharow[e] < params.min_stepsize
                        alpharow[e] = params.min_stepsize * 1.1
                        break
                    end
                end
            end
        end # for e=1:m
        gemm!('T','N',1.0,X,Y,0.0,XY) # Recalculate XY using the new X
        end # inner iteration
# STEP 2: Y update
        for inneri=1:params.inner_iter_Y
        scale!(G, 0)
        for f=1:n
            # compute gradient of L with respect to Yⱼ as follows:
            # ∇{Yⱼ}L = Σⱼ dLⱼ(XᵢYⱼ)/dYⱼ
            for e in glrm.observed_examples[f]
                # but we have no function dLⱼ/dYⱼ, only dLⱼ/d(XᵢYⱼ) aka dLⱼ/du
                # by chain rule, the result is: Σⱼ dLⱼ(XᵢYⱼ)/du * Xᵢ, where dLⱼ/du is our grad() function
                curgrad = grad(losses[f],XY[e,yidxs[f]],A[e,f])
                if isa(curgrad, Number)
                    axpy!(curgrad, ve[e], gf[f])
                else
                    # on v0.4: gemm!('N', 'T', 1.0, ve[e], curgrad, 1.0, gf[f])
                    gemm!('N', 'T', 1.0, ve[e], curgrad, 1.0, gf[f])
                end
            end
            # take a proximal gradient step
            l = length(glrm.observed_examples[f]) + 1
            obj_by_col[f] = col_objective(glrm, f, vf[f])
            while alphacol[f] > params.min_stepsize
                stepsize = alphacol[f]/l
                # newy = prox(ry[f], vf[f] - stepsize*gf[f], stepsize)
                ## gradient step: Yⱼ += -(α/l) * ∇{Yⱼ}L
                axpy!(-stepsize,gf[f],newvf[f])
                ## prox step: Yⱼ = prox_ryⱼ(Yⱼ, α/l)
                prox!(ry[f],newvf[f],stepsize)
                new_obj_by_col = col_objective(glrm, f, newvf[f])
                if new_obj_by_col < obj_by_col[f]
                    copy!(vf[f], newvf[f])
                    alphacol[f] *= 1.05
                    obj_by_col[f] = new_obj_by_col
                    break
                else
                    copy!(newvf[f], vf[f])
                    alphacol[f] *= .7
                    if alphacol[f] < params.min_stepsize
                        alphacol[f] = params.min_stepsize * 1.1
                        break
                    end
                end
            end
        end # for f=1:n
        gemm!('T','N',1.0,X,Y,0.0,XY) # Recalculate XY using the new Y
        end # inner iteration
# STEP 3: Record objective
        obj = sum(obj_by_col)
        t = time() - t
        update_ch!(ch, t, obj)
        t = time()
# STEP 4: Check stopping criterion
        obj_decrease = ch.objective[end-1] - obj
        if i>10 && (obj_decrease < scaled_abs_tol || obj_decrease/obj < params.rel_tol)
            break
        end
        if verbose && i%10==0
            println("Iteration $i: objective value = $(ch.objective[end])")
        end
    end

    return glrm.X, glrm.Y, ch
end


text:
@testset "Points" begin
  @test embeddim(Point(1, 2)) == 2
  @test embeddim(Point(1, 2, 3)) == 3
  @test coordtype(Point(1, 1)) == Int
  @test coordtype(Point(1.,1.)) == Float64
  @test coordtype(Point(1f0, 1f0)) == Float32
  @test coordtype(Point2(1, 1)) == Float64
  @test coordtype(Point3(1, 1, 1)) == Float64
  @test coordtype(Point2f(1, 1)) == Float32
  @test coordtype(Point3f(1, 1, 1)) == Float32

  @test coordtype(Point{2,T}((1, 1))) == T
  @test coordtype(Point{2,T}(1, 1)) == T

  @test coordinates(P2(1, 2)) == [1.0, 2.0]
  @test coordinates(P3(1, 2, 3)) == [1.0, 2.0, 3.0]

  @test P2(1, 2) - P2(1, 1) == T[0, 1]
  @test P3(1, 2, 3) - P3(1, 1, 1) == T[0, 1, 2]
  @test_throws DimensionMismatch P2(1, 2) - P3(1, 2, 3)

  @test P2(1, 2) + V2(0, 0) == P2(1, 2)
  @test P2(2, 3) + V2(2, 1) == P2(4, 4)
  @test P3(1, 2, 3) + V3(0, 0, 0) == P3(1, 2, 3)
  @test P3(2, 3, 4) + V3(2, 1, 0) == P3(4, 4, 4)
  @test_throws DimensionMismatch P2(1, 2) + V3(1, 2, 3)

  @test P2(1, 2) - V2(0, 0) == P2(1, 2)
  @test P2(2, 3) - V2(2, 1) == P2(0, 2)
  @test P3(1, 2, 3) - V3(0, 0, 0) == P3(1, 2, 3)
  @test P3(2, 3, 4) - V3(2, 1, 0) == P3(0, 2, 4)

  @test embeddim(rand(P2)) == 2
  @test embeddim(rand(P3)) == 3
  @test coordtype(rand(P2)) == T
  @test coordtype(rand(P3)) == T

  @test eltype(rand(P2, 3)) == P2
  @test eltype(rand(P3, 3)) == P3

  @test P2(1, 2) ≈ P2(1 + eps(T), T(2))
  @test P3(1, 2, 3) ≈ P3(1 + eps(T), T(2), T(3))

  @test embeddim(Point([1,2])) == 2
  @test coordtype(Point([1,2])) == Int
  @test coordtype(Point([1.,2.])) == Float64

  @test embeddim(Point([1,2,3])) == 3
  @test coordtype(Point([1,2,3])) == Int
  @test coordtype(Point([1.,2.,3.])) == Float64
end


text:
# Autogenerated wrapper script for NOVAS_jll for x86_64-apple-darwin
export libnovas

JLLWrappers.@generate_wrapper_header("NOVAS")
JLLWrappers.@declare_library_product(libnovas, "@rpath/libnovas.3.dylib")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libnovas,
        "lib/libnovas.3.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
import Base.==

@enum(Key,
    BACKSPACE = (@static Sys.iswindows() ? 8 : 127),
    ARROW_LEFT = 1000,
    ARROW_RIGHT,
    ARROW_UP,
    ARROW_DOWN,
    DEL_KEY,
    HOME_KEY,
    END_KEY,
    PAGE_UP,
    PAGE_DOWN,
    S_ARROW_UP,
    S_ARROW_DOWN,
    S_ARROW_LEFT,
    S_ARROW_RIGHT,
    C_ARROW_UP,
    C_ARROW_DOWN,
    C_ARROW_LEFT,
    C_ARROW_RIGHT)

Base.convert(::Type{UInt32}, x::Key) = UInt32(x)
==(c::UInt32, k::Key) = c == UInt32(k)
==(k::Key, c::UInt32) = c == UInt32(k)
==(c::Char, k::Key) = UInt32(c) == UInt32(k)
==(k::Key, c::Char) = UInt32(c) == UInt32(k)

ctrl_key(c::Char)::UInt32 = UInt32(c) & 0x1f

# For debugging
function printNextKey()
	term = REPL.Terminals.TTYTerminal(get(ENV, "TERM", @static Sys.iswindows() ? "" : "dumb"), stdin, stdout, stderr)
	REPL.Terminals.raw!(term, true)
	c = readNextChar()
	print("Code: $(UInt32(c)), Char: $(Char(c))")
	REPL.Terminals.raw!(term, true)
	return nothing
end

readNextChar() = Char(read(stdin,1)[1])

function readKey() ::UInt32
    c = readNextChar()

    # Escape characters
    if c == '\x1b'
        stdin.buffer.size < 3 && return '\x1b'
        esc_a = readNextChar()
        esc_b = readNextChar()

        if esc_a == '['
            if esc_b >= '0' && esc_b <= '9'
                stdin.buffer.size < 4 && return '\x1b'
                esc_c = readNextChar()

                if esc_c == '~'
                    if esc_b == '1'
                        return HOME_KEY
                    elseif esc_b == '4'
                        return END_KEY
                    elseif esc_b == '3'
                        return DEL_KEY
                    elseif esc_b == '5'
                        return PAGE_UP
                    elseif esc_b == '6'
                        return PAGE_DOWN
                    elseif esc_b == '7'
                        return HOME_KEY
                    elseif esc_b == '8'
                        return END_KEY
                    else
                        return '\x1b'
                    end
                elseif esc_c == ';'
                    stdin.buffer.size < 6 && return '\x1b'
                    esc_d = readNextChar()
                    esc_e = readNextChar()

                    if esc_d == '2'
                        # shift + arrorw
                        if esc_e == 'A'
                            return S_ARROW_UP
                        elseif esc_e == 'B'
                            return S_ARROW_DOWN
                        elseif esc_e == 'C'
                            return S_ARROW_RIGHT
                        elseif esc_e == 'D'
                            return S_ARROW_LEFT
                        else
                            return '\x1b'
                        end
                    elseif esc_d == '5'
                        # Ctrl + arrow
                        if esc_e == 'A'
                            return C_ARROW_UP
                        elseif esc_e == 'B'
                            return C_ARROW_DOWN
                        elseif esc_e == 'C'
                            return C_ARROW_RIGHT
                        elseif esc_e == 'D'
                            return C_ARROW_LEFT
                        else
                            return '\x1b'
                        end
                    end
                end
            else
                # Arrow keys
                if esc_b == 'A'
                    return ARROW_UP
                elseif esc_b == 'B'
                    return ARROW_DOWN
                elseif esc_b == 'C'
                    return ARROW_RIGHT
                elseif esc_b == 'D'
                    return ARROW_LEFT
                elseif esc_b == 'H'
                    return HOME_KEY
                elseif esc_b == 'F'
                    return END_KEY
                else
                    return '\x1b'
                end
            end
        elseif esc_a == 'O'
            if esc_a == 'H'
                return HOME_KEY
            elseif esc_a == 'F'
                return END_KEY
            end
        end

        return '\x1b'
    else
        return c;
    end
end


text:
using DoubleExponentialFormulas
using LinearAlgebra: norm
using Test

quadss32 = QuadSS(Float32)
quadss64 = QuadSS(Float64)
quadssBF = QuadSS(BigFloat)

# Test Base.show works without error
let
    io = IOBuffer()
    @test (Base.show(io, MIME("text/plain"), quadss32); true)
    @test (Base.show(io, MIME("text/plain"), quadss64); true)
    @test (Base.show(io, MIME("text/plain"), quadssBF); true)
end

# Test (-∞, ∞) interval (with Sinh-Sinh quadrature)
# Gauss integral
let
    f(x::AbstractFloat) = exp(-x^2)
    expect = sqrt(BigFloat(π))

    I, E = quadss32(f)
    @test I isa Float32
    @test I ≈ expect
    @test E ≤ sqrt(eps(I))*norm(I)

    I, E = quadss64(f)
    @test I isa Float64
    @test I ≈ expect
    @test E ≤ sqrt(eps(I))*norm(I)

    rtol = 1e-30
    I, E = quadssBF(f, rtol=rtol)
    @test I isa BigFloat
    @test isapprox(I, expect, rtol=10rtol)
    @test E ≤ rtol*norm(I)
end

# Test (-∞, ∞) interval (with Sinh-Sinh quadrature)
# Gauss integral
let
    f(x::AbstractFloat) = x^2*exp(-3*x^2)
    expect = sqrt(BigFloat(π)/3)/(2*3)

    I, E = quadss32(f)
    @test I isa Float32
    @test I ≈ expect
    @test E ≤ sqrt(eps(I))*norm(I)

    I, E = quadss64(f)
    @test I isa Float64
    @test I ≈ expect
    @test E ≤ sqrt(eps(I))*norm(I)

    rtol = 1e-30
    I, E = quadssBF(f, rtol=rtol)
    @test I isa BigFloat
    @test isapprox(I, expect, rtol=10rtol)
    @test E ≤ rtol*norm(I)
end

# Test (-∞, ∞) interval (with Sinh-Sinh quadrature)
let
    f(x::AbstractFloat) = 1/(1 + x^2)
    expect = π

    I, E = quadss32(f)
    @test I isa Float32
    @test I ≈ expect
    @test E ≤ sqrt(eps(I))*norm(I)

    I, E = quadss64(f)
    @test I isa Float64
    @test I ≈ expect
    @test E ≤ sqrt(eps(I))*norm(I)

    rtol = 1e-30
    I, E = quadssBF(f, rtol=rtol)
    @test I isa BigFloat
    @test isapprox(I, expect, rtol=10rtol)
    @test E ≤ rtol*norm(I)
end

# Test non-scalar output with QuadSS
let
    f(x::AbstractFloat) = [exp(-x^2), 2*exp(-x^2)]
    expect = [sqrt(BigFloat(π)), 2*sqrt(BigFloat(π))]

    I, E = quadss32(f)
    @test eltype(I) == Float32
    @test I ≈ expect
    @test E ≤ sqrt(eps(eltype(I)))*norm(I)

    I, E = quadss64(f)
    @test eltype(I) == Float64
    @test I ≈ expect
    @test E ≤ sqrt(eps(eltype(I)))*norm(I)

    rtol = 1e-30
    I, E = quadssBF(f, rtol=rtol)
    @test eltype(I) == BigFloat
    @test isapprox(I, expect, rtol=10rtol)
    @test E ≤ rtol*norm(I)
end


text:
module YAXArrays
global const YAXDefaults = (
  workdir = Ref("./"),
  recal   = Ref(false),
  chunksize  = Ref{Any}(:input),
  max_cache  = Ref(1e8),
  cubedir    = Ref(""),
  subsetextensions = [],
)
global const workdir=YAXDefaults.workdir
global const recal=YAXDefaults.recal
function __init__()
  YAXDefaults.workdir[]   = get(ENV,"YAXARRAY_WORKDIR","./")
  YAXDefaults.max_cache[] = parse(Float64,get(ENV,"YAXARRAY_MAX_CACHE","100")) * 1e6
  YAXDefaults.cubedir[]   = get(ENV,"YAXARRAY_CUBEDIR","")
end
YAXdir(x::String)=YAXDefaults.workdir[]=x
recalculate(x::Bool)=YAXDefaults.recal[]=x
recalculate()=YAXDefaults.recal[]
YAXdir()=YAXDefaults.workdir[]
export YAXdir

include("YAXTools.jl")
include("Cubes/Cubes.jl")
include("DatasetAPI/Datasets.jl")
include("DAT/DAT.jl")

using .YAXTools: @reexport
using YAXArrayBase: getattributes

@reexport using Dates: Date, DateTime
@reexport using IntervalSets: (..)
@reexport using .Cubes: concatenatecubes, caxes,
  subsetcube, readcubedata,renameaxis!, YAXArray
@reexport using .Cubes.Axes: CubeAxis, RangeAxis, CategoricalAxis,
  getAxis

@reexport using .DAT: mapCube, getAxis, InDims, OutDims, Dataset,
      CubeTable, cubefittable, fittable, savecube, loadcube, rmcube #From DAT module
@reexport using .Datasets: Dataset, Cube, open_dataset
@reexport using .YAXTools: @loadOrGenerate # from YAXTools

@deprecate saveCube(data, filename) savecube(data,filename)
@deprecate loadCube(filename) loadcube(filename)
@deprecate rmCube(filename) rmcube(filename)
@deprecate exportcube(data, filename; kwargs...) savecube(data, filename; backend=:netcdf, kwargs...)
@deprecate cubeproperties(x) getattributes(x)
@deprecate concantenateCubes(args...) concatenatecubes(args...)

#include("precompile.jl")
#_precompile_()

end # module


text:
k = Knots([1,2,2,3])

@testset "Knots" begin
    @test Knots([1,2,3]) == Knots(1:3) == Knots(1,2,3)
    @test Knots(2,4,5) == Knots([2,4,5])
    @test zero(Knots) == Knots() == Knots([])
    @test Knots(1:3) == Knots([3,2,1])

    @test ♯(k) == length(k) == 4

    @test Knots([-1,2,3]) + 2 * Knots([2,5]) == Knots([-1,2,2,2,3,5,5])
    @test Knots([1,2,3]) + Knots([2,4,5]) == Knots([1,2,2,3,4,5])
    @test 2 * Knots([2,3]) == Knots([2,2,3,3])

    @test unique(k) == Knots([1,2,3])

    @testset "inclusive relation" begin
        @test Knots([1,2,3])     ⊆ Knots([1,2,3])
        @test Knots([1,2,3])     ⊇ Knots([1,2,3])
        @test Knots([1,2,2,3])   ⊆ Knots([1,2,2,3,5])
        @test Knots([1,2,2,3])   ⊈ Knots([1,2,3,5])
        @test Knots([1,2,2,3,5]) ⊇ Knots([1,2,2,3])
        @test Knots([1,2,3,5])   ⊉ Knots([1,2,2,3])
    end

    @test 𝔫(k, 0.3) == 0
    @test 𝔫(k, 1.0) == 1
    @test 𝔫(k, 2.0) == 2
    @test 1 ∈ k
    @test 1.5 ∉ k

    @testset "knotindex for usual case" begin
        @test BasicBSpline._knotindex(k,1.5) == 1
        @test BasicBSpline._knotindex(k,2.5) == 3
        @test BasicBSpline._knotindex₋₀(k,1.5) == 1
        @test BasicBSpline._knotindex₋₀(k,2.5) == 3
        @test BasicBSpline._knotindex₊₀(k,1.5) == 1
        @test BasicBSpline._knotindex₊₀(k,2.5) == 3
    end

    @testset "knotindex for corner case" begin
        @test BasicBSpline._knotindex(k,0) == 0
        @test BasicBSpline._knotindex(k,1) == 1
        @test BasicBSpline._knotindex(k,2) == 3
        @test BasicBSpline._knotindex(k,3) == 3
        @test BasicBSpline._knotindex(k,4) == 4

        @test BasicBSpline._knotindex₋₀(k,0) == 0
        @test BasicBSpline._knotindex₋₀(k,1) == 0
        @test BasicBSpline._knotindex₋₀(k,2) == 1
        @test BasicBSpline._knotindex₋₀(k,3) == 3
        @test BasicBSpline._knotindex₋₀(k,4) == 4

        @test BasicBSpline._knotindex₊₀(k,0) == 0
        @test BasicBSpline._knotindex₊₀(k,1) == 1
        @test BasicBSpline._knotindex₊₀(k,2) == 3
        @test BasicBSpline._knotindex₊₀(k,3) == 4
        @test BasicBSpline._knotindex₊₀(k,4) == 4
    end

    @test string(k) == "Knots([1.0, 2.0, 2.0, 3.0])"
    @test string(Knots()) == "Knots([])"
end


text:
using Test
@time @testset "Mesh core" begin
    include("test_core.jl")
    include("attributes_tests.jl")
end


text:
using Parameters, StatsBase

"""
決定木関連
"""

function branch(x::Matrix,y::Vector,S,rf=0;loss_f=TSS)
    rf==0 && (m=size(x)[2])
    size(x)[2]==0 && return [0,0,0,0,0,0,0]
    best_score = Inf; i_1=1; j_1=1;
    left_1=[];right_1=[]
    left_score_1=Inf; right_score_1=Inf
    for j in 1:m
        for i in S
            left=[k for k in S if x[k,j]<x[i,j]]
            right=setdiff(S,left)
            left_score=loss_f(y[left]); right_score=loss_f(y[right])
            score=left_score + right_score
            if score<best_score
                best_score = score
                i_1=i; j_1=j
                left_1=left;right_1=right
                left_score_1=left_score;right_score_1=right_score
            end
        end
    end
    return i_1,j_1,left_1,right_1,best_score,left_score_1,right_score_1
end

mutable struct Stack
    parent::Union{Int,Vector,AbstractRange} #親ノード番号
    set::Union{Vector,AbstractRange} #教師データの行
    score::Number # 損失関数の値
end

mutable struct Node
    parent::Union{Int,Vector,AbstractRange} #親ノード番号
    j::Int #列番号
    th::Number #閾値
    set::Union{Vector,AbstractRange} #そのノードでの教師データの集合
    left::Int # 分岐後の列番号(左)
    right::Int # 分岐後の列番号(右)
    center::Number # yの平均値
end

function decisiontree(x::Matrix,y::Vector;α=0.0,n_min=1,rf=0,loss_f=TSS)
    rf==0 && (m=size(x)[2])
    stack = [Stack(0,1:size(x)[1],loss_f(y))]
    node=Node[]
    k=0 #Node番号
    while length(stack)>0
        popped=pop!(stack)
        k+=1　#Node番号を増やす。
        i,j,left,right,score,left_score,right_score = branch(x,y,popped.set,rf,;loss_f=loss_f)
        if popped.score-score<α || length(popped.set)<n_min || length(left)==0 || length(right)==0
            push!(node,Node(popped.parent,-1,0,popped.set,0,0,0))
        else
            push!(node,Node(popped.parent,j,x[i,j],popped.set,0,0,0))
            push!(stack,Stack(k,right,right_score))
            push!(stack,Stack(k,left,left_score))
        end
    end
    for h in k:-1:2
        pa=node[h].parent
        if node[pa].right==0
            node[pa].right=h
        else
            node[pa].left = h
        end
    end

    g = loss_f==TSS ? mean : mode_max

    for h in 1:k
        node[h].center =  node[h].j==-1 ? g(y[node[h].set]) : 0
    end
    return node
end

function dtvalue(u,node)
    r=1
    while node[r].j !=-1
        if u[node[r].j] < node[r].th
            r = node[r].left
        else
            r = node[r].right
        end
    end
    return node[r].center
end

function freq(y)
    S = Set(y) |> collect |> sort
    [count(==(i),y) for i in S]
end

mode(y) = length(y)==0 ? 0 : maximum(freq(y))

mis_match(y) = length(y) - mode(y)

function gini(y)
    n = length(y)
    n==0 && (return 0)
    fr = freq(y)
    sum(x->x*(n-x), fr)
end

function entropy(y)
    n = length(y)
    n==0 && (return 0)
    fr = freq(y)
    sum(x->-x*log(x/n), fr)
end

function mode_max(y)
    length(y)==0 && (return -Inf)
    counts(y) |> argmax
end


text:
module TestTables

using Test, Tables, DataFrames

struct NamedTupleIterator{T <: NamedTuple}
    elements::Vector{T}
end
Base.length(n::NamedTupleIterator) = length(n.elements)
Base.IteratorEltype(::Type{<:NamedTupleIterator}) = Base.HasEltype()
Base.eltype(n::NamedTupleIterator{T}) where {T} = T

function Base.iterate(nt::NamedTupleIterator, st=1)
    st > length(nt.elements) && return nothing
    return nt.elements[st], st + 1
end

struct EltypeUnknownIterator
    elements
end
Base.IteratorSize(::Type{EltypeUnknownIterator}) = Base.SizeUnknown()
Base.IteratorEltype(::Type{EltypeUnknownIterator}) = Base.EltypeUnknown()

function Base.iterate(nt::EltypeUnknownIterator, st=1)
    st > length(nt.elements) && return nothing
    return nt.elements[st], st + 1
end

struct DuplicateNamesTable
end
Tables.istable(::Type{DuplicateNamesTable}) = true
Tables.rowaccess(::Type{DuplicateNamesTable}) = true
Tables.rows(x::DuplicateNamesTable) = x
Tables.schema(x::DuplicateNamesTable) = Tables.Schema((:a, :a, :b), Tuple{Float64, Float64, Float64})

Base.length(x::DuplicateNamesTable) = 3
Base.eltype(x::DuplicateNamesTable) = DuplicateRow

function Base.iterate(d::DuplicateNamesTable, st=1)
    st > length(d) && return nothing
    return DuplicateNameRow(), st + 1
end

struct DuplicateNameRow
end
Base.getproperty(d::DuplicateNameRow, nm::Symbol) = 1.0

struct DuplicateNamesColumnTable
end
Tables.istable(::Type{DuplicateNamesColumnTable}) = true
Tables.columnaccess(::Type{DuplicateNamesColumnTable}) = true
Tables.columns(x::DuplicateNamesColumnTable) = x
Tables.schema(x::DuplicateNamesColumnTable) = Tables.Schema((:a, :a, :b), Tuple{Float64, Float64, Float64})
Base.getproperty(d::DuplicateNamesColumnTable, nm::Symbol) = [1.0, 2.0, 3.0]
Base.propertynames(d::DuplicateNamesColumnTable) = (:a, :a, :b)

@testset "Tables" begin
    df = DataFrame(a=Int64[1, 2, 3], b=[:a, :b, :c])

    @testset "basics" begin
        @test Tables.istable(df)
        @test Tables.rowaccess(df)
        @test Tables.columnaccess(df)
        @test Tables.schema(df) === Tables.Schema((:a, :b), Tuple{Int64, Symbol})
        @test Tables.schema(df) == Tables.schema(Tables.rows(df)) == Tables.schema(Tables.columns(df))
        @test @inferred(Tables.materializer(df)(Tables.columns(df))) isa typeof(df)

        row = first(Tables.rows(df))
        @test propertynames(row) == (:a, :b)
        @test getproperty(row, :a) == 1
        @test getproperty(row, :b) == :a
    end

    @testset "Row-style" begin
        bare_rows = Tables.rowtable(df)
        for (actual, expected) in zip(bare_rows, eachrow(df))
            @test actual.a == expected.a
            @test actual.b == expected.b
        end

        and_back = DataFrame(bare_rows)
        @test and_back isa DataFrame
        @test names(and_back) == [:a, :b]
        @test and_back.a == df.a
        @test and_back.b == df.b
    end

    @testset "Column-style" begin
        cols = Tables.columntable(df)
        @test cols.b  ==  df.b
        @test cols.a  ==  df.a

        and_back = DataFrame(cols)
        @test and_back isa DataFrame
        @test names(and_back) == [:a, :b]
        @test and_back.a == df.a
        @test and_back.b == df.b
    end

    @testset "Extras" begin
        # with missing values
        df = DataFrame(a=[1, missing, 3], b=[missing, 'a', "hey"])
        @test isequal(df, DataFrame(Tables.rowtable(df)))
        @test isequal(df, DataFrame(Tables.columntable(df)))

        dn = DuplicateNamesTable()
        @test_throws ErrorException (dn |> DataFrame)

        dn = DuplicateNamesColumnTable()
        @test_throws ArgumentError (dn |> DataFrame)

        # non-Tables.jl constructor fallbacks
        @test DataFrame([(a = 0,), (a = 1,)]) == DataFrame(a = 0:1)

        nt = (a=1, b=:a, c=missing)
        nti = NamedTupleIterator([nt, nt, nt])
        df = DataFrame(nti)
        @test size(df) == (3, 3)
        @test df.a == [1, 1, 1]
        @test df.b == [:a, :a, :a]
        @test isequal(df.c, [missing, missing, missing])

        etu = EltypeUnknownIterator([nt, nt, nt])
        df = DataFrame(etu)
        @test size(df) == (3, 3)
        @test df.a == [1, 1, 1]
        @test df.b == [:a, :a, :a]
        @test isequal(df.c, [missing, missing, missing])

        append!(df, etu)
        @test size(df) == (6, 3)

        # categorical values
        cat = CategoricalVector(["hey", "there", "sailor"])
        cat2 = [c for c in cat] # Vector of CategoricalString
        nt = (a=cat, b=cat2)
        df = DataFrame(nt, copycols=false)
        @test df.a === cat
        @test df.b === cat2
        # test in the unknown schema case that a
        # Vector of CategoricalString is built into CategoricalVector
        ct = Tables.buildcolumns(nothing, Tables.rows(nt))
        @test ct.a !== cat
        @test ct.b !== cat2
        @test ct.a == cat
        @test ct.b == cat == cat2
    end
end

@testset "DataFrame!" begin
    nt = (a=Int64[1, 2, 3], b=[:a, :b, :c])
    df1 = DataFrame!(nt)
    df2 = DataFrame!(df1)
    df3 = DataFrame(nt)
    @test Tables.columntable(df1) === nt
    @test Tables.columntable(df2) === nt
    @test Tables.columntable(df3) == nt
    @test Tables.columntable(df3) !== nt

    v = [(a=1,b=2), (a=3, b=4)]
    df = DataFrame(v)
    @test size(df) == (2, 2)
    @test df.a == [1, 3]
    @test df.b == [2, 4]
    @test DataFrame(v, copycols=false) == DataFrame(v, copycols=true) == df
    @test_throws ArgumentError DataFrame!(v)
end

@testset "columnindex" begin
    df = DataFrame(rand(3,4))
    @test columnindex.(Ref(df), names(df)) == 1:4
    @test columnindex(df, :a) == 0
    @test_throws ErrorException columnindex(df, 1)
    @test_throws ErrorException columnindex(df, "x1")
end

@testset "eachrow and eachcol integration" begin
     df = DataFrame(rand(3,4), [:a, :b, :c, :d])

     df2 = DataFrame(eachrow(df))
     @test df == df2
     @test !any(((a,b),) -> a === b, zip(eachcol(df), eachcol(df2)))

     df2 = DataFrame!(eachrow(df))
     @test df == df2
     @test !any(((a,b),) -> a === b, zip(eachcol(df), eachcol(df2)))

     df2 = DataFrame(eachcol(df, true))
     @test df == df2
     @test !any(((a,b),) -> a === b, zip(eachcol(df), eachcol(df2)))

     df2 = DataFrame!(eachcol(df, true))
     @test df == df2
     @test all(((a,b),) -> a === b, zip(eachcol(df), eachcol(df2)))

     df2 = DataFrame(eachcol(df))
     @test names(df2) == [:x1, :x2, :x3, :x4]
     @test all(((a,b),) -> a == b, zip(eachcol(df), eachcol(df2)))
     @test !any(((a,b),) -> a === b, zip(eachcol(df), eachcol(df2)))

     df2 = DataFrame(eachcol(df))
     @test names(df2) == [:x1, :x2, :x3, :x4]
     @test !any(((a,b),) -> a === b, zip(eachcol(df), eachcol(df2)))
end

end # module