###############################################
# Common refinement of xz- and yz-projections #
###############################################

# The following script computes the Grobner fans of the coefficients of the equation g(y,z) and f(x,z) defining the yz- and xz-projection of the hyperelliptic genus 2 curve with Theta graph as minimal skeleton.

# The input polynomials are obtained from "projectionyz.sing" and "projectionyz.sing", respectively. The variables (y,z) and (x,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")


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

# The following script computes the Grobner fans of the coefficients of the equation g(x,z) defining the xz-projection of the hyperelliptic genus 2 curve with Theta graph as minimal skeleton.

# The input polynomials are obtained from "projectionxz.sing". The variables (x,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 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)



#################################
# Coefficients of xz-projection #
#################################

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

load("coefficientsProjectionMxz.sage")

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

Mx4 = R(expand(oldR(Mx4).substitute(b3=b34+b4)))
Mx3 = R(expand(oldR(Mx3).substitute(b3=b34+b4)))
Mx2 = R(expand(oldR(Mx2).substitute(b3=b34+b4)))


# We factor each one of the above polynomials
# factor(Mx4)
# (-1) * (2*b4^2 + 2*b4*b34 + b34^2 + b2^2)

# factor(Mx3)
# b4^4 + 2*b4^3*b34 - b5^2*b34^2 + b4^2*b34^2 - b5^2*b2^2 + 2*b4^2*b2^2 + 2*b4*b34*b2^2 + b34^2*b2^2

# factor(Mx2)
# b2^2 * (2*b5^2*b4^2 - b4^4 + 2*b5^2*b4*b34 - 2*b4^3*b34 + b5^2*b34^2 - b4^2*b34^2)

#####################
# Common refinement #
#####################

allm = mx4*mx3y*mx2y2*mxy3*my4*mx3*mx2y*mxy2*my3*mx2*mxy*my2*mx*my
# We take the product of all the above fans to compute the common refinement of all Grobner cones:

allmyzxz = allM*allm


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

# We load the relevant dictionary of maximal cones precomputed in "GrobnerConeComputationsCoefficientxz.sage" corresponding to the common refinement of all the Groebner fans of the coefficients of the xy- and yz-projections.

PolyMaximalConesgfmallDict = load("PolyMaximalConesgfCommonRefinementDictyzxz.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 vertexx, 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 otput 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]]



# We save the Pieces of the subdivisions as polyhedra:

# allPiecesTypeIICone = PiecesTypeIICone[1].values() +  PiecesTypeIICone3[1].values() + [Polyhedron(rays=[allWeightsSimplified[-1]], lines=[(1,1,1,1)])]

# save(allPiecesTypeIICone, "allPiecesTypeIICone.sobj")
allPiecesTypeIICone = load("allPiecesTypeIICone.sobj")


##################################
# H-representation of the pieces #
##################################

############
# Chambers #
############

# PiecesTypeIICone[1][0].Hrepresentation()
# (An inequality (1, -1, 0, 0) x + 0 >= 0,
#  An inequality (0, 0, -1, 1) x + 0 >= 0,
#  An inequality (-1, 2, 0, -1) x + 0 >= 0)
# PiecesTypeIICone[1][1].Hrepresentation()
# (An inequality (1, -1, 0, 0) x + 0 >= 0,
#  An inequality (0, 0, 1, -1) x + 0 >= 0,
#  An inequality (-1, 2, -1, 0) x + 0 >= 0)
# PiecesTypeIICone[1][2].Hrepresentation()
# (An inequality (0, 0, -1, 1) x + 0 >= 0,
#  An inequality (1, -2, 0, 1) x + 0 >= 0,
#  An inequality (0, 1, 0, -1) x + 0 >= 0)
# PiecesTypeIICone[1][3].Hrepresentation()
# (An inequality (0, 0, 1, -1) x + 0 >= 0,
#  An inequality (1, -2, 1, 0) x + 0 >= 0,
#  An inequality (0, 1, -1, 0) x + 0 >= 0)

##########
# Facets #
##########

# PiecesTypeIICone3[1][0].Hrepresentation()
# (An equation (0, 0, 1, -1) x + 0 == 0,
#  An inequality (1, -1, 0, 0) x + 0 >= 0,
#  An inequality (-1, 2, -1, 0) x + 0 >= 0)
# PiecesTypeIICone3[1][1].Hrepresentation()
# (An equation (1, -2, 1, 0) x + 0 == 0,
#  An inequality (1, 0, -1, 0) x + 0 >= 0,
#  An inequality (0, 0, 1, -1) x + 0 >= 0)
# PiecesTypeIICone3[1][2].Hrepresentation()
# (An equation (0, 0, 1, -1) x + 0 == 0,
#  An inequality (1, -2, 0, 1) x + 0 >= 0,
#  An inequality (0, 1, 0, -1) x + 0 >= 0)
# PiecesTypeIICone3[1][3].Hrepresentation()
# (An equation (1, -2, 0, 1) x + 0 == 0,
#  An inequality (1, -1, 0, 0) x + 0 >= 0,
#  An inequality (-1, 2, -1, 0) x + 0 >= 0)

#######
# Ray #
#######

# Polyhedron(rays=[allWeightsSimplified[-1]], lines=[(1,1,1,1)]).Hrepresentation()
# (An equation (1, -2, 0, 1) x + 0 == 0,
#  An equation (0, 0, 1, -1) x + 0 == 0,
#  An inequality (1, 0, 0, -1) x + 0 >= 0)



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

#################
# xz-projection #
#################

leadMx4 = {i: factor(LTFromWeights(Mx4,allWeightsSimplified[i],R)) for i in range(len(allWeightsSimplified))}
leadMx3 = {i: factor(LTFromWeights(Mx3,allWeightsSimplified[i],R)) for i in range(len(allWeightsSimplified))}
leadMx2 = {i: factor(LTFromWeights(Mx2,allWeightsSimplified[i],R)) for i in range(len(allWeightsSimplified))}

# # We record the output:
# leadMx4
# {0: (-2) * b4^2,
#  1: (-2) * b4^2,
#  2: (-2) * b4^2,
#  3: (-2) * b4^2,
#  4: (-2) * b4^2,
#  5: (-2) * b4^2,
#  6: (-2) * b4^2,
#  7: (-2) * b4^2,
#  8: (-2) * b4^2}
# leadMx3
# {0: b4^4,
#  1: b4^4,
#  2: (-1) * b2^2 * b5^2,
#  3: (-1) * b34^2 * b5^2,
#  4: b4^4,
#  5: (-1) * (-b4^2 + b5*b34) * (b4^2 + b5*b34),
#  6: (-1) * b5^2 * (b34^2 + b2^2),
#  7: (-1) * (-b4^2 + b5*b2) * (b4^2 + b5*b2),
#  8: (-1) * (-b4^4 + b5^2*b34^2 + b5^2*b2^2)}
# leadMx2
# {0: (2) * b2^2 * b4^2 * b5^2,
#  1: (2) * b2^2 * b4^2 * b5^2,
#  2: (2) * b2^2 * b4^2 * b5^2,
#  3: (2) * b2^2 * b4^2 * b5^2,
#  4: (2) * b2^2 * b4^2 * b5^2,
#  5: (2) * b2^2 * b4^2 * b5^2,
#  6: (2) * b2^2 * b4^2 * b5^2,
#  7: (2) * b2^2 * b4^2 * b5^2,
#  8: (2) * b2^2 * b4^2 * b5^2}



#################
# yz-projection #
#################

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

# # We record the output
# leadmx4
# {0: (2) * b2^2 * b5^3,
#  1: (2) * b34^2 * b5^3,
#  2: (2) * b2^2 * b5^3,
#  3: (2) * b34^2 * b5^3,
#  4: (2) * b5^3 * (b34^2 + b2^2),
#  5: (2) * b34^2 * b5^3,
#  6: (2) * b5^3 * (b34^2 + b2^2),
#  7: (2) * b2^2 * b5^3,
#  8: (2) * b5^3 * (b34^2 + b2^2)}
# leadmx3y
# {0: (-2) * b4^2 * b5^3,
#  1: (-2) * b4^2 * b5^3,
#  2: (-2) * b4^2 * b5^3,
#  3: (-2) * b4^2 * b5^3,
#  4: (-2) * b4^2 * b5^3,
#  5: (-2) * b4^2 * b5^3,
#  6: (-2) * b4^2 * b5^3,
#  7: (-2) * b4^2 * b5^3,
#  8: (-2) * b4^2 * b5^3}
# leadmx2y2
# {0: (4) * b5^5,
#  1: (4) * b5^5,
#  2: (4) * b5^5,
#  3: (4) * b5^5,
#  4: (4) * b5^5,
#  5: (4) * b5^5,
#  6: (4) * b5^5,
#  7: (4) * b5^5,
#  8: (4) * b5^5}
# leadmxy3
# {0: (-4) * b5^5,
#  1: (-4) * b5^5,
#  2: (-4) * b5^5,
#  3: (-4) * b5^5,
#  4: (-4) * b5^5,
#  5: (-4) * b5^5,
#  6: (-4) * b5^5,
#  7: (-4) * b5^5,
#  8: (-4) * b5^5}
# leadmy4
# {0: b5^5,
#  1: b5^5,
#  2: b5^5,
#  3: b5^5,
#  4: b5^5,
#  5: b5^5,
#  6: b5^5,
#  7: b5^5,
#  8: b5^5}
# leadmx3
# {0: (-1) * b2^2 * b4^4 * b5^4,
#  1: (-1) * b34^2 * b4^4 * b5^4,
#  2: b2^4 * b5^6,
#  3: b34^4 * b5^6,
#  4: (-1) * b4^4 * b5^4 * (b34^2 + b2^2),
#  5: b34^2 * b5^4 * (-b4^2 + b5*b34) * (b4^2 + b5*b34),
#  6: b5^6 * (b34^2 + b2^2)^2,
#  7: b2^2 * b5^4 * (-b4^2 + b5*b2) * (b4^2 + b5*b2),
#  8: b5^4 * (b34^2 + b2^2) * (-b4^4 + b5^2*b34^2 + b5^2*b2^2)}
# leadmx2y
# {0: (2) * b2^2 * b4^2 * b5^6,
#  1: (-6) * b34^2 * b4^2 * b5^6,
#  2: (2) * b2^2 * b4^2 * b5^6,
#  3: (-6) * b34^2 * b4^2 * b5^6,
#  4: (-2) * b4^2 * b5^6 * (3*b34^2 - b2^2),
#  5: (-6) * b34^2 * b4^2 * b5^6,
#  6: (-2) * b4^2 * b5^6 * (3*b34^2 - b2^2),
#  7: (2) * b2^2 * b4^2 * b5^6,
#  8: (-2) * b4^2 * b5^6 * (3*b34^2 - b2^2)}
# leadmxy2
# {0: (-3) * b2^2 * b4^2 * b5^6,
#  1: (9) * b34^2 * b4^2 * b5^6,
#  2: (-3) * b2^2 * b4^2 * b5^6,
#  3: (9) * b34^2 * b4^2 * b5^6,
#  4: (3) * b4^2 * b5^6 * (3*b34^2 - b2^2),
#  5: (9) * b34^2 * b4^2 * b5^6,
#  6: (3) * b4^2 * b5^6 * (3*b34^2 - b2^2),
#  7: (-3) * b2^2 * b4^2 * b5^6,
#  8: (3) * b4^2 * b5^6 * (3*b34^2 - b2^2)}
# leadmy3
# {0: b2^2 * b4^2 * b5^6,
#  1: (-3) * b34^2 * b4^2 * b5^6,
#  2: b2^2 * b4^2 * b5^6,
#  3: (-3) * b34^2 * b4^2 * b5^6,
#  4: (-1) * b4^2 * b5^6 * (3*b34^2 - b2^2),
#  5: (-3) * b34^2 * b4^2 * b5^6,
#  6: (-1) * b4^2 * b5^6 * (3*b34^2 - b2^2),
#  7: b2^2 * b4^2 * b5^6,
#  8: (-1) * b4^2 * b5^6 * (3*b34^2 - b2^2)}
# leadmx2
# {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,
#  4: (-4) * b2^2 * b34^2 * b4^4 * b5^7,
#  5: (-4) * b2^2 * b34^2 * b4^4 * b5^7,
#  6: (-4) * b2^2 * b34^2 * b4^4 * b5^7,
#  7: (-4) * b2^2 * b34^2 * b4^4 * b5^7,
#  8: (-4) * b2^2 * b34^2 * b4^4 * b5^7}
# leadmxy
# {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,
#  4: (2) * b34^2 * b4^6 * b5^7,
#  5: (2) * b34^2 * b4^6 * b5^7,
#  6: (2) * b34^2 * b4^6 * b5^7,
#  7: (2) * b34^2 * b4^6 * b5^7,
#  8: (2) * b34^2 * b4^6 * b5^7}
# leadmy2
# {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,
#  4: (-1) * b34^2 * b4^6 * b5^7,
#  5: (-1) * b34^2 * b4^6 * b5^7,
#  6: (-1) * b34^2 * b4^6 * b5^7,
#  7: (-1) * b34^2 * b4^6 * b5^7,
#  8: (-1) * b34^2 * b4^6 * b5^7}
# leadmx
# {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,
#  4: b2^2 * b34^2 * b4^8 * b5^8,
#  5: b2^2 * b34^2 * b4^8 * b5^8,
#  6: b2^2 * b34^2 * b4^8 * b5^8,
#  7: b2^2 * b34^2 * b4^8 * b5^8,
#  8: b2^2 * b34^2 * b4^8 * b5^8}
# leadmxy
# {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,
#  4: (2) * b34^2 * b4^6 * b5^7,
#  5: (2) * b34^2 * b4^6 * b5^7,
#  6: (2) * b34^2 * b4^6 * b5^7,
#  7: (2) * b34^2 * b4^6 * b5^7,
#  8: (2) * b34^2 * b4^6 * b5^7}


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

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

#################
# xz-projection #
#################

labelCoefficientsxz = [Mx4, Mx3, Mx2]

allLeadingTermsxz ={}
for p in range(len(labelCoefficientsxz)):
    LTsp = {}
    for k in range(len(allWeightsSimplified)):
    	LTsp[k] = LTFromWeights(labelCoefficientsxz[p],allWeightsSimplified[k],R)
    allLeadingTermsxz[p] =  LTsp

# # We save the output:
# save(allLeadingTermsxz, "allLeadingTermsxz.sobj")

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

# # We record the number of leading terms per coefficient:
# [len(allLeadingTermsNoRepetitionsxz[p].keys()) for p in range(len(labelCoefficientsxz))]
# [1, 7, 1]

# labelCoefficientsxz = [Mx4, Mx3, Mx2]

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

#################
# yz-projection #
#################

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

allLeadingTermsyz ={}
for p in range(len(labelCoefficientsyz)):
    LTsp = {}
    for k in range(len(allWeightsSimplified)):
    	LTsp[k] = LTFromWeights(labelCoefficientsyz[p],allWeightsSimplified[k],R)
    allLeadingTermsyz[p] =  LTsp

# # We save the output:
# save(allLeadingTermsyz, "allLeadingTermsyz.sobj")


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

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

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

# # We record the possible leading terms:
# allLeadingTermsNoRepetitionsyz
# {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]}}