[FTheoryTools] Implement method for well-quantized G4-fluxes

experimental/FTheoryTools/src/FTheoryTools.jl

@@ -30,6 +30,7 @@ include("G4Fluxes/constructors.jl")
 include("G4Fluxes/attributes.jl")
 include("G4Fluxes/properties.jl")
 include("G4Fluxes/special_attributes.jl")
+include("G4Fluxes/special-intersection-theory.jl")
 
 include("Serialization/tate_models.jl")
 include("Serialization/weierstrass_models.jl")

experimental/FTheoryTools/src/G4Fluxes/special-intersection-theory.jl

@@ -0,0 +1,269 @@
+# ---------------------------------------------------------------------------------------------------------
+# (1) Compute the intersection product of an algebraic cycle with a hypersurface.
+# ---------------------------------------------------------------------------------------------------------
+
+function sophisticated_intersection_product(v::NormalToricVariety, indices::NTuple{4, Int64}, hypersurface_equation::MPolyDecRingElem{QQFieldElem, QQMPolyRingElem}, intersection_dict::Dict{NTuple{4, Int64}, ZZRingElem}, special_intersection_dict::Dict{String, ZZRingElem})
+
+  # (A) Have we computed this intersection number in the past? If so, just use that result...
+  if haskey(intersection_dict, indices)
+    return intersection_dict[indices]
+  end
+
+  # (B) Get the indices of the variables that we try to intersect and, by virtue of the SR-ideal, intersect trivially
+  mnf = Oscar._minimal_nonfaces(v)
+  ignored_sets = Set([Tuple(sort(Vector{Int}(Polymake.row(mnf, i)))) for i in 1:Polymake.nrows(mnf)])
+  for sr_set in ignored_sets
+    if is_subset(sr_set, indices)
+      intersection_dict[indices] = ZZ(0)
+      return ZZ(0)
+    end
+  end
+
+  # (C) Deal with self-intersection and should-never-happen case.
+  variable_pos = Set(indices)
+  if length(variable_pos) < 4 && length(variable_pos) >= 1
+    return intersection_from_equivalent_cycle(v, indices, hypersurface_equation, intersection_dict, special_intersection_dict)
+  end
+  if length(variable_pos) == 0
+    println("WEIRD! THIS SHOULD NEVER HAPPEN! INFORM THE AUTHORS!")
+    println("")
+  end
+
+
+  # (D) Deal with transverse intersection...
+
+  # D.1 Work out the intersection locus in detail.
+  pt_reduced, gs_reduced, remaining_vars, reduced_scaling_relations = Oscar._reduce_hypersurface_equation(v, hypersurface_equation, indices)
+
+  # D.2 If pt == 0, then we are not looking at a transverse intersection. So take an equivalent cycle and try again...
+  if is_zero(pt_reduced)
+    return intersection_from_equivalent_cycle(v, indices, hypersurface_equation, intersection_dict, special_intersection_dict)
+  end
+
+  # D.3 If pt is constant and non-zero, then the intersection is trivial.
+  if is_constant(pt_reduced) && is_zero(pt_reduced) == false
+    intersection_dict[indices] = ZZ(0)
+    return ZZ(0)
+  end
+
+  # D.4 Helper function for the cases below
+  function has_one_and_rest_zero(vec)
+    return count(==(1), vec) == 1 && all(x -> x == 0 || x == 1, vec)
+  end
+
+  # C.5 Cover a case that seems to appear frequently for our investigation:
+  # pt_reduced of the form a * x + b * y for non-zero number a,b and remaining variables x, y subject to a reduced SR generator x * y and scaling relation [1,1].
+  # This will thus always give exactly one solution (x = 1, y = -a/b), and so the intersection number is one.
+  if length(gs_reduced) == 1 && length(remaining_vars) == 2
+    mons_list = collect(monomials(pt_reduced))
+    if length(mons_list) == 2
+      if all(x -> x != 0, collect(coefficients(pt_reduced)))
+        exps_list = [collect(exponents(k))[1] for k in mons_list]
+        if has_one_and_rest_zero(exps_list[1]) && has_one_and_rest_zero(exps_list[2])
+          if gs_reduced[1] == remaining_vars[1] * remaining_vars[2]
+            if reduced_scaling_relations == matrix(ZZ, [[1,1]])
+              intersection_dict[indices] = ZZ(1)
+              return ZZ(1)
+            end
+          end
+        end
+      end
+    end
+  end
+
+  # C.6 Cover a case that seems to appear frequently for our investigation:
+  # pt_reduced of the form a * x for non-zero number a and remaining variables x, y subject to a reduced SR generator x * y and scaling relation [*, != 0].
+  # This only gives the solution [0:1], so one intersection point.
+  if length(gs_reduced) == 1 && length(remaining_vars) == 2
+    mons_list = collect(monomials(pt_reduced))
+    if length(mons_list) == 1 && collect(coefficients(pt_reduced))[1] != 0
+      list_of_exps = collect(exponents(mons_list[1]))[1]
+      number_of_zeros = count(==(0), list_of_exps)
+      if number_of_zeros == length(list_of_exps) - 1
+        if gs_reduced[1] == remaining_vars[1] * remaining_vars[2]
+          if reduced_scaling_relations[1,1] != 0 && reduced_scaling_relations[1,2] != 0
+            highest_power = list_of_exps[findfirst(x -> x > 0, list_of_exps)]
+            if highest_power == 1
+              intersection_dict[indices] = highest_power
+              return highest_power
+            end
+          end
+        end
+      end
+    end
+  end
+
+  # C.7 Cover a case that seems to appear frequently for our investigation. It looks as follows:
+  # pt_reduced = -5700*w8*w10
+  # remaining_vars = MPolyDecRingElem{QQFieldElem, QQMPolyRingElem}[w8, w10]
+  # gs_reduced = MPolyDecRingElem{QQFieldElem, QQMPolyRingElem}[w8*w10]
+  # reduced_scaling_relations = [1 1]
+  # This gives exactly two solutions, namely [0:1] and [1:0].
+  if length(gs_reduced) == 1 && length(remaining_vars) == 2
+    mons_list = collect(monomials(pt_reduced))
+    if length(mons_list) == 1 && mons_list[1] == remaining_vars[1] * remaining_vars[2]
+      if gs_reduced[1] == remaining_vars[1] * remaining_vars[2]
+        if reduced_scaling_relations == matrix(ZZ, [[1,1]])
+          intersection_dict[indices] = 2
+          return 2
+        end
+      end
+    end
+  end
+
+  # C.8 Check if this was covered in our special cases
+  if haskey(special_intersection_dict, string([pt_reduced, gs_reduced, remaining_vars, reduced_scaling_relations]))
+    numb = special_intersection_dict[string([pt_reduced, gs_reduced, remaining_vars, reduced_scaling_relations])]
+    intersection_dict[indices] = numb
+    return numb
+  end
+
+  # C.9 In all other cases, proceed via a rationally equivalent cycle
+  println("")
+  println("FOUND CASE THAT CANNOT YET BE DECIDED!")
+  println("$pt_reduced")
+  println("$remaining_vars")
+  println("$gs_reduced")
+  println("$indices")
+  println("$reduced_scaling_relations")
+  println("TRYING WITH EQUIVALENT CYCLE")
+  println("")
+  numb = intersection_from_equivalent_cycle(v, indices, hypersurface_equation, intersection_dict, special_intersection_dict)
+  special_intersection_dict[string([pt_reduced, gs_reduced, remaining_vars, reduced_scaling_relations])] = numb
+  return numb
+
+end
+
+
+
+# ---------------------------------------------------------------------------------------------------------
+# (2) Compute the intersection product from a rationally equivalent cycle.
+# ---------------------------------------------------------------------------------------------------------
+
+function intersection_from_equivalent_cycle(v::NormalToricVariety, indices::NTuple{4, Int64}, hypersurface_equation::MPolyDecRingElem{QQFieldElem, QQMPolyRingElem}, intersection_dict::Dict{NTuple{4, Int64}, ZZRingElem}, special_intersection_dict::Dict{String, ZZRingElem})
+  coeffs_list, tuple_list = Oscar._rationally_equivalent_cycle(v, indices)
+  intersect_numb = 0
+  for k in 1:length(tuple_list)
+    intersect_numb += coeffs_list[k] * sophisticated_intersection_product(v, tuple_list[k], hypersurface_equation, intersection_dict, special_intersection_dict)
+  end
+  @req is_integer(intersect_numb) "Should have expected to find only integer intersection numbers..."
+  intersection_dict[indices] = ZZ(intersect_numb)
+  return ZZ(intersect_numb)
+end
+
+
+
+# ---------------------------------------------------------------------------------------------------------
+# (3) A function to reduce the hypersurface polynomial to {xi = 0} with i in indices
+# ---------------------------------------------------------------------------------------------------------
+
+function _reduce_hypersurface_equation(v::NormalToricVariety, p_hyper::MPolyRingElem, indices::NTuple{4, Int64})
+
+  # Set variables to zero in the hypersurface equation
+  vanishing_vars_pos = unique(indices)  
+  S = cox_ring(v)
+  gS = gens(cox_ring(v))
+  new_p_hyper = divrem(p_hyper, gS[vanishing_vars_pos[1]])[2]
+  for m in 2:length(vanishing_vars_pos)
+    new_p_hyper = divrem(new_p_hyper, gS[vanishing_vars_pos[m]])[2]
+  end
+
+  # Is the resulting polynomial constant?
+  if is_constant(new_p_hyper)
+    return [new_p_hyper, [], [], zero_matrix(ZZ, 0, 0)]
+  end
+
+  # Identify the remaining variables
+  mnf = Oscar._minimal_nonfaces(v)
+  sr_ideal_pos = [Vector{Int}(Polymake.row(mnf, i)) for i in 1:Polymake.nrows(mnf)]
+  remaining_vars_pos = Set(1:length(gS))
+  for my_exps in sr_ideal_pos
+    len_my_exps = length(my_exps)
+    inter_len = count(idx -> idx in vanishing_vars_pos, my_exps)
+    if len_my_exps == inter_len + 1
+      delete!(remaining_vars_pos, my_exps[findfirst(idx -> !(idx in vanishing_vars_pos), my_exps)])
+    end
+  end
+  set_to_one_list = sort([k for k in 1:length(gS) if k ∉ remaining_vars_pos])
+  remaining_vars_pos = setdiff(collect(remaining_vars_pos), vanishing_vars_pos)
+  remaining_vars = [gS[k] for k in remaining_vars_pos]
+
+  # Extract remaining Stanley-Reisner ideal relations
+  sr_reduced = Vector{MPolyDecRingElem{QQFieldElem, QQMPolyRingElem}}()
+  for k in 1:length(sr_ideal_pos)
+    if isdisjoint(set_to_one_list, sr_ideal_pos[k])
+      push!(sr_reduced, prod(gS[m] for m in sr_ideal_pos[k]))
+    end
+  end
+
+  # Identify the remaining scaling relations
+  prepared_scaling_relations = zero_matrix(ZZ, torsion_free_rank(grading_group(S)), length(set_to_one_list) + length(remaining_vars_pos))
+  for k in 1:(length(set_to_one_list) + length(remaining_vars_pos))
+    col = k <= length(set_to_one_list) ? S.d[set_to_one_list[k]].coeff : S.d[remaining_vars_pos[k - length(set_to_one_list)]].coeff
+    for l in 1:length(col)
+      prepared_scaling_relations[l, k] = col[l]
+    end
+  end
+  prepared_scaling_relations = hnf(prepared_scaling_relations)
+  reduced_scaling_relations = prepared_scaling_relations[length(set_to_one_list) + 1: nrows(prepared_scaling_relations), length(set_to_one_list) + 1 : ncols(prepared_scaling_relations)]
+
+  # Identify the final form of the reduced hypersurface equation, by setting all variables to one that we can
+  images = [k in remaining_vars_pos ? gS[k] : one(S) for k in 1:length(gS)]
+  pt_reduced = evaluate(new_p_hyper, images)
+
+  # Return the result
+  return [pt_reduced, sr_reduced, remaining_vars, reduced_scaling_relations]
+end
+
+
+
+# ---------------------------------------------------------------------------------------------------------
+# (4) A function to find a rationally equivalent algebraic cycle.
+# ---------------------------------------------------------------------------------------------------------
+
+function _rationally_equivalent_cycle(v::NormalToricVariety, indices::NTuple{4, Int64})
+
+  # Identify positions of the single and triple variable
+  power_variable = nothing
+  for k in Set(indices)
+    if count(==(k), indices) > 1
+      power_variable = k
+      break
+    end
+  end
+  if power_variable === nothing
+    index = rand(1:length(indices))
+    power_variable = indices[index]
+  end
+  other_variables = [k for k in Set(indices) if k != power_variable]
+  @req length(other_variables) + 1 <= 5 "Found too many variables -- will likely not find a suitable relation!"
+
+  # Let us simplify the problem by extracting the entries in the columns of single_variables and double_variables of the linear relation matrix
+  linear_relations_matrix = matrix(QQ, rays(v))
+  simpler_matrix = linear_relations_matrix[vcat(other_variables, power_variable), :]
+  b = zero_matrix(QQ, length(other_variables) + 1, 1)
+  b[nrows(b), 1] = 1
+  A = solve(simpler_matrix, b; side =:right)
+
+  # Now form the relation in case...
+  employed_relation = -sum((linear_relations_matrix[:, k] .* A[k]) for k in 1:5)
+  employed_relation[power_variable] = 0
+
+  # Generate coefficients and tuples
+  coeffs = Vector{QQFieldElem}()
+  tuples = Vector{NTuple{4, Int64}}()
+  prepared_list = collect(indices)
+  pos_power_variable = findfirst(==(power_variable), prepared_list)
+
+  # Populate `coeffs` and `tuples`
+  for k in 1:length(employed_relation)
+    if employed_relation[k] != 0
+      push!(coeffs, employed_relation[k])
+      new_tuple = copy(prepared_list)
+      new_tuple[pos_power_variable] = k
+      push!(tuples, Tuple(new_tuple))
+    end
+  end
+  return [coeffs, tuples]
+
+end