Initial refactoring

sage/toric_topdata.sage

@@ -1,309 +1,703 @@
-### 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
+import logging
+import os
+import json
+from datetime import datetime, timezone
+
+logger = logging.getLogger(__name__)
+logging.basicConfig(level=logging.DEBUG)  # or logging.DEBUG to see both
+
 # For grep-commands
 import re 
+find_zs = re.compile('z\d+')
+find_ls = re.compile('l\d+')
+find_index_z = re.compile('\d+')
 
 
-class ToricTopData:
+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
+    are excluded from the convex hull and stored in a separate list.
+    """
 
-    findzs = re.compile('z\d+')
+    #############################
-    findls = re.compile('l\d+')
+    # Initialisation
-    findindexz = re.compile('\d+')
+    #############################
 
-    def pointstopoly(self, points):
+    def _validate_points_and_get_dimension(self) -> int:
-        # Takes a list of points and returns the points inside its convex hull while omitting points inside faces of co-dimension one
+        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)
+        points_in_convex_hull.remove(self.origin)
-        polytope.append(zeros)
+        points_in_convex_hull.append(self.origin)
-        return polytope
 
-    def disc(self, polytope, only_sc=False, no_triangulation=0):
+        return lattice_polytope, points_in_convex_hull, points_in_codim1_faces
-        # 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
 
-        self.pc = LatticePolytope(polytope)
+    def __init__(self, points):
-        dimp = len(polytope[0])
+        self.defining_points = points
-        zeros = [0 for i in range(dimp)]
+        self.dimension = self._validate_points_and_get_dimension()
-        # Compute all relations in all faces
+        self.origin = [0] * self.dimension
-        self.globalks = []
+        self.lattice_polytope, self.relevant_lattice_points, self.points_in_codim1_faces = self._points_to_polytope()
-        if only_sc:
+
-            pcfaces = self.pc.faces()[:-1]
+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.
+    """
+
+    ##############################
+    # 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 = [range(len(self.relevant_lattice_points) - 1)]  # Default to a trivial partition if none is provided
 
-        for facesofdim in pcfaces: # for each face-dimension
+    def set_Mori_cone(self, lvec = None):
-            atdim = []
+        """
-            for face in facesofdim: # for each face of fixed dimension
+        Sets the Mori cone of the toric variety associated with the polytope.
-                atface = []
+        """
-                facepoints = [list(i) for i in matrix(face.points())]
+        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(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.")
+            logger.info("Using the first %d linearly independent vectors.", len(self.Mori_cone.rays()))
+            self.Mori_cone = Cone(
+                                transpose(
+                                    transpose(self.Mori_cone.rays())
+                                    *
+                                    matrix(
+                                            transpose(matrix(self.Mori_cone.rays()).kernel().gens())
+                                    ).kernel().gens()
+                                ).transpose()
+                             )
+        
+        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.fan = self.triangulation.fan(self.origin)
+        self.ToricVariety = ToricVariety(self.fan)
+        self.set_Mori_cone(lvec = 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
+
+        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():
+        return polynomials
-                    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
+    def _append_unique_discriminants(self, discriminants, seen_discriminants, polynomials, index, total_relation_blocks, only_strong_coupling):
-                kernel = matrix(facepoints).kernel().gens()
+        for polynomial in polynomials:
-                for k in kernel:                                     # just formatting
+            if polynomial == 0:
-                    globalk = [0 for i in range(len(polytope)-1)]    # MAY HAVE TO BE ADAPTED FOR CASES WITH MORE VERTICES
+                continue
-                    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))]
+            polynomial_key = str(polynomial)
-        z = var('z', n=len(ls)+1, latex_name='z') # z[0] is superfluous
+            if polynomial_key in seen_discriminants:
-        l = var('l', n=len(ls)+1, latex_name='l') # l[0] is superfluous
+                continue
-        #self.globalks.reverse()
+
-        self.globalks = [k for k in self.globalks if k not in ([[]],)]
+            seen_discriminants.add(polynomial_key)
-        discs = []
+            if only_strong_coupling:
-        for i,atdim in enumerate(self.globalks):
+                codim = total_relation_blocks - index
-            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
+                codim = total_relation_blocks - 1 - index
-            except:  # this is just for 1-parameter cases where there is nothing to solve
+            discriminants.append([polynomial, codim])
-                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
+    def disc(self, only_strong_coupling=False, no_triangulation=0):
-            discstmp = discs
+        """
-            discstmp.reverse()
+        Computes the A-discriminant for the toric polytope (following Aspinwall, Plesser, Wang),
-            discstmp.append([])
+        which gives the singular loci of the moduli space in Batyrev coordinates.
-            discstmp.reverse()
-            discs = discstmp
-        return discs
 
-    def topdata(self, polytope, nef_partition=0, lvec=0, no_triangulation=0, returnvars=False):
+        Parameters:
-        # Computes topological data for hypersurfaces and CICYs in toric ambient spaces.
+        - only_strong_coupling: If True, only gives loci arising from edges (one-dimensional faces).
-        # For hypersurfaces, nef_partition should remain untouched (=0).
+        - no_triangulation: Index of the triangulation to use.
-        # 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:
+        Returns:
-        #         - elliptic:           if t_i^n = 0 and t_i^{n-1} != 0
+        A list of tuples [disc_i, codim_i] of discriminant factors disc_i
-        #         - K3 (only 3-folds):  if c2.t_i = 24
+        coming from a relation inside a face of codimension codim_i.
+        """
 
-        # input:  - array of points in polytope
+        all_linear_dependencies_among_points = self._collect_linear_dependencies(only_strong_coupling)
-        #         - (for CICYs:) nef-partition
+        self._append_origin_weight(all_linear_dependencies_among_points)
-        #         - (optional:) set of l-vectors to be used
+        a_vars, z_vars, lambda_vars, a_row, mori_matrix = self._setup_discriminant_symbols()
-        #         - (optional:) index of triangulation
 
-        # output: - Mori-cone generators
+        all_linear_dependencies_among_points = [relations for relations in all_linear_dependencies_among_points if relations]
-        #         - Intersection rings of CY and ambient space
+        discriminants = []
-        #         - Topological data
+        seen_discriminants = set()
-        #         - Possible fibrations
+        total_relation_blocks = len(all_linear_dependencies_among_points)
-        #         - instance attributes: triangulation and Mori-cone generators
 
-        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
+            equation_system = self._build_equation_system(
-        if( nef_partition == 0):
+                relations_at_given_dimension,
-            no_polys = 1
+                a_vars,
-            nef_partition = [[i for i in range(len(polytope)-1)]] # trivial partition
+                z_vars,
-        else:
+                lambda_vars,
-            if( sorted(flatten(nef_partition)) != [ i for i in range(len(polytope)-1)] ):
+                a_row,
-                print("Nef-partition not valid!")
+                mori_matrix,
-            no_polys = len(nef_partition)
+            )
-        zeros=zero_vector(dimp)
+            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)
+        if only_strong_coupling:
-        pc_star_and_fine=pc_star.restrict_to_fine_triangulations()
+            # Insert an empty codimension-0 entry for compatibility with prior behavior.
-        if(len(pc_star_and_fine.triangulations_list())==0):
+            discriminants = [[]] + discriminants
-            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)
 
-        # Change l-vectors if given:
+        return discriminants
-        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
 
-        HH=self.X.cohomology_ring()
+    ##############################
-        self.D = [HH(self.X.divisor(i)) for i in [0..len(polytope)-2]]
+    # Topological Data Computation
-        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))))])
+    def _compute_Kahler_generators(self):
-        if (Bsinv.det()==0):
+        """
-            print("l-vectors not independent in divisor basis. Pick them manually with ``lvec=...''")
+        Computes the Kähler cone generators J_a dual to the Mori cone generators l^a.
-            return
+
+        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))))]))])
+        return Bs * vector([self.D[i] for i in z_indices])
-        # Hyperplane divisor class:
-        self.Ydualform = product([sum([HH(self.D[p]) for p in part]) for part in nef_partition])
 
-        # Tuple lists for computations below
+    def _compute_intersection_ring(self):
-        TupleListCY = UnorderedTuples(range(nodivs),dimp-no_polys)
+        """
-        TupleListCYordered = Tuples(range(nodivs),dimp-no_polys) # for intersection ring multiplicities must be included
+        Builds the intersection ring polynomial by adding up the intersection numbers
-        TupleListAm = UnorderedTuples(range(nodivs),dimp)
+        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:
+    def _intersection_substitution_rules(self):
-        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]
+        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:
+    def _replace_intersection_monomials(self, expression):
-        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])
+        Replaces top-degree Kähler monomials in an expression by CY intersection numbers.
-        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:
+        expanded_expression = expression.expand()
-        var(self.findzs.findall(str([self.J[i] for i in range(nodivs)])))                    # Introduces all z in Kähler-forms as variables...
+        replaced_expression = expanded_expression
-        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
 
-        print('--- Toric divisors (ambient space) -----------')
+        # Replace each degree-n monomial by its corresponding intersection number,
-        print(self.D)
+        # preserving the polynomial coefficient in front of that monomial.
-        print('\n--- Kähler cone generators (ambient space) --- ')
+        for idx_tuple, value in self.intersection_numbers_CY.items():
-        print(self.J)
+            monomial = product([self.t[i] for i in idx_tuple]).expand()
-        print('\n--- Mori-cone-generators (ambient space) ------')
+            coeff = expanded_expression.coefficient(monomial)
-        for j in range(nodivs):
+            if coeff != 0:
-            print(self.MoriMatrix[0][j],self.MoriMatrix[1][j])
+                replaced_expression -= coeff * monomial
+                replaced_expression += coeff * value
 
-        print('\n--- Intersection on CY ------------------------')
+        return replaced_expression.expand()
-        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)
 
-        print('\n--- Intersection in ambient space -------------')
+    def _compute_Chern_character(self):
-        Cam=[product([t[iden[0][i]] for i in range(dimp)])==iden[1] for iden in intAm]
+        """
-        print(Cam)
+        Computes the total Chern class of the CY via the adjunction formula:
 
-        print('\n--- Topological data --------------------------')
+            c(Y) = c(TX) / c(N_{Y/X})
-        #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)
 
+        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:
+        Chern_character_taylored = (Chern_character).subs({t: c * t for t in self.t}).taylor(c, 0, self.cy_dimension)
-        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))
 
-        msg = ""
+        Chern_polynomials = [
-        # elliptic:
+            [
-        ## Indicates whether J_i^n==0 with J_i^(n-1)!=0 for n-folds and some J_i
+                (prod(self.t[i] for i in tple) * Chern_character_taylored.coefficient(c, k)).expand()
-        for i in range(nodivs):
+                for tple in UnorderedTuples(range(self.no_divs), self.cy_dimension - k)
-            for intn in intCY:
+            ]
-                if( intn[0] == [i]*(dimp-no_polys) and intn[1] == 0):
+            for k in range(1, self.cy_dimension + 1) # for each degree
-                    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"
 
-        # K3 for 3-folds
+        integrated_Chern_classes = [
-        ## Indicates whether c2.J_i =24 for some J_i
+            [
-        if( dimp+1-no_polys == 4):
+                self._replace_intersection_monomials(poly) for poly in polys
-            TupleList2 = UnorderedTuples(range(nodivs),2)
+            ]
-            for j in range(nodivs):
+            for polys in Chern_polynomials
-                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"
 
-        if( msg != "" ):
+        return Chern_polynomials, integrated_Chern_classes
-            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
 
+    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):
+        for i in range(self.no_divs):
-        # Input: a list l with entries l[i,j] that give the weight in the
+            self_intersection = self.intersection_numbers_CY[tuple([i] * self.cy_dimension)]
-        #        ambient projective space i of polynomial j.
+            if self_intersection != 0:
-        #        E.g. (P^3| 3 1)
+                continue
-        #             (P^2| 2 0)
+            has_nonzero_neighbor = any(
-        #        corresponds to the list ((3,1),(2,0))
+                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" + "\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 (relative to the repository root), if a
+        model_name is given. Plain JSON, so it is directly loadable in Python
+        (json.load), C++ (e.g. nlohmann::json), or any other language without any
+        custom parsing code; lvec loads back as a plain list of lists, e.g. lvec[0].
+        """
+        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.
+        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:  - (for CICYs:) nef-partition
+                - (optional:) set of l-vectors to be used
+                - (optional:) index of triangulation
+        """
+        no_polys = len(self.nef_partition)
+        self.cy_dimension = self.dimension - no_polys
+        # 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 +726,17 @@ class ToricTopData:
                 counters[j]+=len(part[j])
             partition.append(flatten(part))
         
-        if justdata:
+        return list(polytope), partition
-            return [list(polytope),partition]
+
-        elif justmori:
+
-            moricone((list(polytope)))
+    def __init__(self, CICY, no_triangulation=0, model_name=None, lvec=None):
-        else:
+        points, nef_partition = self._CICY_to_points(CICY)
-            self.topdata(list(polytope),nef_partition=partition,lvec=lvec)
+        print(points, nef_partition)
+        super().__init__(
+            points,
+            no_triangulation=no_triangulation,
+            model_name=model_name,
+            nef_partition=nef_partition,
+            lvec=lvec
+        )
+