##########################################################
# Computing potential heights of yz-projections of Theta #
##########################################################

# In order to find all possible combinatorial types of resulting modified 3D-tropical curves with Theta graphs as minimal skeletons, we must compute the Groebner fans of all the coefficients. Afterwards, we intersect the fan with the TypeIICone, check for subdivisions (in "ComputingSubdivisionsCoefficientsyz.sage"), and compute the leading terms for each piece in the subdivision.
# The following script computes the Groebner fans of the coefficients of the equation g(y,z) defining the yz-projection of the hyperelliptic genus 2 curve with Theta graph as minimal skeleton.

# The input polynomials are obtained from "projectionyz.sing". The variables (y,z) are substituted for (x,y).

# Our first task is to rewrite the coefficients (by replacing a3 by a4 + a34, where x34 = a3-a4.)

# We load the auxiliary files for computation of Leading terms for weight orders and Grobner fans:

load("macroslocalJInvariantComputationsAndPrimesWithBadReduction.sage")

#########################################
# Generate the ambient ring and Weights #
#########################################

# We take two of the branch points to be 0 and 'Infinity', so we need only consider 4 variables b5, b4, b3, and b2, giving the branch points: a5 = -b5^2, a4 = b4^2, a3 = b3^2, a2 = b2^2.

As = [var("b5"),var("b4"),var("b34"), var("b2")]
R = PolynomialRing(QQ,As)
K = FractionField(R); 

# Our first task is to rewrite the coefficients, by replacing a3 by a4 + a34, where a34 = a3-a4. Since the initials of a3 and a4 agree, we may assume the same happens for b3 and b4 (otherwise, we would have replaced b4 with -b4). Then, a3-a4 = (b3-b4)(b3+b4). If we write b3-b4 = b34, we get a34 = b34*(b4+2b4). 

oldAs = [var("b5"),var("b4"),var("b3"),var("b34"), var("b2")]
oldR = PolynomialRing(QQ,oldAs)
oldK = FractionField(oldR); 

# K
# Fraction Field of Multivariate Polynomial Ring in a5, a4, a3, x34, a2 over Rational Field



################
# Type II Cone #
################

# Recall: the variables of R are ordered a5,a3,a34,a2

# The following script constructs the Type II cone, with respect to the variables a5,a4,a34,a2:
# Its defining inequalities are:
# a5 - a4 >= 0
# a4 - a34>=0
# a4 - a2 >= 0
TypeIICone = Polyhedron(ieqs=[(0,1,-1,0,0),(0,0,1,-1,0),(0,0,1,0,-1)])

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

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

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

# dim(TypeIICone)
# 4


#################################
# Coefficients of yz-projection #
#################################

# We load the relevant coefficients of the yz-projection, computed with Singular in "projectionyz.sing". WARNING: we label our variables as x, y (refering to y, z).


load("coefficientsProjectionyz.sage")

# We turn the expressions into polynomials and perform the substitution b3 = b3+b4

mx4 = R(expand(oldR(mx4).substitute(b3=b34+b4)))
mx3y = R(expand(oldR(mx3y).substitute(b3=b34+b4)))
mx2y2 = R(expand(oldR(mx2y2).substitute(b3=b34+b4)))
mxy3 = R(expand(oldR(mxy3).substitute(b3=b34+b4)))
my4 = R(expand(oldR(my4).substitute(b3=b34+b4)))
mx3 = R(expand(oldR(mx3).substitute(b3=b34+b4)))
mx2y = R(expand(oldR(mx2y).substitute(b3=b34+b4)))
mxy2 = R(expand(oldR(mxy2).substitute(b3=b34+b4)))
my3 = R(expand(oldR(my3).substitute(b3=b34+b4)))
mx2 = R(expand(oldR(mx2).substitute(b3=b34+b4)))
mxy = R(expand(oldR(mxy).substitute(b3=b34+b4)))
my2 = R(expand(oldR(my2).substitute(b3=b34+b4)))
mx = R(expand(oldR(mx).substitute(b3=b34+b4)))
my = R(expand(oldR(my).substitute(b3=b34+b4)))

# We take the product of all the above fans to compute the common refinement of all Grobner cones:

allm = mx4*mx3y*mx2y2*mxy3*my4*mx3*mx2y*mxy2*my3*mx2*mxy*my2*mx*my

# # We factor each one of the above polynomials
# factor(mx4)
# b5 * (2*b5^2*b34^2 + 3*b4^2*b34^2 + 3*b4*b34^3 + b34^4 + 2*b5^2*b2^2 - b4^2*b2^2 - b4*b34*b2^2 + b2^4)
# factor(mx3y)
# (-2) * b5 * (b5^2*b4^2 + b5^2*b4*b34 + 2*b5^2*b34^2 + 6*b4^2*b34^2 + 6*b4*b34^3 + 2*b34^4 + 2*b5^2*b2^2 - 2*b4^2*b2^2 - 2*b4*b34*b2^2 + 2*b2^4)
# factor(mx2y2)
# b5 * (4*b5^4 + 5*b5^2*b4^2 + 5*b5^2*b4*b34 + 2*b5^2*b34^2 + 18*b4^2*b34^2 + 18*b4*b34^3 + 6*b34^4 + 2*b5^2*b2^2 - 6*b4^2*b2^2 - 6*b4*b34*b2^2 + 6*b2^4)
# factor(mxy3)
# (-4) * b5 * (b5^4 + b5^2*b4^2 + b5^2*b4*b34 + 3*b4^2*b34^2 + 3*b4*b34^3 + b34^4 - b4^2*b2^2 - b4*b34*b2^2 + b2^4)
# factor(my4)
# b5 * (b5^4 + b5^2*b4^2 + b5^2*b4*b34 + 3*b4^2*b34^2 + 3*b4*b34^3 + b34^4 - b4^2*b2^2 - b4*b34*b2^2 + b2^4)
# factor(mx3)
# b5^2 * (-b5^2*b4^4*b34^2 - b4^6*b34^2 - 2*b5^2*b4^3*b34^3 - 3*b4^5*b34^3 + b5^4*b34^4 - 3*b4^4*b34^4 + b5^2*b4*b34^5 - b4^3*b34^5 - b5^2*b4^4*b2^2 - 2*b5^2*b4^3*b34*b2^2 + 2*b5^4*b34^2*b2^2 + 2*b5^2*b4^2*b34^2*b2^2 - 3*b4^4*b34^2*b2^2 + 3*b5^2*b4*b34^3*b2^2 - 6*b4^3*b34^3*b2^2 - 4*b4^2*b34^4*b2^2 - b4*b34^5*b2^2 + b5^4*b2^4 + b5^2*b4^2*b2^4 + b5^2*b4*b34*b2^4 + 3*b4^2*b34^2*b2^4 + 3*b4*b34^3*b2^4 + b34^4*b2^4)
# factor(mx2y)
# (-1) * b5^2 * (6*b5^4*b4^2*b34^2 + 5*b5^2*b4^4*b34^2 - 3*b4^6*b34^2 + 6*b5^4*b4*b34^3 + 10*b5^2*b4^3*b34^3 - 9*b4^5*b34^3 + 3*b5^4*b34^4 + 8*b5^2*b4^2*b34^4 - 9*b4^4*b34^4 + 3*b5^2*b4*b34^5 - 3*b4^3*b34^5 - 2*b5^4*b4^2*b2^2 - 3*b5^2*b4^4*b2^2 - 2*b5^4*b4*b34*b2^2 - 6*b5^2*b4^3*b34*b2^2 + 2*b5^4*b34^2*b2^2 - 9*b4^4*b34^2*b2^2 + 3*b5^2*b4*b34^3*b2^2 - 18*b4^3*b34^3*b2^2 - 12*b4^2*b34^4*b2^2 - 3*b4*b34^5*b2^2 + 3*b5^4*b2^4 + 3*b5^2*b4^2*b2^4 + 3*b5^2*b4*b34*b2^4 + 9*b4^2*b34^2*b2^4 + 9*b4*b34^3*b2^4 + 3*b34^4*b2^4)
# factor(mxy2)
# (3) * b5^2 * (3*b5^4*b4^2*b34^2 + 3*b5^2*b4^4*b34^2 - b4^6*b34^2 + 3*b5^4*b4*b34^3 + 6*b5^2*b4^3*b34^3 - 3*b4^5*b34^3 + b5^4*b34^4 + 4*b5^2*b4^2*b34^4 - 3*b4^4*b34^4 + b5^2*b4*b34^5 - b4^3*b34^5 - b5^4*b4^2*b2^2 - b5^2*b4^4*b2^2 - b5^4*b4*b34*b2^2 - 2*b5^2*b4^3*b34*b2^2 - b5^2*b4^2*b34^2*b2^2 - 3*b4^4*b34^2*b2^2 - 6*b4^3*b34^3*b2^2 - 4*b4^2*b34^4*b2^2 - b4*b34^5*b2^2 + b5^4*b2^4 + b5^2*b4^2*b2^4 + b5^2*b4*b34*b2^4 + 3*b4^2*b34^2*b2^4 + 3*b4*b34^3*b2^4 + b34^4*b2^4)
# factor(my3)
# (-1) * b5^2 * (3*b5^4*b4^2*b34^2 + 3*b5^2*b4^4*b34^2 - b4^6*b34^2 + 3*b5^4*b4*b34^3 + 6*b5^2*b4^3*b34^3 - 3*b4^5*b34^3 + b5^4*b34^4 + 4*b5^2*b4^2*b34^4 - 3*b4^4*b34^4 + b5^2*b4*b34^5 - b4^3*b34^5 - b5^4*b4^2*b2^2 - b5^2*b4^4*b2^2 - b5^4*b4*b34*b2^2 - 2*b5^2*b4^3*b34*b2^2 - b5^2*b4^2*b34^2*b2^2 - 3*b4^4*b34^2*b2^2 - 6*b4^3*b34^3*b2^2 - 4*b4^2*b34^4*b2^2 - b4*b34^5*b2^2 + b5^4*b2^4 + b5^2*b4^2*b2^4 + b5^2*b4*b34*b2^4 + 3*b4^2*b34^2*b2^4 + 3*b4*b34^3*b2^4 + b34^4*b2^4)
# factor(mx2)
# b2^2 * b34^2 * b5^3 * (-4*b5^4*b4^4 - 4*b5^2*b4^6 + b4^8 - 8*b5^4*b4^3*b34 - 12*b5^2*b4^5*b34 + 4*b4^7*b34 - 5*b5^4*b4^2*b34^2 - 13*b5^2*b4^4*b34^2 + 6*b4^6*b34^2 - b5^4*b4*b34^3 - 6*b5^2*b4^3*b34^3 + 4*b4^5*b34^3 - b5^2*b4^2*b34^4 + b4^4*b34^4 + 3*b5^4*b4^2*b2^2 + 3*b5^2*b4^4*b2^2 - b4^6*b2^2 + 3*b5^4*b4*b34*b2^2 + 6*b5^2*b4^3*b34*b2^2 - 3*b4^5*b34*b2^2 + b5^4*b34^2*b2^2 + 4*b5^2*b4^2*b34^2*b2^2 - 3*b4^4*b34^2*b2^2 + b5^2*b4*b34^3*b2^2 - b4^3*b34^3*b2^2)
# factor(mxy)
# (-2) * b34^2 * b5^3 * (-b5^4*b4^6 - b5^2*b4^8 - 3*b5^4*b4^5*b34 - 4*b5^2*b4^7*b34 - 3*b5^4*b4^4*b34^2 - 6*b5^2*b4^6*b34^2 - b5^4*b4^3*b34^3 - 4*b5^2*b4^5*b34^3 - b5^2*b4^4*b34^4 - 3*b5^4*b4^4*b2^2 - 3*b5^2*b4^6*b2^2 + b4^8*b2^2 - 6*b5^4*b4^3*b34*b2^2 - 9*b5^2*b4^5*b34*b2^2 + 4*b4^7*b34*b2^2 - 4*b5^4*b4^2*b34^2*b2^2 - 10*b5^2*b4^4*b34^2*b2^2 + 6*b4^6*b34^2*b2^2 - b5^4*b4*b34^3*b2^2 - 5*b5^2*b4^3*b34^3*b2^2 + 4*b4^5*b34^3*b2^2 - b5^2*b4^2*b34^4*b2^2 + b4^4*b34^4*b2^2 + 3*b5^4*b4^2*b2^4 + 3*b5^2*b4^4*b2^4 - b4^6*b2^4 + 3*b5^4*b4*b34*b2^4 + 6*b5^2*b4^3*b34*b2^4 - 3*b4^5*b34*b2^4 + b5^4*b34^2*b2^4 + 4*b5^2*b4^2*b34^2*b2^4 - 3*b4^4*b34^2*b2^4 + b5^2*b4*b34^3*b2^4 - b4^3*b34^3*b2^4)
# factor(my2)
# b34^2 * b5^3 * (-b5^4*b4^6 - b5^2*b4^8 - 3*b5^4*b4^5*b34 - 4*b5^2*b4^7*b34 - 3*b5^4*b4^4*b34^2 - 6*b5^2*b4^6*b34^2 - b5^4*b4^3*b34^3 - 4*b5^2*b4^5*b34^3 - b5^2*b4^4*b34^4 - 3*b5^4*b4^4*b2^2 - 3*b5^2*b4^6*b2^2 + b4^8*b2^2 - 6*b5^4*b4^3*b34*b2^2 - 9*b5^2*b4^5*b34*b2^2 + 4*b4^7*b34*b2^2 - 4*b5^4*b4^2*b34^2*b2^2 - 10*b5^2*b4^4*b34^2*b2^2 + 6*b4^6*b34^2*b2^2 - b5^4*b4*b34^3*b2^2 - 5*b5^2*b4^3*b34^3*b2^2 + 4*b4^5*b34^3*b2^2 - b5^2*b4^2*b34^4*b2^2 + b4^4*b34^4*b2^2 + 3*b5^4*b4^2*b2^4 + 3*b5^2*b4^4*b2^4 - b4^6*b2^4 + 3*b5^4*b4*b34*b2^4 + 6*b5^2*b4^3*b34*b2^4 - 3*b4^5*b34*b2^4 + b5^4*b34^2*b2^4 + 4*b5^2*b4^2*b34^2*b2^4 - 3*b4^4*b34^2*b2^4 + b5^2*b4*b34^3*b2^4 - b4^3*b34^3*b2^4)
# factor(mx)
# (-1) * b2^2 * b34^2 * b4^3 * (b4 + b34)^3 * b5^6 * (-b4^2 - b4*b34 + b2^2) * (b5^2 + b4^2 + b4*b34)
# factor(my)
# b2^2 * b34^2 * b4^3 * (b4 + b34)^3 * b5^6 * (-b4^2 - b4*b34 + b2^2) * (b5^2 + b4^2 + b4*b34)


###############################
# Computation of Grobner fans #
###############################

# We load the relevant dictionary of maximal cones precomputed in "GrobnerConeComputationsCoefficientyz.sage" corresponding to the common refinement of all the Groebner fans of the coefficients.

PolyMaximalConesgfmallDict = load("coefficientsyz/PolyMaximalConesgfCommonRefinementDictyz.sobj")

###############################
# Computation of subdivisions #
###############################

# Our calculations are valid when the residue field has characteristic 0.

# We confirm that indeed the TypeIICone is subdivided by the Grobner fan of allm.

# PolyMaximalConesAllm = PolyMaximalConesgfmallDict
# checkNoSubdivisions(TypeIICone,PolyMaximalConesAllm)
# ([31, 32, 33, 34], False)

# Next, we compute the induced subdivisions:

Subdivisions = computeSubdivisionsOfTypeConeByGrobnerFan(TypeIICone,PolyMaximalConesAllm)

# # We record the output
# Subdivisions
# ([31, 32, 33, 34],
#  [A 4-dimensional polyhedron in QQ^4 defined as the convex hull of 1 vertex, 3 rays, 1 line,
#   A 4-dimensional polyhedron in QQ^4 defined as the convex hull of 1 vertex, 3 rays, 1 line,
#   A 4-dimensional polyhedron in QQ^4 defined as the convex hull of 1 vertex, 3 rays, 1 line,
#   A 4-dimensional polyhedron in QQ^4 defined as the convex hull of 1 vertex, 3 rays, 1 line])


# We compute the pieces of the subdivision and sample points to compute the leading terms:

PiecesTypeIICone = findRaysPiecesSamplePoint(Subdivisions[1])
# # We record the output
# PiecesTypeIICone
# ({0: [[1, 1, 0, 1], [2, 1, 0, 0], [1, 1, 0, 0]],
#   1: [[1, 1, 1, 0], [2, 1, 0, 0], [1, 1, 0, 0]],
#   2: [[0, 0, -1, 0], [0, -1, -2, -2], [0, -1, -1, -1]],
#   3: [[0, 0, 0, -1], [0, -1, -1, -1], [0, -1, -2, -2]]},
#  {0: A 4-dimensional polyhedron in QQ^4 defined as the convex hull of 1 vertex, 3 rays, 1 line,
#   1: A 4-dimensional polyhedron in QQ^4 defined as the convex hull of 1 vertex, 3 rays, 1 line,
#   2: A 4-dimensional polyhedron in QQ^4 defined as the convex hull of 1 vertex, 3 rays, 1 line,
#   3: A 4-dimensional polyhedron in QQ^4 defined as the convex hull of 1 vertex, 3 rays, 1 line},
#  {0: [8, 6, 0, 2], 1: [8, 6, 2, 0], 2: [0, -4, -8, -6], 3: [0, -4, -6, -8]})
 

# We record the sample points to compute the leading terms:

allWeightsMxl = PiecesTypeIICone[-1]


###########################
# Lower dimensional Cones #
###########################

Subdivisions3 = computeFacets(Subdivisions)

# # We record the output
# Subdivisions3
# ([(0, 1), (1, 3), (2, 3), (0, 2)],
#  [A 3-dimensional polyhedron in QQ^4 defined as the convex hull of 1 vertex, 2 rays, 1 line,
#   A 3-dimensional polyhedron in QQ^4 defined as the convex hull of 1 vertex, 2 rays, 1 line,
#   A 3-dimensional polyhedron in QQ^4 defined as the convex hull of 1 vertex, 2 rays, 1 line,
#   A 3-dimensional polyhedron in QQ^4 defined as the convex hull of 1 vertex, 2 rays, 1 line])

PiecesTypeIICone3 = findRaysPiecesSamplePoint(Subdivisions3[1])

# # We record the output:
# PiecesTypeIICone3
# ({0: [[2, 1, 0, 0], [1, 1, 0, 0]],
#   1: [[1, 1, 1, 0], [2, 1, 0, 0]],
#   2: [[0, -1, -2, -2], [0, -1, -1, -1]],
#   3: [[1, 1, 0, 1], [2, 1, 0, 0]]},
#  {0: A 3-dimensional polyhedron in QQ^4 defined as the convex hull of 1 vertex, 2 rays, 1 line,
#   1: A 3-dimensional polyhedron in QQ^4 defined as the convex hull of 1 vertex, 2 rays, 1 line,
#   2: A 3-dimensional polyhedron in QQ^4 defined as the convex hull of 1 vertex, 2 rays, 1 line,
#   3: A 3-dimensional polyhedron in QQ^4 defined as the convex hull of 1 vertex, 2 rays, 1 line},
#  {0: [6, 4, 0, 0], 1: [6, 4, 2, 0], 2: [0, -4, -6, -6], 3: [6, 4, 0, 2]})


allWeightsFacets = PiecesTypeIICone3[-1]
		 
# Since all the rays except for one ([17, 10, 3, 0]) lie in the the boundary of the TypeIICone, we need only consider one ray for our Leading term computations:

allWeightsRays = {0: [0,-1,-2,-2]}

###################################
# Simplify Weight vectors to Test #
###################################

allWeights = allWeightsMxl.values() + allWeightsFacets.values() + allWeightsRays.values()

# We compute primitive vectors with non-positive coefficients equivalent to each one of the sample vectors above, by substracting appropriate multiplies of the all-ones vector  and taking a primitive vector thereafter. This is implemented in the function "makeNegativeAndPrimitive"

allWeightsSimplified = []
for k in range(len(allWeights)):
    allWeightsSimplified.append(makeNegativeAndPrimitive(allWeights[k]))

# # We record the output vectors: we will use these to compute examples with Singular
# allWeightsSimplified
# [[0, -1, -4, -3],
#  [0, -1, -3, -4],
#  [0, -2, -4, -3],
#  [0, -2, -3, -4],
#  [0, -1, -3, -3],
#  [0, -1, -2, -3],
#  [0, -2, -3, -3],
#  [0, -1, -3, -2],
#  [0, -1, -2, -2]]



################################
# Computation of Leading Terms #
################################

leadmx4 = {i: LTFromWeights(mx4,allWeightsMxl[i],R) for i in range(len(allWeightsMxl))}
leadmx3y = {i: LTFromWeights(mx3y,allWeightsMxl[i],R) for i in range(len(allWeightsMxl))}
leadmx2y2 = {i: LTFromWeights(mx2y2,allWeightsMxl[i],R) for i in range(len(allWeightsMxl))}
leadmxy3 = {i: LTFromWeights(mxy3,allWeightsMxl[i],R) for i in range(len(allWeightsMxl))}
leadmy4 = {i: LTFromWeights(my4,allWeightsMxl[i],R) for i in range(len(allWeightsMxl))}
leadmx3 = {i: LTFromWeights(mx3,allWeightsMxl[i],R) for i in range(len(allWeightsMxl))}
leadmx2y = {i: LTFromWeights(mx2y,allWeightsMxl[i],R) for i in range(len(allWeightsMxl))}
leadmxy2 = {i: LTFromWeights(mxy2,allWeightsMxl[i],R) for i in range(len(allWeightsMxl))}
leadmy3 = {i: LTFromWeights(my3,allWeightsMxl[i],R) for i in range(len(allWeightsMxl))}
leadmx2 = {i: LTFromWeights(mx2,allWeightsMxl[i],R) for i in range(len(allWeightsMxl))}
leadmxy = {i: LTFromWeights(mxy,allWeightsMxl[i],R) for i in range(len(allWeightsMxl))}
leadmy2 = {i: LTFromWeights(my2,allWeightsMxl[i],R) for i in range(len(allWeightsMxl))}
leadmx = {i: LTFromWeights(mx,allWeightsMxl[i],R) for i in range(len(allWeightsMxl))}
leadmy = {i: LTFromWeights(my,allWeightsMxl[i],R) for i in range(len(allWeightsMxl))}


# # We record the output:
# leadmx4
# {0: 2*b5^3*b2^2, 1: 2*b5^3*b34^2, 2: 2*b5^3*b2^2, 3: 2*b5^3*b34^2}
# leadmx3y
# {0: -2*b5^3*b4^2, 1: -2*b5^3*b4^2, 2: -2*b5^3*b4^2, 3: -2*b5^3*b4^2}
# leadmx2y2
# {0: 4*b5^5, 1: 4*b5^5, 2: 4*b5^5, 3: 4*b5^5}
# leadmxy3
# {0: -4*b5^5, 1: -4*b5^5, 2: -4*b5^5, 3: -4*b5^5}
# leadmy4
# {0: b5^5, 1: b5^5, 2: b5^5, 3: b5^5}
# leadmx3
# {0: -b5^4*b4^4*b2^2, 1: -b5^4*b4^4*b34^2, 2: b5^6*b2^4, 3: b5^6*b34^4}
# leadmx2y
# {0: 2*b5^6*b4^2*b2^2,
#  1: -6*b5^6*b4^2*b34^2,
#  2: 2*b5^6*b4^2*b2^2,
#  3: -6*b5^6*b4^2*b34^2}
# leadmxy2
# {0: -3*b5^6*b4^2*b2^2,
#  1: 9*b5^6*b4^2*b34^2,
#  2: -3*b5^6*b4^2*b2^2,
#  3: 9*b5^6*b4^2*b34^2}
# leadmy3
# {0: b5^6*b4^2*b2^2,
#  1: -3*b5^6*b4^2*b34^2,
#  2: b5^6*b4^2*b2^2,
#  3: -3*b5^6*b4^2*b34^2}
# leadmx2
# {0: -4*b5^7*b4^4*b34^2*b2^2,
#  1: -4*b5^7*b4^4*b34^2*b2^2,
#  2: -4*b5^7*b4^4*b34^2*b2^2,
#  3: -4*b5^7*b4^4*b34^2*b2^2}
# leadmxy
# {0: 2*b5^7*b4^6*b34^2,
#  1: 2*b5^7*b4^6*b34^2,
#  2: 2*b5^7*b4^6*b34^2,
#  3: 2*b5^7*b4^6*b34^2}
# leadmy2
# {0: -b5^7*b4^6*b34^2,
#  1: -b5^7*b4^6*b34^2,
#  2: -b5^7*b4^6*b34^2,
#  3: -b5^7*b4^6*b34^2}
# leadmx
# {0: b5^8*b4^8*b34^2*b2^2,
#  1: b5^8*b4^8*b34^2*b2^2,
#  2: b5^8*b4^8*b34^2*b2^2,
#  3: b5^8*b4^8*b34^2*b2^2}
# leadmy
# {0: -b5^8*b4^8*b34^2*b2^2,
#  1: -b5^8*b4^8*b34^2*b2^2,
#  2: -b5^8*b4^8*b34^2*b2^2,
#  3: -b5^8*b4^8*b34^2*b2^2}



#######################################
# Computation of Leading Terms Facets #
#######################################

# Labels of facets = Subdivisions3[0] = [(0, 1), (1, 3), (2, 3), (0, 2)]

leadFacetsmx4 = {i: factor(LTFromWeights(mx4,allWeightsFacets[i],R)) for i in range(len(allWeightsFacets))}
leadFacetsmx3y = {i: factor(LTFromWeights(mx3y,allWeightsFacets[i],R)) for i in range(len(allWeightsFacets))}
leadFacetsmx2y2 = {i: factor(LTFromWeights(mx2y2,allWeightsFacets[i],R)) for i in range(len(allWeightsFacets))}
leadFacetsmxy3 = {i: factor(LTFromWeights(mxy3,allWeightsFacets[i],R)) for i in range(len(allWeightsFacets))}
leadFacetsmy4 = {i: factor(LTFromWeights(my4,allWeightsFacets[i],R)) for i in range(len(allWeightsFacets))}
leadFacetsmx3 = {i: factor(LTFromWeights(mx3,allWeightsFacets[i],R)) for i in range(len(allWeightsFacets))}
leadFacetsmx2y = {i: factor(LTFromWeights(mx2y,allWeightsFacets[i],R)) for i in range(len(allWeightsFacets))}
leadFacetsmxy2 = {i: factor(LTFromWeights(mxy2,allWeightsFacets[i],R)) for i in range(len(allWeightsFacets))}
leadFacetsmy3 = {i: factor(LTFromWeights(my3,allWeightsFacets[i],R)) for i in range(len(allWeightsFacets))}
leadFacetsmx2 = {i: factor(LTFromWeights(mx2,allWeightsFacets[i],R)) for i in range(len(allWeightsFacets))}
leadFacetsmxy = {i: factor(LTFromWeights(mxy,allWeightsFacets[i],R)) for i in range(len(allWeightsFacets))}
leadFacetsmy2 = {i: factor(LTFromWeights(my2,allWeightsFacets[i],R)) for i in range(len(allWeightsFacets))}
leadFacetsmx = {i: factor(LTFromWeights(mx,allWeightsFacets[i],R)) for i in range(len(allWeightsFacets))}
leadFacetsmy = {i: factor(LTFromWeights(my,allWeightsFacets[i],R)) for i in range(len(allWeightsFacets))}


# We record the outputs:
# leadFacetsmx4
# {0: (2) * b5^3 * (b34^2 + b2^2),
#  1: (2) * b34^2 * b5^3,
#  2: (2) * b5^3 * (b34^2 + b2^2),
#  3: (2) * b2^2 * b5^3}
# leadFacetsmx3y
# {0: (-2) * b4^2 * b5^3,
#  1: (-2) * b4^2 * b5^3,
#  2: (-2) * b4^2 * b5^3,
#  3: (-2) * b4^2 * b5^3}
# leadFacetsmx2y2
# {0: (4) * b5^5, 1: (4) * b5^5, 2: (4) * b5^5, 3: (4) * b5^5}
# leadFacetsmxy3
# {0: (-4) * b5^5, 1: (-4) * b5^5, 2: (-4) * b5^5, 3: (-4) * b5^5}
# leadFacetsmy4
# {0: b5^5, 1: b5^5, 2: b5^5, 3: b5^5}
# leadFacetsmx3
# {0: (-1) * b4^4 * b5^4 * (b34^2 + b2^2),
#  1: b34^2 * b5^4 * (-b4^2 + b5*b34) * (b4^2 + b5*b34),
#  2: b5^6 * (b34^2 + b2^2)^2,
#  3: b2^2 * b5^4 * (-b4^2 + b5*b2) * (b4^2 + b5*b2)}
# leadFacetsmx2y
# {0: (-2) * b4^2 * b5^6 * (3*b34^2 - b2^2),
#  1: (-6) * b34^2 * b4^2 * b5^6,
#  2: (-2) * b4^2 * b5^6 * (3*b34^2 - b2^2),
#  3: (2) * b2^2 * b4^2 * b5^6}
# leadFacetsmxy2
# {0: (3) * b4^2 * b5^6 * (3*b34^2 - b2^2),
#  1: (9) * b34^2 * b4^2 * b5^6,
#  2: (3) * b4^2 * b5^6 * (3*b34^2 - b2^2),
#  3: (-3) * b2^2 * b4^2 * b5^6}
# leadFacetsmy3
# {0: (-1) * b4^2 * b5^6 * (3*b34^2 - b2^2),
#  1: (-3) * b34^2 * b4^2 * b5^6,
#  2: (-1) * b4^2 * b5^6 * (3*b34^2 - b2^2),
#  3: b2^2 * b4^2 * b5^6}
# leadFacetsmx2
# {0: (-4) * b2^2 * b34^2 * b4^4 * b5^7,
#  1: (-4) * b2^2 * b34^2 * b4^4 * b5^7,
#  2: (-4) * b2^2 * b34^2 * b4^4 * b5^7,
#  3: (-4) * b2^2 * b34^2 * b4^4 * b5^7}
# leadFacetsmxy
# {0: (2) * b34^2 * b4^6 * b5^7,
#  1: (2) * b34^2 * b4^6 * b5^7,
#  2: (2) * b34^2 * b4^6 * b5^7,
#  3: (2) * b34^2 * b4^6 * b5^7}
# leadFacetsmy2
# {0: (-1) * b34^2 * b4^6 * b5^7,
#  1: (-1) * b34^2 * b4^6 * b5^7,
#  2: (-1) * b34^2 * b4^6 * b5^7,
#  3: (-1) * b34^2 * b4^6 * b5^7}
# leadFacetsmx
# {0: b2^2 * b34^2 * b4^8 * b5^8,
#  1: b2^2 * b34^2 * b4^8 * b5^8,
#  2: b2^2 * b34^2 * b4^8 * b5^8,
#  3: b2^2 * b34^2 * b4^8 * b5^8}
# leadFacetsmy
# {0: (-1) * b2^2 * b34^2 * b4^8 * b5^8,
#  1: (-1) * b2^2 * b34^2 * b4^8 * b5^8,
#  2: (-1) * b2^2 * b34^2 * b4^8 * b5^8,
#  3: (-1) * b2^2 * b34^2 * b4^8 * b5^8}




#####################################
# Computation of Leading Terms Rays #
#####################################

leadRaysmx4 = {i: factor(LTFromWeights(mx4,allWeightsRays[i],R)) for i in range(len(allWeightsRays))}
leadRaysmx3y = {i: factor(LTFromWeights(mx3y,allWeightsRays[i],R)) for i in range(len(allWeightsRays))}
leadRaysmx2y2 = {i: factor(LTFromWeights(mx2y2,allWeightsRays[i],R)) for i in range(len(allWeightsRays))}
leadRaysmxy3 = {i: factor(LTFromWeights(mxy3,allWeightsRays[i],R)) for i in range(len(allWeightsRays))}
leadRaysmy4 = {i: factor(LTFromWeights(my4,allWeightsRays[i],R)) for i in range(len(allWeightsRays))}
leadRaysmx3 = {i: factor(LTFromWeights(mx3,allWeightsRays[i],R)) for i in range(len(allWeightsRays))}
leadRaysmx2y = {i: factor(LTFromWeights(mx2y,allWeightsRays[i],R)) for i in range(len(allWeightsRays))}
leadRaysmxy2 = {i: factor(LTFromWeights(mxy2,allWeightsRays[i],R)) for i in range(len(allWeightsRays))}
leadRaysmy3 = {i: factor(LTFromWeights(my3,allWeightsRays[i],R)) for i in range(len(allWeightsRays))}
leadRaysmx2 = {i: factor(LTFromWeights(mx2,allWeightsRays[i],R)) for i in range(len(allWeightsRays))}
leadRaysmxy = {i: factor(LTFromWeights(mxy,allWeightsRays[i],R)) for i in range(len(allWeightsRays))}
leadRaysmy2 = {i: factor(LTFromWeights(my2,allWeightsRays[i],R)) for i in range(len(allWeightsRays))}
leadRaysmx = {i: factor(LTFromWeights(mx,allWeightsRays[i],R)) for i in range(len(allWeightsRays))}
leadRaysmy = {i: factor(LTFromWeights(my,allWeightsRays[i],R)) for i in range(len(allWeightsRays))}


# We record the outputs:
# leadRaysmx4[0]
# (2) * b5^3 * (b34^2 + b2^2)
# leadRaysmx4[0]
# (2) * b5^3 * (b34^2 + b2^2)
# leadRaysmx3y[0]
# (-2) * b4^2 * b5^3
# leadRaysmx2y2[0]
# (4) * b5^5
# leadRaysmxy3[0]
# (-4) * b5^5
# leadRaysmy4[0]
# b5^5
# leadRaysmx3[0]
# b5^4 * (b34^2 + b2^2) * (-b4^4 + b5^2*b34^2 + b5^2*b2^2)
# leadRaysmx2y[0]
# (-2) * b4^2 * b5^6 * (3*b34^2 - b2^2)
# leadRaysmxy2[0]
# (3) * b4^2 * b5^6 * (3*b34^2 - b2^2)
# leadRaysmy3[0]
# (-1) * b4^2 * b5^6 * (3*b34^2 - b2^2)
# leadRaysmx2[0]
# (-4) * b2^2 * b34^2 * b4^4 * b5^7
# leadRaysmxy[0]
# (2) * b34^2 * b4^6 * b5^7
# leadRaysmy2[0]
# (-1) * b34^2 * b4^6 * b5^7
# leadRaysmx[0]
# b2^2 * b34^2 * b4^8 * b5^8
# leadRaysmy[0]
# (-1) * b2^2 * b34^2 * b4^8 * b5^8


###########################
# Values of Leading Terms #
###########################

# Finally, we put together the indices giving the same leadingTerm values for each

lw = len(allWeightsMxl)
lf = len(allWeightsFacets)

numCells = lw + lf + 1
# numCells
# 9

labelCoefficients = [mx4, mx3y, mx2y2, mxy3, my4, mx3, mx2y, mxy2, my3, mx2, mxy, my2, mx,my]

allLeadingTerms ={}
for p in range(len(labelCoefficients)):
    LTsp = {}
    for k in range(lw):
    	LTsp[k] = LTFromWeights(labelCoefficients[p],allWeightsMxl[k],R)
    for k in range(lf):
    	LTsp[k+lw] = LTFromWeights(labelCoefficients[p],allWeightsFacets[k],R)
    LTsp[lw+lf]= LTFromWeights(labelCoefficients[p],allWeightsRays[0],R)
    allLeadingTerms[p] =  LTsp

# We compute the leading terms without repetitions and all the indices with the given leading terms
allLeadingTermsNoRepetitions = {}
for p in range(len(labelCoefficients)):
    norepetitions=list(set(allLeadingTerms[p].values()))
    numCells = len(allLeadingTerms[p].keys())
    indicesNoRep = {}
    for x in norepetitions:
    	indicesNoRep[x] = sorted([j for j in allLeadingTerms[p].keys() if allLeadingTerms[p][j] == x])
    # print norepetitions
    allLeadingTermsNoRepetitions[p] = {factor(x):indicesNoRep[x] for x in norepetitions}

# # We record the number of leading terms per coefficient:
# [len(allLeadingTermsNoRepetitions[p].keys()) for p in range(len(labelCoefficients))]
# [3, 1, 1, 1, 1, 9, 3, 3, 3, 1, 1, 1, 1, 1]


# # We record the possible leading terms:
# allLeadingTermsNoRepetitions
# {0: {(2) * b2^2 * b5^3: [0, 2, 7],
#   (2) * b34^2 * b5^3: [1, 3, 5],
#   (2) * b5^3 * (b34^2 + b2^2): [4, 6, 8]},
#  1: {(-2) * b4^2 * b5^3: [0, 1, 2, 3, 4, 5, 6, 7, 8]},
#  2: {(4) * b5^5: [0, 1, 2, 3, 4, 5, 6, 7, 8]},
#  3: {(-4) * b5^5: [0, 1, 2, 3, 4, 5, 6, 7, 8]},
#  4: {b5^5: [0, 1, 2, 3, 4, 5, 6, 7, 8]},
#  5: {b2^4 * b5^6: [2],
#   (-1) * b2^2 * b4^4 * b5^4: [0],
#   b2^2 * b5^4 * (-b4^2 + b5*b2) * (b4^2 + b5*b2): [7],
#   b34^4 * b5^6: [3],
#   b5^6 * (b34^2 + b2^2)^2: [6],
#   (-1) * b34^2 * b4^4 * b5^4: [1],
#   (-1) * b4^4 * b5^4 * (b34^2 + b2^2): [4],
#   b34^2 * b5^4 * (-b4^2 + b5*b34) * (b4^2 + b5*b34): [5],
#   b5^4 * (b34^2 + b2^2) * (-b4^4 + b5^2*b34^2 + b5^2*b2^2): [8]},
#  6: {(2) * b2^2 * b4^2 * b5^6: [0, 2, 7],
#   (-6) * b34^2 * b4^2 * b5^6: [1, 3, 5],
#   (-2) * b4^2 * b5^6 * (3*b34^2 - b2^2): [4, 6, 8]},
#  7: {(-3) * b2^2 * b4^2 * b5^6: [0, 2, 7],
#   (9) * b34^2 * b4^2 * b5^6: [1, 3, 5],
#   (3) * b4^2 * b5^6 * (3*b34^2 - b2^2): [4, 6, 8]},
#  8: {b2^2 * b4^2 * b5^6: [0, 2, 7],
#   (-3) * b34^2 * b4^2 * b5^6: [1, 3, 5],
#   (-1) * b4^2 * b5^6 * (3*b34^2 - b2^2): [4, 6, 8]},
#  9: {(-4) * b2^2 * b34^2 * b4^4 * b5^7: [0, 1, 2, 3, 4, 5, 6, 7, 8]},
#  10: {(2) * b34^2 * b4^6 * b5^7: [0, 1, 2, 3, 4, 5, 6, 7, 8]},
#  11: {(-1) * b34^2 * b4^6 * b5^7: [0, 1, 2, 3, 4, 5, 6, 7, 8]},
#  12: {b2^2 * b34^2 * b4^8 * b5^8: [0, 1, 2, 3, 4, 5, 6, 7, 8]},
#  13: {(-1) * b2^2 * b34^2 * b4^8 * b5^8: [0, 1, 2, 3, 4, 5, 6, 7, 8]}}