### 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)