################################################################
# Groebner cone, subdivision and initial forms computations Q4 #
################################################################


# In order to certify the calculations in "Qrefinementwithx34Noa4BadPrimes.sage" does not depend on the chosen weight vector, we need to check the  Dumbbell Cone and its faces are not subdivided by the Grobner complex of  Q4. This script computes the Grobner Cones for this polynomial and checks for potential subdivisions. In the presence of subdivisions (which happens for the Theta Cone) it computes sample points on each cell.

# Our calculations are valid when the residue field has characteristic 0. Based on the results, we later refine the calculation to include positive characteristic, based on the initial forms we obtain.

# The initial forms are computed in the script "initialFormsABCDQ4.sage"

# The necessary data (including the polynomials) and macros are obtained by loading the following file:

load("macroslocalJInvariantComputationsAndPrimesWithBadReduction.sage")

load("Qrefinementwithx34Noa4.sage")

# We create the fraction field of the ambient polynomial rings
K = FractionField(S); 
KK = FractionField(SS);



# The computation of the Grobner fan of Q4 was not done in "Qrefinementwithx34Noa4.sage" (it was done for Q4')




####################################### 
# Groebner complex computation for Q4 #
#######################################

# Recall: the variables of S are ordered a6,a5,a4,a3,a2,a1

Q = Q4

# We record the monomials and coefficients of Q, as well about the exponent vectors.
monQ = Q.monomials()
coeffQ = Q.coefficients()
exponentVectorsQ = [x.exponents() for x in monQ]

# We convert each the tuple of each entry (a list) into a list of integers.
exponentVectorsQInZZ = []
m = len(exponentVectorsQ[0][0])
for x in exponentVectorsQ:
      v = [ZZ(x[0][j]) for j in range(0,m)]
      exponentVectorsQInZZ.append(v)



#######################################
# Construct the Newton Polytope of Q4 #
#######################################

# The calculations of Grobner fans only require the vertices of the Newton polytope of Q. The number of vertices of this polytope is 366 compared to the 14111 monomials in Q. The computation of the Grobner fan of the full Q does not terminate, as opposed to that of the alternative polynomial AltQ.

NPQ = Polyhedron(vertices=  exponentVectorsQInZZ)
# NPQ
# A 5-dimensional polyhedron in ZZ^6 defined as the convex hull of 180 vertices


# we record the extremal monomials of Q.

extremalExponents = [tuple(list(x)) for x in NPQ.vertices()]
# print len(extremalExponents)
# 180

# We construct AltQ, a polynomial obtained as the sum of the terms of Q involving extremal monomials.

AltQMonomials = [exponentVectorsQInZZ.index(list(x)) for x in extremalExponents]

# len(AltQMonomials)
# 180

AltQ = S(0)
for j in AltQMonomials:
    AltQ = AltQ + S(coeffQ[j]*monQ[j])

# AltQ
# -8*a6^6*a5^6*a4^4*a3^2 - 8*a6^6*a5^4*a4^6*a3^2 - 8*a6^4*a5^6*a4^6*a3^2 - 8*a6^6*a5^6*a4^2*a3^4 - 8*a6^6*a5^2*a4^6*a3^4 - 8*a6^2*a5^6*a4^6*a3^4 - 8*a6^6*a5^4*a4^2*a3^6 - 8*a6^4*a5^6*a4^2*a3^6 - 8*a6^6*a5^2*a4^4*a3^6 - 8*a6^2*a5^6*a4^4*a3^6 - 8*a6^4*a5^2*a4^6*a3^6 - 8*a6^2*a5^4*a4^6*a3^6 - 8*a6^6*a5^6*a4^4*a2^2 - 8*a6^6*a5^4*a4^6*a2^2 - 8*a6^4*a5^6*a4^6*a2^2 - 8*a6^6*a5^6*a3^4*a2^2 - 8*a6^6*a4^6*a3^4*a2^2 - 8*a5^6*a4^6*a3^4*a2^2 - 8*a6^6*a5^4*a3^6*a2^2 - 8*a6^4*a5^6*a3^6*a2^2 - 8*a6^6*a4^4*a3^6*a2^2 - 8*a5^6*a4^4*a3^6*a2^2 - 8*a6^4*a4^6*a3^6*a2^2 - 8*a5^4*a4^6*a3^6*a2^2 - 8*a6^6*a5^6*a4^2*a2^4 - 8*a6^6*a5^2*a4^6*a2^4 - 8*a6^2*a5^6*a4^6*a2^4 - 8*a6^6*a5^6*a3^2*a2^4 - 8*a6^6*a4^6*a3^2*a2^4 - 8*a5^6*a4^6*a3^2*a2^4 - 8*a6^6*a5^2*a3^6*a2^4 - 8*a6^2*a5^6*a3^6*a2^4 - 8*a6^6*a4^2*a3^6*a2^4 - 8*a5^6*a4^2*a3^6*a2^4 - 8*a6^2*a4^6*a3^6*a2^4 - 8*a5^2*a4^6*a3^6*a2^4 - 8*a6^6*a5^4*a4^2*a2^6 - 8*a6^4*a5^6*a4^2*a2^6 - 8*a6^6*a5^2*a4^4*a2^6 - 8*a6^2*a5^6*a4^4*a2^6 - 8*a6^4*a5^2*a4^6*a2^6 - 8*a6^2*a5^4*a4^6*a2^6 - 8*a6^6*a5^4*a3^2*a2^6 - 8*a6^4*a5^6*a3^2*a2^6 - 8*a6^6*a4^4*a3^2*a2^6 - 8*a5^6*a4^4*a3^2*a2^6 - 8*a6^4*a4^6*a3^2*a2^6 - 8*a5^4*a4^6*a3^2*a2^6 - 8*a6^6*a5^2*a3^4*a2^6 - 8*a6^2*a5^6*a3^4*a2^6 - 8*a6^6*a4^2*a3^4*a2^6 - 8*a5^6*a4^2*a3^4*a2^6 - 8*a6^2*a4^6*a3^4*a2^6 - 8*a5^2*a4^6*a3^4*a2^6 - 8*a6^4*a5^2*a3^6*a2^6 - 8*a6^2*a5^4*a3^6*a2^6 - 8*a6^4*a4^2*a3^6*a2^6 - 8*a5^4*a4^2*a3^6*a2^6 - 8*a6^2*a4^4*a3^6*a2^6 - 8*a5^2*a4^4*a3^6*a2^6 - 8*a6^6*a5^6*a4^4*a1^2 - 8*a6^6*a5^4*a4^6*a1^2 - 8*a6^4*a5^6*a4^6*a1^2 - 8*a6^6*a5^6*a3^4*a1^2 - 8*a6^6*a4^6*a3^4*a1^2 - 8*a5^6*a4^6*a3^4*a1^2 - 8*a6^6*a5^4*a3^6*a1^2 - 8*a6^4*a5^6*a3^6*a1^2 - 8*a6^6*a4^4*a3^6*a1^2 - 8*a5^6*a4^4*a3^6*a1^2 - 8*a6^4*a4^6*a3^6*a1^2 - 8*a5^4*a4^6*a3^6*a1^2 - 8*a6^6*a5^6*a2^4*a1^2 - 8*a6^6*a4^6*a2^4*a1^2 - 8*a5^6*a4^6*a2^4*a1^2 - 8*a6^6*a3^6*a2^4*a1^2 - 8*a5^6*a3^6*a2^4*a1^2 - 8*a4^6*a3^6*a2^4*a1^2 - 8*a6^6*a5^4*a2^6*a1^2 - 8*a6^4*a5^6*a2^6*a1^2 - 8*a6^6*a4^4*a2^6*a1^2 - 8*a5^6*a4^4*a2^6*a1^2 - 8*a6^4*a4^6*a2^6*a1^2 - 8*a5^4*a4^6*a2^6*a1^2 - 8*a6^6*a3^4*a2^6*a1^2 - 8*a5^6*a3^4*a2^6*a1^2 - 8*a4^6*a3^4*a2^6*a1^2 - 8*a6^4*a3^6*a2^6*a1^2 - 8*a5^4*a3^6*a2^6*a1^2 - 8*a4^4*a3^6*a2^6*a1^2 - 8*a6^6*a5^6*a4^2*a1^4 - 8*a6^6*a5^2*a4^6*a1^4 - 8*a6^2*a5^6*a4^6*a1^4 - 8*a6^6*a5^6*a3^2*a1^4 - 8*a6^6*a4^6*a3^2*a1^4 - 8*a5^6*a4^6*a3^2*a1^4 - 8*a6^6*a5^2*a3^6*a1^4 - 8*a6^2*a5^6*a3^6*a1^4 - 8*a6^6*a4^2*a3^6*a1^4 - 8*a5^6*a4^2*a3^6*a1^4 - 8*a6^2*a4^6*a3^6*a1^4 - 8*a5^2*a4^6*a3^6*a1^4 - 8*a6^6*a5^6*a2^2*a1^4 - 8*a6^6*a4^6*a2^2*a1^4 - 8*a5^6*a4^6*a2^2*a1^4 - 8*a6^6*a3^6*a2^2*a1^4 - 8*a5^6*a3^6*a2^2*a1^4 - 8*a4^6*a3^6*a2^2*a1^4 - 8*a6^6*a5^2*a2^6*a1^4 - 8*a6^2*a5^6*a2^6*a1^4 - 8*a6^6*a4^2*a2^6*a1^4 - 8*a5^6*a4^2*a2^6*a1^4 - 8*a6^2*a4^6*a2^6*a1^4 - 8*a5^2*a4^6*a2^6*a1^4 - 8*a6^6*a3^2*a2^6*a1^4 - 8*a5^6*a3^2*a2^6*a1^4 - 8*a4^6*a3^2*a2^6*a1^4 - 8*a6^2*a3^6*a2^6*a1^4 - 8*a5^2*a3^6*a2^6*a1^4 - 8*a4^2*a3^6*a2^6*a1^4 - 8*a6^6*a5^4*a4^2*a1^6 - 8*a6^4*a5^6*a4^2*a1^6 - 8*a6^6*a5^2*a4^4*a1^6 - 8*a6^2*a5^6*a4^4*a1^6 - 8*a6^4*a5^2*a4^6*a1^6 - 8*a6^2*a5^4*a4^6*a1^6 - 8*a6^6*a5^4*a3^2*a1^6 - 8*a6^4*a5^6*a3^2*a1^6 - 8*a6^6*a4^4*a3^2*a1^6 - 8*a5^6*a4^4*a3^2*a1^6 - 8*a6^4*a4^6*a3^2*a1^6 - 8*a5^4*a4^6*a3^2*a1^6 - 8*a6^6*a5^2*a3^4*a1^6 - 8*a6^2*a5^6*a3^4*a1^6 - 8*a6^6*a4^2*a3^4*a1^6 - 8*a5^6*a4^2*a3^4*a1^6 - 8*a6^2*a4^6*a3^4*a1^6 - 8*a5^2*a4^6*a3^4*a1^6 - 8*a6^4*a5^2*a3^6*a1^6 - 8*a6^2*a5^4*a3^6*a1^6 - 8*a6^4*a4^2*a3^6*a1^6 - 8*a5^4*a4^2*a3^6*a1^6 - 8*a6^2*a4^4*a3^6*a1^6 - 8*a5^2*a4^4*a3^6*a1^6 - 8*a6^6*a5^4*a2^2*a1^6 - 8*a6^4*a5^6*a2^2*a1^6 - 8*a6^6*a4^4*a2^2*a1^6 - 8*a5^6*a4^4*a2^2*a1^6 - 8*a6^4*a4^6*a2^2*a1^6 - 8*a5^4*a4^6*a2^2*a1^6 - 8*a6^6*a3^4*a2^2*a1^6 - 8*a5^6*a3^4*a2^2*a1^6 - 8*a4^6*a3^4*a2^2*a1^6 - 8*a6^4*a3^6*a2^2*a1^6 - 8*a5^4*a3^6*a2^2*a1^6 - 8*a4^4*a3^6*a2^2*a1^6 - 8*a6^6*a5^2*a2^4*a1^6 - 8*a6^2*a5^6*a2^4*a1^6 - 8*a6^6*a4^2*a2^4*a1^6 - 8*a5^6*a4^2*a2^4*a1^6 - 8*a6^2*a4^6*a2^4*a1^6 - 8*a5^2*a4^6*a2^4*a1^6 - 8*a6^6*a3^2*a2^4*a1^6 - 8*a5^6*a3^2*a2^4*a1^6 - 8*a4^6*a3^2*a2^4*a1^6 - 8*a6^2*a3^6*a2^4*a1^6 - 8*a5^2*a3^6*a2^4*a1^6 - 8*a4^2*a3^6*a2^4*a1^6 - 8*a6^4*a5^2*a2^6*a1^6 - 8*a6^2*a5^4*a2^6*a1^6 - 8*a6^4*a4^2*a2^6*a1^6 - 8*a5^4*a4^2*a2^6*a1^6 - 8*a6^2*a4^4*a2^6*a1^6 - 8*a5^2*a4^4*a2^6*a1^6 - 8*a6^4*a3^2*a2^6*a1^6 - 8*a5^4*a3^2*a2^6*a1^6 - 8*a4^4*a3^2*a2^6*a1^6 - 8*a6^2*a3^4*a2^6*a1^6 - 8*a5^2*a3^4*a2^6*a1^6 - 8*a4^2*a3^4*a2^6*a1^6




#########################################
# Calculation of the Groebner fan of Q4 #
#########################################

# We use AltQ to compute the Grobner fan of Q4:

# Iq = ideal(AltQ)
# GFq = Iq.groebner_fan()

# We pre-compute the Grobner fan:
# GFqFan = GFq.polyhedralfan()
# save(GFqFan,"GrobnerFanAltQ4.sobj")

# We convert the data to Rational Polyhedral fan:
# newGF = GFqFan.to_RationalPolyhedralFan()
# save(newGF, "GFq4Fan_as_RationalPolyhedra.sobj")
# newGF=load("GFq4Fan_as_RationalPolyhedra.sobj")

# We load the pre-computed the Grober fan.
GFqFan = load("GrobnerFanAltQ4.sobj")


# We record the f-vector:
# print GFqFan.f_vector()
# [1, 32, 180, 420, 450, 180]

# We confirm the dimension of the fan is the expected one:
# print GFqFan.dim()
# 6

# We record the coordinates of the rays of the Grobner fan:
RaysGFqfan = GFqFan.rays()
# print RaysGFqfan
# [[-5, 1, 1, 1, 1, 1], [-1, -1, -1, -1, -1, 5], [-1, -1, -1, -1, 5, -1], [-1, -1, -1, 1, 1, 1], [-1, -1, -1, 5, -1, -1], [-1, -1, 1, -1, 1, 1], [-1, -1, 1, 1, -1, 1], [-1, -1, 1, 1, 1, -1], [-1, -1, 5, -1, -1, -1], [-1, 1, -1, -1, 1, 1], [-1, 1, -1, 1, -1, 1], [-1, 1, -1, 1, 1, -1], [-1, 1, 1, -1, -1, 1], [-1, 1, 1, -1, 1, -1], [-1, 1, 1, 1, -1, -1], [-1, 5, -1, -1, -1, -1], [1, -5, 1, 1, 1, 1], [1, -1, -1, -1, 1, 1], [1, -1, -1, 1, -1, 1], [1, -1, -1, 1, 1, -1], [1, -1, 1, -1, -1, 1], [1, -1, 1, -1, 1, -1], [1, -1, 1, 1, -1, -1], [1, 1, -5, 1, 1, 1], [1, 1, -1, -1, -1, 1], [1, 1, -1, -1, 1, -1], [1, 1, -1, 1, -1, -1], [1, 1, 1, -5, 1, 1], [1, 1, 1, -1, -1, -1], [1, 1, 1, 1, -5, 1], [1, 1, 1, 1, 1, -5], [5, -1, -1, -1, -1, -1]]


# We record the maximal cones of the fan. Each cone is described by the indices of the extremal rays spanning these cones. Such indices refer to the list of 32 rays:
MaximalConesQAlt = GFqFan.maximal_cones().values()[0]
# print MaximalConesQAlt
# [[0, 1, 2, 3, 16], [0, 1, 3, 4, 16], [0, 2, 3, 4, 16], [0, 1, 2, 5, 16], [0, 1, 4, 6, 16], [0, 2, 4, 7, 16], [0, 1, 5, 8, 16], [0, 2, 5, 8, 16], [0, 1, 6, 8, 16], [0, 2, 7, 8, 16], [0, 4, 6, 8, 16], [0, 4, 7, 8, 16], [0, 1, 2, 3, 23], [0, 1, 3, 4, 23], [0, 2, 3, 4, 23], [0, 1, 2, 5, 27], [0, 1, 4, 6, 29], [0, 2, 4, 7, 30], [0, 1, 5, 8, 27], [0, 2, 5, 8, 27], [0, 1, 6, 8, 29], [0, 2, 7, 8, 30], [0, 4, 6, 8, 29], [0, 4, 7, 8, 30], [0, 1, 2, 9, 23], [0, 1, 4, 10, 23], [0, 2, 4, 11, 23], [0, 1, 2, 9, 27], [0, 1, 4, 10, 29], [0, 2, 4, 11, 30], [0, 1, 8, 12, 27], [0, 2, 8, 13, 27], [0, 1, 8, 12, 29], [0, 2, 8, 13, 30], [0, 4, 8, 14, 29], [0, 4, 8, 14, 30], [0, 1, 9, 15, 23], [0, 2, 9, 15, 23], [0, 1, 10, 15, 23], [0, 2, 11, 15, 23], [0, 4, 10, 15, 23], [0, 4, 11, 15, 23], [0, 1, 9, 15, 27], [0, 2, 9, 15, 27], [0, 1, 10, 15, 29], [0, 2, 11, 15, 30], [0, 4, 10, 15, 29], [0, 4, 11, 15, 30], [0, 1, 12, 15, 27], [0, 2, 13, 15, 27], [0, 1, 12, 15, 29], [0, 2, 13, 15, 30], [0, 4, 14, 15, 29], [0, 4, 14, 15, 30], [0, 8, 12, 15, 27], [0, 8, 13, 15, 27], [0, 8, 12, 15, 29], [0, 8, 13, 15, 30], [0, 8, 14, 15, 29], [0, 8, 14, 15, 30], [1, 2, 3, 16, 23], [1, 3, 4, 16, 23], [2, 3, 4, 16, 23], [1, 2, 5, 16, 27], [1, 4, 6, 16, 29], [2, 4, 7, 16, 30], [1, 5, 8, 16, 27], [2, 5, 8, 16, 27], [1, 6, 8, 16, 29], [2, 7, 8, 16, 30], [4, 6, 8, 16, 29], [4, 7, 8, 16, 30], [1, 2, 9, 23, 27], [1, 4, 10, 23, 29], [2, 4, 11, 23, 30], [1, 8, 12, 27, 29], [2, 8, 13, 27, 30], [4, 8, 14, 29, 30], [1, 9, 15, 23, 27], [2, 9, 15, 23, 27], [1, 10, 15, 23, 29], [2, 11, 15, 23, 30], [4, 10, 15, 23, 29], [4, 11, 15, 23, 30], [1, 12, 15, 27, 29], [2, 13, 15, 27, 30], [4, 14, 15, 29, 30], [8, 12, 15, 27, 29], [8, 13, 15, 27, 30], [8, 14, 15, 29, 30], [1, 2, 16, 17, 23], [1, 4, 16, 18, 23], [2, 4, 16, 19, 23], [1, 2, 16, 17, 27], [1, 4, 16, 18, 29], [2, 4, 16, 19, 30], [1, 8, 16, 20, 27], [2, 8, 16, 21, 27], [1, 8, 16, 20, 29], [2, 8, 16, 21, 30], [4, 8, 16, 22, 29], [4, 8, 16, 22, 30], [1, 2, 17, 23, 27], [1, 4, 18, 23, 29], [2, 4, 19, 23, 30], [1, 8, 20, 27, 29], [2, 8, 21, 27, 30], [4, 8, 22, 29, 30], [1, 15, 23, 24, 27], [2, 15, 23, 25, 27], [1, 15, 23, 24, 29], [2, 15, 23, 25, 30], [4, 15, 23, 26, 29], [4, 15, 23, 26, 30], [1, 15, 24, 27, 29], [2, 15, 25, 27, 30], [4, 15, 26, 29, 30], [8, 15, 27, 28, 29], [8, 15, 27, 28, 30], [8, 15, 28, 29, 30], [1, 16, 17, 23, 31], [2, 16, 17, 23, 31], [1, 16, 18, 23, 31], [2, 16, 19, 23, 31], [4, 16, 18, 23, 31], [4, 16, 19, 23, 31], [1, 16, 17, 27, 31], [2, 16, 17, 27, 31], [1, 16, 18, 29, 31], [2, 16, 19, 30, 31], [4, 16, 18, 29, 31], [4, 16, 19, 30, 31], [1, 16, 20, 27, 31], [2, 16, 21, 27, 31], [1, 16, 20, 29, 31], [2, 16, 21, 30, 31], [4, 16, 22, 29, 31], [4, 16, 22, 30, 31], [8, 16, 20, 27, 31], [8, 16, 21, 27, 31], [8, 16, 20, 29, 31], [8, 16, 21, 30, 31], [8, 16, 22, 29, 31], [8, 16, 22, 30, 31], [1, 17, 23, 27, 31], [2, 17, 23, 27, 31], [1, 18, 23, 29, 31], [2, 19, 23, 30, 31], [4, 18, 23, 29, 31], [4, 19, 23, 30, 31], [1, 20, 27, 29, 31], [2, 21, 27, 30, 31], [4, 22, 29, 30, 31], [8, 20, 27, 29, 31], [8, 21, 27, 30, 31], [8, 22, 29, 30, 31], [1, 23, 24, 27, 31], [2, 23, 25, 27, 31], [1, 23, 24, 29, 31], [2, 23, 25, 30, 31], [4, 23, 26, 29, 31], [4, 23, 26, 30, 31], [1, 24, 27, 29, 31], [2, 25, 27, 30, 31], [4, 26, 29, 30, 31], [8, 27, 28, 29, 31], [8, 27, 28, 30, 31], [8, 28, 29, 30, 31], [15, 23, 24, 27, 31], [15, 23, 25, 27, 31], [15, 23, 24, 29, 31], [15, 23, 25, 30, 31], [15, 23, 26, 29, 31], [15, 23, 26, 30, 31], [15, 24, 27, 29, 31], [15, 25, 27, 30, 31], [15, 26, 29, 30, 31], [15, 27, 28, 29, 31], [15, 27, 28, 30, 31], [15, 28, 29, 30, 31]]


# We record the maximal cones in the Grobner Fan of Q as a dictionary of polyhedral objects. The keys coincide with the indices of the maximal cones in the output of the Gfan computation.

PolyMaximalCones = MaximalCones(GFqFan)

# We check that they have the expected dimension:
# set([dim(x) for x in PolyMaximalCones.values()])
# {6}



###############################################
# Subdivisions induced on each non-Theta Cone #
###############################################

# We start by loading the cones:
load("allConesGenus2.sage")



############################################################
# Subdivision of Dumbbell induced by the Grobner fan of Q4 #
############################################################

# The following script computes the intersection of the Dumbbell cone with all maximal cells in the Grobner fan of Q4.

IntersectionsQDumbbellCone = {}
for i in range(0,len(PolyMaximalCones.keys())):
    # print str(i)
    polyCone =PolyMaximalCones[i]
    IntersectionQPolyCone = DumbbellCone.intersection(polyCone)
    IntersectionsQDumbbellCone[i] = IntersectionQPolyCone

set([dim(x) for x in IntersectionsQDumbbellCone.values()])
# {1, 2, 3, 4, 5, 6}


# We record those intersection of maximal dimension, i.e. the dimension of the Theta Cone.
maxIntersections = [i for i in IntersectionsQDumbbellCone.keys() if dim(IntersectionsQDumbbellCone[i])==6]

# print maxIntersections
# [179]


SubdivisionsOfDumbbellByQ = [IntersectionsQDumbbellCone[i] for i in maxIntersections]

# print [len(SubdivisionsOfDumbbellByQ[i].rays()) for i in range(0,len(SubdivisionsOfDumbbellByQ))]
# [5]

# print [dim(SubdivisionsOfDumbbellByQ[i]) for i in range(0,len(SubdivisionsOfDumbbellByQ))]
# [6]

# No subdivision is produced.





################################################################
# Subdivision of Figure Eight induced by the Grobner fan of Q4 #
################################################################

# The following script computes the intersection of the Figure Eight cone with all maximal cells in the Grobner fan of Q4.

IntersectionsQFigEightCone = {}
for i in range(0,len(PolyMaximalCones.keys())):
    # print str(i)
    polyCone =PolyMaximalCones[i]
    IntersectionQPolyCone = FigEightCone.intersection(polyCone)
    IntersectionsQFigEightCone[i] = IntersectionQPolyCone

set([dim(x) for x in IntersectionsQFigEightCone.values()])
# {1, 2, 3, 4, 5}


# We record those intersection of maximal dimension, i.e. the dimension of the Theta Cone.
maxIntersections = [i for i in IntersectionsQFigEightCone.keys() if dim(IntersectionsQFigEightCone[i])==5]

# print maxIntersections
# [176, 179]


SubdivisionsOfFigEightByQ = [IntersectionsQFigEightCone[i] for i in maxIntersections]

# print [len(SubdivisionsOfFigEightByQ[i].rays()) for i in range(0,len(SubdivisionsOfFigEightByQ))]
# [4, 4]

# print [dim(SubdivisionsOfFigEightByQ[i]) for i in range(0,len(SubdivisionsOfFigEightByQ))]
# [5, 5]

# We check the two cells agree:

# SubdivisionsOfFigEightByQ[0] == SubdivisionsOfFigEightByQ[1]
# True

# No subdivision is produced.

###########################################################
# Subdivision of Type IV induced by the Grobner fan of Q4 #
###########################################################

# The following script computes the intersection of the Type IV cone with all maximal cells in the Grobner fan of Q4.

IntersectionsQTypeIVCone = {}
for i in range(0,len(PolyMaximalCones.keys())):
    # print str(i)
    polyCone =PolyMaximalCones[i]
    IntersectionQPolyCone = TypeIVCone.intersection(polyCone)
    IntersectionsQTypeIVCone[i] = IntersectionQPolyCone

set([dim(x) for x in IntersectionsQTypeIVCone.values()])
# {1, 2, 3, 4, 5}


# We record those intersection of maximal dimension, i.e. the dimension of the Theta Cone.
maxIntersections = [i for i in IntersectionsQTypeIVCone.keys() if dim(IntersectionsQTypeIVCone[i])==5]

# print maxIntersections
# [167, 179]

SubdivisionsOfTypeIVByQ = [IntersectionsQTypeIVCone[i] for i in maxIntersections]
			
# print [len(SubdivisionsOfTypeIVByQ[i].rays()) for i in range(0,len(SubdivisionsOfTypeIVByQ))]
# [4, 4]

# print [dim(SubdivisionsOfTypeIVByQ[i]) for i in range(0,len(SubdivisionsOfTypeIVByQ))]
# [5, 5]

# We check that these intersections are the same. We conclude there is no subdivision.

# print [(i,j) for i in range(0,len(maxIntersections)-1) for j in range(i+1,len(maxIntersections)) if SubdivisionsOfTypeIVByQ[i] == SubdivisionsOfTypeIVByQ[j]]
# [(0,1)]

# No subdivision is produced.


##########################################################
# Subdivision of Type V induced by the Grobner fan of Q4 #
##########################################################

# The following script computes the intersection of the Type V cone with all maximal cells in the Grobner fan of Q4.

IntersectionsQTypeVCone = {}
for i in range(0,len(PolyMaximalCones.keys())):
    # print str(i)
    polyCone =PolyMaximalCones[i]
    IntersectionQPolyCone = TypeVCone.intersection(polyCone)
    IntersectionsQTypeVCone[i] = IntersectionQPolyCone

set([dim(x) for x in IntersectionsQTypeVCone.values()])
# {1, 2, 3, 4}


# We record those intersection of maximal dimension, i.e. the dimension of the Theta Cone.
maxIntersections = [i for i in IntersectionsQTypeVCone.keys() if dim(IntersectionsQTypeVCone[i])==4]

# print maxIntersections
# [166, 167, 178, 179]

SubdivisionsOfTypeVByQ = [IntersectionsQTypeVCone[i] for i in maxIntersections]

# print [len(SubdivisionsOfTypeVByQ[i].rays()) for i in range(0,len(SubdivisionsOfTypeVByQ))]
# [3, 3, 3, 3]

# print [dim(SubdivisionsOfTypeVByQ[i]) for i in range(0,len(SubdivisionsOfTypeVByQ))]
# [4, 4, 4, 4]

# We check that these intersections are the same. We conclude there is no subdivision.

# print [(i,j) for i in range(0,len(maxIntersections)-1) for j in range(i+1,len(maxIntersections)) if SubdivisionsOfTypeVByQ[i] == SubdivisionsOfTypeVByQ[j]]
# [(0, 1), (0, 2), (0, 3), (1, 2), (1, 3), (2, 3)]

# No subdivision is produced.


###########################################################
# Subdivision of Type VI induced by the Grobner fan of Q4 #
###########################################################

# The following script computes the intersection of the Type VI cone with all maximal cells in the Grobner fan of Q4.

IntersectionsQTypeVICone = {}
for i in range(0,len(PolyMaximalCones.keys())):
    # print str(i)
    polyCone =PolyMaximalCones[i]
    IntersectionQPolyCone = TypeVICone.intersection(polyCone)
    IntersectionsQTypeVICone[i] = IntersectionQPolyCone

set([dim(x) for x in IntersectionsQTypeVICone.values()])
# {1, 2, 3, 4}


# We record those intersection of maximal dimension, i.e. the dimension of the Theta Cone.
maxIntersections = [i for i in IntersectionsQTypeVICone.keys() if dim(IntersectionsQTypeVICone[i])==4]

# print maxIntersections
# [152, 155, 164, 167, 176, 179]



SubdivisionsOfTypeVIByQ = [IntersectionsQTypeVICone[i] for i in maxIntersections]

# print [len(SubdivisionsOfTypeVIByQ[i].rays()) for i in range(0,len(SubdivisionsOfTypeVIByQ))]
# [3, 3, 3, 3, 3, 3]


# print [dim(SubdivisionsOfTypeVIByQ[i]) for i in range(0,len(SubdivisionsOfTypeVIByQ))]
# [4, 4, 4, 4, 4, 4]


# We check that these intersections are the same. We conclude there is no subdivision.

# all( [SubdivisionsOfTypeVIByQ[i] == SubdivisionsOfTypeVIByQ[j] for i in range(0,len(maxIntersections)-1) for j in range(i+1,len(maxIntersections))])
# True


# No subdivision is produced.

############################################################
# Subdivision of Type VII induced by the Grobner fan of Q4 #
############################################################

# The following script computes the intersection of the Type VII cone with all maximal cells in the Grobner fan of Q4.

IntersectionsQTypeVIICone = {}
for i in range(0,len(PolyMaximalCones.keys())):
    # print str(i)
    polyCone =PolyMaximalCones[i]
    IntersectionQPolyCone = TypeVIICone.intersection(polyCone)
    IntersectionsQTypeVIICone[i] = IntersectionQPolyCone

set([dim(x) for x in IntersectionsQTypeVIICone.values()])
# {1, 2, 3}


# We record those intersection of maximal dimension, i.e. the dimension of the Theta Cone.
maxIntersections = [i for i in IntersectionsQTypeVIICone.keys() if dim(IntersectionsQTypeVIICone[i])==3]

# print maxIntersections
# [129, 131, 135, 137, 141, 143, 147, 149, 151, 152, 154, 155, 159, 161, 163, 164, 166, 167, 171, 173, 175, 176, 178, 179]


SubdivisionsOfTypeVIIByQ = [IntersectionsQTypeVIICone[i] for i in maxIntersections]

# print [len(SubdivisionsOfTypeVIIByQ[i].rays()) for i in range(0,len(SubdivisionsOfTypeVIIByQ))]
# [2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2]


# print [dim(SubdivisionsOfTypeVIIByQ[i]) for i in range(0,len(SubdivisionsOfTypeVIIByQ))]
# [3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3]


# We check that these intersections are the same. We conclude there is no subdivision.

# all( [SubdivisionsOfTypeVIIByQ[i] == SubdivisionsOfTypeVIIByQ[j] for i in range(0,len(maxIntersections)-1) for j in range(i+1,len(maxIntersections))])
# True

# No subdivision is produced.


#############################################################################
# CONCLUSION: No subdivisions are produced and we can use any sample point to compute initial forms. 
##############################################################################








######################################
# Subdivisions induced on Theta Cone #
######################################


##############
# Theta Cone #
##############

# Recall: the variables of SS are ordered a6,a5,x34,a3,a2,a1

# The following script constructs the Theta cone, with respect to the variables a6,a5,x34,a3,a2,a1:
# Its defining inequalities are:
# a6 - a5 >= 0
# a5 - a3 >= 0
# a3 - a2 >= 0
# a2 - a1 >= 0
# a3 - x34 >= 0

ThetaCone = Polyhedron(ieqs=[(0,1,-1,0,0,0,0),(0,0,1,0,-1,0,0), (0,0,0,0,1,-1,0),(0,0,0,0,0,1,-1), (0,0,0,-1,1,0,0)])

# ThetaCone.Hrepresentation()
# (An inequality (0, 0, -1, 1, 0, 0) x + 0 >= 0,
#  An inequality (0, 0, 0, 0, 1, -1) x + 0 >= 0,
#  An inequality (0, 0, 0, 1, -1, 0) x + 0 >= 0,
#  An inequality (0, 1, 0, -1, 0, 0) x + 0 >= 0,
#  An inequality (1, -1, 0, 0, 0, 0) x + 0 >= 0)

# ThetaCone.rays()
# (A ray in the direction (0, -1, -1, -1, -1, -1),
#  A ray in the direction (0, 0, -1, -1, -1, -1),
#  A ray in the direction (0, 0, 0, 0, -1, -1),
#  A ray in the direction (0, 0, 0, 0, 0, -1),
#  A ray in the direction (0, 0, -1, 0, 0, 0))

# ThetaCone.lines()
# (A line in the direction (1, 1, 1, 1, 1, 1),)


# print dim(ThetaCone)
# 6


#########################################################
# Subdivision of Theta induced by the Grobner fan of Q4 #
#########################################################

# We load the Groebner fan and compute the dictionary of maximal cones. This fan was computed in "Qrefinementwithx34Noa4.sage"

GFqFan = load("GrobnerFanAltQ4x34.sobj")
PolyMaximalCones = MaximalCones(GFqFan)

# The following script computes the intersection of the Theta cone with all maximal cells in the Grobner fan of Q4.

IntersectionsQThetaCone = {}
for i in range(0,len(PolyMaximalCones.keys())):
    # print str(i)
    polyCone =PolyMaximalCones[i]
    IntersectionQPolyCone = ThetaCone.intersection(polyCone)
    IntersectionsQThetaCone[i] = IntersectionQPolyCone

set([dim(x) for x in IntersectionsQThetaCone.values()])
# {1, 2, 3, 4, 5, 6}


# We record those intersection of maximal dimension, i.e. the dimension of the Theta Cone.
maxIntersections = [i for i in IntersectionsQThetaCone.keys() if dim(IntersectionsQThetaCone[i])==6]

# print maxIntersections
# [163, 164, 167]


SubdivisionsOfThetaByQ = [IntersectionsQThetaCone[i] for i in maxIntersections]

# print [len(SubdivisionsOfThetaByQ[i].rays()) for i in range(0,len(SubdivisionsOfThetaByQ))]
# [5, 5, 5]

# print [dim(SubdivisionsOfThetaByQ[i]) for i in range(0,len(SubdivisionsOfThetaByQ))]
# [6, 6, 6]



# We check that these 3 intersections are distinct. We know there are no containments by construction since they are maximal dimension and the maximal cells of Q have disjoint relative interiors. Thus, there are precisely 3 cells in the subdivision of Theta induced by the Grobner fan of Q.

# SubdivisionsOfThetaByQ[0] == SubdivisionsOfThetaByQ[1]
# # False
# SubdivisionsOfThetaByQ[0] == SubdivisionsOfThetaByQ[2]
# # False
# SubdivisionsOfThetaByQ[1] == SubdivisionsOfThetaByQ[2]
# # False




########################################
# Subdivision pieces and sample points #
########################################

# We record the list of rays of each of the 3 cells of the Theta cone, the polyhedron they generate, and a point in the relative interior of each cone. The latter will allow us to compute the initial form of Q induced by the weight given by each of these points.

raysOfThetaGfan = {}
ThetaSubdividedGfan = {}
pointThetaGfan = {}

CellIndices = [0,1,2]
for i in CellIndices:
    raysOfThetaGfan[i] = [list(x) for x in SubdivisionsOfThetaByQ[i].rays()]
    ThetaSubdividedGfan[i] = Polyhedron(rays=raysOfThetaGfan[i])
    pointThetaGfan[i]= list(2*sum([vector(r) for r in raysOfThetaGfan[i]]))


# We record the Rays of the Subdivisions:

# for i in CellIndices:
#     print 'Cell ' + str(i)
#     for j in SubdivisionsOfThetaByQ[i].rays():
#     	print j
#     print SubdivisionsOfThetaByQ[i].lines()
#     print ''

# Cell 0 -- Theta 2
# A ray in the direction (0, 0, 0, 0, 0, -1)
# A ray in the direction (0, 0, -1, 0, 0, 0)
# A ray in the direction (0, 0, -2, -1, -2, -2)
# A ray in the direction (0, -1, -1, -1, -1, -1)
# A ray in the direction (0, 0, -1, -1, -1, -1)
# (A line in the direction (1, 1, 1, 1, 1, 1),)

# Cell 1 -- Theta 0
# A ray in the direction (0, 0, 0, 0, 0, -1)
# A ray in the direction (0, 0, 0, 0, -1, -1)
# A ray in the direction (0, 0, -2, -1, -2, -2)
# A ray in the direction (0, 0, -1, -1, -1, -1)
# A ray in the direction (0, -1, -1, -1, -1, -1)
# (A line in the direction (1, 1, 1, 1, 1, 1),)

# Cell 2 -- Theta 1
# A ray in the direction (0, -1, -1, -1, -1, -1)
# A ray in the direction (0, 0, 0, 0, 0, -1)
# A ray in the direction (0, 0, 0, 0, -1, -1)
# A ray in the direction (0, 0, -1, 0, 0, 0)
# A ray in the direction (0, 0, -2, -1, -2, -2)
# (A line in the direction (1, 1, 1, 1, 1, 1),)



# We record the H-representation of each piece: this gives us the inequalities satisfied by the three edges of the Theta graph.
ThetaCellPiecesHRep = {}
for i in CellIndices:
    ThetaCellPiecesHRep[i] =  SubdivisionsOfThetaByQ[i].Hrepresentation()    

# for i in CellIndices:
#     print 'Cell ' + str(i)
#     ThetaCellPiecesHRep[i]

# Cell 0 -- Theta 2
# (An inequality (0, 0, -1, 0, 1, 0) x + 0 >= 0,
#  An inequality (1, -1, 0, 0, 0, 0) x + 0 >= 0,
#  An inequality (0, 0, 0, 0, 1, -1) x + 0 >= 0,
#  An inequality (0, 0, 0, 1, -1, 0) x + 0 >= 0,
#  An inequality (0, 1, 0, -2, 1, 0) x + 0 >= 0)
# Cell 1 -- Theta 0
# (An inequality (0, 0, -1, 1, 0, 0) x + 0 >= 0,
#  An inequality (0, 0, 0, 0, 1, -1) x + 0 >= 0,
#  An inequality (1, -1, 0, 0, 0, 0) x + 0 >= 0,
#  An inequality (0, 0, 1, 0, -1, 0) x + 0 >= 0,
#  An inequality (0, 1, 1, -2, 0, 0) x + 0 >= 0)
# Cell 2 -- Theta 1
# (An inequality (0, -1, -1, 2, 0, 0) x + 0 >= 0,
#  An inequality (0, -1, 0, 2, -1, 0) x + 0 >= 0,
#  An inequality (1, -1, 0, 0, 0, 0) x + 0 >= 0,
#  An inequality (0, 1, 0, -1, 0, 0) x + 0 >= 0,
#  An inequality (0, 0, 0, 0, 1, -1) x + 0 >= 0)


findingCones = {}
for i in CellIndices:
    testPoint = pointThetaGfan[i]
    findingCones[i]=[]
    for j in range(0,len(PolyMaximalCones.keys())):
    	# print str(j)
    	polyCone =PolyMaximalCones[j]
	if testPoint in polyCone:
	   findingCones[i].append(j)

# findingCones
# {0: [163], 1: [164], 2: [167]}


# As sample points to compute the initial forms, we take the sum of the rays on each piece.

pointThetaGfan = {i: sum([vector(raysOfThetaGfan[i][j]) for j in range(0, len(raysOfThetaGfan[i]))]) for i in CellIndices}


# We record the sample points:
# pointThetaGfan
# {0: (0, -1, -5, -3, -4, -5),
#  1: (0, -1, -4, -3, -5, -6),
#  2: (0, -1, -4, -2, -4, -5)}

# We take appropriate shifts by the all-ones vector to ensure that the entries of the points are positive.

shiftedPointThetaGfan = {}
for k in pointThetaGfan.keys():
    shiftedPointThetaGfan[k] = [pointThetaGfan[k][j]+6 for j in range(0,6)]

shiftedPointThetaGfan
# {0: [6, 5, 1, 3, 2, 1],
#  1: [6, 5, 2, 3, 1, 0],
#  2: [6, 5, 2, 4, 2, 1]}




###############################
# Computing the intersections #
###############################

Cell01 = SubdivisionsOfThetaByQ[0].intersection(SubdivisionsOfThetaByQ[1])
# print dim(Cell01)
# 5
# Cell01.Hrepresentation()    
# (An equation (0, 0, 1, 0, -1, 0) x + 0 == 0,
#  An inequality (1, -1, 0, 0, 0, 0) x + 0 >= 0,
#  An inequality (0, 0, 1, 0, 0, -1) x + 0 >= 0,
#  An inequality (0, 0, -1, 1, 0, 0) x + 0 >= 0,
#  An inequality (0, 1, 1, -2, 0, 0) x + 0 >= 0)

Cell02 =  SubdivisionsOfThetaByQ[0].intersection(SubdivisionsOfThetaByQ[2])
# print dim(Cell02)
# 5
# Cell02.Hrepresentation()    
# (An equation (0, 1, 0, -2, 1, 0) x + 0 == 0,
#  An inequality (1, -1, 0, 0, 0, 0) x + 0 >= 0,
#  An inequality (0, -1, -1, 2, 0, 0) x + 0 >= 0,
#  An inequality (0, -1, 0, 2, 0, -1) x + 0 >= 0,
#  An inequality (0, 1, 0, -1, 0, 0) x + 0 >= 0)

Cell12 =  SubdivisionsOfThetaByQ[1].intersection(SubdivisionsOfThetaByQ[2])
# print dim(Cell12)
# 5
# Cell12.Hrepresentation()    
# (An equation (0, 1, 1, -2, 0, 0) x + 0 == 0,
#  An inequality (1, -1, 0, 0, 0, 0) x + 0 >= 0,
#  An inequality (0, 1, -1, 0, 0, 0) x + 0 >= 0,
#  An inequality (0, 0, 0, 0, 1, -1) x + 0 >= 0,
#  An inequality (0, 0, 1, 0, -1, 0) x + 0 >= 0)

Cell012 =  Cell01.intersection(SubdivisionsOfThetaByQ[2])
print dim(Cell012)
4
Cell012.Hrepresentation()    
(An equation (0, 0, 1, 0, -1, 0) x + 0 == 0,
 An equation (0, 1, 1, -2, 0, 0) x + 0 == 0,
 An inequality (1, -1, 0, 0, 0, 0) x + 0 >= 0,
 An inequality (0, 0, 1, 0, 0, -1) x + 0 >= 0,
 An inequality (0, 1, -1, 0, 0, 0) x + 0 >= 0)


#####################################################
# Checking subdivision on the low dimensional cells #
#####################################################

# We check that the low dimensional cells Cell01 throught Cell012 don't get subdivided by the Grobner fan of QT4. This will allow us to pick any sample point to compute initial forms on these cells.

# We load the Groebner fan and compute the dictionary of maximal cones. This fan was computed in "Qrefinementwithx34Noa4.sage"

GFqFan = load("GrobnerFanAltQ4x34.sobj")
PolyMaximalCones = MaximalCones(GFqFan)


###########
# Cell 01 #
###########

# The following script computes the intersection of the Cell cone with all maximal cells in the Grobner fan of QT4.

CellCone = Cell01
IntersectionsQCellCone = {}
for i in range(0,len(PolyMaximalCones.keys())):
    # print str(i)
    polyCone =PolyMaximalCones[i]
    IntersectionQPolyCone = CellCone.intersection(polyCone)
    IntersectionsQCellCone[i] = IntersectionQPolyCone

set([dim(x) for x in IntersectionsQCellCone.values()])
# {1, 2, 3, 4, 5}


# We record those intersection of maximal dimension, i.e. the dimension of the Cell Cone.
maxIntersections = [i for i in IntersectionsQCellCone.keys() if dim(IntersectionsQCellCone[i])==5]

# print maxIntersections
# [163, 164]


SubdivisionsOfCellByQ = [IntersectionsQCellCone[i] for i in maxIntersections]

# print [len(SubdivisionsOfCellByQ[i].rays()) for i in range(0,len(SubdivisionsOfCellByQ))]
# [4, 4]

# print [dim(SubdivisionsOfCellByQ[i]) for i in range(0,len(SubdivisionsOfCellByQ))]
# [5, 5]

# all( [SubdivisionsOfCellByQ[i] == SubdivisionsOfCellByQ[j] for i in range(0,len(maxIntersections)-1) for j in range(i+1,len(maxIntersections))])
# True

# No subdivision needed.



###########
# Cell 02 #
###########

# The following script computes the intersection of the Cell cone with all maximal cells in the Grobner fan of QT4.

CellCone = Cell02
IntersectionsQCellCone = {}
for i in range(0,len(PolyMaximalCones.keys())):
    # print str(i)
    polyCone =PolyMaximalCones[i]
    IntersectionQPolyCone = CellCone.intersection(polyCone)
    IntersectionsQCellCone[i] = IntersectionQPolyCone

set([dim(x) for x in IntersectionsQCellCone.values()])
# {1, 2, 3, 4, 5}


# We record those intersection of maximal dimension, i.e. the dimension of the Cell Cone.
maxIntersections = [i for i in IntersectionsQCellCone.keys() if dim(IntersectionsQCellCone[i])==5]

# print maxIntersections
# [163, 167]


SubdivisionsOfCellByQ = [IntersectionsQCellCone[i] for i in maxIntersections]

# print [len(SubdivisionsOfCellByQ[i].rays()) for i in range(0,len(SubdivisionsOfCellByQ))]
# [4, 4]

# print [dim(SubdivisionsOfCellByQ[i]) for i in range(0,len(SubdivisionsOfCellByQ))]
# [5, 5]

# all( [SubdivisionsOfCellByQ[i] == SubdivisionsOfCellByQ[j] for i in range(0,len(maxIntersections)-1) for j in range(i+1,len(maxIntersections))])
# True

# No subdivision needed.


###########
# Cell 12 #
###########

# The following script computes the intersection of the Cell cone with all maximal cells in the Grobner fan of QT4.

CellCone = Cell12
IntersectionsQCellCone = {}
for i in range(0,len(PolyMaximalCones.keys())):
    # print str(i)
    polyCone =PolyMaximalCones[i]
    IntersectionQPolyCone = CellCone.intersection(polyCone)
    IntersectionsQCellCone[i] = IntersectionQPolyCone

set([dim(x) for x in IntersectionsQCellCone.values()])
# {1, 2, 3, 4, 5}


# We record those intersection of maximal dimension, i.e. the dimension of the Cell Cone.
maxIntersections = [i for i in IntersectionsQCellCone.keys() if dim(IntersectionsQCellCone[i])==5]

# print maxIntersections
# [164, 167]


SubdivisionsOfCellByQ = [IntersectionsQCellCone[i] for i in maxIntersections]

# print [len(SubdivisionsOfCellByQ[i].rays()) for i in range(0,len(SubdivisionsOfCellByQ))]
# [4, 4]

# print [dim(SubdivisionsOfCellByQ[i]) for i in range(0,len(SubdivisionsOfCellByQ))]
# [5, 5]

# all( [SubdivisionsOfCellByQ[i] == SubdivisionsOfCellByQ[j] for i in range(0,len(maxIntersections)-1) for j in range(i+1,len(maxIntersections))])
# True

# No subdivision needed.


############
# Cell 012 #
############

# The following script computes the intersection of the Cell cone with all maximal cells in the Grobner fan of QT4.

CellCone = Cell012
IntersectionsQCellCone = {}
for i in range(0,len(PolyMaximalCones.keys())):
    # print str(i)
    polyCone =PolyMaximalCones[i]
    IntersectionQPolyCone = CellCone.intersection(polyCone)
    IntersectionsQCellCone[i] = IntersectionQPolyCone

set([dim(x) for x in IntersectionsQCellCone.values()])
# {1, 2, 3, 4}


# We record those intersection of maximal dimension, i.e. the dimension of the Cell Cone.
maxIntersections = [i for i in IntersectionsQCellCone.keys() if dim(IntersectionsQCellCone[i])==4]

# print maxIntersections
# [163, 164, 167]


SubdivisionsOfCellByQ = [IntersectionsQCellCone[i] for i in maxIntersections]

# print [len(SubdivisionsOfCellByQ[i].rays()) for i in range(0,len(SubdivisionsOfCellByQ))]
# [3, 3, 3]

# print [dim(SubdivisionsOfCellByQ[i]) for i in range(0,len(SubdivisionsOfCellByQ))]
# [4, 4, 4]

# all( [SubdivisionsOfCellByQ[i] == SubdivisionsOfCellByQ[j] for i in range(0,len(maxIntersections)-1) for j in range(i+1,len(maxIntersections))])
# True

# No subdivision needed.



##################################################
# Computing sample points for the  intersections #
##################################################

subdivisionsOfTypeCone = Cell01
samplePoint = findRaysPiecesSamplePoint([subdivisionsOfTypeCone])[-1]
pointCell01 = makeNegativeAndPrimitive(samplePoint[0])
# print pointCell01
# [0, -1, -4, -3, -4, -5]


subdivisionsOfTypeCone = Cell02
samplePoint = findRaysPiecesSamplePoint([subdivisionsOfTypeCone])[-1]
pointCell02 = makeNegativeAndPrimitive(samplePoint[0])
# print pointCell02
# [0, -1, -4, -2, -3, -4]


subdivisionsOfTypeCone = Cell12
samplePoint = findRaysPiecesSamplePoint([subdivisionsOfTypeCone])[-1]
pointCell12 = makeNegativeAndPrimitive(samplePoint[0])
# print pointCell12
# [0, -1, -3, -2, -4, -5]

subdivisionsOfTypeCone = Cell012
samplePoint = findRaysPiecesSamplePoint([subdivisionsOfTypeCone])[-1]
pointCell012 = makeNegativeAndPrimitive(samplePoint[0])
# print pointCell012
# [0, -1, -3, -2, -3, -4]