Merge pull request 'Initial commit of topdata' (#2) from 1-topdata into master

Reviewed-on: #2

playground/.gitignore

@@ -0,0 +1,4 @@
+# Exclude testing files
+
+*
+!.gitignore

sage/toric_topdata.sage

@@ -0,0 +1,340 @@
+### 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 re
+
+
+class ToricTopData:
+
+    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
+
+        zeros = [0 for i in range(len(points[0]))] #origin
+        pc = LatticePolytope(points) # all points
+        polytope = [list(m) for m in pc.points()]
+        self.pcodim1 = [] #points in co-dimension one
+        for f in pc.facets():
+            for i in f.interior_point_indices():
+                try:
+                    polytope.remove([list(m) for m in f.points(i)][0]) # remove those inside codim1 faces
+                    self.pcodim1.append([list(m) for m in f.points(i)][0]) # save omitted points in pcodim1
+                except:
+                    pass
+        # Moving zeros to the end
+        polytope.remove(zeros)
+        polytope.append(zeros)
+        return polytope
+
+    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
+
+        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]
+        else:
+            pcfaces = self.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 = [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
+        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)
+
+        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)
+
+        # 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
+
+        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
+        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])
+
+        # 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,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]
+
+        # 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
+
+        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])
+
+        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)
+
+        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)
+
+        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 = 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))
+
+        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"
+
+        # 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"
+
+        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
+
+
+    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))
+
+        ambients=(np.array(CICY).transpose()).sum(axis=0)-1
+        if ambients in ZZ:
+            ambients=[ambients]
+        dimp=sum(ambients)
+        polytope=identity_matrix(int(dimp))
+        zeros=[0]*dimp
+        polytope=polytope.insert_row(0,zeros)
+
+        offset=0
+        for i in range(len(ambients)):
+            newrow=[0]*dimp
+            entry=[]
+            for j in range(ambients[i]):
+                newrow[offset]=-1
+                offset+=1
+                entry.append(offset)
+            polytope=polytope.insert_row(sum([ambients[n]+1 for n in [0..i]]),newrow)
+
+        partition=[]
+        eqs=matrix(CICY).transpose()
+        counters=[0]*len(ambients)
+        for eq in eqs:
+            part=[]
+            for j in range(len(eq)):
+                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))
+
+        if justdata:
+            return [list(polytope),partition]
+        elif justmori:
+            moricone((list(polytope)))
+        else:
+            self.topdata(list(polytope),nef_partition=partition,lvec=lvec)

sage/toric_topdata.sage.py

@@ -0,0 +1,348 @@
+
+
+# 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 )]))
+