{{{id=18|
PointConfiguration.set_engine('internal')
///
}}}

{{{id=17|
# this constructs a fan from a polytope by performing a fine star triangulation
# need TOPCOM for regularity

def P_to_F(P):
    Pint = P.integral_points()
    m = P.dim()
    origin = [0 for i in range(m)]
    PC = PointConfiguration(Pint).restrict_to_fine_triangulations(fine=True).restrict_to_star_triangulations(origin)

        
    tr = PC.triangulate()    
    F = Fan(tr,Pint)
    return F
///
}}}

{{{id=6|
# this prints (as many of) the 16 reflexive pairs in 2 dimensions + the FaceFan(blue) = NormalFan(red) polytope

for i in range(6,7):
    P = Polyhedron(vertices=ReflexivePolytope(2, i).vertices_pc())
    PD = P.polar()
    F = FaceFan(P)
    P.plot(fill=False) + F.plot(ray_color='black')+ PD.plot(fill=False, color = 'red')
///
}}}

{{{id=9|
# relation between monomials and points for P123 is
#[[(-1, -1), z^6],
# [(-1, 0), z^4*x],
# [(-1, 1), z^2*x^2],
# [(-1, 2), x^3],
# [(0, -1), z^3*y],
# [(0, 0), z*x*y],
# [(1, -1), y^2]]

# exercise: write a code that prints this output
///
}}}

{{{id=1|
from subprocess import Popen, PIPE
import subprocess
///
}}}

{{{id=2|
# this will use palp to compute the vertices of a polyhedron given the corresponding weight system; 
# input is a string, output is a list
# format for input
# a a0 a1 a2 a3 a4 ... an b b0 b1 b2 b3 b4 ... bn ...
# where a = sum a_i, b = sum b_i 


def palp_vertices(str):
    q1 = Popen(["echo",str], stdout=PIPE)
    q2 = Popen(["poly.x", "-ef"],stdin=q1.stdout,stdout=PIPE)
    A=q2.communicate()[0]
    print A
    for i in range(len(A)):
        if A[i]=='\n':
            B = A[-(len(A)-i+1):] 
            break
    Result = B.replace("P\n    ","[[").replace("P\n   ","[[").replace("P\n  ","[[").replace("P\n ","[[").replace("   ", ",").replace("  ", ",").replace(" ", ",").replace("\n, ","],[").replace("\n,","],[").replace("\n", "]]")
    return eval(Result)
///
}}}

{{{id=3|
# this computes the hodge numbers h0i i = 0,1,2 for a divisor of a CY threefold which descends from a toric divisor of the 
# ambient space; input is the polytope and the lattice point

def h0i(P,pt):
    
    # first we need to know in which face pt is sitting and if it is in the polytope at all

    if vector(pt) not in P.integral_points():
        print 'not contained in polytope'  

        return 0

    else:

        FF = FaceFan(P)
        d = P.dim()
    
        for i in range(d+1):
            for cone in FF.cones(i):
                if cone.relative_interior_contains(pt):
                    Pfaceverts = cone.rays()
                    k = i-1
        #print Pfaceverts
        
        
    # compute which are the vertices of the dual face on the M-lattice polytope    
        
        PD=P.polar()
        for face in PD.faces(d-k-1):

            intersections = [vector(dualpt)*vector(pt) for dualpt in list(face.as_polyhedron().vertices_matrix().transpose())]
            # print intersections 
            if all(int == -1 for int in intersections):
                Pdualfacevert = list(face.as_polyhedron().vertices_matrix().transpose())
                #print Pdualfacevert
                break   
                
    # now we need to know how many interior points this face has 
    
        dualfacep = Polyhedron(vertices=Pdualfacevert)     
        PCdF = PointConfiguration(dualfacep.integral_points())
        numintpts = len(PCdF.face_interior(dim=d-k-1))
        
        if d-k == d:
            h00=1
            h10=0
            h20=numintpts 
            return (h00,h10,h20)  
        if d-k == d-1:
            h00=1
            h10=numintpts
            h20=0 
            return (h00,h10,h20)  
        if d-k == d-2:
            h00= numintpts + 1
            h10=0
            h20=0  
            return (h00,h10,h20)     
        else:
            #print 'higher codimension than 2'
            return (0,0,0)
///
}}}

{{{id=4|
# this computes h11 using Batyrev's formula for Calabi-Yau threefolds

def H11(P):
    Pint = P.integral_points()
    h11 = 0
    for pt in Pint:
       
        if pt != vector((0,0,0,0)):
            #print pt
            hodgediv = h0i(P,pt)
            h11 = h11 + hodgediv[0]
    h11 = h11 -4
    return h11
///
}}}

{{{id=14|
from sage.geometry.polyhedron.plot import Projection
Projection(P).schlegel([1,0,0,0]).plot()
///
}}}

{{{id=13|

///
}}}