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