3 refactor topdata (#8)

---------

Co-authored-by: Julian Piribauer <julian.piribauer@gmail.com>
Reviewed-on: #8
This commit was merged in pull request #8.
This commit is contained in:
2026-07-24 20:01:14 +02:00
co-authored by Julian Piribauer
parent 28aada2bfe
commit 1efe3065ef
2 changed files with 681 additions and 627 deletions
+668 -266
View File
@@ -1,309 +1,704 @@
### CLASS: ToricTopData ####
# Methods:
# - pointstopoly: converts list of points to CY-polytope
# - disc: computes discriminant for polytope
# - topdata: computes all sorts of topological data for polytope
# - CICYtopdata: computes topdata in the simple CICY format,
# e.g. [[3,3]] for two cubics in P5 or [[3,0,1],[0,3,1]] for the Tian--Yau manifold
import numpy as np
# For grep-commands
import logging
import os
import json
from datetime import datetime, timezone
import re
# Logger
logger = logging.getLogger(__name__)
logging.basicConfig(level=logging.DEBUG)
class ToricTopData:
# For grep-commands
find_zs = re.compile(r'z\d+')
find_ls = re.compile(r'l\d+')
find_index_z = re.compile(r'\d+')
findzs = re.compile('z\d+')
findls = re.compile('l\d+')
findindexz = re.compile('\d+')
def pointstopoly(self, points):
# Takes a list of points and returns the points inside its convex hull while omitting points inside faces of co-dimension one
class Polytope:
"""
Class representing a lattice polytope. Objects contain a list of defining points
and a list of points in the convex hull of the polytope. Points inside codimension
one are excluded from the convex hull and stored in a separate list.
"""
#############################
# Initialisation
#############################
def _validate_points_and_get_dimension(self) -> int:
if len(self.defining_points) < 1:
raise ValueError("At least one point is required.")
elif len(set([len(p) for p in self.defining_points])) != 1:
raise ValueError("All points must have the same dimension.")
elif len(self.defining_points) <= len(self.defining_points[0]):
raise ValueError("The number of points must be greater than the dimension of the points.")
return len(self.defining_points[0])
def _points_to_polytope(self) -> tuple:
"""
Takes the defining points of the polytope and returns the points inside its convex hull
while removing and returning points inside faces of co-dimension one
as a second return value.
"""
lattice_polytope = LatticePolytope(self.defining_points)
points_in_convex_hull = [list(m) for m in lattice_polytope.points()]
points_in_codim1_faces = [] # is populated below
for facet in lattice_polytope.facets():
for i in facet.interior_point_indices():
points_inside_facet = [list(m) for m in facet.points(i)][0]
if points_inside_facet in points_in_convex_hull:
points_in_convex_hull.remove(points_inside_facet) # remove those inside codim1 faces
points_in_codim1_faces.append(points_inside_facet) # save omitted points in pcodim1
zeros = [0 for i in range(len(points[0]))] #origin
pc = LatticePolytope(points) # all points
polytope = [list(m) for m in pc.points()]
self.pcodim1 = [] #points in co-dimension one
for f in pc.facets():
for i in f.interior_point_indices():
try:
polytope.remove([list(m) for m in f.points(i)][0]) # remove those inside codim1 faces
self.pcodim1.append([list(m) for m in f.points(i)][0]) # save omitted points in pcodim1
except:
pass
# Moving zeros to the end
polytope.remove(zeros)
polytope.append(zeros)
return polytope
points_in_convex_hull.remove(self.origin)
points_in_convex_hull.append(self.origin)
def disc(self, polytope, only_sc=False, no_triangulation=0):
# Takes a poltyope and returns a list of tuples [disc_i,codim_i] of discriminant factors disc_i coming from a relation inside a face of codimension codim_i
return lattice_polytope, points_in_convex_hull, points_in_codim1_faces
self.pc = LatticePolytope(polytope)
dimp = len(polytope[0])
zeros = [0 for i in range(dimp)]
# Compute all relations in all faces
self.globalks = []
if only_sc:
pcfaces = self.pc.faces()[:-1]
def __init__(self, points):
self.defining_points = points
self.dimension = self._validate_points_and_get_dimension()
self.origin = [0] * self.dimension
self.lattice_polytope, self.relevant_lattice_points, self.points_in_codim1_faces = self._points_to_polytope()
class ToricPolytope(Polytope):
"""
Used to represent the anti-canonical Weyl divisor -K of a mirror family.
By definition, it is also given by the polar dual polytope to the
toric polytope describing the original family's ambient space.
Example: For the mirror quintic, the vertices are e_1, ..., e_4, -e_1-...-e_4.
For CICYs, a nef-partition must be provided, defaulting to hypersurface if none given.
"""
##############################
# Initialisation
##############################
def _set_fine_star_triangulations(self):
"""
Returns a list of fine star triangulations of the polytope.
"""
self.point_configuration = PointConfiguration(self.relevant_lattice_points)
star_triangulations = self.point_configuration.restrict_to_star_triangulations(self.origin)
fine_star_triangulations = star_triangulations.restrict_to_fine_triangulations()
if len(fine_star_triangulations.triangulations_list()) == 0:
raise ValueError("No fine star triangulation found for the given polytope.")
if len(fine_star_triangulations.triangulations_list()) > 1:
logger.warning("More than one fine star triangulation! (%d)", len(fine_star_triangulations.triangulations_list()))
self.triangulations = fine_star_triangulations.triangulations_list()
def set_triangulation(self, no_triangulation):
"""
Sets the triangulation of the polytope based on the provided index,
ensuring that the index is valid.
"""
if no_triangulation < 0 or no_triangulation >= len(self.triangulations):
raise IndexError("Invalid triangulation index {} (available: 0..{}).".format(no_triangulation, len(self.triangulations) - 1))
logger.info("Using triangulation index %d of [0..%d].", no_triangulation, len(self.triangulations) - 1)
self.triangulation = self.triangulations[no_triangulation]
def set_nef_partition(self, nef_partition):
"""
Sets the nef partition for the polytope, ensuring that it is valid.
"""
if nef_partition is not None:
if sorted(flatten(nef_partition)) != list(range(len(self.relevant_lattice_points) - 1)):
raise ValueError("Invalid nef partition: {}".format(nef_partition))
self.nef_partition = nef_partition
else:
pcfaces = self.pc.faces()
self.nef_partition = [list(range(len(self.relevant_lattice_points) - 1))] # Default to a trivial partition if none is provided
for facesofdim in pcfaces: # for each face-dimension
atdim = []
for face in facesofdim: # for each face of fixed dimension
atface = []
facepoints = [list(i) for i in matrix(face.points())]
def set_Mori_cone(self, lvec = None):
"""
Sets the Mori cone of the toric variety associated with the polytope.
"""
default_Mori_cone = self.ToricVariety.Mori_cone()
if lvec is None:
self.Mori_cone = default_Mori_cone
else:
if not isinstance(lvec, (list, tuple, Matrix)):
raise TypeError("lvec must be a list, tuple, or Matrix.")
elif len(set([len(row) for row in lvec])) != 1:
raise ValueError("All rows in lvec must have the same length.")
elif len(lvec[0]) != len(self.relevant_lattice_points):
raise ValueError("The length of the l-vectors is invalid.")
elif matrix(default_Mori_cone.rays()).stack(matrix(lvec)).rank() != matrix(default_Mori_cone.rays()).rank():
raise ValueError("The provided l-vectors are not in the linear span of the original Mori cone.")
else:
self.Mori_cone = Cone(lvec)
if len(matrix(self.Mori_cone.rays()).kernel().gens()) > 0:
logger.warning("Non-simplicial Mori-cone encountered.")
self.Mori_cone = Cone(
transpose(
transpose(self.Mori_cone.rays())
*
matrix(
transpose(matrix(self.Mori_cone.rays()).kernel().gens())
).kernel().gens()
).transpose()
)
logger.info("Using the first %d linearly independent vectors.", len(self.Mori_cone.rays()))
self.lvec = matrix(self.Mori_cone.rays())
def __init__(
self,
points,
no_triangulation = 0,
model_name = None,
nef_partition = None,
lvec = None, # Note that the inner point is the last entry
):
super().__init__(points)
self.model_name = model_name
self._set_fine_star_triangulations()
self.set_triangulation(no_triangulation)
self.set_nef_partition(nef_partition)
self.cy_dimension = self.dimension - len(self.nef_partition)
self.fan = self.triangulation.fan(self.origin)
self.ToricVariety = ToricVariety(self.fan)
self.set_Mori_cone(lvec = lvec)
self.model_name = model_name if model_name else "Unnamed Model"
logger.info(
"Initialised %s: %s%d with h21 = %d.",
self.model_name,
"CY" if len(self.nef_partition) == 1 else "CICY",
self.cy_dimension,
len(list(self.lvec)),
)
###############################
# Discriminant Computation
###############################
def _collect_linear_dependencies(self, only_strong_coupling):
# Compute linear relations in all faces and map them to global coordinates.
all_linear_dependencies_among_points = []
polytope_faces = self.lattice_polytope.faces()
relevant_faces = polytope_faces[:-1] if only_strong_coupling else polytope_faces
for faces_of_given_dimension in relevant_faces:
relations_at_given_dimension = []
for face in faces_of_given_dimension:
points_in_face = [list(point) for point in matrix(face.points())]
if self.origin in points_in_face:
points_in_face.remove(self.origin)
# Keep only points that survived codim-1 face filtering.
points_in_face = [point for point in points_in_face if point in self.relevant_lattice_points]
if len(points_in_face) == 0:
continue
# In principle, the relations are just the kernel of the matrix of points in the face, but
# to express the discriminants in Batyrev coordinates, we need to map them to the global coordinates
# of the polytope.
point_index_dic = {str(point): self.relevant_lattice_points.index(point) for point in points_in_face}
linear_dependencies_in_face = matrix(points_in_face).kernel().gens()
for generator in linear_dependencies_in_face:
# Exclude origin here; its weight is appended later as minus the sum.
generator_as_global_relation = [0] * (len(self.relevant_lattice_points) - 1)
for idx in range(len(generator)):
generator_as_global_relation[point_index_dic[str(points_in_face[idx])]] = generator[idx]
relations_at_given_dimension.append(generator_as_global_relation)
if len(relations_at_given_dimension) == 0:
unique_relations = []
else:
unique_relations = np.unique(np.matrix(relations_at_given_dimension), axis=0).tolist()
all_linear_dependencies_among_points.append(unique_relations)
return all_linear_dependencies_among_points
def _append_origin_weight(self, all_linear_dependencies_among_points):
# Adding the origin weight to each relation.
for relations_at_dim in all_linear_dependencies_among_points:
for relation in relations_at_dim:
if len(relation) > 0:
relation.append(-sum(relation))
def _setup_discriminant_symbols(self):
mori_rays = self.Mori_cone.rays()
a_vars = [var("a_{}".format(u), latex_name="a_{{}}".format(u)) for u in (1..len(mori_rays))]
z_vars = var('z', n=len(mori_rays)+1, latex_name='z') # z[0] is superfluous
lambda_vars = var('l', n=len(mori_rays)+1, latex_name='l') # l[0] is superfluous
a_row = matrix(a_vars)
mori_matrix = matrix(mori_rays)
return a_vars, z_vars, lambda_vars, a_row, mori_matrix
def _build_equation_system(self, relations_at_given_dimension, a_vars, z_vars, lambda_vars, a_row, mori_matrix):
equation_system = []
number_of_points = len(relations_at_given_dimension[0])
# Eq. (8)-style building blocks: linear forms in lambda_a for each point index i.
linear_forms = [
sum([
relations_at_given_dimension[relation_id][point_id] * lambda_vars[relation_id]
for relation_id in range(len(relations_at_given_dimension))
])
for point_id in range(number_of_points)
]
for relation in relations_at_given_dimension:
# This determines the Batyrev coordinate the dependence represents:
solution = solve((a_row * mori_matrix - matrix(relation)).list(), a_row.list())
z_index = [abs(ai.subs(solution)) for ai in a_vars].index(1) + 1
equation_system.append(
z_vars[z_index] - prod([linear_forms[point_id] ** relation[point_id] for point_id in range(number_of_points)])
)
return equation_system
def _eliminate_and_normalize(self, equation_system, lambda_vars):
lambda_names = sorted(list(set(find_ls.findall(str(equation_system)))))
lambda_symbols = [lambda_vars[int(name[1:])] for name in lambda_names]
reverse_solve = False
try:
facepoints.remove(zeros)
polynomials = maxima.eliminate(equation_system, lambda_symbols[1:]).sage()
if polynomials[0] == 0:
polynomials = maxima.eliminate(equation_system, lambda_symbols[:-1]).sage()
reverse_solve = True
except:
# In one-parameter cases there may be nothing to eliminate.
polynomials = equation_system
logger.debug("No elimination needed for the equation system: %s", equation_system)
try:
if len(lambda_symbols) == 0:
raise ValueError("No lambda variables found in elimination system.")
if not reverse_solve:
lambda0 = lambda_symbols[0]
polynomials = [poly / lambda0 ** (poly.degree(lambda0)) for poly in polynomials]
else:
last_lambda = lambda_symbols[-1]
polynomials = [poly / last_lambda ** (poly.degree(last_lambda)) for poly in polynomials]
except:
pass
for p in facepoints.copy():
try:
if p not in polytope:
facepoints.remove(p) # remove points that were omitted for polytope
except:
pass
return polynomials
rels = {str(i):polytope.index(i) for i in facepoints} # dictionary for getting positions in relations correct
kernel = matrix(facepoints).kernel().gens()
for k in kernel: # just formatting
globalk = [0 for i in range(len(polytope)-1)] # MAY HAVE TO BE ADAPTED FOR CASES WITH MORE VERTICES
for i in range(len(k)):
globalk[rels[str(facepoints[i])]] = k[i]
atdim.append(globalk)
self.globalks.append([[i[j] for j in range(len(i))] for i in np.unique(np.matrix(atdim),axis=0)])
for atdim in self.globalks: # inserting weight for innter point
for k in atdim:
if len(k)>0:
k.append(-sum(k))
# Computing discriminants
self.pc = PointConfiguration(polytope)
pc_star = self.pc.restrict_to_star_triangulations(zeros)
pc_star_and_fine=pc_star.restrict_to_fine_triangulations()
if(len(pc_star_and_fine.triangulations_list())==0):
print("No fine star triangulation!")
return 0
if(len(pc_star_and_fine.triangulations_list())>1):
print("More than one fine star triangulation! (",len(pc_star_and_fine.triangulations_list()),")")
triangulation=pc_star_and_fine.triangulations_list()[no_triangulation]
fan=triangulation.fan(zeros)
X=ToricVariety(fan)
ls = [l for l in X.Mori_cone().rays()]
if (len(matrix(ls).kernel().gens())>0):
print("Non-simplicical Mori-cone")
def _append_unique_discriminants(self, discriminants, seen_discriminants, polynomials, index, total_relation_blocks, only_strong_coupling):
for polynomial in polynomials:
if polynomial == 0:
continue
As = [var("a_{}".format(u), latex_name="a_{{}}".format(u)) for u in (1..len(ls))]
z = var('z', n=len(ls)+1, latex_name='z') # z[0] is superfluous
l = var('l', n=len(ls)+1, latex_name='l') # l[0] is superfluous
#self.globalks.reverse()
self.globalks = [k for k in self.globalks if k not in ([[]],)]
discs = []
for i,atdim in enumerate(self.globalks):
solsys = []
for kindex in range(len(atdim)):
sol=solve((matrix(As)*matrix(ls)-matrix(atdim[kindex])).list(),matrix(As).list())
zindex = [abs(i.subs(sol)) for i in As].index(1) +1
solsys.append(z[zindex]-prod([sum([atdim[j][i]*l[j] for j in range(len(atdim))])**(atdim[kindex][i]) for i in range(len(atdim[0]))]))
lambdas = list(set(self.findls.findall(str(solsys))))
lambdas.sort()
try:
pol = maxima.eliminate(solsys,[eval(i) for i in lambdas[1:]]).sage()
if pol[0]==0:
pol = maxima.eliminate(solsys,[eval(i) for i in lambdas[:-1]]).sage()
reversesolve = True
polynomial_key = str(polynomial)
if polynomial_key in seen_discriminants:
continue
seen_discriminants.add(polynomial_key)
if only_strong_coupling:
codim = total_relation_blocks - index
else:
reversesolve = False
except: # this is just for 1-parameter cases where there is nothing to solve
pol = solsys
pass
try:
if not reversesolve:
pol = [i/l0**(i.degree(l0)) for i in pol]
else:
pol = [i/eval(lambdas[-1])**(i.degree(eval(lambdas[-1]))) for i in pol]
except:
pass
for poli in pol:
if (poli not in [a[0] for a in discs] and poli!=0):
if only_sc:
discs.append([poli,len(self.globalks)-i]) # account for offset
else:
discs.append([poli,len(self.globalks)-1-i])
codim = total_relation_blocks - 1 - index
discriminants.append([polynomial, codim])
if only_sc: #insert an empty list for the codimension 0 face discriminant
discstmp = discs
discstmp.reverse()
discstmp.append([])
discstmp.reverse()
discs = discstmp
return discs
def disc(self, only_strong_coupling=False):
"""
Computes the A-discriminant for the toric polytope (following Aspinwall, Plesser, Wang),
which gives the singular loci of the moduli space in Batyrev coordinates. There is an option
to only compute strong coupling discriminant factors, so those coming from dependencies
in one-dimensional faces.
def topdata(self, polytope, nef_partition=0, lvec=0, no_triangulation=0, returnvars=False):
# Computes topological data for hypersurfaces and CICYs in toric ambient spaces.
# For hypersurfaces, nef_partition should remain untouched (=0).
# For CICYs with d polynomials, a nef_partition has to be supplied:
# its format should be a list of d lists as in the examples below
# which is a decomposition of the N points of $(polytope); the number
# should corresponds to an enumeration of the points of $(polytope) with
# the inner point omitted.
Returns:
A list of tuples [disc_i, codim_i] of discriminant factors disc_i
coming from a relation inside a face of codimension codim_i.
"""
# Conditions for possible fibrations are:
# - elliptic: if t_i^n = 0 and t_i^{n-1} != 0
# - K3 (only 3-folds): if c2.t_i = 24
logger.info("Computing discriminant for %s", self.model_name)
if only_strong_coupling:
logger.info("Restricting to strong-coupling loci (codimension-1 faces).")
# input: - array of points in polytope
# - (for CICYs:) nef-partition
# - (optional:) set of l-vectors to be used
# - (optional:) index of triangulation
all_linear_dependencies_among_points = self._collect_linear_dependencies(only_strong_coupling)
self._append_origin_weight(all_linear_dependencies_among_points)
a_vars, z_vars, lambda_vars, a_row, mori_matrix = self._setup_discriminant_symbols()
# output: - Mori-cone generators
# - Intersection rings of CY and ambient space
# - Topological data
# - Possible fibrations
# - instance attributes: triangulation and Mori-cone generators
all_linear_dependencies_among_points = [relations for relations in all_linear_dependencies_among_points if relations]
discriminants = []
seen_discriminants = set()
total_relation_blocks = len(all_linear_dependencies_among_points)
pc = PointConfiguration(polytope)
for index, relations_at_given_dimension in enumerate(all_linear_dependencies_among_points):
if len(relations_at_given_dimension) == 0:
continue
dimp=len(polytope[0]) # dimension of polytope
if( nef_partition == 0):
no_polys = 1
nef_partition = [[i for i in range(len(polytope)-1)]] # trivial partition
else:
if( sorted(flatten(nef_partition)) != [ i for i in range(len(polytope)-1)] ):
print("Nef-partition not valid!")
no_polys = len(nef_partition)
zeros=zero_vector(dimp)
equation_system = self._build_equation_system(
relations_at_given_dimension,
a_vars,
z_vars,
lambda_vars,
a_row,
mori_matrix,
)
polynomials = self._eliminate_and_normalize(equation_system, lambda_vars)
self._append_unique_discriminants(
discriminants,
seen_discriminants,
polynomials,
index,
total_relation_blocks,
only_strong_coupling,
)
pc_star = pc.restrict_to_star_triangulations(zeros)
pc_star_and_fine=pc_star.restrict_to_fine_triangulations()
if(len(pc_star_and_fine.triangulations_list())==0):
print("No fine star triangulation!")
return 0
if(len(pc_star_and_fine.triangulations_list())>1):
print("More than one fine star triangulation! (",len(pc_star_and_fine.triangulations_list()),")")
self.triangulation=pc_star_and_fine.triangulations_list()[no_triangulation]
self.fan=self.triangulation.fan(zeros)
self.X=ToricVariety(self.fan)
if only_strong_coupling:
# Insert an empty codimension-0 entry for compatibility with prior behavior.
discriminants = [[]] + discriminants
# Change l-vectors if given:
if lvec==0:
self.lall=matrix(self.X.Mori_cone().rays())
else:
self.lall=lvec
nodivs = len(polytope)-dimp-1
self.MoriMatrix=matrix([[-sum([self.lall[j][i] for i in part]) for part in nef_partition] for j in range(nodivs)]), self.lall[:,:-1]
# Check whether Mori-cone is simplicial and continue with first $(nodivs) vectors
if ( self.lall.dimensions()[0] > nodivs ):
print("Non-simplicial Kähler-cone! (",self.lall.dimensions()[0]," > ",nodivs,")")
print("Picking ",nodivs," linearly independent vectors.")
l = transpose(transpose(self.lall)*(matrix(transpose(matrix(self.lall.kernel().gens())).kernel().gens()).transpose()))
self.MoriMatrix=matrix([[-sum([l[j][i] for i in part]) for part in nef_partition] for j in range(nodivs)]), l[:,:-1]
else:
l = self.lall
logger.info("Found %d discriminant factors.", len(discriminants))
logger.info("Discriminant factors: %s", discriminants)
HH=self.X.cohomology_ring()
self.D = [HH(self.X.divisor(i)) for i in [0..len(polytope)-2]]
zs = list(set(self.findzs.findall(str(self.D)))) # gives the $(nodivs) z-variables present in divisor classes
# Find Kähler-cone generators as duals to l-vectors:
Bsinv=l.matrix_from_columns([eval(m) for m in sorted(list(self.findindexz.findall(str(zs))))])
if (Bsinv.det()==0):
print("l-vectors not independent in divisor basis. Pick them manually with ``lvec=...''")
return
return discriminants
##############################
# Topological Data Computation
##############################
def _compute_Kahler_generators(self):
"""
Computes the Kähler cone generators J_a dual to the Mori cone generators l^a.
The toric divisors D_i satisfy ∑_i l^a_i D_i = 0 in cohomology, so only a
subset of them are linearly independent. We identify the free generators via
the z-variables appearing in the cohomology ring, extract the corresponding
columns of the Mori matrix l, invert it, and express J = Bs * D_basis.
"""
z_names = list(set(find_zs.findall(str(self.D))))
z_indices = sorted([int(find_index_z.search(name).group()) for name in z_names])
Bsinv = self.lvec.matrix_from_columns(z_indices)
if Bsinv.det() == 0:
raise ValueError("l-vectors are not independent in the divisor basis. Supply them manually via lvec=...")
Bs = transpose(Bsinv**(-1))
self.J = Bs*vector([self.D[i] for i in list(set([eval(i) for i in self.findindexz.findall(str(self.findzs.findall(str(self.D))))]))])
# Hyperplane divisor class:
self.Ydualform = product([sum([HH(self.D[p]) for p in part]) for part in nef_partition])
return Bs * vector([self.D[i] for i in z_indices])
# Tuple lists for computations below
TupleListCY = UnorderedTuples(range(nodivs),dimp-no_polys)
TupleListCYordered = Tuples(range(nodivs),dimp-no_polys) # for intersection ring multiplicities must be included
TupleListAm = UnorderedTuples(range(nodivs),dimp)
def _compute_intersection_ring(self):
"""
Builds the intersection ring polynomial by adding up the intersection numbers
in self.intersection_numbers_CY.
"""
# Full ring: ordered tuples, so t_i*t_j and t_j*t_i both contribute.
intring = sum(
product([self.t[i] for i in tl]) * self.intersection_numbers_CY[tuple(sorted(tl))]
for tl in Tuples(range(self.no_divs), self.cy_dimension)
)
# Ring without multiplicities: each distinct monomial once.
intringnomults = sum(
product([self.t[i] for i in tl]) * val
for tl, val in self.intersection_numbers_CY.items()
)
return intring, intringnomults
# Find intersection numbers on CY and on ambient space:
intCY = [list((tl,self.X.integrate(self.Ydualform*product([self.J[tl[j]] for j in range(dimp-no_polys)])))) for tl in TupleListCY]
intAm = [list((tl,self.X.integrate(product([self.J[tl[j]] for j in range(dimp)])))) for tl in TupleListAm]
def _intersection_substitution_rules(self):
"""
Builds monomial substitution rules from intersection numbers, e.g.
t0^3 -> kappa_(0,0,0), t0*t1^2 -> kappa_(0,1,1).
"""
rules = {}
for idx_tuple, value in self.intersection_numbers_CY.items():
monomial = product([self.t[i] for i in idx_tuple]).expand()
rules[monomial] = value
return rules
# Form intersection ring on CY:
t = var('t', n=nodivs, latex_name='t')
intring = sum([product([t[index] for index in tl])*self.X.integrate(self.Ydualform*product([self.J[tl[j]] for j in range(dimp-no_polys)])) for tl in TupleListCYordered])
intringnomults = sum([product([t[index] for index in tl])*self.X.integrate(self.Ydualform*product([self.J[tl[j]] for j in range(dimp-no_polys)])) for tl in set([tuple(sorted(i)) for i in TupleListCYordered])])
# Compute Chern-character for topological data:
var(self.findzs.findall(str([self.J[i] for i in range(nodivs)]))) # Introduces all z in Kähler-forms as variables...
zs = [eval(z) for z in self.findzs.findall(str([self.J[i] for i in range(nodivs)]))] # ... and puts them into a vector.
tsubs=solve([lift(self.J[i])==t[i] for i in range(nodivs)],zs) # Finds substitution rule for zs in terms Kähler-cone generators
cc = (product([1+eval(str((lift(d)))).subs(tsubs[0]) for d in self.D]))/(product([1+sum([eval(str(lift(HH(self.D[p])))) for p in part]) for part in nef_partition]).subs(tsubs[0])) # Adjunction-formula
def _replace_intersection_monomials(self, expression):
"""
Replaces top-degree Kähler monomials in an expression by CY intersection numbers.
"""
expanded_expression = expression.expand()
replaced_expression = expanded_expression
print('--- Toric divisors (ambient space) -----------')
print(self.D)
print('\n--- Kähler cone generators (ambient space) --- ')
print(self.J)
print('\n--- Mori-cone-generators (ambient space) ------')
for j in range(nodivs):
print(self.MoriMatrix[0][j],self.MoriMatrix[1][j])
# Replace each degree-n monomial by its corresponding intersection number,
# preserving the polynomial coefficient in front of that monomial.
for idx_tuple, value in self.intersection_numbers_CY.items():
monomial = product([self.t[i] for i in idx_tuple]).expand()
coeff = expanded_expression.coefficient(monomial)
if coeff != 0:
replaced_expression -= coeff * monomial
replaced_expression += coeff * value
print('\n--- Intersection on CY ------------------------')
Ccy=[product([t[iden[0][i]] for i in range(dimp-no_polys)])==iden[1] for iden in intCY]
print(Ccy)
print('\nR = ',intring)
print('\nR (no multiplicities) = ',intringnomults)
return replaced_expression.expand()
print('\n--- Intersection in ambient space -------------')
Cam=[product([t[iden[0][i]] for i in range(dimp)])==iden[1] for iden in intAm]
print(Cam)
def _compute_Chern_character(self):
"""
Computes the total Chern class of the CY via the adjunction formula:
print('\n--- Topological data --------------------------')
#print((product([1+lift(d) for d in self.D]))/(product([1+sum([lift(HH(self.D[p])) for p in part]) for part in nef_partition])))
chi = self.X.integrate(self.Ydualform*(product([1+lift(d) for d in self.D]))/(product([1+sum([lift(HH(self.D[p])) for p in part]) for part in nef_partition])))
print('chi = ',chi)
c(Y) = c(TX) / c(N_{Y/X})
where c(TX) = ∏_i (1 + D_i) and c(N_{Y/X}) = ∏_k (1 + Y_k*) with Y_k*
the nef partition divisor classes. Everything is expressed in the Kähler
basis by solving lift(J_i) = t_i for the cohomology ring z-variables.
"""
z_names = find_zs.findall(str([self.J[i] for i in range(self.no_divs)]))
var(z_names)
zs = [eval(z) for z in z_names]
tsubs = solve([lift(self.J[i]) == self.t[i] for i in range(self.no_divs)], zs)
numerator = product([1 + eval(str(lift(d))).subs(tsubs[0]) for d in self.D])
denominator = product([
1 + sum([eval(str(lift(self.HH(self.D[p])))) for p in part])
for part in self.nef_partition
]).subs(tsubs[0])
return numerator / denominator
def _compute_Chern_polynomials(self, Chern_character):
"""
Computes the Chern polynomials of the CY by expanding the total Chern character
and substituting intersection numbers.
Returns both a Chern polynomial for each degree restricted to suitable intersections
of Kähler divisors and their integrated versions.
"""
c = var('c')
if( round((cc).subs({t:c*t for t in t}).taylor(c,0,1).coefficient(c,1).subs({t[i]:1 for i in range(nodivs)}),8) != 0):
return "First Chern-class not zero!"
# Print Chern classes:
cJ = [[0]]*(dimp-no_polys)
for i in range(2,dimp-no_polys):
print('\nc',i,' = ',(cc).subs({t:c*t for t in t}).taylor(c,0,i).coefficient(c,i))
TupleListi = UnorderedTuples(range(nodivs),i)
TupleListz = UnorderedTuples(range(nodivs),dimp-no_polys-i)
cJ[i] = [sum([((cc).subs({t:c*t for t in t}).taylor(c,0,i).coefficient(c,i)).coefficient(product([t[i] for i in tupl]))*((product([t[j] for j in tuplz])*product([t[i] for i in tupl])).subs(Ccy)) for tupl in TupleListi]) for tuplz in TupleListz]
print('integrated: '+str([cJ[i][j]*product([t[k] for k in TupleListz[j]]) for j in range(len(TupleListz))]))
print('\nc',dimp-no_polys,' = ',(cc).subs({t:c*t for t in t}).taylor(c,0,dimp-no_polys).coefficient(c,dimp-no_polys))
Chern_character_taylored = (Chern_character).subs({t: c * t for t in self.t}).taylor(c, 0, self.cy_dimension)
msg = ""
# elliptic:
## Indicates whether J_i^n==0 with J_i^(n-1)!=0 for n-folds and some J_i
for i in range(nodivs):
for intn in intCY:
if( intn[0] == [i]*(dimp-no_polys) and intn[1] == 0):
if not(all([round((t[i]**(dimp-no_polys-1)*t[j]).subs(Ccy),10) == 0 for j in range(nodivs) if j != i]) ):
msg += "Possible elliptic fibration in cycle dual to t"+str(i)+".\n"
Chern_polynomials = [
[
(prod(self.t[i] for i in tple) * Chern_character_taylored.coefficient(c, k)).expand()
for tple in UnorderedTuples(range(self.no_divs), self.cy_dimension - k)
]
for k in range(1, self.cy_dimension + 1) # for each degree
]
# K3 for 3-folds
## Indicates whether c2.J_i =24 for some J_i
if( dimp+1-no_polys == 4):
TupleList2 = UnorderedTuples(range(nodivs),2)
for j in range(nodivs):
if( sum([((cc).subs({t:c*t for t in t}).taylor(c,0,2).coefficient(c,2)).coefficient(product([t[i] for i in tupl]))*((product([t[i] for i in tupl])*t[j]).subs(Ccy)) for tupl in TupleList2]) == 24 ):
msg += "Possible K3-fibration in divisor t"+str(j)+".\n"
integrated_Chern_classes = [
[
self._replace_intersection_monomials(poly) for poly in polys
]
for polys in Chern_polynomials
]
if( msg != "" ):
print('\n--- Fibrations --------------------------------')
print(msg)
if returnvars:
return str([(self.MoriMatrix[0][j],self.MoriMatrix[1][j]) for j in range(nodivs)]).replace("(","{").replace(")","}").replace("[","{").replace("]","}") # might return more in the future if necessary
return Chern_polynomials, integrated_Chern_classes
def _check_for_fibrations(self):
"""
Checks for possible elliptic fibrations via the condition
elliptic in direction i <=> J_i^n == 0 but J_i^(n-1)*J_j != 0 for some j != i.
"""
messages = []
def CICYtopdata(self, CICY, justdata=False, justmori=False, lvec=0):
# Input: a list l with entries l[i,j] that give the weight in the
# ambient projective space i of polynomial j.
# E.g. (P^3| 3 1)
# (P^2| 2 0)
# corresponds to the list ((3,1),(2,0))
for i in range(self.no_divs):
self_intersection = self.intersection_numbers_CY[tuple([i] * self.cy_dimension)]
if self_intersection != 0:
continue
has_nonzero_neighbor = any(
self.intersection_numbers_CY[tuple(sorted([i] * (self.cy_dimension - 1) + [j]))] != 0
for j in range(self.no_divs) if j != i
)
if has_nonzero_neighbor:
messages.append("Possible elliptic fibration in the divisor class dual to t{}.".format(i))
if messages:
logger.info("\n\n--- Elliptic fibrations --------------------------\n%s", "\n".join(messages))
return messages
def _get_polytope_and_glsm_table(self):
"""
Builds a list of lines showing the polytope points as columns of a
(dimension x n) matrix (one column per point), with the GLSM charge
l-vectors listed directly underneath: one row per Mori cone
generator, with each l-vector entry aligned below the point
coordinate it belongs to. Returned as a list of lines (rather than a
single newline-joined string) so that each line becomes its own JSON
array entry and stays on its own line when written to file.
To-do: add inner point and CICY functionality.
"""
points = list(self.relevant_lattice_points) # create a copy
points.remove(self.origin) # Remove the origin point
n_outer_points = len(points)
nef_unit_vectors = [tuple([1 if i == j else 0 for j in range(len(self.nef_partition))]) for i in range(len(self.nef_partition))]
for i, partition in enumerate(self.nef_partition):
for p in partition:
points[p] = tuple(list(nef_unit_vectors[i]) + list(points[p])) # Append nef partition index to each point
points.append(tuple(list(nef_unit_vectors[i]) + [0] * self.dimension)) # Add a point for the nef partition itself
point_rows = [[str(point[row]) for point in points] for row in range(self.dimension + len(self.nef_partition))]
lvec_rows = self._get_lvec_rows()
lvec_rows_string = [[str(entry) for entry in row] for row in lvec_rows] # Exclude the last entry of each l-vector (origin point)
column_widths = [
max(len(row[col]) for row in point_rows + lvec_rows_string)
for col in range(n_outer_points + len(self.nef_partition))
]
def format_row(row_entries, separator):
# Double the separator right before the trailing nef-partition-marker columns.
row_string = row_entries[0].rjust(column_widths[0])
for col in range(1, len(row_entries)):
sep = separator * 2 if col == n_outer_points else separator
row_string += sep + row_entries[col].rjust(column_widths[col])
return row_string
point_lines = [format_row(row, "|") for row in point_rows]
lvec_lines = [format_row(row, ",") for row in lvec_rows_string]
lvec_lines = [line.replace(',,', '; ') for line in lvec_lines]
# Blank line after the leading nef-partition-indicator rows.
num_nef_rows = len(self.nef_partition)
point_lines = point_lines[:num_nef_rows] + [""] + point_lines[num_nef_rows:]
return point_lines + ["-" * len(point_lines[-1])] + lvec_lines
def _get_lvec_rows(self):
"""
Returns the GLSM charge l-vectors (one row per Mori cone generator), with the
nef-partition weight(s) appended, as a plain list of lists of ints - accessible
as lvec[0], lvec[1], ... once loaded from the JSON output.
"""
lvec_rows = [list(row[:-1]) for row in self.lvec.rows()] # strip the total weight entry
for row in lvec_rows:
for partition in self.nef_partition:
row.append(-sum([row[j] for j in partition])) # Append the weight of the NEF partition
return lvec_rows
def _get_output(self):
"""
Returns a dictionary containing all relevant topological data of the CY.
"""
output = {
"polytope_and_glsm_table": self._get_polytope_and_glsm_table(),
"lvec": self._get_lvec_rows(),
"kahler_cone_generators": [str(j) for j in self.J],
"intersection_numbers_CY": dict(self.intersection_numbers_CY),
"intersection_numbers_ambient_space": dict(self.intersection_numbers_ambient_space),
"intersection_ring": self.intersection_ring,
"intersection_ring_no_multiplicities": self.intersection_ring_no_multiplicities,
"chern_polynomials": self.Chern_polynomials,
"integrated_Chern_classes": self.integrated_Chern_classes,
}
debug_string = "\n\n--- Kähler cone generators (ambient space) ---\n" + str(output["kahler_cone_generators"]) + \
"\n\n--- Intersection numbers ambient space -------\n" + str(output["intersection_numbers_ambient_space"])
info_string = "\n\n--- Polytope and GLSM table ------------------\n" + "\n".join(output["polytope_and_glsm_table"]) + \
"\n\n--- L-vectors (GLSM charges) -----------------\n" + str(output["lvec"]) + \
"\n\n--- Intersection numbers CY ------------------\n" + str(output["intersection_numbers_CY"]) + \
"\n\n--- Intersection ring ------------------------\n" + str(output["intersection_ring"]) + \
"\n\n--- Intersection ring no multiplicities ------\n" + str(output["intersection_ring_no_multiplicities"]) + \
"\n\n--- Chern polynomials ------------------------\n" + str(output["chern_polynomials"]) + \
"\n\n--- Integrated Chern classes -----------------\n" + str(output["integrated_Chern_classes"])
logger.info(info_string)
logger.debug(debug_string)
return output
@staticmethod
def _to_json_safe(value):
"""
Recursively converts a value into JSON-safe primitives: dicts/lists are
recursed into (tuple dict keys become strings), numeric Sage objects
(Integer, Rational, ...) become plain Python int/float, and anything else
(e.g. symbolic expressions) is stringified.
"""
if isinstance(value, dict):
return {str(k): ToricPolytope._to_json_safe(v) for k, v in value.items()}
if isinstance(value, (list, tuple)):
return [ToricPolytope._to_json_safe(v) for v in value]
if value is None or isinstance(value, (bool, str, int, float)):
return value
try:
as_int = int(value)
if as_int == value:
return as_int
except (TypeError, ValueError):
pass
try:
return float(value)
except (TypeError, ValueError):
pass
return str(value)
def _write_output_json(self, output):
"""
Writes the given output dictionary (from _get_output, i.e. exactly what is
logged) to data/<model_name>.json."""
if not self.model_name:
logger.info("No model_name given, skipping writing topological data to JSON.")
return None
repo_root = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
out_dir = os.path.join(repo_root, "data/topdata")
os.makedirs(out_dir, exist_ok=True)
out_path = os.path.join(out_dir, "{}.json".format(self.model_name))
with open(out_path, "w") as f:
json.dump(self._to_json_safe(output), f, indent=2)
logger.info("Wrote topological data to %s", out_path)
return out_path
def topdata(self):
"""
Computes topological data for hypersurfaces and CICYs in toric ambient spaces
and saves the result to /data/topdata/<model_name>.json if a model_name is provided.
"""
# Each point gives a toric divisor, each "column" gives one linear relation
self.no_divs = len(self.relevant_lattice_points) - self.dimension - 1
# Cohomology ring and toric divisors D_0, ..., D_{N-1}.
self.HH = self.ToricVariety.cohomology_ring()
self.D = [self.HH(self.ToricVariety.divisor(i)) for i in range(len(self.relevant_lattice_points) - 1)]
# Kähler cone generators J_a (dual to Mori cone generators l^a).
self.J = self._compute_Kahler_generators()
# Poincaré dual of the CY: product over nef partition parts.
self.Ydualform = product([sum([self.D[p] for p in part]) for part in self.nef_partition])
# Kähler parameter symbols t_0, ..., t_{nodivs-1}.
self.t = var('t', n=self.no_divs, latex_name='t')
# Intersection numbers κ_{i_1...i_n} = ∫_Y J_{i_1}∧...∧J_{i_n}.
self.intersection_numbers_CY = {
tuple(tl): self.ToricVariety.integrate(self.Ydualform * product([self.J[tl[j]] for j in range(self.cy_dimension)]))
for tl in UnorderedTuples(range(self.no_divs), self.cy_dimension) # Unordered means only increasing tuples
}
self.intersection_numbers_ambient_space = {
tuple(tl): self.ToricVariety.integrate(product([self.J[tl[j]] for j in range(self.dimension)]))
for tl in UnorderedTuples(range(self.no_divs), self.dimension)
}
# Intersection ring polynomial and version without multiplicities.
self.intersection_ring, self.intersection_ring_no_multiplicities = self._compute_intersection_ring()
Chern_character = self._compute_Chern_character()
self.Chern_polynomials, self.integrated_Chern_classes = self._compute_Chern_polynomials(Chern_character)
output = self._get_output()
self._write_output_json(output)
self._check_for_fibrations()
class ToricPolytopeCICY(ToricPolytope):
"""
Computes topological data for CICYs in toric ambient spaces.
"""
def _CICY_to_points(self, CICY):
"""
Converts a CICY configuration matrix to a list of points in the ambient space.
Each row of the CICY corresponds to a projective space, and each column corresponds to a polynomial.
The entries give the weights of the polynomials in the respective projective spaces.
"""
ambients=(np.array(CICY).transpose()).sum(axis=0)-1
if ambients in ZZ:
ambients=[ambients]
@@ -332,9 +727,16 @@ class ToricTopData:
counters[j]+=len(part[j])
partition.append(flatten(part))
if justdata:
return [list(polytope),partition]
elif justmori:
moricone((list(polytope)))
else:
self.topdata(list(polytope),nef_partition=partition,lvec=lvec)
return list(polytope), partition
def __init__(self, CICY, no_triangulation=0, model_name=None, lvec=None):
points, nef_partition = self._CICY_to_points(CICY)
super().__init__(
points,
no_triangulation=no_triangulation,
model_name=model_name,
nef_partition=nef_partition,
lvec=lvec
)
-348
View File
@@ -1,348 +0,0 @@
# This file was *autogenerated* from the file sage/toric_topdata.sage
from sage.all_cmdline import * # import sage library
_sage_const_0 = Integer(0); _sage_const_1 = Integer(1); _sage_const_2 = Integer(2); _sage_const_8 = Integer(8); _sage_const_10 = Integer(10); _sage_const_4 = Integer(4); _sage_const_24 = Integer(24)### FUNCTIONS: ####
# - pointstopoly: converts list of points to CY-polytope
# - disc: computes discriminant for polytope
# - topdata: computes all sorts of topological data for polytope
# - CICYtopdata: computes topdata in the simple CICY format, e.g. [[3,3]] for two cubics in P5 or [[3,0,1],[0,3,1]] for the Tian--Yau manifold
import numpy as np
# For grep-commands
import re
findzs = re.compile('z\d+')
findls = re.compile('l\d+')
findindexz = re.compile('\d+')
def pointstopoly(points):
# Takes a list of points and returns the points inside its convex hull while omitting points inside faces of co-dimension one
global pcodim1 #points in co-dimension one
zeros = [_sage_const_0 for i in range(len(points[_sage_const_0 ]))] #origin
pc = LatticePolytope(points) # all points
polytope = [list(m) for m in pc.points()]
pcodim1 = []
for f in pc.facets():
for i in f.interior_point_indices():
try:
polytope.remove([list(m) for m in f.points(i)][_sage_const_0 ]) # remove those inside codim1 faces
pcodim1.append([list(m) for m in f.points(i)][_sage_const_0 ]) # save omitted points in pcodim1
except:
pass
# Moving zeros to the end
polytope.remove(zeros)
polytope.append(zeros)
return polytope
def disc(polytope, only_sc=False, no_triangulation=_sage_const_0 ):
# Takes a poltyope and returns a list of tuples [disc_i,codim_i] of discriminant factors disc_i coming from a relation inside a face of codimension codim_i
global globalks, pc
pc = LatticePolytope(polytope)
dimp = len(polytope[_sage_const_0 ])
zeros = [_sage_const_0 for i in range(dimp)]
# Compute all relations in all faces
globalks = []
if only_sc:
pcfaces = pc.faces()[:-_sage_const_1 ]
else:
pcfaces = pc.faces()
for facesofdim in pcfaces: # for each face-dimension
atdim = []
for face in facesofdim: # for each face of fixed dimension
atface = []
facepoints = [list(i) for i in matrix(face.points())]
try:
facepoints.remove(zeros)
except:
pass
for p in facepoints.copy():
try:
if p not in polytope:
facepoints.remove(p) # remove points that were omitted for polytope
except:
pass
rels = {str(i):polytope.index(i) for i in facepoints} # dictionary for getting positions in relations correct
kernel = matrix(facepoints).kernel().gens()
for k in kernel: # just formatting
globalk = [_sage_const_0 for i in range(len(polytope)-_sage_const_1 )] # MAY HAVE TO BE ADAPTED FOR CASES WITH MORE VERTICES
for i in range(len(k)):
globalk[rels[str(facepoints[i])]] = k[i]
atdim.append(globalk)
globalks.append([[i[j] for j in range(len(i))] for i in np.unique(np.matrix(atdim),axis=_sage_const_0 )])
for atdim in globalks: # inserting weight for innter point
for k in atdim:
if len(k)>_sage_const_0 :
k.append(-sum(k))
# Computing discriminants
pc = PointConfiguration(polytope)
pc_star = pc.restrict_to_star_triangulations(zeros)
pc_star_and_fine=pc_star.restrict_to_fine_triangulations()
if(len(pc_star_and_fine.triangulations_list())==_sage_const_0 ):
print("No fine star triangulation!")
return _sage_const_0
if(len(pc_star_and_fine.triangulations_list())>_sage_const_1 ):
print("More than one fine star triangulation! (",len(pc_star_and_fine.triangulations_list()),")")
triangulation=pc_star_and_fine.triangulations_list()[no_triangulation]
fan=triangulation.fan(zeros)
X=ToricVariety(fan)
ls = [l for l in X.Mori_cone().rays()]
if (len(matrix(ls).kernel().gens())>_sage_const_0 ):
print("Non-simplicical Mori-cone")
As = [var("a_{}".format(u), latex_name="a_{{}}".format(u)) for u in (ellipsis_iter(_sage_const_1 ,Ellipsis,len(ls)))]
z = var('z', n=len(ls)+_sage_const_1 , latex_name='z') # z[0] is superfluous
l = var('l', n=len(ls)+_sage_const_1 , latex_name='l') # l[0] is superfluous
#globalks.reverse()
globalks = [k for k in globalks if k not in ([[]],)]
discs = []
for i,atdim in enumerate(globalks):
solsys = []
for kindex in range(len(atdim)):
sol=solve((matrix(As)*matrix(ls)-matrix(atdim[kindex])).list(),matrix(As).list())
zindex = [abs(i.subs(sol)) for i in As].index(_sage_const_1 ) +_sage_const_1
solsys.append(z[zindex]-prod([sum([atdim[j][i]*l[j] for j in range(len(atdim))])**(atdim[kindex][i]) for i in range(len(atdim[_sage_const_0 ]))]))
lambdas = list(set(findls.findall(str(solsys))))
lambdas.sort()
try:
pol = maxima.eliminate(solsys,[eval(i) for i in lambdas[_sage_const_1 :]]).sage()
if pol[_sage_const_0 ]==_sage_const_0 :
pol = maxima.eliminate(solsys,[eval(i) for i in lambdas[:-_sage_const_1 ]]).sage()
reversesolve = True
else:
reversesolve = False
except: # this is just for 1-parameter cases where there is nothing to solve
pol = solsys
pass
try:
if not reversesolve:
pol = [i/l0**(i.degree(l0)) for i in pol]
else:
pol = [i/eval(lambdas[-_sage_const_1 ])**(i.degree(eval(lambdas[-_sage_const_1 ]))) for i in pol]
except:
pass
for poli in pol:
if (poli not in [a[_sage_const_0 ] for a in discs] and poli!=_sage_const_0 ):
if only_sc:
discs.append([poli,len(globalks)-i]) # account for offset
else:
discs.append([poli,len(globalks)-_sage_const_1 -i])
if only_sc: #insert an empty list for the codimension 0 face discriminant
discstmp = discs
discstmp.reverse()
discstmp.append([])
discstmp.reverse()
discs = discstmp
return discs
def topdata(polytope, nef_partition=_sage_const_0 , lvec=_sage_const_0 , no_triangulation=_sage_const_0 , returnvars=False):
# Computes topological data for hypersurfaces and CICYs in toric ambient spaces.
# For hypersurfaces, nef_partition should remain untouched (=0).
# For CICYs with d polynomials, a nef_partition has to be supplied:
# its format should be a list of d lists as in the examples below
# which is a decomposition of the N points of $(polytope); the number
# should corresponds to an enumeration of the points of $(polytope) with
# the inner point omitted.
# Conditions for possible fibrations are:
# - elliptic: if t_i^n = 0 and t_i^{n-1} != 0
# - K3 (only 3-folds): if c2.t_i = 24
# input: - array of points in polytope
# - (for CICYs:) nef-partition
# - (optional:) set of l-vectors to be used
# - (optional:) index of triangulation
# output: - Mori-cone generators
# - Intersection rings of CY and ambient space
# - Topological data
# - Possible fibrations
# - GLOBALS: triangulation and Mori-cone generators
global triangulation, lall, MoriMatrix, fan, X, Ydualform, J, D
pc = PointConfiguration(polytope)
dimp=len(polytope[_sage_const_0 ]) # dimension of polytope
if( nef_partition == _sage_const_0 ):
no_polys = _sage_const_1
nef_partition = [[i for i in range(len(polytope)-_sage_const_1 )]] # trivial partition
else:
if( sorted(flatten(nef_partition)) != [ i for i in range(len(polytope)-_sage_const_1 )] ):
print("Nef-partition not valid!")
no_polys = len(nef_partition)
zeros=zero_vector(dimp)
pc_star = pc.restrict_to_star_triangulations(zeros)
pc_star_and_fine=pc_star.restrict_to_fine_triangulations()
if(len(pc_star_and_fine.triangulations_list())==_sage_const_0 ):
print("No fine star triangulation!")
return _sage_const_0
if(len(pc_star_and_fine.triangulations_list())>_sage_const_1 ):
print("More than one fine star triangulation! (",len(pc_star_and_fine.triangulations_list()),")")
triangulation=pc_star_and_fine.triangulations_list()[no_triangulation]
fan=triangulation.fan(zeros)
X=ToricVariety(fan)
# Change l-vectors if given:
if lvec==_sage_const_0 :
lall=matrix(X.Mori_cone().rays())
else:
lall=lvec
nodivs = len(polytope)-dimp-_sage_const_1
MoriMatrix=matrix([[-sum([lall[j][i] for i in part]) for part in nef_partition] for j in range(nodivs)]), lall[:,:-_sage_const_1 ]
# Check whether Mori-cone is simplicial and continue with first $(nodivs) vectors
if ( lall.dimensions()[_sage_const_0 ] > nodivs ):
print("Non-simplicial Kähler-cone! (",lall.dimensions()[_sage_const_0 ]," > ",nodivs,")")
print("Picking ",nodivs," linearly independent vectors.")
l = transpose(transpose(lall)*(matrix(transpose(matrix(lall.kernel().gens())).kernel().gens()).transpose()))
MoriMatrix=matrix([[-sum([l[j][i] for i in part]) for part in nef_partition] for j in range(nodivs)]), l[:,:-_sage_const_1 ]
else:
l = lall
HH=X.cohomology_ring()
D = [HH(X.divisor(i)) for i in (ellipsis_range(_sage_const_0 ,Ellipsis,len(polytope)-_sage_const_2 ))]
zs = list(set(findzs.findall(str(D)))) # gives the $(nodivs) z-variables present in divisor classes
# Find Kähler-cone generators as duals to l-vectors:
Bsinv=l.matrix_from_columns([eval(m) for m in sorted(list(findindexz.findall(str(zs))))])
if (Bsinv.det()==_sage_const_0 ):
print("l-vectors not independent in divisor basis. Pick them manually with ``lvec=...''")
return
Bs = transpose(Bsinv**(-_sage_const_1 ))
J = Bs*vector([D[i] for i in list(set([eval(i) for i in findindexz.findall(str(findzs.findall(str(D))))]))])
# Hyperplane divisor class:
Ydualform = product([sum([HH(D[p]) for p in part]) for part in nef_partition])
# Tuple lists for computations below
TupleListCY = UnorderedTuples(range(nodivs),dimp-no_polys)
TupleListCYordered = Tuples(range(nodivs),dimp-no_polys) # for intersection ring multiplicities must be included
TupleListAm = UnorderedTuples(range(nodivs),dimp)
# Find intersection numbers on CY and on ambient space:
intCY = [list((tl,X.integrate(Ydualform*product([J[tl[j]] for j in range(dimp-no_polys)])))) for tl in TupleListCY]
intAm = [list((tl,X.integrate(product([J[tl[j]] for j in range(dimp)])))) for tl in TupleListAm]
# Form intersection ring on CY:
t = var('t', n=nodivs, latex_name='t')
intring = sum([product([t[index] for index in tl])*X.integrate(Ydualform*product([J[tl[j]] for j in range(dimp-no_polys)])) for tl in TupleListCYordered])
intringnomults = sum([product([t[index] for index in tl])*X.integrate(Ydualform*product([J[tl[j]] for j in range(dimp-no_polys)])) for tl in set([tuple(sorted(i)) for i in TupleListCYordered])])
# Compute Chern-character for topological data:
var(findzs.findall(str([J[i] for i in range(nodivs)]))) # Introduces all z in Kähler-forms as variables...
zs = [eval(z) for z in findzs.findall(str([J[i] for i in range(nodivs)]))] # ... and puts them into a vector.
tsubs=solve([lift(J[i])==t[i] for i in range(nodivs)],zs) # Finds substitution rule for zs in terms Kähler-cone generators
cc = (product([_sage_const_1 +eval(str((lift(d)))).subs(tsubs[_sage_const_0 ]) for d in D]))/(product([_sage_const_1 +sum([eval(str(lift(HH(D[p])))) for p in part]) for part in nef_partition]).subs(tsubs[_sage_const_0 ])) # Adjunction-formula
print('--- Toric divisors (ambient space) -----------')
print(D)
print('\n--- Kähler cone generators (ambient space) --- ')
print(J)
print('\n--- Mori-cone-generators (ambient space) ------')
for j in range(nodivs):
print(MoriMatrix[_sage_const_0 ][j],MoriMatrix[_sage_const_1 ][j])
print('\n--- Intersection on CY ------------------------')
Ccy=[product([t[iden[_sage_const_0 ][i]] for i in range(dimp-no_polys)])==iden[_sage_const_1 ] for iden in intCY]
print(Ccy)
print('\nR = ',intring)
print('\nR (no multiplicities) = ',intringnomults)
print('\n--- Intersection in ambient space -------------')
Cam=[product([t[iden[_sage_const_0 ][i]] for i in range(dimp)])==iden[_sage_const_1 ] for iden in intAm]
print(Cam)
print('\n--- Topological data --------------------------')
#print((product([1+lift(d) for d in D]))/(product([1+sum([lift(HH(D[p])) for p in part]) for part in nef_partition])))
chi = X.integrate(Ydualform*(product([_sage_const_1 +lift(d) for d in D]))/(product([_sage_const_1 +sum([lift(HH(D[p])) for p in part]) for part in nef_partition])))
print('chi = ',chi)
c = var('c')
if( round((cc).subs({t:c*t for t in t}).taylor(c,_sage_const_0 ,_sage_const_1 ).coefficient(c,_sage_const_1 ).subs({t[i]:_sage_const_1 for i in range(nodivs)}),_sage_const_8 ) != _sage_const_0 ):
return "First Chern-class not zero!"
# Print Chern classes:
cJ = [[_sage_const_0 ]]*(dimp-no_polys)
for i in range(_sage_const_2 ,dimp-no_polys):
print('\nc',i,' = ',(cc).subs({t:c*t for t in t}).taylor(c,_sage_const_0 ,i).coefficient(c,i))
TupleListi = UnorderedTuples(range(nodivs),i)
TupleListz = UnorderedTuples(range(nodivs),dimp-no_polys-i)
cJ[i] = [sum([((cc).subs({t:c*t for t in t}).taylor(c,_sage_const_0 ,i).coefficient(c,i)).coefficient(product([t[i] for i in tupl]))*((product([t[j] for j in tuplz])*product([t[i] for i in tupl])).subs(Ccy)) for tupl in TupleListi]) for tuplz in TupleListz]
print('integrated: '+str([cJ[i][j]*product([t[k] for k in TupleListz[j]]) for j in range(len(TupleListz))]))
print('\nc',dimp-no_polys,' = ',(cc).subs({t:c*t for t in t}).taylor(c,_sage_const_0 ,dimp-no_polys).coefficient(c,dimp-no_polys))
msg = ""
# elliptic:
## Indicates whether J_i^n==0 with J_i^(n-1)!=0 for n-folds and some J_i
for i in range(nodivs):
for intn in intCY:
if( intn[_sage_const_0 ] == [i]*(dimp-no_polys) and intn[_sage_const_1 ] == _sage_const_0 ):
if not(all([round((t[i]**(dimp-no_polys-_sage_const_1 )*t[j]).subs(Ccy),_sage_const_10 ) == _sage_const_0 for j in range(nodivs) if j != i]) ):
msg += "Possible elliptic fibration in cycle dual to t"+str(i)+".\n"
# K3 for 3-folds
## Indicates whether c2.J_i =24 for some J_i
if( dimp+_sage_const_1 -no_polys == _sage_const_4 ):
TupleList2 = UnorderedTuples(range(nodivs),_sage_const_2 )
for j in range(nodivs):
if( sum([((cc).subs({t:c*t for t in t}).taylor(c,_sage_const_0 ,_sage_const_2 ).coefficient(c,_sage_const_2 )).coefficient(product([t[i] for i in tupl]))*((product([t[i] for i in tupl])*t[j]).subs(Ccy)) for tupl in TupleList2]) == _sage_const_24 ):
msg += "Possible K3-fibration in divisor t"+str(j)+".\n"
if( msg != "" ):
print('\n--- Fibrations --------------------------------')
print(msg)
if returnvars:
return str([(MoriMatrix[_sage_const_0 ][j],MoriMatrix[_sage_const_1 ][j]) for j in range(nodivs)]).replace("(","{").replace(")","}").replace("[","{").replace("]","}") # might return more in the future if necessary
def CICYtopdata(CICY,justdata=False,justmori=False, lvec=_sage_const_0 ):
# Input: a list l with entries l[i,j] that give the weight in the
# ambient projective space i of polynomial j.
# E.g. (P^3| 3 1)
# (P^2| 2 0)
# corresponds to the list ((3,1),(2,0))
ambients=(np.array(CICY).transpose()).sum(axis=_sage_const_0 )-_sage_const_1
if ambients in ZZ:
ambients=[ambients]
dimp=sum(ambients)
polytope=identity_matrix(int(dimp))
zeros=[_sage_const_0 ]*dimp
polytope=polytope.insert_row(_sage_const_0 ,zeros)
offset=_sage_const_0
for i in range(len(ambients)):
newrow=[_sage_const_0 ]*dimp
entry=[]
for j in range(ambients[i]):
newrow[offset]=-_sage_const_1
offset+=_sage_const_1
entry.append(offset)
polytope=polytope.insert_row(sum([ambients[n]+_sage_const_1 for n in (ellipsis_range(_sage_const_0 ,Ellipsis,i))]),newrow)
partition=[]
eqs=matrix(CICY).transpose()
counters=[_sage_const_0 ]*len(ambients)
for eq in eqs:
part=[]
for j in range(len(eq)):
part.append([sum([ambients[m]+_sage_const_1 for m in (ellipsis_range(_sage_const_0 ,Ellipsis,(j-_sage_const_1 )))])+counters[j]+a for a in (ellipsis_range(_sage_const_0 ,Ellipsis,(eq[j]-_sage_const_1 )))])
counters[j]+=len(part[j])
partition.append(flatten(part))
if justdata:
return [list(polytope),partition]
elif justmori:
moricone((list(polytope)))
else:
topdata(list(polytope),nef_partition=partition,lvec=lvec)
topdata(pointstopoly([(_sage_const_1 ,_sage_const_0 ),(_sage_const_0 ,_sage_const_1 ),(-_sage_const_1 ,-_sage_const_1 )]))