From 1efe3065ef6306390520314b6350cbaaddcbeecd Mon Sep 17 00:00:00 2001 From: julian Date: Fri, 24 Jul 2026 20:01:14 +0200 Subject: [PATCH] 3 refactor topdata (#8) --------- Co-authored-by: Julian Piribauer Reviewed-on: https://gitea.piribauer.ch/julian/Calabi-Yau-Period-Geometry/pulls/8 --- sage/toric_topdata.sage | 960 ++++++++++++++++++++++++++----------- sage/toric_topdata.sage.py | 348 -------------- 2 files changed, 681 insertions(+), 627 deletions(-) delete mode 100644 sage/toric_topdata.sage.py diff --git a/sage/toric_topdata.sage b/sage/toric_topdata.sage index 9ee0b96..c5fe0a3 100644 --- a/sage/toric_topdata.sage +++ b/sage/toric_topdata.sage @@ -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. + """ - 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 + ############################# + # 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 + # 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())] - 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 = [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") - - 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 - 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]) - - 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(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. - - # 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 - # - instance attributes: triangulation and Mori-cone generators - - pc = PointConfiguration(polytope) - - 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 + 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( 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) + 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())) - 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) + self.lvec = matrix(self.Mori_cone.rays()) - # 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 + 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() - 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 + 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: + 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 + + return polynomials + + def _append_unique_discriminants(self, discriminants, seen_discriminants, polynomials, index, total_relation_blocks, only_strong_coupling): + for polynomial in polynomials: + if polynomial == 0: + continue + + 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: + codim = total_relation_blocks - 1 - index + discriminants.append([polynomial, codim]) + + 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. + + 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. + """ + + logger.info("Computing discriminant for %s", self.model_name) + if only_strong_coupling: + logger.info("Restricting to strong-coupling loci (codimension-1 faces).") + + 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() + + 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) + + for index, relations_at_given_dimension in enumerate(all_linear_dependencies_among_points): + if len(relations_at_given_dimension) == 0: + continue + + 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, + ) + + if only_strong_coupling: + # Insert an empty codimension-0 entry for compatibility with prior behavior. + discriminants = [[]] + discriminants + + logger.info("Found %d discriminant factors.", len(discriminants)) + logger.info("Discriminant factors: %s", discriminants) + + 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/.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/.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] @@ -331,10 +726,17 @@ class ToricTopData: part.append([sum([ambients[m]+1 for m in [0..(j-1)]])+counters[j]+a for a in [0..(eq[j]-1)]]) counters[j]+=len(part[j]) partition.append(flatten(part)) + + return list(polytope), partition - if justdata: - return [list(polytope),partition] - elif justmori: - moricone((list(polytope))) - else: - self.topdata(list(polytope),nef_partition=partition,lvec=lvec) + + 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 + ) + \ No newline at end of file diff --git a/sage/toric_topdata.sage.py b/sage/toric_topdata.sage.py deleted file mode 100644 index e39c79c..0000000 --- a/sage/toric_topdata.sage.py +++ /dev/null @@ -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 )])) -