text:
################################################################################
#
#  (Generalized) norm relations
#
################################################################################

# Example
#
# julia> f = x^8+8*x^6+64*x^4-192*x^2+576
# x^8+8*x^6+64*x^4-192*x^2+576
#
# julia> K, a = number_field(f);
#
# julia> N  = Hecke._norm_relation_setup(K)
# Norm relation of
#   Number field over Rational Field with defining polynomial x^8+8*x^6+64*x^4-192*x^2+576
# of index 4 and the subfields
#   Number field over Rational Field with defining polynomial x^2+1
#   Number field over Rational Field with defining polynomial x^2-6
#   Number field over Rational Field with defining polynomial x^2-2
#   Number field over Rational Field with defining polynomial x^2+2
#   Number field over Rational Field with defining polynomial x^2+6
#   Number field over Rational Field with defining polynomial x^2-3
#   Number field over Rational Field with defining polynomial x^2-x+1

mutable struct NormRelation{T}
  K::AnticNumberField
  subfields::Vector{Tuple{AnticNumberField, NfToNfMor}}
  denominator::Int
  coefficients::Vector{Vector{Tuple{NfToNfMor, Int}}}
  ispure::Bool
  pure_coefficients::Vector{Int}
  isabelian::Bool
  isgeneric::Bool
  coefficients_gen::Vector{Vector{Tuple{Int, NfToNfMor, NfToNfMor}}}
  embed_cache::Dict{Tuple{Int, Int}, Dict{nf_elem, nf_elem}}
  mor_cache::Dict{NfToNfMor, Dict{nf_elem, nf_elem}}
  induced::Dict{NfToNfMor, Perm{Int}}
  embed_cache_triv::Vector{Dict{nf_elem, nf_elem}}
  nonredundant::Vector{Int}
  isoptimal::Int
  isnormal::BitArray{1} # if i-th subfield is normal

  function NormRelation{T}() where {T}
    z = new{T}()
    z.ispure = false
    z.isabelian = false
    z.isgeneric = false
    z.denominator = 0
    z.embed_cache = Dict{Tuple{Int, Int}, Dict{nf_elem, nf_elem}}()
    z.mor_cache = Dict{NfToNfMor, Dict{nf_elem, nf_elem}}()
    z.induced = Dict{NfToNfMor, Perm{Int}}()
    z.isoptimal = -1
    return z
  end
end

function Base.show(io::IO, N::NormRelation)
  if !isdefined(N, :K)
    print(io, "Invalid norm relation")
    return
  end

  print(io, "Norm relation of\n  ", N.K, "\nof index ", index(N), " and the subfields")
  for i in 1:length(N)
    print(io, "\n  ", subfields(N)[i][1])
    N.isnormal[i] ? print(io, " (normal)") : nothing
  end
end

function _norm_relation_setup_abelian(K::AnticNumberField; small_degree::Bool = true, pure::Bool = true, index::fmpz = zero(fmpz))
  G = automorphisms(K)
  A, GtoA, AtoG = Hecke.find_isomorphism_with_abelian_group(G);
  if iszero(index)
    subs = [f for f in subgroups(A) if order(f[1]) > 1]
  else
    subs = [f for f in subgroups(A) if order(f[1]) > 1 && order(f[1]) == index || order(f[1]) == index^2]
  end

  b, den, ls = Hecke._has_norm_relation_abstract(A, subs, large_index = !small_degree, pure = pure)
  @assert b
  n = length(ls)

  z = NormRelation{Int}()
  z.K = K
  z.subfields = Vector{Tuple{AnticNumberField, NfToNfMor}}(undef, n)
  z.denominator = den
  z.ispure = false

  for i in 1:n
    F, mF = Hecke.fixed_field1(K, NfToNfMor[AtoG[f] for f in ls[i][2]])
    S, mS = simplify(F, cached = false)
    L = S
    mL = mS * mF
    z.subfields[i] = L, mL
  end

  if pure
    z.ispure = true
    z.pure_coefficients = Vector{Int}(undef, n)
    for i in 1:n
      z.pure_coefficients[i] = ls[i][1]
    end
  else
    z.coefficients = Vector{Vector{Tuple{NfToNfMor, Int}}}(undef, n)
    for i in 1:n
      w = Vector{Tuple{NfToNfMor, Int}}(undef, length(ls[i][1]))
      for (j, (expo, auto)) in enumerate(ls[i][1])
        w[j] = (AtoG[auto], expo)
      end
      z.coefficients[i] = w
    end

  end

  z.isabelian = true

  return z
end

function _norm_relation_for_sunits(K::AnticNumberField; small_degree::Bool = true,  pure::Bool = false, target_den::fmpz = zero(fmpz), max_degree::Int = degree(K))
  @vprint :NormRelation 1 "Computing automorphisms\n"
  A = automorphisms(K)
  G, AtoG, GtoA = generic_group(A, *)
  if iszero(target_den)
     b, den, ls = _has_norm_relation_abstract(G, [f for f in subgroups(G, conjugacy_classes = false) if order(f[1]) > 1 && div(order(G), order(f[1])) <= max_degree], pure = pure, large_index = small_degree)
  else
    b, den, ls = _has_norm_relation_abstract(G, [f for f in subgroups(G, conjugacy_classes = false) if order(f[1]) > 1 && div(order(G), order(f[1])) <= max_degree], target_den = target_den, pure = pure, large_index = small_degree)
  end

  if !b
    @show find_small_group(G)
    throw(error("Galois group does not admit Brauer relation"))
  end

  nonredundant = trues(length(ls))

  for i in 1:length(ls)
    if !nonredundant[i]
      continue
    end
    for j = i+1:length(ls)
      if Base.issubset(ls[j][1], ls[i][1])
        nonredundant[i] = false
        break
      elseif Base.issubset(ls[i][1], ls[j][1])
        nonredundant[j] = false
      end
    end
  end
  nsubs = findall(nonredundant)
  n = length(nsubs)

  z = NormRelation{Int}()
  z.K = K
  z.isnormal = falses(n)
  z.subfields = Vector{Tuple{AnticNumberField, NfToNfMor}}(undef, n)
  z.denominator = den
  z.ispure = pure
  z.embed_cache_triv = Vector{Dict{nf_elem, nf_elem}}(undef, n)
  z.nonredundant = Int[i for i = 1:n]

  for i in 1:n
    z.embed_cache_triv[i] = Dict{nf_elem, nf_elem}()
  end

  @vprint :NormRelation 1 "Computing subfields\n"
  for i in 1:n
    @vprint :NormRelation 3 "Computing subfield $i / $n \n"
		auts = NfToNfMor[GtoA[f] for f in ls[nsubs[i]][1]]
    @vtime :NormRelation 3 F, mF = Hecke.fixed_field1(K, auts)
    @vprint :NormRelation 1 "Simplifying \n $F\n"
    @vtime :NormRelation 3 S, mS = simplify(F, cached = false)
    L = S
    mL = mS * mF
    z.subfields[i] = L, mL
    z.isnormal[i] = Hecke._isnormal(ls[i][1])
  end
  z.isabelian = false

  return z
end

function _norm_relation_setup_generic(K::AnticNumberField; small_degree::Bool = true, pure::Bool = false, target_den::fmpz = zero(fmpz), max_degree::Int = degree(K))
  @vprint :NormRelation 1 "Computing automorphisms\n"
  A = automorphisms(K)
  G, AtoG, GtoA = generic_group(A, *)
  if iszero(target_den)
     b, den, ls = _has_norm_relation_abstract(G, [f for f in subgroups(G, conjugacy_classes = false) if order(f[1]) > 1 && div(order(G), order(f[1])) <= max_degree], pure = pure, large_index = small_degree)
  else
    b, den, ls = _has_norm_relation_abstract(G, [f for f in subgroups(G, conjugacy_classes = false) if order(f[1]) > 1 && div(order(G), order(f[1])) <= max_degree], target_den = target_den, pure = pure, large_index = small_degree)
  end

  if !b
    @show find_small_group(G)
    throw(error("Galois group does not admit Brauer relation"))
  end
  n = length(ls)

  nonredundant = Dict{Int, Bool}()

  for i in 1:length(ls)
    red = false
    for j in keys(nonredundant)
      if Base.issubset(ls[j][1], ls[i][1])
        red = true
        break
      elseif Base.issubset(ls[i][1], ls[j][1])
        delete!(nonredundant, j)
        nonredundant[i] = true
      end
    end

    if !red
      nonredundant[i] = true
    end
  end

  z = NormRelation{Int}()
  z.K = K
  z.isnormal = falses(n)
  z.subfields = Vector{Tuple{AnticNumberField, NfToNfMor}}(undef, n)
  z.denominator = den
  z.ispure = pure
  z.embed_cache_triv = Vector{Dict{nf_elem, nf_elem}}(undef, n)
  z.nonredundant = collect(keys(nonredundant))

  for i in 1:n
    z.embed_cache_triv[i] = Dict{nf_elem, nf_elem}()
  end

  @vprint :NormRelation 1 "Computing subfields\n"
  for i in 1:n
    @vprint :NormRelation 3 "Computing subfield $i / $n \n"
		auts = NfToNfMor[GtoA[f] for f in ls[i][1]]
    @vtime :NormRelation 3 F, mF = Hecke.fixed_field1(K, auts)
    @vprint :NormRelation 1 "Simplifying \n $F\n"
    @vtime :NormRelation 3 S, mS = simplify(F, cached = false)
    L = S
    mL = mS * mF
    z.subfields[i] = L, mL
    z.isnormal[i] = Hecke._isnormal(ls[i][1])
  end

  z.coefficients_gen = Vector{Vector{Tuple{Int, NfToNfMor, NfToNfMor}}}(undef, n)

  for i in 1:n
    w = Vector{Tuple{Int, NfToNfMor, NfToNfMor}}(undef, length(ls[i][2]))
    for (j, (expo, auto_pre, auto_post)) in enumerate(ls[i][2])
      w[j] = (expo, GtoA[auto_pre], GtoA[auto_post])
    end
    z.coefficients_gen[i] = w
  end

  z.isabelian = false

  return z
end

Base.length(N::NormRelation) = length(N.subfields)

field(N::NormRelation) = N.K

subfields(N::NormRelation) = N.subfields

subfield(N::NormRelation, i::Int) = N.subfields[i]

embedding(N::NormRelation, i::Int) = N.subfields[i][2]

index(N::NormRelation) = N.denominator

norm(N::NormRelation, i::Int, a) = norm(N.embeddings[i], a)

ispure(N::NormRelation) = N.ispure

function check_relation(N::NormRelation, a::nf_elem)
  @assert !iszero(a)
  if false#ispure(N)
    z = one(field(N))
    for i in 1:length(N)
      z = z * N(norm(embedding(N, i), a), i)
    end
    return a^index(N) == z
  else
    # I applied the exponent when embedding
    z = one(field(N))
    for i in 1:length(N)
      z = z * N([Hecke.norm(embedding(N, i), auto(a)^exp) for  (exp, auto, _) in N.coefficients_gen[i]], i)
    end
    return a^index(N) == z
  end
end

function _apply_auto(N::NormRelation, auto::NfToNfMor, x::nf_elem)
  if haskey(N.mor_cache, auto)
    d = N.mor_cache[auto]
    if haskey(d, x)
      return d[x]
    else
      y = auto(x)
      d[x] = y
      return y
    end
  else
    d = Dict{nf_elem, nf_elem}()
    y = auto(x)
    d[x] = y
    N.mor_cache[auto] = d
    return y
  end
end

function (N::NormRelation)(x::Union{nf_elem, FacElem{nf_elem, AnticNumberField}}, i::Int, j::Int)
  @assert (1 <= j) && (j <= length(N.coefficients_gen[i]))
  if haskey(N.embed_cache, (i, j)) && haskey(N.embed_cache[i, j], x)
    return N.embed_cache[i, j][x]
  end

  _, mk = subfield(N, i)
  y = _apply_auto(N, mk, x)

  exp, _, auto = N.coefficients_gen[i][j]

  z = _apply_auto(N, auto, y)
  return z
end

function (N::NormRelation)(x::Vector{<:Union{nf_elem, FacElem{nf_elem, AnticNumberField}}}, i::Int)
  @assert length(x) == length(N.coefficients_gen[i])
  z = one(N.K)
  for j in 1:length(N.coefficients_gen[i])
    z = z * N(x[j], i, j)
  end
  return z
end

################################################################################
#
#  Cheap induce action
#
################################################################################

function induce_action_from_subfield(N::NormRelation, i, s, FB, cache)
  S = FB.ideals
  ZK = order(S[1])

  z = SMat{fmpz}[zero_matrix(SMat, FlintZZ, 0, length(S)) for i in 1:degree(field(N))]

  mk = embedding(N, i)
  zk = order(s[1])

  if length(cache) == 0
    cache = resize!(cache, length(s))
    cached = false
  else
    cached = true
  end

  autos = automorphisms(field(N), copy = false)

  for auto in autos
    if haskey(N.induced, auto)
      p = N.induced[auto]
    else
      p = induce(FB, auto)
      N.induced[auto] = p
    end
  end

  @assert mod(degree(ZK), degree(zk)) == 0
  reldeg = divexact(degree(ZK), degree(zk))

  for l in 1:length(s)
    if cached
      v = cache[l]
    else
      v = Tuple{Int, fmpz}[]
      P = s[l]
      genofsl = elem_in_nf(_apply_auto(N, mk, P.gen_two.elem_in_nf))
      pmin = minimum(P, copy = false)
      # Compute the number of ideals
      # Assume that L/K and L/Q are normal
      rele = divexact(ramification_index((FB.fb[pmin].lp)[1][2]), ramification_index(P))
      relf = divexact(degree((FB.fb[pmin].lp)[1][2]), degree(P))
      @assert mod(reldeg, rele * relf) == 0
      numb_ideal = divexact(reldeg, rele * relf)
      found = 0
      for k in 1:length(S)
        Q = S[k]
        if minimum(Q, copy = false) == pmin
          ant = anti_uniformizer(Q)
          if (genofsl * ant) in ZK
            found += 1
            @assert mod(ramification_index(Q), ramification_index(s[l])) == 0
            push!(v, (k, divexact(ramification_index(Q), ramification_index(s[l]))))
          end
        end
        if found == numb_ideal
          break
        end
      end
      cache[l] = v
    end

    # We still need to sort the positions of the non-zero entries
    #@show l, v
    sort!(v, by = x -> x[1])
    #@show l, v
    # Now apply the automorphism
    for i in 1:degree(field(N))
      p = N.induced[autos[i]]
      r = [(p[i], e) for (i, e) in v]
      sort!(r, lt = (a,b)->a[1]<b[1])
    #@show l, r
      push!(z[i], sparse_row(FlintZZ, r, sort = false))
    end
  end

  #for l in 1:length(s)
  #  @show l
  #  a = uniformizer(s[l])
  #  @show [ valuation(a, Q) for Q in s]
  #  @show [ valuation(auto(mk(a.elem_in_nf)), Q) for Q in FB.ideals]
  #  @show [ z.rows[l][m] for m in 1:length(FB.ideals) ]
  #  @assert [ valuation(auto(mk(a.elem_in_nf)), Q) for Q in FB.ideals] == [ z.rows[l][m] for m in 1:length(FB.ideals) ]
  #end
  return z
end

function induce_action(N::NormRelation, i, j, s, FB, cache)
  S = FB.ideals
  ZK = order(S[1])
  z = zero_matrix(SMat, FlintZZ, 0, length(S))
  mk = embedding(N, i)
  zk = order(s[1])

  _ , _, auto = N.coefficients_gen[i][j]
  #@show auto
  if haskey(N.induced, auto)
    p = N.induced[auto]
  else
    p = induce(FB, auto)
    N.induced[auto] = p
  end
  #@show p

  if length(cache) == 0
    cache = resize!(cache, length(s))
    cached = false
  else
    cached = true
  end

  @assert mod(degree(ZK), degree(zk)) == 0
  reldeg = divexact(degree(ZK), degree(zk))

  for l in 1:length(s)
    if cached
      v = cache[l]
    else
      v = Tuple{Int, fmpz}[]
      P = s[l]
      genofsl = elem_in_nf(_apply_auto(N, mk, P.gen_two.elem_in_nf))
      pmin = minimum(P, copy = false)
      # Compute the number of ideals
      # Assume that L/K and L/Q are normal
      rele = divexact(ramification_index((FB.fb[pmin].lp)[1][2]), ramification_index(P))
      relf = divexact(degree((FB.fb[pmin].lp)[1][2]), degree(P))
      @assert mod(reldeg, rele * relf) == 0
      numb_ideal = divexact(reldeg, rele * relf)
      found = 0
      for k in 1:length(S)
        Q = S[k]
        if minimum(Q, copy = false) == pmin
          ant = anti_uniformizer(Q)
          if (genofsl * ant) in ZK
            found += 1
            @assert mod(ramification_index(Q), ramification_index(s[l])) == 0
            push!(v, (k, divexact(ramification_index(Q), ramification_index(s[l]))))
          end
        end
        if found == numb_ideal
          break
        end
      end
      cache[l] = v
    end

    # We still need to sort the positions of the non-zero entries
    @show l, v
    sort!(v, by = x -> x[1])
    #@show l, v
    # Now apply the automorphism
    r = [(p[i], e) for (i, e) in v]
    sort!(r, lt = (a,b)->a[1]<b[1])
    #@show l, r
    push!(z, sparse_row(FlintZZ, r, sort = false))
  end

  #for l in 1:length(s)
  #  @show l
  #  a = uniformizer(s[l])
  #  @show [ valuation(a, Q) for Q in s]
  #  @show [ valuation(auto(mk(a.elem_in_nf)), Q) for Q in FB.ideals]
  #  @show [ z.rows[l][m] for m in 1:length(FB.ideals) ]
  #  @assert [ valuation(auto(mk(a.elem_in_nf)), Q) for Q in FB.ideals] == [ z.rows[l][m] for m in 1:length(FB.ideals) ]
  #end
  return z
end

function induce_action(N::NormRelation, i, s::Vector, S::Vector)
  if !ispure(N)
    throw(error("Not implemented yet"))
  end
  ZK = order(S[1])
  z = zero_matrix(SMat, FlintZZ, 0, length(S))
  mk = embedding(N, i)
  for j in 1:length(s)
    v = Tuple{Int, fmpz}[]
    for k in 1:length(S)
      Q = S[k]
      if minimum(Q, copy = false) == minimum(s[j], copy = false)
        ant = anti_uniformizer(Q)
        if (elem_in_nf(mk(s[j].gen_two.elem_in_nf)) * ant) in ZK
          push!(v, (k, divexact(ramification_index(Q), ramification_index(s[j]))))
        end
      end
    end

    # We still need to sort the positions of the non-zero entries
    sort!(v, by = x -> x[1])

    push!(z, N.pure_coefficients[i] * sparse_row(FlintZZ, v, sort = false))
  end
  return z
end

one(T::FacElemMon{AnticNumberField}) = T()

function Hecke.simplify(c::Hecke.ClassGrpCtx)
  d = Hecke.class_group_init(c.FB, SMat{fmpz}, add_rels = false)
  U = Hecke.UnitGrpCtx{FacElem{nf_elem, AnticNumberField}}(order(d))

  Hecke.module_trafo_assure(c.M)
  trafos = c.M.trafo

  for i=1:length(c.FB.ideals)
    x = zeros(fmpz, length(c.R_gen) + length(c.R_rel))
    x[i] = 1
    for j in length(trafos):-1:1
      Hecke.apply_right!(x, trafos[j])
    end
    y = FacElem(vcat(c.R_gen, c.R_rel), x)
    fl = Hecke.class_group_add_relation(d, y, deepcopy(c.M.basis.rows[i]))
    @assert fl
  end
  for i=1:nrows(c.M.rel_gens)
    if iszero(c.M.rel_gens.rows[i])
      Hecke._add_unit(U, c.R_rel[i])
    end
  end
  for i=1:length(U.units)
    Hecke.class_group_add_relation(d, U.units[i], SRow(FlintZZ))
  end
  return d, U
end

function units(c::Hecke.ClassGrpCtx)
  d = Hecke.class_group_init(c.FB, SMat{fmpz}, add_rels = false)
  U = Hecke.UnitGrpCtx{FacElem{nf_elem, AnticNumberField}}(order(d))

  Hecke.module_trafo_assure(c.M)
  trafos = c.M.trafo

  for i=1:nrows(c.M.rel_gens)
    if iszero(c.M.rel_gens.rows[i])
      Hecke._add_unit(U, c.R_rel[i])
    end
  end

  U.units = Hecke.reduce(U.units, U.tors_prec)
  U.tentative_regulator = Hecke.regulator(U.units, 64)

  return U
end

################################################################################
#
#  Residue computation via Brauer relations
#
################################################################################

function zeta_log_residue(O::NfOrd, N::NormRelation, abs_error::Float64)
  degree(O) == 1 && error("Number field must be of degree > 1")
  !ispure(N) && error("Norm relation must be a Brauer relation")
  @show index(N)
  target_prec = Int(floor(log(abs_error)))
  @show target_prec
  residues = arb[]
  for i in 1:length(N)
    v = N.coefficients_gen[i]
    @assert length(v) == 1
    c = (v[1][1] * divexact(degree(N.K), degree(N.subfields[i][1]))) // index(N)
    push!(residues, c * zeta_log_residue(maximal_order(N.subfields[i][1]), abs_error/(2*length(N))))
  end
  z = sum(residues)
  @assert radiuslttwopower(z, target_prec)
  return z
end

function _lift_to_normalized_brauer_relation(N)
  @assert ispure(N)
  rel = fmpq[]
  for i in 1:length(N)
    v = N.coefficients_gen[i]
    push!(rel, fmpq(v[1][1] * divexact(degree(N.K), degree(N.subfields[i][1]))) // index(N))
  end
  return rel
end

function hR_from_subfields(N, abs_prec::Int = 64)
  @assert ispure(N)
  rel = _lift_to_normalized_brauer_relation(N)
  vals = arb[]
  for i in 1:length(N)
    zk = lll(maximal_order(subfield(N, i)[1]))
    c, _ = class_group(zk)
    u, mu = unit_group_fac_elem(zk)
    push!(vals, (1//order(torsion_unit_group(zk)[1]) * order(c) * regulator([mu(u[i]) for i in 2:rank(u) + 1], abs_prec))^rel[i])
  end
  return prod(vals) * order(torsion_unit_group(maximal_order(field(N)))[1])
end

################################################################################
#
#  Code for abstract norm relations
#
################################################################################

# TODO: Make this great again
function _has_norm_relation_abstract(G::GrpGen, H::Vector{Tuple{GrpGen, GrpGenToGrpGenMor}};
                                     primitive::Bool = false,
                                     greedy::Bool = false,
                                     large_index::Bool = false,
                                     pure::Bool = false,
                                     index_bound::Int = -1,
                                     target_den::fmpz = zero(fmpz))
  if index_bound != -1
    H = [h for h in H if order(G) <= order(h[1]) * index_bound]
  end

  sort!(H, by = x -> order(x[1]))

  if large_index
    reverse!(H)
  end

  QG = AlgGrp(FlintQQ, G)
  norms_rev = Dict{elem_type(QG), Int}()
  norms = Vector{elem_type(QG)}(undef, length(H))
  for i in 1:length(H)
    norms_rev[sum([QG(H[i][2](h)) for h in H[i][1]])] = i
    norms[i] = sum([QG(H[i][2](h)) for h in H[i][1]])
  end

  n = Int(order(G))

  if pure
    if iszero(target_den)
      m = zero_matrix(FlintQQ, length(H), n)
      for i in 1:length(H)
        for j in 1:n
          m[i, j] = norms[i].coeffs[j]
        end
      end

      onee = matrix(FlintQQ, 1, n, coefficients(one(QG)))

      b, v, K = can_solve_with_kernel(m, onee, side = :left)
    else
      m = zero_matrix(FlintZZ, length(H), n)
      for i in 1:length(H)
        for j in 1:n
          m[i, j] = FlintZZ(norms[i].coeffs[j])
        end
      end

      onee = matrix(FlintZZ, 1, n, coefficients(one(QG)))

      b, w, K = can_solve_with_kernel(m, target_den * onee, side = :left)
      v = 1//target_den * change_base_ring(FlintQQ, w)
    end

    if !b
      return false, zero(FlintZZ), Vector{Tuple{Vector{Tuple{fmpz, GrpAbFinGenElem}}, Vector{GrpAbFinGenElem}}}()
    end

    @assert b

    if iszero(target_den)
      v = _reduce_modulo(v, K)
    end

    subgroups_needed = Int[ i for i in 1:length(H) if !iszero(v[1, i])]

    den = denominator(v[1, 1])
    for i in 2:length(H)
      den = lcm!(den, den, denominator(v[1, i]))
    end

    vvv = Vector{Vector{Tuple{fmpz, GrpGenElem, GrpGenElem}}}(undef, length(subgroups_needed))
    for i in 1:length(vvv)
      vvv[i] = Vector{Tuple{fmpz, GrpGenElem, GrpGenElem}}()
    end

    solutions = Vector{Tuple{Vector{GrpGenElem}, Vector{Tuple{fmpz, GrpGenElem, GrpGenElem}}}}()

    for i in 1:length(subgroups_needed)
      push!(vvv[i], (numerator(den * v[1, subgroups_needed[i]]), id(G), id(G)))
      sgroup = [H[subgroups_needed[i]][2](h) for h in H[subgroups_needed[i]][1]]
      push!(solutions, (sgroup, vvv[i]))
    end

    z = zero(QG)
    for cc in 1:length(subgroups_needed)
      z = z + v[1, subgroups_needed[cc]] * norms[subgroups_needed[cc]]
    end

    @assert isone(z)

    return true, den, solutions
  end

  wd = decompose(QG)
  idempotents = [ f(one(A)) for (A, f) in wd]

  nonannihilating = Vector{Vector{Int}}(undef, length(idempotents))

  for j in 1:length(idempotents)
    nonannihilating_ej = Vector{Int}()
    for i in 1:length(norms)
      if !iszero(idempotents[j] * norms[i])
        push!(nonannihilating_ej, i)
      end
    end
    nonannihilating[j] = nonannihilating_ej
  end

  l = length(idempotents)

  subgroups_needed = Int[]

  if any(isempty, nonannihilating)
    return false, zero(FlintZZ), Vector{Tuple{Vector{Tuple{fmpz, GrpAbFinGenElem}}, Vector{GrpAbFinGenElem}}}()
  end

  subgroups_needed = Int[]

  for i in 1:length(nonannihilating)
    if greedy || !large_index
      sort!(nonannihilating[i], by = x -> order(H[x][1]), rev = true)
    elseif large_index
      sort!(nonannihilating[i], by = x -> order(H[x][1]))
    end
  end

  k = nonannihilating[1][1]
  push!(subgroups_needed, k)
  i = 2
  for i in 2:l
    if length(intersect(subgroups_needed, nonannihilating[i])) > 0
      continue
    else
      push!(subgroups_needed, nonannihilating[i][1])
    end
  end

  if !iszero(target_den)
    subgroups_needed = collect(1:length(H))
  end

  if primitive
    # Check if we can replace a subgroup by bigger subgroups
    for i in 1:length(subgroups_needed)
      Q, mQ = quo(G, H[subgroups_needed[i]][1])
      if length(Q) == 1
        continue
      end
      b, den, ls = _has_norm_relation_abstract(Q, [f for f in subgroups(Q) if order(f[1]) > 1])
      if b
        return _has_norm_relation_abstract(G, deleteat!(H, subgroups_needed[i]))
      end
    end
  end

  # Now solve the equation to find a representation 1 = sum_{H} x_H N_H y_H
  n = Int(order(G))

  # We use the following trick to reduce the number of equations
  # Instead of working with the Q-basis (g N_H h)_(H, g in G, h in H), we take
  # (g N_H h)_(H, g in G/H, h in G\N_G(H))

  left_cosets_for_sub = Vector{GrpGenElem}[]
  right_cosets_for_normalizer = Vector{GrpGenElem}[]

  for i in 1:length(subgroups_needed)
    lc = Hecke.left_cosets(G, H[subgroups_needed[i]][2])
    push!(left_cosets_for_sub, lc)
    NGH = Hecke.normalizer(G, H[subgroups_needed[i]][2])
    push!(right_cosets_for_normalizer, Hecke.right_cosets(G, NGH[2]))
  end

  if iszero(target_den)
    onee = matrix(FlintQQ, 1, n, coefficients(one(QG)))

    m = zero_matrix(FlintQQ, dot(length.(left_cosets_for_sub), length.(right_cosets_for_normalizer)), n)

    B = basis(QG)

    help_with_indicies = Vector{Tuple{Int, Int, Int}}(undef, nrows(m))
    cc = 1
    for i in 1:length(subgroups_needed)
      lcosub = length(left_cosets_for_sub[i])
      lconor = length(right_cosets_for_normalizer[i])
      for (k, g) in enumerate(left_cosets_for_sub[i])
        for (l, h) in enumerate(right_cosets_for_normalizer[i])
          for j in 1:n
            m[cc, j] = (QG(g) * norms[subgroups_needed[i]] * QG(h)).coeffs[j]
          end
          help_with_indicies[cc] = (i, k, l)
          cc += 1
        end
      end
    end

    b, w, K = can_solve_with_kernel(m, onee, side = :left)
    v = w
  elseif true
    onee = matrix(FlintZZ, 1, n, coefficients(target_den * one(QG)))

    m = zero_matrix(FlintZZ, dot(length.(left_cosets_for_sub), length.(right_cosets_for_normalizer)), n)

    help_with_indicies = Vector{Tuple{Int, Int, Int}}(undef, nrows(m))
    cc = 1
    for i in 1:length(subgroups_needed)
      for (k, g) in enumerate(left_cosets_for_sub[i])
        for (l, h) in enumerate(right_cosets_for_normalizer[i])
          for j in 1:n
            m[cc, j] = FlintZZ((QG(g) * norms[subgroups_needed[i]] * QG(h)).coeffs[j])
          end
          help_with_indicies[cc] = (i, k, l)
          cc += 1
        end
      end
    end

    b, w, K = can_solve_with_kernel(m, onee, side = :left)

    v = 1//target_den * change_base_ring(FlintQQ, w)
  end

  @assert b

  z = zero(QG)
  for cc in 1:ncols(v)
    if iszero(v[1, cc])
      continue
    else
      i, k, l = help_with_indicies[cc]
      z = z + v[1, cc] *  QG(left_cosets_for_sub[i][k]) * norms[subgroups_needed[i]] * QG(right_cosets_for_normalizer[i][l])
    end
  end

  @assert isone(z)

  den = one(FlintZZ)

  for i in 1:ncols(v)
    den = lcm(den, denominator(v[1, i]))
  end

  solutions = Vector{Tuple{Vector{GrpGenElem}, Vector{Tuple{fmpz, GrpGenElem, GrpGenElem}}}}()

  vvv = Vector{Vector{Tuple{fmpz, GrpGenElem, GrpGenElem}}}(undef, length(subgroups_needed))
  for i in 1:length(vvv)
    vvv[i] = Vector{Tuple{fmpz, GrpGenElem, GrpGenElem}}()
  end

  for cc in 1:ncols(v)
    if iszero(v[1, cc])
      continue
    else
      i, k, l = help_with_indicies[cc]
      push!(vvv[i], (FlintZZ(den * v[1, cc]), left_cosets_for_sub[i][k], right_cosets_for_normalizer[i][l]))
    end
  end

  for i in 1:length(vvv)
    sgroup = [H[subgroups_needed[i]][2](h) for h in H[subgroups_needed[i]][1]]
    if length(vvv[i]) > 0
      push!(solutions, (sgroup, vvv[i]))
    end
  end

  if !iszero(target_den)
    @assert iszero(mod(target_den, den))
  end

  return true, den, solutions
end

################################################################################
#
#  Existence of norm relations
#
################################################################################

function has_useful_brauer_relation(G)
  subs = subgroups(G, conjugacy_classes = true)
  for (H,_) in subs
    if iscyclic(H)
      continue
    end
    fac = factor(order(H))
    if length(fac) > 2
      continue
    end
    if length(fac) == 2
      e = prod(e for (p, e) in fac)
      if e == 1
        return true
      end
    end
    if length(fac) == 1
      for (p, e) in fac
        if e == 2
          return true
        end
      end
    end
  end
  return false
end

const _fermat_primes = Int[3, 5, 17, 257, 65537]

function has_useful_generalized_norm_relation(G)
  subs = subgroups(G, conjugacy_classes = true)
  for (H,_) in subs
    if iscyclic(H)
      continue
    end

    # This is the check for SL(2, 3)
    # SL(2, 5) has order 120
    if order(H) == 24 && find_small_group(H)[1] == (24, 3)
      return true
    end

    fac = factor(order(H))
    if length(fac) > 2
      continue
    end
    if length(fac) == 2
      e = prod(e for (p, e) in fac)
      if e == 1
        return true
      end
    end
    if length(fac) == 1
      for (p, e) in fac
        if e == 2
          return true
        end
      end
    end
  end
  return false
end

################################################################################
#
#  Get smallest norm relation coprime to something
#
################################################################################

function _smallest_scalar_norm_relation_coprime(G::GrpGen, m::fmpz)

  primes = fmpz[ p for (p, _) in factor(m)]

  S = Localization(FlintZZ, primes)

  all_non_trivial_subs = [ (H, mH) for (H, mH) in subgroups(G) if order(H) > 1]

  sort!(all_non_trivial_subs, by = x -> order(x[1]))

  reverse!(all_non_trivial_subs)

  QG = AlgGrp(FlintQQ, G)
  norms_rev = Dict{elem_type(QG), Int}()
  norms = Vector{elem_type(QG)}(undef, length(all_non_trivial_subs))
  for i in 1:length(all_non_trivial_subs)
    norms_rev[sum([QG(all_non_trivial_subs[i][2](h)) for h in all_non_trivial_subs[i][1]])] = i
    norms[i] = sum([QG(all_non_trivial_subs[i][2](h)) for h in all_non_trivial_subs[i][1]])
  end

  n = Int(order(G))

  M = zero_matrix(S, length(all_non_trivial_subs), n)

  for i in 1:length(all_non_trivial_subs)
    for j in 1:n
      M[i, j] = norms[i].coeffs[j]
    end
  end

  onee = matrix(S, 1, n, coefficients(one(QG)))

  fl, v = can_solve_with_solution(M, onee, side = :left)

  if !fl
    solutions = Vector{Tuple{Vector{GrpGenElem}, Vector{Tuple{fmpz, GrpGenElem, GrpGenElem}}}}()
    return false, fmpz(0), solutions
  end

  k = 0

  for i in 1:length(all_non_trivial_subs)
    fl, v = can_solve_with_solution(view(M, 1:i, 1:ncols(M)), onee, side = :left)
    if fl
      k = i
      break
    end
  end

  _,_, K = can_solve_with_kernel(view(M, 1:k, 1:ncols(M)), onee, side = :left)

  v = _reduce_modulo(v, K)

  subgroups_needed = Int[ i for i in 1:k if !iszero(v[1, i])]

  den = denominator(v[1, 1])
  for i in 2:k
    den = lcm!(den, den, denominator(v[1, i]))
  end

  @assert iscoprime(den, m)

  vvv = Vector{Vector{Tuple{fmpz, GrpGenElem, GrpGenElem}}}(undef, length(subgroups_needed))
  for i in 1:length(vvv)
    vvv[i] = Vector{Tuple{fmpz, GrpGenElem, GrpGenElem}}()
  end

  solutions = Vector{Tuple{Vector{GrpGenElem}, Vector{Tuple{fmpz, GrpGenElem, GrpGenElem}}}}()

  for i in 1:length(subgroups_needed)
    push!(vvv[i], (numerator(den * v[1, subgroups_needed[i]]), id(G), id(G)))
    sgroup = [all_non_trivial_subs[subgroups_needed[i]][2](h) for h in all_non_trivial_subs[subgroups_needed[i]][1]]
    push!(solutions, (sgroup, vvv[i]))
  end

  z = zero(QG)
  for cc in 1:length(subgroups_needed)
    z = z + v[1, subgroups_needed[cc]].data * norms[subgroups_needed[cc]]
  end

  @assert isone(z)

  return true, den, solutions
end

function has_coprime_norm_relation(K::AnticNumberField, m::fmpz)
  G, mG = automorphism_group(K)
  b, den, ls = _smallest_scalar_norm_relation_coprime(G, m)

  if !b
    return false, NormRelation{Int}()
  end

  nonredundant = Dict{Int, Bool}()

  for i in 1:length(ls)
    red = false
    for j in keys(nonredundant)
      if Base.issubset(ls[j][1], ls[i][1])
        red = true
        break
      elseif Base.issubset(ls[i][1], ls[j][1])
        delete!(nonredundant, j)
        nonredundant[i] = true
      end
    end

    if !red
      nonredundant[i] = true
    end
  end

  n = length(ls)

  z = NormRelation{Int}()
  z.K = K
  z.isnormal = falses(n)
  z.subfields = Vector{Tuple{AnticNumberField, NfToNfMor}}(undef, n)
  z.denominator = den
  z.ispure = true
  z.embed_cache_triv = Vector{Dict{nf_elem, nf_elem}}(undef, n)
  z.nonredundant = collect(keys(nonredundant))

  for i in 1:n
    z.embed_cache_triv[i] = Dict{nf_elem, nf_elem}()
  end

  for i in 1:n
    F, mF = Hecke.fixed_field1(K, NfToNfMor[mG(f) for f in ls[i][1]])
    S, mS = simplify(F, cached = false)
    L = S
    mL = mS * mF
    z.subfields[i] = L, mL
    z.isnormal[i] = Hecke._isnormal(ls[i][1])
  end

  z.coefficients_gen = Vector{Vector{Tuple{Int, NfToNfMor, NfToNfMor}}}(undef, n)

  for i in 1:n
    w = Vector{Tuple{Int, NfToNfMor, NfToNfMor}}(undef, length(ls[i][2]))
    for (j, (expo, auto_pre, auto_post)) in enumerate(ls[i][2])
      w[j] = (expo, mG(auto_pre), mG(auto_post))
    end
    z.coefficients_gen[i] = w
  end

  z.isabelian = false
  z.isoptimal = m

  return true, z
end

function _reduce_modulo(v, K)
  H = hnf(K)
  w = deepcopy(v)
  for i in nrows(H):-1:1
    k = ncols(H)
    while iszero(H[i, k])
      k = k - 1
    end
    if !iszero(w[1, k])
      fl, c = divides(w[1, k], H[i, k])
      if fl
        for j in 1:ncols(H)
          w[1, j] = w[1, j] - c * H[i, j]
        end
      end
    end
  end
  return w
end


text:
# This file is a part of TypeDBClient.  License is MIT: https://github.com/Humans-of-Julia/TypeDBClient.jl/blob/main/LICENSE

#
# package typedb.client.cluster;
#
# import typedb.client.api.TypeDBClient;
# import typedb.client.api.TypeDBOptions;
# import typedb.client.api.TypeDBSession;
# import typedb.client.common.exception.TypeDBClientException;
# import typedb.client.common.rpc.TypeDBStub;
# import typedb.client.core.CoreClient;
# import typedb.common.collection.Pair;
# import typedb.protocol.ClusterServerProto;
# import io.grpc.StatusRuntimeException;
# import org.slf4j.Logger;
# import org.slf4j.LoggerFactory;
#
# import java.util.Map;
# import java.util.Set;
# import java.util.concurrent.ConcurrentHashMap;
# import java.util.concurrent.ConcurrentMap;
#
# import static typedb.client.common.exception.ErrorMessage.Client.CLUSTER_UNABLE_TO_CONNECT;
# import static typedb.client.common.rpc.RequestBuilder.Cluster.Server.allReq;
# import static typedb.common.collection.Collections.pair;
# import static java.util.stream.Collectors.toMap;
# import static java.util.stream.Collectors.toSet;
#
# public class ClusterClient implements TypeDBClient.Cluster {
#
#     private static final Logger LOG = LoggerFactory.getLogger(ClusterClient.class);
#     private final Map<String, CoreClient> coreClients;
#     private final Map<String, TypeDBStub.Cluster> stubs;
#     private final ClusterDatabaseManager databaseMgrs;
#     private final ConcurrentMap<String, ClusterDatabase> clusterDatabases;
#     private boolean isOpen;
#
#     public ClusterClient(Set<String> addresses) {
#         this(addresses, CoreClient.calculateParallelisation());
#     }
#
#     public ClusterClient(Set<String> addresses, int parallelisation) {
#         coreClients = fetchServerAddresses(addresses).stream()
#                 .map(address -> pair(address, new CoreClient(address, parallelisation)))
#                 .collect(toMap(Pair::first, Pair::second));
#         stubs = coreClients.entrySet().stream()
#                 .map(client -> pair(client.getKey(), TypeDBStub.cluster(client.getValue().channel())))
#                 .collect(toMap(Pair::first, Pair::second));
#         databaseMgrs = new ClusterDatabaseManager(this);
#         clusterDatabases = new ConcurrentHashMap<>();
#         isOpen = true;
#     }
#
#     private Set<String> fetchServerAddresses(Set<String> addresses) {
#         for (String address : addresses) {
#             try (CoreClient client = new CoreClient(address)) {
#                 LOG.debug("Fetching list of cluster servers from {}...", address);
#                 TypeDBStub.Cluster stub = TypeDBStub.cluster(client.channel());
#                 ClusterServerProto.ServerManager.All.Res res = stub.serversAll(allReq());
#                 Set<String> members = res.getServersList().stream().map(ClusterServerProto.Server::getAddress).collect(toSet());
#                 LOG.debug("The cluster servers are {}", members);
#                 return members;
#             } catch (StatusRuntimeException e) {
#                 LOG.error("Fetching cluster servers from {} failed.", address);
#             }
#         }
#         throw new TypeDBClientException(CLUSTER_UNABLE_TO_CONNECT, String.join(",", addresses));
#     }
#
#     @Override
#     public boolean isOpen() {
#         return isOpen;
#     }
#
#     @Override
#     public ClusterDatabaseManager databases() {
#         return databaseMgrs;
#     }
#
#     @Override
#     public ClusterSession session(String database, TypeDBSession.Type type) {
#         return session(database, type, TypeDBOptions.cluster());
#     }
#
#     @Override
#     public ClusterSession session(String database, TypeDBSession.Type type, TypeDBOptions options) {
#         TypeDBOptions.Cluster clusterOptions = options.asCluster();
#         if (clusterOptions.readAnyReplica().isPresent() && clusterOptions.readAnyReplica().get()) {
#             return sessionAnyReplica(database, type, clusterOptions);
#         } else {
#             return sessionPrimaryReplica(database, type, clusterOptions);
#         }
#     }
#
#     private ClusterSession sessionPrimaryReplica(String database, TypeDBSession.Type type, TypeDBOptions.Cluster options) {
#         return openSessionFailsafeTask(database, type, options, this).runPrimaryReplica();
#     }
#
#     private ClusterSession sessionAnyReplica(String database, TypeDBSession.Type type, TypeDBOptions.Cluster options) {
#         return openSessionFailsafeTask(database, type, options, this).runAnyReplica();
#     }
#
#     private FailsafeTask<ClusterSession> openSessionFailsafeTask(
#             String database, TypeDBSession.Type type, TypeDBOptions.Cluster options, ClusterClient client) {
#         return new FailsafeTask<ClusterSession>(this, database) {
#             @Override
#             ClusterSession run(ClusterDatabase.Replica replica) {
#                 return new ClusterSession(client, replica.address(), database, type, options);
#             }
#         };
#     }
#
#     // TODO: this is not good - we should not pass an internal object to be modified outside of this class
#     ConcurrentMap<String, ClusterDatabase> databaseByName() {
#         return clusterDatabases;
#     }
#
#     Map<String, CoreClient> coreClients() {
#         return coreClients;
#     }
#
#     Set<String> clusterMembers() {
#         return coreClients.keySet();
#     }
#
#     CoreClient coreClient(String address) {
#         return coreClients.get(address);
#     }
#
#     TypeDBStub.Cluster stub(String address) {
#         return stubs.get(address);
#     }
#
#     @Override
#     public boolean isCluster() {
#         return true;
#     }
#
#     @Override
#     public Cluster asCluster() {
#         return this;
#     }
#
#     @Override
#     public void close() {
#         coreClients.values().forEach(CoreClient::close);
#         isOpen = false;
#     }
# }

ext: jl


text:
module ProgressTimer

export TimingInit, TimingProgress
using Dates

const sec_len=Dates.Millisecond(Dates.Second(1))
const min_len=Dates.Millisecond(Dates.Minute(1))
const hr_len=Dates.Millisecond(Dates.Hour(1))

function TimingInit()
    ###Construcor. Set the clock and such things
    StartWallTime=now()
    NextWallTime=now()+sec_len ###Wait for a second before showing the clock
    return (StartWallTime,NextWallTime)
end

function get_time_and_unit(Time;Decending=false)
    if Time >= hr_len ###More than an hour
        #println("$Time is bigger than an hour")
        TimeValue=round(Time.value/hr_len.value,digits=1)
        if Decending && Time<(2*hr_len)
            return (TimeValue," hr",Time-hr_len)
        else
            return (TimeValue," hr",hr_len)
        end
    elseif Time >= min_len ###More than a minute
        #println("$Time is bigger than a minute")
        TimeValue=round(Time.value/min_len.value,digits=1)
        if Decending && Time<(2*min_len)
            return (TimeValue," min",Time-min_len)
        else
            return (TimeValue," min",min_len)
        end
    else
        #println("$Time is on the order of seconds")
        return (Int(round(Time.value/sec_len.value))," sec",sec_len)
    end
    return TimeValue, TimeStr, TimeUnit
end

function TimingProgress(TimingObject,Nr,NSamples;Message=nothing)
    (StartWallTime,NextTime)=TimingObject
    NowWallTime=now()
    #println("Nr=$Nr NSamples=$NSamples")
    #println("Wait Time: ",WaitTime)
    #println("Time for next write: ",NextTime)
    #println("Clock Time now:      ",NowWallTime)
    if NowWallTime > NextTime
        #println(".................")
        #println("Time to write!")
        TimeWallElapsed=Dates.Millisecond(NowWallTime-StartWallTime)
        #println("Wall Time Elapsed: ",TimeWallElapsed)
        ##Computes samples left
        SamplesLeft = NSamples - Nr
        #println("SamplesLeft=$SamplesLeft")
        ##Compute time left
        TimeWallLeft = Dates.Millisecond(ceil((SamplesLeft/(1.0*Nr))*TimeWallElapsed.value))
        #println("TimeWallLeft=$TimeWallLeft")
        (t_el,dim_el,step_el)=get_time_and_unit(TimeWallElapsed)
        (t_rem,dim_rem,step_rem)=get_time_and_unit(TimeWallLeft,Decending=true)
        (t_tot,dim_tot,step_tot)=get_time_and_unit(TimeWallElapsed+TimeWallLeft)
        WaitTime=maximum([minimum([step_el,step_rem,step_tot]),sec_len])
        #println(Nr," done ", SamplesLeft," to go. Time: ",TimeWallElapsed,
        #        " Time left: ",TimeWallLeft,
        #        " of ",TimeWallElapsed+TimeWallLeft)
        if Message!=nothing
            println(Message)
        end
        println(Nr," done ", SamplesLeft," to go. Time: ",t_el,dim_el,
                ". Time left: ",t_rem,dim_rem,
                " of ",t_tot,dim_tot,".")
        #println("Wait Time: ",WaitTime)
        NextTime=now()+WaitTime
        #println("Time for new next write: ",NextTime)
        return (StartWallTime,NextTime)
    else
        #println("Skip writing")
        return TimingObject
    end
end

end #module

max_stars_repo_name: fremling/ProgressTimer

max_issues_repo_name: fremling/ProgressTimer

max_forks_repo_name: fremling/ProgressTimer


text:
include("src/base.jl")
include("src/envs/toytext_mdp.jl")
include("src/policies/crossentropy_policy.jl")
include("src/solvers/crossentropy_policy_solver.jl")

#gym_env = ToyTextMDP(gym.make("Taxi-v2"))
gym_env = ToyTextMDP(gym.make("FrozenLake8x8-v0"))

solver = CrossentropyPolicySolver(100, 500, 0.75, 0.1; max_frame_iterations = 1000, print_every_n = 100, min_score = 0.0, min_score_step = 0.0001)
policy = solve(solver, gym_env, verbose = true)

println(policy.action_map)

max_stars_repo_name: dmitrijsc/practical-rl

max_issues_repo_name: dmitrijsc/practical-rl

max_forks_repo_name: dmitrijsc/practical-rl


text:
@testset "`get_package_file`" begin
    using Pkg, JET
    using JET: get_package_file

    target = pathof(JET)

    @test get_package_file("JET") == target
    @test_throws ErrorException get_package_file("unknown")

    @test get_package_file(JET) == target
    @test_throws ErrorException get_package_file(Module())

    # suppress logs from Pkg.jl if possible
    function pkg_activate(args...; io = IOBuffer(), kwargs...)
        @static if VERSION ≥ v"1.7-DEV"
            Pkg.activate(args...; io, kwargs...)
        else
            Pkg.activate(args...; kwargs...)
        end
    end
    old = Pkg.project().path
    try
        pkg_activate(pkgdir(JET))
        @test get_package_file(nothing) == target

        pkg_activate(; temp = true)
        @test_throws ErrorException get_package_file(nothing)
    catch err
        rethrow(err)
    finally
        pkg_activate(old)
    end
end

using Base.TOML
using JET: process_config_dict!
macro toml_str(s); TOML.parse(TOML.Parser(s)); end

@testset "`process_config_dict`" begin
    let
        config_dict = toml"""
        # usual
        analyze_from_definitions = true

        # will be `parse`d or `eval`ed
        context = "Base"
        concretization_patterns = ["const x_ = y_"]
        toplevel_logger = "stdout"
        inference_logger = "stdout"
        """

        config = process_config_dict!(config_dict)
        @test (:context => Base) in config
        @test (:concretization_patterns => [:(const x_ = y_)]) in config
        @test (:toplevel_logger => stdout) in config
        @test (:inference_logger => stdout) in config
    end

    # error when invalid expression given
    let
        config_dict = toml"""
        concretization_patterns = ["const x_ = end"]
        """
        @test_throws Meta.ParseError process_config_dict!(config_dict)
    end

    # error when incomplete expression given
    let
        config_dict = toml"""
        concretization_patterns = ["const x_ = "]
        """
        @test_throws ErrorException process_config_dict!(config_dict)
    end

    # should be whitespece/newline insensitive
    let
        config_dict = toml"""
        concretization_patterns = [
            \"\"\"
            const x_ = y_
            \"\"\"
        ]
        """

        config = process_config_dict!(config_dict)
        @test (:concretization_patterns => [:(const x_ = y_)]) in config
    end
end


text:
function vertCongruenceAA(W)
	err, i, todelete, vclasses = 10^-6, 1, [], []
	verts = convert(Lar.Points, W')
	kdtree = NearestNeighbors.KDTree(verts)
	newverts = zeros(Int, size(verts,2))
	for vi in 1:size(verts,2)
		if !(vi in todelete)
			nearvs = NearestNeighbors.inrange(kdtree, verts[:,vi], err)
			push!(vclasses,nearvs)
			newverts[nearvs] .= i
			nearvs = setdiff(nearvs, vi)
			todelete = union(todelete, nearvs)
			i += 1
		end
	end
	V = zeros(3,length(vclasses))
	for (k,class) in enumerate(vclasses)
			V[:,k] = sum(W[class,:],dims=1)/length(class)
	end
	return V, vclasses
end

function cellCongruenceAA(Delta,inclasses)
	cellarray = Lar.cop2lar(Delta)
	new_e = Array{Int64,1}(undef,size(Delta,2))
	for (k,class) in enumerate(inclasses)
		for e in class
			new_e[e] = k
		end
	end
  cells = [map(e->new_e[e], face) for face in cellarray]
  outclasses = DefaultOrderedDict{Array{Int64,1},Array{Int64,1}}([])
  for (k,face) in enumerate(cells)
    if outclasses[face] == []
    	outclasses[face] = [k]
    else
    	append!(outclasses[face],[k])
    end
  end
	FEnew = sort(collect(keys(outclasses)))
  outclasses = sort(collect(values(outclasses)))
  return FEnew,outclasses
end

function chainCongruenceAA(W, T)
	V, vclasses = vertCongruenceAA(W)
	EV, eclasses = cellCongruenceAA(T[1],vclasses)
	FE, fclasses = cellCongruenceAA(T[2],eclasses)
	#@show size(V,2) - size(EV,1) + size(FE,1)
	return V,EV,FE
end


text:
function msd!(positions::AbstractArray, ref_pos::AbstractArray)
    result = 0.0f0

    @inbounds for i = axes(positions, 2)
        δx = positions[1, i] - ref_pos[1, i]
        δy = positions[2, i] - ref_pos[2, i]
        δz = positions[3, i] - ref_pos[3, i]

        result += δx*δx + δy*δy + δz*δz
    end

    return result
end

function savemsd(positions::AbstractArray, path::String)

    @inbounds for i ∈ axes(positions, 2)
        filename = joinpath(path, "msd_$(i).csv")

        open(filename, "a") do io
            println(io, "$(positions[1, i]),$(positions[2, i]),$(positions[3, i])")
        end
    end

    return nothing
end


text:
function backend()
    if STATE.backend == :NotSet
        backend = get(ENV, "JULIA_CONDAPKG_BACKEND", "")
        exe = get(ENV, "JULIA_CONDAPKG_EXE", "")
        if backend == ""
            backend = exe == "" ? "MicroMamba" : "System"
        end
        if backend == "MicroMamba"
            STATE.backend = :MicroMamba
        elseif backend == "Null"
            STATE.backend = :Null
        elseif backend == "System"
            ok = false
            for exe in (exe == "" ? ["micromamba", "mamba", "conda"] : [exe])
                exe2 = Sys.which(exe)
                if exe2 !== nothing
                    STATE.backend = :System
                    STATE.condaexe = exe2
                    ok = true
                    break
                end
            end
            if !ok
                if exe == ""
                    error("could not find a conda, mamba or micromamba executable")
                else
                    error("not an executable: JULIA_CONDAPKG_EXE=$exe")
                end
            end
        else
            error("invalid backend: JULIA_CONDAPKG_BACKEND=$backend")
        end
    end
    STATE.backend
end

function conda_cmd(args=``; io::IO=stderr)
    b = backend()
    if b == :MicroMamba
        MicroMamba.cmd(args, io=io)
    elseif b == :System
        `$(STATE.condaexe) $args`
    elseif b == :Null
        error("Can not run conda command when backend is Null. Manage conda actions outside of julia.")
    else
        @assert false
    end
end


text:
import Downloads
using GitHub
using Pkg.Artifacts
using Base.BinaryPlatforms
using Tar
using SHA

const gh_auth = GitHub.AnonymousAuth()

version(release::Release) = try
    VersionNumber(release.tag_name)
catch
    v"0.0.0"
end

latest_release(repo; auth) = reduce(releases(repo; auth)[1]) do releaseA, releaseB
    version(releaseA) > version(releaseB) ? releaseA : releaseB
end

function make_artifacts(dir)
    release = latest_release("wasmerio/wasmer"; auth=gh_auth)
    release_version = VersionNumber(release.tag_name)

    platforms = Dict([
        "wasmer-linux-amd64.tar.gz" => Platform("x86_64", "linux"; libc="glibc"),
        "wasmer-linux-musl-amd64.tar.gz" => Platform("x86_64", "linux"; libc="musl"),
        "wasmer-darwin-amd64.tar.gz" => Platform("x86_64", "macos"),
        "wasmer-darwin-arm64.tar.gz" => Platform("aarch64", "macos"),
   ])

    assets = filter(asset -> haskey(platforms, asset["name"]), release.assets)
    artifacts_toml = joinpath(@__DIR__, "Artifacts.toml")

    for asset in assets
        platform = platforms[asset["name"]]

        @info "Downloading $(asset["browser_download_url"]) for $platform"
        archive_location = joinpath(dir, asset["name"])
        download_url = asset["browser_download_url"]
        Downloads.download(download_url, archive_location;
            progress=(t,n) -> print("$(floor(100*n/t))%\r"))
        println()

        artifact_hash = create_artifact() do artifact_dir
            run(`tar -xvf $archive_location -C $artifact_dir`)
        end

        download_hash = open(archive_location, "r") do f
            bytes2hex(sha256(f))
        end
        bind_artifact!(artifacts_toml, "libwasmer", artifact_hash; platform, force=true, download_info=[
            (download_url, download_hash)
        ])
        @info "done $platform"
    end
end

main() = mktempdir(make_artifacts)


text:
export satellite_plan_mdp_mcts
function satellite_plan_mdp_mcts(problem::SatPlanningProblem; enable_resources::Bool=false)

    # Sort opportunities
    sort!(problem.opportunities, by = x -> x.t_start)

    # Set initial state
    init_opp = problem.opportunities[1]
    state = initialstate(problem, Random.MersenneTwister(4))
    states = SatMDPState[state]
    plan = Opportunity[state.last_action]
    plan_rewards = Real[]
    action = nothing
    total_reward = 0.0
    r  = 0.0

    solver = MCTSSolver(n_iterations=problem.mcts_n_iterations, depth=problem.solve_depth, exploration_constant=problem.mcts_exploration_constant)
    # @POMDPs.requirements_info(solver, problem)
    policy = POMDPs.solve(solver, problem);

    state = initialstate(problem, Random.MersenneTwister(4))
    push!(states, state)

    # Compute initial state reward
    if typeof(state.last_action) == Collect
        r += state.last_action.location.reward
    end
    push!(plan_rewards, r)
    total_reward += r

    while !isterminal(problem, state)
        # Compute optimal action
        action = POMDPs.action(policy, state)

        # Advance state with next action
        next_state = POMDPs.gen(problem, state, action, Random.MersenneTwister(4)).sp
        push!(states, state)

        # Compute action reward
        r = POMDPs.reward(problem, state, action)

        # # Check next-state is safe, otherwise override with safe action
        # if next_state.power <= 0.0
        #     # Get Sunpointed action
        #     candidate_actions = problem.lt_feasible_actions[(state.last_cdo_action.id, state.last_action.id)]
        #     action = candidate_actions[findfirst(x -> typeof(x) == Sunpoint, candidate_actions)]

        #     next_state = POMDPs.gen(problem, state, action, Random.MersenneTwister(4)).sp
        #     r = POMDPs.reward(problem, state, action)
        # end

        # Alternate method for regenerating action
        if r <= 0.0 && enable_resources == true
            # Get Sunpointed action
            candidate_actions = problem.lt_feasible_actions[(state.last_cdo_action.id, state.last_action.id)]
            action = candidate_actions[findfirst(x -> typeof(x) == Sunpoint, candidate_actions)]

            next_state = POMDPs.gen(problem, state, action, Random.MersenneTwister(4)).sp

            # Compute estimated reward action
            r = POMDPs.reward(problem, state, action)
        end

        # Update State, action, reward
        state = next_state
        push!(plan, action)
        push!(states, state)
        push!(plan_rewards, r)
        total_reward += r
    end

    return states, plan, plan_rewards, total_reward
end


text:
# ---
# title: 207. Course Schedule
# id: problem207
# author: Indigo
# date: 2021-01-20
# difficulty: Medium
# categories: Depth-first Search, Breadth-first Search, Graph, Topological Sort
# link: <https://leetcode.com/problems/course-schedule/description/>
# hidden: true
# ---
#
# There are a total of `numCourses` courses you have to take, labeled from `0`
# to `numCourses-1`.
#
# Some courses may have prerequisites, for example to take course 0 you have to
# first take course 1, which is expressed as a pair: `[0,1]`
#
# Given the total number of courses and a list of prerequisite **pairs** , is it
# possible for you to finish all courses?
#
#
#
# **Example 1:**
#
#
#
#     Input: numCourses = 2, prerequisites = [[1,0]]
#     Output: true
#     Explanation:  There are a total of 2 courses to take.
#                  To take course 1 you should have finished course 0. So it is possible.
#
#
# **Example 2:**
#
#
#
#     Input: numCourses = 2, prerequisites = [[1,0],[0,1]]
#     Output: false
#     Explanation:  There are a total of 2 courses to take.
#                  To take course 1 you should have finished course 0, and to take course 0 you should
#                  also have finished course 1. So it is impossible.
#
#
#
#
# **Constraints:**
#
#   * The input prerequisites is a graph represented by **a list of edges** , not adjacency matrices. Read more about [how a graph is represented](https://www.khanacademy.org/computing/computer-science/algorithms/graph-representation/a/representing-graphs).
#   * You may assume that there are no duplicate edges in the input prerequisites.
#   * `1 <= numCourses <= 10^5`
#
#
## @lc code=start
using LeetCode
function can_finish_courses(num_courses, prerequisites::Vector{Vector{Int}})
    in_deg = fill(0, num_courses)
    point_to = [Int[] for _ in 1:num_courses] # fill(Int[], num_courses) is not proper!
    for edge in prerequisites
        in_deg[edge[1]] += 1
        push!(point_to[edge[2]], edge[1])
    end
    q = Queue{Int}()
    for idx in 1:length(in_deg)
        (in_deg[idx] == 0) && (enqueue!(q, idx))
    end
    while !isempty(q)
        frt = dequeue!(q)
        for pt in point_to[frt]
            (in_deg[pt] -= 1) == 0 && enqueue!(q, pt)
        end
    end
    return all(==(0), in_deg)
end
## @lc code=end


text:
#!/usr/bin/env julia

##########################################################################
# run_performancetest.jl
#
# Compares runtimes of different methods.
##########################################################################

using CoordinatedPrecoding
using Compat

##########################################################################
# Set up logging
Lumberjack.configure(Lumberjack._lumber_mill.timber_trucks["console"]; mode = "warn")
Lumberjack.add_truck(Lumberjack.LumberjackTruck("debug.log", "debug"), "debug")

##########################################################################
# General settings
srand(8071232234)

##########################################################################
# Performance test
simulation_params = @compat Dict(
    "I" => 3, "Kc" => 1, "N" => 2, "M" => 2,
    "P_dBm" => 20.,
    "d" => 1,
    "Ntest" => 100,
    "precoding_methods" => [
        Shi2011_WMMSE,
        Gomadam2008_MaxSINR,
        Komulainen2013_WMMSE,
        # Razaviyayn2013_MinMaxWMMSE,
        Eigenprecoding
    ],
    "aux_precoding_params" => Dict(
        "initial_precoders" => "dft",
        "stop_crit" => 0.,
        "max_iters" => 100,
    )
)
network =
    setup_interfering_broadcast_channel(simulation_params["I"],
        simulation_params["Kc"], simulation_params["N"], simulation_params["M"],
        transmit_power=10^(simulation_params["P_dBm"]/10),
        num_streams=simulation_params["d"])

timing(network, simulation_params)


text:
#> This code implements the MLJ model interface for models in the
#> DecisionTree.jl package. It is annotated so that it may serve as a
#> template for other supervised models of type `Probabilistic`. The
#> annotations, which begin with "#>", should be removed (but copy
#> this file first!). See also the model interface specification at
#> "doc/adding_new_models.md".

#> Note that all models need to "register" their location by setting
#> `load_path(<:ModelType)` appropriately.

module DecisionTree_

#> export the new models you're going to define (and nothing else):
export DecisionTreeClassifier, DecisionTreeRegressor

import MLJBase

#> needed for all classifiers:
using CategoricalArrays

#> import package:
import ..DecisionTree # strange syntax b/s we are lazy-loading

# here T is target type, and the `Vector{T}` is for storing target levels:
const DecisionTreeClassifierFitResultType{T} =
    Tuple{Union{DecisionTree.Node{Float64,T}, DecisionTree.Leaf{T}}, Vector{T}}

"""
    DecisionTreeClassifer(; kwargs...)

CART decision tree classifier from
[https://github.com/bensadeghi/DecisionTree.jl/blob/master/README.md](https://github.com/bensadeghi/DecisionTree.jl/blob/master/README.md). Predictions
are probabilistic.

For post-fit pruning, set `post-prune=true` and set
`min_purity_threshold` appropriately. Other hyperparameters as per
package documentation cited above.

"""
mutable struct DecisionTreeClassifier{T} <: MLJBase.Probabilistic{DecisionTreeClassifierFitResultType{T}}
    target_type::Type{T}  # target is CategoricalArray{target_type}
    pruning_purity::Float64
    max_depth::Int
    min_samples_leaf::Int
    min_samples_split::Int
    min_purity_increase::Float64
    n_subfeatures::Float64
    display_depth::Int
    post_prune::Bool
    merge_purity_threshold::Float64
end

# constructor:
#> all arguments are kwargs with a default value
function DecisionTreeClassifier(
    ; target_type=Int
    , pruning_purity=1.0
    , max_depth=-1
    , min_samples_leaf=1
    , min_samples_split=2
    , min_purity_increase=0.0
    , n_subfeatures=0
    , display_depth=5
    , post_prune=false
    , merge_purity_threshold=0.9)

    model = DecisionTreeClassifier{target_type}(
        target_type
        , pruning_purity
        , max_depth
        , min_samples_leaf
        , min_samples_split
        , min_purity_increase
        , n_subfeatures
        , display_depth
        , post_prune
        , merge_purity_threshold)

    message = MLJBase.clean!(model)       #> future proof by including these
    isempty(message) || @warn message #> two lines even if no clean! defined below

    return model
end

#> The following optional method (the fallback does nothing, returns
#> empty warning) is called by the constructor above but also by the
#> fit methods below:
function MLJBase.clean!(model::DecisionTreeClassifier)
    warning = ""
    if  model.pruning_purity > 1
        warning *= "Need pruning_purity < 1. Resetting pruning_purity=1.0.\n"
        model.pruning_purity = 1.0
    end
    if model.min_samples_split < 2
        warning *= "Need min_samples_split < 2. Resetting min_samples_slit=2.\n"
        model.min_samples_split = 2
    end
    return warning
end

#> A required `fit` method returns `fitresult, cache, report`. (Return
#> `cache=nothing` unless you are overloading `update`)
function MLJBase.fit(model::DecisionTreeClassifier{T2}
             , verbosity::Int   #> must be here (and typed!!) even if not used (as here)
             , X
             , y::CategoricalVector{T}) where {T,T2}

    T == T2 || throw(ErrorException("Type, $T, of target incompatible "*
                                    "with type, $T2, of $model."))

    Xmatrix = MLJBase.matrix(X)

    classes = levels(y) # *all* levels in pool of y, not just observed ones
    decoder = MLJBase.CategoricalDecoder(y)
    y_plain = MLJBase.transform(decoder, y)

    tree = DecisionTree.build_tree(y_plain
                                   , Xmatrix
                                   , model.n_subfeatures
                                   , model.max_depth
                                   , model.min_samples_leaf
                                   , model.min_samples_split
                                   , model.min_purity_increase)
    if model.post_prune
        tree = DecisionTree.prune_tree(tree, model.merge_purity_threshold)
    end

    verbosity < 3 || DecisionTree.print_tree(tree, model.display_depth)

    fitresult = (tree, classes)

    #> return package-specific statistics (eg, feature rankings,
    #> internal estimates of generalization error) in `report`, which
    #> should be `nothing` or a dictionary keyed on symbols.

    cache = nothing
    report = nothing

    return fitresult, cache, report

end

function MLJBase.predict(model::DecisionTreeClassifier{T}
                     , fitresult
                     , Xnew) where T
    Xmatrix = MLJBase.matrix(Xnew)
    tree, classes = fitresult

    # apply_tree_proba returns zero probabilities on levels unseen in
    # train, so we can give it all the levels in the pool of the
    # training vector:
    y_probabilities = DecisionTree.apply_tree_proba(tree, Xmatrix, classes)
    return [MLJBase.UnivariateNominal(classes, y_probabilities[i,:])
            for i in 1:size(y_probabilities, 1)]
end

# metadata:
MLJBase.load_path(::Type{<:DecisionTreeClassifier}) = "MLJModels.DecisionTree_.DecisionTreeClassifier"
MLJBase.package_name(::Type{<:DecisionTreeClassifier}) = "DecisionTree"
MLJBase.package_uuid(::Type{<:DecisionTreeClassifier}) = "7806a523-6efd-50cb-b5f6-3fa6f1930dbb"
MLJBase.package_url(::Type{<:DecisionTreeClassifier}) = "https://github.com/bensadeghi/DecisionTree.jl"
MLJBase.is_pure_julia(::Type{<:DecisionTreeClassifier}) = true
MLJBase.input_scitypes(::Type{<:DecisionTreeClassifier}) = MLJBase.Continuous
MLJBase.target_scitype(::Type{<:DecisionTreeClassifier}) = MLJBase.Multiclass
MLJBase.input_is_multivariate(::Type{<:DecisionTreeClassifier}) = true

end # module


text:
"""
Author: Dylan Morris
Date: 10 Feb 2022 at 6:15:09 pm
Description:
Use this script to quickly plot diagnostics of a chain. This could involve scatter-plots
to assess correlations between parameters, traceplots to check for divergent behaviours and
also just assessing the overall reasonableness of the posterior samples.
"""

let
    # for dataframe handling
    using DataFrames
    using CSV
    using Chain

    # nice plotting
    using StatsPlots
    using PlotThemes
    using Plots
    using LaTeXStrings

    # analysis/commonly used
    using Statistics
    using StatsBase
    using Dates
    using KernelDensity
    using Random
    using Distributions
    using BenchmarkTools
    using ProgressBars
end

gr()

default(
	linewidth=1.5,
	label=nothing,
	framestyle=:box,
)
scalefontsizes(0.95)

##

samples = CSV.read(
    "results/2022-04-05/posterior_sample_2022-04-05.csv",
    DataFrame,
)

par = "phi[1]"
par = "prop_md_third[391]"
par = "prop_md_third[391]"

num_samples = nrow(samples)
num_chains = 4
m = num_samples ÷ num_chains

par = "mu_hat_third[392]"
par = "md_third[392]"
par = "theta_md"

fig = plot(layout = (2, 1))
for i in 1:4
    idx = 1 + m * (i - 1):m * i
    plot!(fig, subplot = 1, samples[idx, par])
    plot!(fig, subplot = 2, kde(samples[idx, par]))
end
display(fig)

plot(samples[1:1000,par])
plot!(samples[1001:2000,par])
plot!(samples[2001:3000,par])
plot!(samples[3001:4000,par])

plot(kde(samples[1:1000,par]))
plot!(kde(samples[1001:2000,par]))
plot!(kde(samples[2001:3000,par]))
plot!(kde(samples[3001:4000,par]))

scatter(samples[!, "m1[7]"], samples[!, "r[7]"])

x = [samples[!, "phi[" * string(i) * "]"] for i in range(1, 4)]

fig = plot(layout = 12, legend = false)
for i in 369:380
    plot!(fig, subplot = i - 368, samples[!, "prop_md_third_wave." * string(i)])
end
fig

fig = plot(layout = 12, legend = false)
for i in 369:380
    plot!(fig, subplot = i - 368, kde(samples[!, "prop_md_third_wave." * string(i)]))
end
fig

plot(samples[1:1000, "prop_md_third_wave.370"])
plot!(samples[1001:2000, "prop_md_third_wave.370"])
plot!(samples[2001:3000, "prop_md_third_wave.370"])
plot!(samples[3001:4000, "prop_md_third_wave.370"])

plot(kde(samples[1:1000, "prop_md_third_wave.370"]))
plot!(kde(samples[1001:2000, "prop_md_third_wave.370"]))
plot!(kde(samples[2001:3000, "prop_md_third_wave.370"]))
plot!(kde(samples[3001:4000, "prop_md_third_wave.370"]))

plot(samples[1:1000, "bet.5"])
plot!(samples[1001:2000, "bet.5"])
plot!(samples[2001:3000, "bet.5"])
plot!(samples[3001:4000, "bet.5"])

plot((cumsum(samples[1:1000, "bet.5"])) ./ (1:1000))
plot!((cumsum(samples[1001:2000, "bet.5"])) ./ (1:1000))
plot!((cumsum(samples[2001:3000, "bet.5"])) ./ (1:1000))
plot!((cumsum(samples[3001:4000, "bet.5"])) ./ (1:1000))

plot(kde(samples[1:1000, "bet.4"]))

plot!(kde(samples[1001:2000, "bet.4"]))
plot!(kde(samples[2001:3000, "bet.4"]))
plot!(kde(samples[3001:4000, "bet.4"]))

plot(samples[!, "theta_md"])
plot(samples[1:2000, "theta_md"])
plot!(samples[2001:4000, "theta_md"])
plot!(samples[4001:6000, "theta_md"])
plot!(samples[6001:8000, "theta_md"])

plot(kde(samples[!, "theta_md"]))

# scatter(samples[!, "tau.8"], samples[!, "r.8"])
# xlabel!("days after 15/11/2021")
# ylabel!("rate of increase for the sigmoid")

# f(x, m0, m1, r, tau) = m0 + (m1 - m0) / (1 + exp(-r * (x - tau)))
# (m01, m11, r1, tau1) = (0.05, 0.95, 0.025, 25)
# f1(x) = f(x, m01, m11, r1, tau1)
# (m02, m12, r2, tau2) = (0.05, 0.95, 0.6, 17)
# f2(x) = f(x, m02, m12, r2, tau2)

# plot(0:50, f1.(0:50))
# plot!(0:50, f2.(0:50))

using MCMCDiagnostics

function ess(X)

    # get number of samples and parameters
    (samp, par) = size(X)

    # vector for storing ESS
    n_eff = zeros(par)

    # lags to evaluate the correlations at
    lags = convert.(Int, round.(0:ceil(samp / 2), digits = 0))

    # now evaluate ess for each parameter
    for ii = 1:par
        ac1 = autocor(X[:,ii], lags)
        cut = findfirst(x -> x < 0, ac1)
        pt_total = sum(ac1[2:cut])
        n_eff[ii] = samp / (1 + 2 * pt_total)
    end

    return n_eff

end

cut = findfirst(x -> x < 0, ac1)
pt_total = sum(ac1[2:cut])
n_eff[ii] = samp / (1 + 2 * pt_total)

samples_only = samples[1:1000, 9:end]
samples_mat = Matrix(samples_only)

# get number of samples and parameters
(samp, par) = size(samples_mat)

# vector for storing ESS
n_eff = zeros(par)

# lags to evaluate the correlations at
lags = convert.(Int, round.(0:ceil(samp / 2), digits = 0))

# now evaluate ess for each parameter
ii = 5
ac1 = autocor(samples_mat[:,ii], lags)
plot(ac1[1:10])

ess_samples_only = ess(samples_mat)
# some negative ess???
bad_samples = findall(ess_samples_only .< 0)

plot(samples_mat[:, bad_samples[1:5]], layout = 5)

samples_only[!,9136:9138]

plot_kde = true
plot_traceplots = false

# filter by divergent samples
# samples_good = filter(:divergent__ => ==(0.0), samples)
# samples_divergent = filter(:divergent__ => ==(1.0), samples)

# plot traceplots for ALL parameters in samples (this will take a few minutes and
# will produce quite a large pdf).
samples = CSV.read(
    "results/2022-03-22/50_case_ascertainment/posterior_sample_2022-03-22.csv",
    DataFrame,
)

let

    num_samples = nrow(samples)
    num_each_chain = num_samples ÷ 4

    sp = 1
    page = 1

    names_to_plot = names(samples)

    fig = plot(
        layout = (8,5),
        dpi=200,
        size=(1000,1200),
        # link=:x,
        framestyle=:minimal,
        legend = false,
    )

    for name in ProgressBar(names_to_plot)
        if sp == 1
            fig = plot(
                layout = (8,5),
                dpi=200,
                size=(1000,1200),
                # link=:x,
                framestyle=:minimal,
                legend = false,
            )
        end
        # plot!(fig, subplot = sp, samples[!,name])
        for i in 1:4
            plot!(fig, subplot = sp, kde(samples[(1 + (i - 1) * num_each_chain):(i * num_each_chain),name]))
        end
        xlabel!(fig, subplot = sp, name)
        sp += 1
        if sp == 8*5 + 1 || sp == length(names_to_plot) + 1
            sp = 1
            savefig(fig, "tmp_plots/page"*string(page)*".pdf")
            page += 1
        end
    end

    using PDFmerger

    pdf_filenames = [
        "tmp_plots/page"*string(p)*".pdf" for p in 1:page-1
    ]

    merge_pdfs(
        pdf_filenames,
        "tmp_plots/kde.pdf",
        cleanup=true,
    )

    sp = 1
    page = 1

    names_to_plot = names(samples)

    fig = plot(
        layout = (8,5),
        dpi=200,
        size=(1000,1200),
        # link=:x,
        framestyle=:minimal,
        legend = false,
    )

    for name in ProgressBar(names_to_plot)
        if sp == 1
            fig = plot(
                layout = (8,5),
                dpi=200,
                size=(1000,1200),
                # link=:x,
                framestyle=:minimal,
                legend = false,
            )
        end
        # plot!(fig, subplot = sp, samples[!,name])
        for i in 1:4
            plot!(fig, subplot = sp, samples[(1 + (i - 1) * num_each_chain):(i * num_each_chain),name])
        end
        xlabel!(fig, subplot = sp, name)
        sp += 1
        if sp == 8*5 + 1 || sp == length(names_to_plot) + 1
            sp = 1
            savefig(fig, "tmp_plots/page"*string(page)*".pdf")
            page += 1
        end
    end

    using PDFmerger

    pdf_filenames = [
        "tmp_plots/page"*string(p)*".pdf" for p in 1:page-1
    ]

    merge_pdfs(
        pdf_filenames,
        "tmp_plots/traceplots.pdf",
        cleanup=true,
    )

end

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

# Code for assessing the fitted TP (μ(t)) coming out of stan against the forecasted
# estimates. More of a sanity check than anything else and just helps to see that the
# posterior sampled μ(t) and recombined/forecasted estimates align.

state = "TAS"

# mu_hat_filtered_old = CSV.read(
#     "results/UoA_forecast_output/2022-02-15/" * state * "_2022-02-15_TP.csv",
#     DataFrame,
# )

mu_hat_local = CSV.read(
    "mu_hat_Australian Capital Territory.csv",
    DataFrame,
    drop = [1],
)

mu_hat_all = CSV.read(
    "results/soc_mob_R2022-02-22.csv",
    DataFrame,
    drop = [1],
)

mu_hat_forecast = @chain mu_hat_all begin
    subset(
        :type => ByRow(==("R_L")),
        :state => ByRow(==(state)),
    )
end

mu_hat_filtered = CSV.read(
    "results/UoA_forecast_output/2022-02-22/" * state * "_2022-02-22_TP.csv",
    DataFrame,
)

Reff_delta = CSV.read(
    "results/EpyReff/Reff_delta2022-03-08tau_5.csv",
    DataFrame
)
Reff_omicron = CSV.read(
    "results/EpyReff/Reff_omicron2022-03-08tau_5.csv",
    DataFrame
)

fig = plot(layout = (4, 2), size = (800, 800))
plot!(fig, Reff_delta.INFECTION_DATES, Reff_delta.mean, group = Reff_delta.STATE, lc = 1)
plot!(fig, Reff_delta.INFECTION_DATES, Reff_delta.bottom, group = Reff_delta.STATE, lc = 1, ls = :dash)
plot!(fig, Reff_delta.INFECTION_DATES, Reff_delta.top, group = Reff_delta.STATE, lc = 1, ls = :dash)
plot!(fig, Reff_omicron.INFECTION_DATES, Reff_omicron.mean, group = Reff_omicron.STATE, lc = 2)
plot!(fig, Reff_omicron.INFECTION_DATES, Reff_omicron.bottom, group = Reff_omicron.STATE, lc = 2, ls = :dash)
plot!(fig, Reff_omicron.INFECTION_DATES, Reff_omicron.top, group = Reff_omicron.STATE, lc = 2, ls = :dash)
xlims!((Dates.Date("2021-12-01"), Dates.Date("2022-03-05")))
ylims!(0.5, 3)

Reff_delta = filter(:STATE => ==("NSW"), Reff_delta)
Reff_omicron = filter(:STATE => ==("NSW"), Reff_omicron)

fig = plot(size = (800, 800))
plot!(fig, Reff_delta.INFECTION_DATES, Reff_delta.mean, lc = 1)
plot!(fig, Reff_delta.INFECTION_DATES, Reff_delta.bottom, lc = 1, ls = :dash)
plot!(fig, Reff_delta.INFECTION_DATES, Reff_delta.top, lc = 1, ls = :dash)
plot!(fig, Reff_omicron.INFECTION_DATES, Reff_omicron.mean, lc = 2)
plot!(fig, Reff_omicron.INFECTION_DATES, Reff_omicron.bottom, lc = 2, ls = :dash)
plot!(fig, Reff_omicron.INFECTION_DATES, Reff_omicron.top, lc = 2, ls = :dash)
vline!(fig, [Dates.Date("2021-12-13")], lc = :black, ls = :dash, lw = 2)
vline!(fig, [Dates.Date("2021-11-15")], lc = :green, ls = :dash, lw = 2)
xlims!((Dates.Date("2021-11-10"), Dates.Date("2022-03-05")))
ylims!(0.5, 1.7)

Reff_state = filter(:STATE => ==(state), Reff)
plot(
    mu_hat_forecast.date,
    mean(Matrix(mu_hat_forecast[!,9:end]), dims = 2),
    label = "mu"
)
mu_hat_filtered_median = median(Matrix(mu_hat_filtered[:,3:end-1]), dims = 2)
plot!(
    mu_hat_filtered[:,"onset date"],
    mu_hat_filtered_median,
    label = "filtered mu"
)
mu_hat_local_median = median(Matrix(mu_hat_local[:,3:end-1]), dims = 1)[:]
t0 = Dates.Date("2021-12-01")
Δt = Dates.Day(1)
t1 = t0 + Dates.Day(length(mu_hat_local_median)-1)
dr = t0:Δt:t1
plot!(
    dr,
    mu_hat_local_median,
    label = "fitted mu"
)
plot!(
    Reff_state[:,"INFECTION_DATES"],
    Reff_state[:,"median"],
    label = "Reff"
)
vline!(
    [Dates.Date("2021-12-19")]
)
ylims!((0,2))
# xlims!((NSW_date_range[1], NSW_date_range[end] + Dates.Day(35)))
xlims!((Dates.Date("2021-12-15"), Dates.Date("2021-12-15") + Dates.Day(60)))

Reff_state = filter(:STATE => ==(state), Reff)
plot(
    mu_hat_forecast.date,
    Matrix(mu_hat_forecast[!,9:end]),
    label = "mu",
    lc = 1,
    leg = false,
)
# mu_hat_filtered_median = median(Matrix(mu_hat_filtered[:,3:end-1]), dims = 2)
# plot!(
#     mu_hat_filtered[:,"onset date"],
#     mu_hat_filtered_median,
#     label = "filtered mu"
# )
mu_hat_local_mat = Matrix(mu_hat_local[:,3:end-1])'
t0 = Dates.Date("2021-12-01")
Δt = Dates.Day(1)
t1 = t0 + Dates.Day(length(mu_hat_local_median)-1)
dr = t0:Δt:t1
plot!(
    dr,
    mu_hat_local_mat,
    lc = 2,
)
vline!(
    [Dates.Date("2021-12-19")]
)
ylims!((0,2))
# xlims!((NSW_date_range[1], NSW_date_range[end] + Dates.Day(35)))
xlims!((Dates.Date("2021-12-15"), Dates.Date("2021-12-15") + Dates.Day(60)))

mu_hat_VIC = CSV.read(
    "./mu_hat_VIC.csv",
    DataFrame,
    drop = [1],
)

mu_hat_all = CSV.read(
    "results/soc_mob_R2022-02-01.csv",
    DataFrame,
    drop = [1],
)

mu_hat_VIC_filtered = CSV.read(
    "results/UoA_forecast_output/2022-02-09/VIC_2022-02-09_TP.csv",
    DataFrame,
)

mu_hat_VIC_forecast = @chain mu_hat_all begin
    subset(
        :type => ByRow(==("R_L")),
        :state => ByRow(==("VIC")),
    )
end

VIC_date_range = Dates.Date("2021-08-15"):Dates.Day(1):Dates.Date("2021-06-15")+Dates.Day(nrow(mu_hat_VIC)-1)
plot(
    VIC_date_range,
    median(Matrix(mu_hat_VIC), dims = 2),
)
plot!(
    mu_hat_VIC_forecast.date,
    mu_hat_VIC_forecast.median,
)

mu_hat_VIC_filtered_median = median(Matrix(mu_hat_VIC_filtered[:,3:end-1]), dims = 2)
plot!(
    mu_hat_VIC_filtered[:,"onset date"],
    mu_hat_VIC_filtered_median,
)
xlims!((VIC_date_range[1], VIC_date_range[end]))

# pars = [
#     "bet.1",
#     "bet.2",
#     "bet.3",
#     "bet.4",
#     "bet.5",
#     "theta_masks",
#     "voc_effect_delta",
#     "voc_effect_omicron",
#     "susceptible_depletion_factor",
# ]

# sp = 1
# page = 1

# fig = plot(
#     layout = 9,
#     dpi=200,
#     size=(1000,1000),
#     framestyle=:minimal,
#     legend = false,
# )

# for x in ProgressBar(pars)
#     for y in pars
#         if sp == 1
#             fig = plot(
#                 layout = 9,
#                 dpi=200,
#                 size=(1000,1000),
#                 framestyle=:minimal,
#                 legend = false,
#             )
#         end
#         scatter!(
#             fig, subplot = sp, samples[:,x], samples[:,y]
#         )
#         xlabel!(fig, subplot = sp, x)
#         ylabel!(fig, subplot = sp, y)
#         sp += 1
#     end
#     sp = 1
#     savefig(fig, "tmp_plots/page"*string(page)*".pdf")
#     page += 1
# end

# pdf_filenames = [
#     "tmp_plots/page"*string(p)*".pdf" for p in 1:page-1
# ]

# merge_pdfs(
#     pdf_filenames,
#     cleanup=true,
# )

# f = plot(layout = 2)
# scatter!(f, subplot = 1, samples[:,"voc_effect_delta"], samples[:, "theta_md"])
# xlabel!(f, subplot = 1, "voc_effect_delta")
# ylabel!(f, subplot = 1, "theta_md")
# scatter!(f, subplot = 2, samples[:,"voc_effect_omicron"], samples[:, "theta_md"])
# xlabel!(f, subplot = 2, "voc_effect_omicron")
# ylabel!(f, subplot = 2, "theta_md")
# savefig(f, "correlations.pdf")

mu_hat_filtered_old = CSV.read(
    "results/UoA_forecast_output/2022-02-02/NT_2022-02-02_TP.csv",
    DataFrame,
)

mu_hat_filtered = CSV.read(
    "results/UoA_forecast_output/2022-02-09/NSW_2022-02-09_TP.csv",
    DataFrame,
)

# nice hack for datetime range
start_date = Dates.Date("2021-06-15")
Δt = Dates.Day(1)
end_date = start_date + Dates.Day(nrow(mu_hat)-1)
date_range = start_date:Δt:end_date

mu_hat_filtered_old_median = median(Matrix(mu_hat_filtered_old[:,3:end-1]), dims = 2)
mu_hat_filtered_median = median(Matrix(mu_hat_filtered[:,3:end-1]), dims = 2)
plot(
    mu_hat_filtered_old[:,"onset date"],
    mu_hat_filtered_old_median,
)
plot!(
    mu_hat_filtered[:,"onset date"],
    mu_hat_filtered_median,
)
# xlims!((NSW_date_range[1], NSW_date_range[end] + Dates.Day(35)))

D = rand(100, 3, 1000)
case_counts = zeros(100)
case_counts1 = zeros(100)
@btime $case_counts1 .= sum(@view($D[:, 1:2, 1]), dims = 2)

function count_cases!(case_counts, D, sim)

    for i in 1:length(case_counts)
        case_counts[i] = D[i,1,sim] + D[i,2,sim]
    end

    return nothing
end

@btime count_cases!($case_counts, $D, 1)


text:
# =================================================
# DummySet: a set to test default methods
# =================================================
D = DummySet()
@test space_dimension(D) == nothing
@test dimension(D) == nothing

# =================================================
# Singleton: a singleton is a set with one element
# =================================================
s = Singleton([1.0, 2.0, 3.0])

@test space_dimension(s) == 3
@test dimension(s) == 0

# =================================================
# Ball: a ball in n-dimensional space
# =================================================
B = Ball(zeros(3), 1.0)

@test space_dimension(B) == 3
@test dimension(B) == 3


text:
module QiskitQuantumInfo

import QuantumOps

export Pauli, symplectic
export PauliList, SparsePauliOp
export hellinger_distance!, hellinger_distance

include("pauli.jl")
include("pauli_list.jl")
include("distance.jl")
include("util.jl")

end


text:
# A tolerance graph is a graph whose vertices are represented by pairs
# (J,t) where J is a closed interval and t is a real number. Vertices
# are adjacent if and only if length(J * JJ) >= min(t,tt).

export ToleranceGraph

"""
`ToleranceGraph(Jlist, tlist)` creates a tolerance graph where `Jlist`
is a list of intervals and `tlist` is a list of tolerances.
"""
function ToleranceGraph(
    Jlist::Vector{ClosedInterval{S}},
    tlist::Vector{T},
) where {S<:Real,T<:Real}
    n = length(Jlist)
    G = IntGraph(n)

    for i = 1:n-1
        J = Jlist[i]
        t = tlist[i]
        for j = i+1:n
            JJ = Jlist[j]
            tt = tlist[j]

            K = J * JJ
            if !isempty(K) && length(K) >= min(t, tt)
                add!(G, i, j)
            end
        end
    end
    return G
end

"""
`ToleranceGraph(f)` creates a tolerance graph where `f` is a `Dict`
mapping vertices to pairs `(J,t)` where `J` is a `ClosedInterval` and
`t` is a `Real` tolerance.
"""
function ToleranceGraph(f::Dict{X,Tuple{ClosedInterval{S},T}}) where {S<:Real,T<:Real,X}
    vtcs = collect(keys(f))
    G = SimpleGraph{X}()
    for v in vtcs
        add!(G, v)
    end

    n = length(vtcs)
    for i = 1:n-1
        v = vtcs[i]
        J, t = f[v]
        for j = i+1:n
            w = vtcs[j]
            JJ, tt = f[w]

            K = J * JJ
            if !isempty(K) && length(K) >= min(t, tt)
                add!(G, v, w)
            end
        end
    end
    return G
end


text:
using Distributed;
addprocs(20);
@everywhere using PowerModels, Gurobi, Ipopt, Random, Distributions;
@everywhere using CSV;
@everywhere include("parse_dict.jl");
@everywhere using DataFrames;

list_file = readdir("./data/")
perturb_file = readdir("./data/perturb/")
gen_data_file = readdir("./data/perturb/perturb_time_data/")
stop_time = Float64[100, 200, 400, 600, 800, 1000] # the time to stop solver
round = [i for i in 1:1000]

case_name = AbstractString[]
objective_value_100 = []
objective_value_200 = []
objective_value_400 = []
objective_value_600 = []
objective_value_800 = []
objective_value_1000 = []
dual_inf_100 = []
dual_inf_200 = []
dual_inf_400 = []
dual_inf_600 = []
dual_inf_800 = []
dual_inf_1000 = []
cv_100 = []
cv_200 = []
cv_400 = []
cv_600 = []
cv_800 = []
cv_1000 = []
nlp_err_100 = []
nlp_err_200 = []
nlp_err_400 = []
nlp_err_600 = []
nlp_err_800 = []
nlp_err_1000 = []

function exist_in_datafile(gen_data_file, itema)
    for time in stop_time
        data_file_name = itema[1:end-2] * "_" * string(time) * ".csv"
        if (data_file_name in gen_data_file) == false
            return false
        end
    end
    return true

end

function exist_perturb_file(perturb_file, itema)
    perturb_data_file = itema[1:end-2] * "_" * "perturb" * ".csv"
    if perturb_data_file in perturb_file
        return true
    end
    return false
end

for itema in list_file
    if occursin(".m", itema) && (itema == "case5.m")
        if exist_perturb_file(perturb_file, itema)
            if (exist_in_datafile(gen_data_file, itema)) == false

                network_data = parse_file("./data/" * itema)

                pd_pro_ls = []

                df = CSV.read("./data/perturb/" * string(itema[1:end-2]) * "_" * "perturb.csv", DataFrame)

                for pd_profile in df[!, "pd"]
                    push!(pd_pro_ls, pd_profile[16:end-1])
                end

                pd_ls = []
                for pd_pro in pd_pro_ls
                    push!(pd_ls, parse2dict(pd_pro))

                end

                for iter in round

                    network_data_new = deepcopy(network_data)

                    gen_keys = [i for i in keys(network_data["gen"])]

                    #relax the inf term
                    for i in 1:length(gen_keys)
                        if network_data_new["gen"][gen_keys[i]]["qmax"] == Inf
                            network_data_new["gen"][gen_keys[i]]["qmax"] = 10^20
                        end

                        if network_data_new["gen"][gen_keys[i]]["qmax"] == -Inf
                            network_data_new["gen"][gen_keys[i]]["qmax"] = -10^20
                        end

                        if network_data_new["gen"][gen_keys[i]]["qmin"] == Inf
                            network_data_new["gen"][gen_keys[i]]["qmin"] = 10^20
                        end

                        if network_data_new["gen"][gen_keys[i]]["qmin"] == -Inf
                            network_data_new["gen"][gen_keys[i]]["qmin"] = -10^20
                        end
                    end

                    if haskey(network_data_new, "dcline")
                        dcline_keys = [i for i in keys(network_data["dcline"])]

                        for i in 1:length(dcline_keys)
                            if network_data_new["dcline"][dcline_keys[i]]["qmaxf"] == Inf
                                network_data_new["dcline"][dcline_keys[i]]["qmaxf"] = 10^20
                            end

                            if network_data_new["dcline"][dcline_keys[i]]["qmaxf"] == -Inf
                                network_data_new["dcline"][dcline_keys[i]]["qmaxf"] = -10^20
                            end

                            if network_data_new["dcline"][dcline_keys[i]]["qmint"] == Inf
                                network_data_new["dcline"][dcline_keys[i]]["qmint"] = 10^20
                            end

                            if network_data_new["dcline"][dcline_keys[i]]["qmint"] == -Inf
                                network_data_new["dcline"][dcline_keys[i]]["qmint"] = -10^20
                            end

                            if network_data_new["dcline"][dcline_keys[i]]["qminf"] == Inf
                                network_data_new["dcline"][dcline_keys[i]]["qminf"] = 10^20
                            end

                            if network_data_new["dcline"][dcline_keys[i]]["qminf"] == -Inf
                                network_data_new["dcline"][dcline_keys[i]]["qminf"] = -10^20
                            end

                            if network_data_new["dcline"][dcline_keys[i]]["qmaxt"] == Inf
                                network_data_new["dcline"][dcline_keys[i]]["qmaxt"] = 10^20
                            end

                            if network_data_new["dcline"][dcline_keys[i]]["qmaxt"] == -Inf
                                network_data_new["dcline"][dcline_keys[i]]["qmaxt"] = -10^20
                            end
                        end
                    end
                    #modify the demand
                    pd_keys = [i for i in keys(network_data_new["load"])]
                    for key in pd_keys
                        network_data_new["load"][key]["pd"] = pd_ls[iter][parse(Float64, key)]
                    end

                    for time in stop_time
                        push!(case_name, itema)
                        out_file_name = "./data/perturb/perturb_time_data/" * itema[1:end-2] * "_" * string(time) *"_"* string(iter) * ".csv"
                        pmap(ω -> run_ac_opf(network_data_new, optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0, "max_wall_time" => time, "file_print_level" => 3,
                                "output_file" => out_file_name)), 1:1000)

                        open(out_file_name, "r") do file_object
                            for line in eachline(file_object)
                                if occursin("Objective...............", line)
                                    if time == 100
                                        push!(objective_value_100, line[29:50])
                                    end

                                    if time == 200
                                        push!(objective_value_200, line[29:50])
                                    end

                                    if time == 400
                                        push!(objective_value_400, line[29:50])
                                    end

                                    if time == 600
                                        push!(objective_value_600, line[29:50])
                                    end

                                    if time == 800
                                        push!(objective_value_800, line[29:50])
                                    end

                                    if time == 1000
                                        push!(objective_value_1000, line[29:50])
                                    end
                                end

                                if occursin("Dual infeasibility....", line)
                                    if time == 100
                                        push!(dual_inf_100, line[29:50])
                                    end

                                    if time == 200
                                        push!(dual_inf_200, line[29:50])
                                    end

                                    if time == 400
                                        push!(dual_inf_400, line[29:50])
                                    end

                                    if time == 600
                                        push!(dual_inf_600, line[29:50])
                                    end

                                    if time == 800
                                        push!(dual_inf_800, line[29:50])
                                    end

                                    if time == 1000
                                        push!(dual_inf_1000, line[29:50])
                                    end
                                end

                                if occursin("Constraint violation....", line)
                                    if time == 100
                                        push!(cv_100, line[29:50])
                                    end

                                    if time == 200
                                        push!(cv_200, line[29:50])
                                    end

                                    if time == 400
                                        push!(cv_400, line[29:50])
                                    end

                                    if time == 600
                                        push!(cv_600, line[29:50])
                                    end

                                    if time == 800
                                        push!(cv_800, line[29:50])
                                    end

                                    if time == 1000
                                        push!(cv_1000, line[29:50])
                                    end
                                end

                                if occursin("Overall NLP error....", line)
                                    if time == 100
                                        push!(nlp_err_100, line[29:50])
                                    end

                                    if time == 200
                                        push!(nlp_err_200, line[29:50])
                                    end

                                    if time == 400
                                        push!(nlp_err_400, line[29:50])
                                    end

                                    if time == 600
                                        push!(nlp_err_600, line[29:50])
                                    end

                                    if time == 800
                                        push!(nlp_err_800, line[29:50])
                                    end

                                    if time == 1000
                                        push!(nlp_err_1000, line[29:50])
                                    end
                                end

                            end
                        end

                    end
                end
            end
        end
    end
end

dict_cat = Dict("case_name" => case_name, "obj_100" => objective_value_100, "obj_200" => objective_value_200, "obj_400" => objective_value_400,
    "obj_600" => objective_value_600, "obj_800" => objective_value_800, "obj_1000" => objective_value_1000, "dual_inf_100" => dual_inf_100, "dual_inf_200" => dual_inf_200,
    "dual_inf_400" => dual_inf_400, "dual_inf_600" => dual_inf_600, "dual_inf_800" => dual_inf_800, "dual_inf_1000" => dual_inf_1000, "cv_100" => cv_100, "cv_200" => cv_200,
    "cv_400" => cv_400, "cv_600" => cv_600, "cv_800" => cv_800, "cv_1000" => cv_1000, "nlp_err_100" => nlp_err_100, "nlp_err_200" => nlp_err_200,
    "nlp_err_400" => nlp_err_400, "nlp_err_600" => nlp_err_600, "nlp_err_800" => nlp_err_800, "nlp_err_1000" => nlp_err_1000)

df_result = DataFrame(dict_cat)
CSV.write("result_perturb_time.csv", df_result, bufsize = 2^30)

max_stars_repo_name: zhouliang-yu/Grid-ML

max_issues_repo_name: zhouliang-yu/Grid-ML

max_forks_repo_name: zhouliang-yu/Grid-ML


text:
using LightGraphs

@testset "Graph creation" begin
    g = StorageGraph()
    @test get_prop(g) == 1
    @test nextid(g, (A=1,)=>(D=0.4,)) == 1
    add_nodes!(g, (A=1,)=>(D=0.4,)=>(B=0.5,))
    @test length(g.data) == length(g.index)
    @test nv(g) == 3
    @test ne(g) == 2
    @test g.data == Dict(1=>(D=0.4,), 2=>(B=0.5,), 3=>(A=1,))
    @test g.index == Dict((D=0.4,)=>1, (B=0.5,)=>2, (A=1,)=>3)
    @test g.paths == Dict(Edge(1,2)=>Set(1), Edge(3,1)=>Set(1))
    @test get_prop(g) == 2
    dep1 = (A=1,)=>(D=0.4,)=>(B=0.5,)
    dep2 = (A=1,)=>(D=0.4,)=>(B=0.6,)
    dep3 = (A=1,)=>(D=0.5,)=>(B=0.5,)
    dep4 = (A=2,)=>(D=0.4,)
    @test nextid(g, dep1) == 1
    @test nextid(g, dep2) == 2
    @test nextid(g, dep3) == 2
    @test nextid(g, dep3) == 2
end

@testset "Adding data proggresively" begin
    g = StorageGraph()

    @testset "add_nodes! simple" begin
        add_nodes!(g, (a=1,)=>(b=1,))
        @test nv(g) == 2
        @test ne(g) == 1
        @test get_prop(g, 1) == (a=1,)
        @test get_prop(g, 2) == (b=1,)
        @test get_prop(g, 1, 2) == Set(1)
        @test get_prop(g) == 2
    end

    @testset "add_nodes! continuing" begin
        add_nodes!(g, (a=1,)=>(b=1,)=>(c=1,))
        @test nv(g) == 3
        @test ne(g) == 2
        @test get_prop(g, 1) == (a=1,)
        @test get_prop(g, 2) == (b=1,)
        @test get_prop(g, 3) == (c=1,)
        @test get_prop(g, 1, 2) == get_prop(g, 2, 3) == Set(1)
        @test get_prop(g) == 2

        add_nodes!(g, (a=2,)=>(b=1,)=>(c=1,))
        @test nv(g) == 4
        @test ne(g) == 3
        @test get_prop(g, 1) == (a=1,)
        @test get_prop(g, 2) == (b=1,)
        @test get_prop(g, 3) == (c=1,)
        @test get_prop(g, 4) == (a=2,)
        @test get_prop(g, 1, 2) == Set(1)
        @test get_prop(g, 2, 3) == Set([1,2])
        @test get_prop(g, 4, 2) == Set([2])
        @test get_prop(g) == 3

        add_nodes!(g, (a=1,)=>(b=2,)=>(c=1,))
        @test nv(g) == 5
        @test ne(g) == 5
        @test get_prop(g, 1) == (a=1,)
        @test get_prop(g, 2) == (b=1,)
        @test get_prop(g, 3) == (c=1,)
        @test get_prop(g, 4) == (a=2,)
        @test get_prop(g, 5) == (b=2,)
        @test get_prop(g, 1, 2) == Set(1)
        @test get_prop(g, 2, 3) == Set([1,2])
        @test get_prop(g, 4, 2) == Set([2])
        @test get_prop(g, 1, 5) == Set(3)
        @test get_prop(g, 5, 3) == Set(3)
        @test get_prop(g) == 4

        add_nodes!(g, (a=1,)=>(b=1,)=>(c=1,)=>(d=1,))
        @test nv(g) == 6
        @test ne(g) == 6
        @test get_prop(g, 1) == (a=1,)
        @test get_prop(g, 2) == (b=1,)
        @test get_prop(g, 3) == (c=1,)
        @test get_prop(g, 4) == (a=2,)
        @test get_prop(g, 5) == (b=2,)
        @test get_prop(g, 6) == (d=1,)
        @test get_prop(g, 1, 2) == Set(1)
        @test get_prop(g, 2, 3) == Set([1,2])
        @test get_prop(g, 4, 2) == Set([2])
        @test get_prop(g, 1, 5) == Set(3)
        @test get_prop(g, 5, 3) == Set(3)
        @test get_prop(g, 3, 6) == Set(1)
        @test get_prop(g) == 4

        add_nodes!(g, (a=2,)=>(b=2,))
        @test nv(g) == 6
        @test ne(g) == 7
        @test get_prop(g, 1) == (a=1,)
        @test get_prop(g, 2) == (b=1,)
        @test get_prop(g, 3) == (c=1,)
        @test get_prop(g, 4) == (a=2,)
        @test get_prop(g, 5) == (b=2,)
        @test get_prop(g, 6) == (d=1,)
        @test get_prop(g, 1, 2) == Set(1)
        @test get_prop(g, 2, 3) == Set([1,2])
        @test get_prop(g, 4, 2) == Set([2])
        @test get_prop(g, 1, 5) == Set(3)
        @test get_prop(g, 5, 3) == Set(3)
        @test get_prop(g, 3, 6) == Set(1)
        @test get_prop(g, 4, 5) == Set(4)
        @test get_prop(g) == 5
    end

# FIXME!
    # val1 = (q₀=[[0,1],[0,2],[0,3]], q₂=[[2,1],[2,2],[2,3]])
    # @testset "add_bulk!" begin
    #     add_bulk!(g, (A=1,)=>(B=0.4,), val1)
    #     @test nv(g) == 5
    #     @test ne(g) == 4
    #     @test get_prop(g, 1) == (A=1,)
    #     @test get_prop(g, 2) == (B=0.4,)
    #     @test get_prop(g, 3) == (q₀=[0,1],q₂=[2,1])
    #     @test get_prop(g, 4) == (q₀=[0,2],q₂=[2,2])
    #     @test get_prop(g, 5) == (q₀=[0,3],q₂=[2,3])
    #     @test get_prop(g, 1, 2) == Set([1,2,3])
    #     @test get_prop(g, 2, 3) == Set(1)
    #     @test get_prop(g, 2, 4) == Set(2)
    #     @test get_prop(g, 2, 5) == Set(3)
    #
    #     @test nextid(g, (A=1,)=>(B=0.5,)=>(x=1,)) == 4
    #     add_bulk!(g, (A=1,)=>(B=0.5,), (q₀=[[0,-1],[0,-2]], q₂=[[2,-1],[2,-2]]))
    #     @test nv(g) == 8
    #     @test ne(g) == 7
    #     @test get_prop(g, 6) == (B=0.5,)
    #     @test get_prop(g, 7) == (q₀=[0,-1],q₂=[2,-1])
    #     @test get_prop(g, 8) == (q₀=[0,-2],q₂=[2,-2])
    #     @test get_prop(g, 1, 6) == Set([4,5])
    #     @test get_prop(g, 6, 7) == Set(4)
    #     @test get_prop(g, 6, 8) == Set(5)
    # end

end

@testset "Adding data from scratch" begin
    g = StorageGraph()

    @testset "add_nodes!" begin
        add_nodes!(g, (A=1,)=>(B=0.4,))
        @test nv(g) == 2
        @test ne(g) == 1
    end

    g = StorageGraph()

    @testset "add_bulk!" begin
        dep = (q₀=[[0,1],[0,2],[0,3]], q₂=[[2,1],[2,2],[2,3]])
        add_bulk!(g, (A=1,)=>(B=0.4,), dep)
        @test nv(g) == 5
        @test ne(g) == 4
    end

end


text:
#***********************************************************************
#  Copyright 2006 John A. Trangenstein
#
#  This software is made available for research and instructional use
#  only.
#  You may copy and use this software without charge for these
#  non-commercial purposes, provided that the copyright notice and
#  associated text is reproduced on all copies.
#  For all other uses (including distribution of modified versions),
#  please contact the author at
#    John A. Trangenstein
#    Department of Mathematics
#    Duke University
#    Durham, NC 27708-0320
#    USA
#  or
#    johnt@math.duke.edu
#
#  This software is made available "as is" without any assurance that it
#  is completely correct, or that it will work for your purposes.
#  Use the software at your own risk.
#***********************************************************************

module Consts

    const roundoff  =  1.0-14
    const small     =  1.0-20
    const huge      =  1.0e300
    const undefind  =  1.0e300
    const lnrndoff  =  14.0e0
    const lnsmall   =  20.0e0

    # problem-specific parameters:
    const jump      =  0.0
    const x_left    = -0.2
    const x_right   =  1.0
    const statelft  =  2.0
    const statergt  =  0.0
    const velocity  =  1.0

    export roundoff, huge, jump, x_left, x_right, statelft, statergt, velocity

end