####################################################################################
# Groebner cone computations common refinement coefficients xz- and yz-projections #
####################################################################################

# In order to find all possible combinatorial types of resulting modified 3D-tropical curves with Theta graphs as minimal skeletons, we compute the Groebner fans of the product of all non-monomial coefficients of the (x,z) and (y,z)-projections. Afterwards, we intersect the fan with the TypeIICone, check for subdivisions (in "ComputingSubdivisionsCoefficientsxzyz.sage"), and compute the leading terms for each piece in the subdivision.

# This script computes the Grobner Cones for each of the polynomials.

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

# The necessary date (including the polynomials) is obtained by loading the following file:

load("macroslocalJInvariantComputationsAndPrimesWithBadReduction.sage")

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

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

###################
# XZ Coefficients #
###################

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 record the output as plain text file:
# f = file("RelevantCoefficientsXZThetaCone.txt",'w')
# f.write('Mx4 = ' + str(Mx4) +'\n'+'\n')
# f.write('Mx3 = ' + str(Mx3) +'\n'+'\n')
# f.write('Mx2 = ' + str(Mx2) +'\n'+'\n')
# f.close()

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

allM = Mx4*Mx3*Mx2

###################
# YZ Coefficients #
###################

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 record the output as plain text file:
# f = file("RelevantCoefficientsYZThetaCone.txt",'w')
# f.write('my4 = ' + str(mx4) +'\n'+'\n')
# f.write('my3z = ' + str(mx3y) +'\n'+'\n')
# f.write('my2z2 = ' + str(mx2y2) +'\n'+'\n')
# f.write('myz3 = ' + str(mxy3) +'\n'+'\n')
# f.write('mz4 = ' + str(my4) +'\n'+'\n')
# f.write('my3 = ' + str(mx3) +'\n'+'\n')
# f.write('my2z = ' + str(mx2y) +'\n'+'\n')
# f.write('myz2 = ' + str(mxy2) +'\n'+'\n')
# f.write('mz3 = ' + str(my3) +'\n'+'\n')
# f.write('my2 = ' + str(mx2) +'\n'+'\n')
# f.write('my2 = ' + str(mxy) +'\n'+'\n')
# f.write('mz2 = ' + str(my2) +'\n'+'\n')
# f.write('my = ' + str(mx) +'\n'+'\n')
# f.write('mz = ' + str(my) +'\n'+'\n')
# f.close()

# 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

################
# Combination  #
################


allmyzxz = allM*allm


# gfm = GrobnerFanComputations(fromPolynomialsToExtremalMonomials(allmyzxz,R))

# We save it:
# save(gfm, "GroebnerFansCommonRefinementyzxz.sobj")

gfm = load("GroebnerFansCommonRefinementyzxz.sobj")

# We record the f-vector:
# print gfm.f_vector()
# [1, 21, 54, 35]

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

# We record the coordinates of the rays of the Grobner fan:
RaysGFm = gfm.rays()
# print RaysGFm
# [[-3, -3, 5, 1], [-3, 1, 1, 1], [-1, -5, 3, 3], [-1, -1, -1, 3], [-1, -1, 1, 1], [-1, -1, 3, -1], [-1, 1, -1, 1], [-1, 1, 1, -1], [-1, 3, -1, -1], [0, -1, 1, 0], [0, 1, -1, 0], [1, -3, 1, 1], [1, -3, 5, -3], [1, -1, -1, 1], [1, -1, 1, -1], [1, 1, -3, 1], [1, 1, -1, -1], [1, 1, 1, -3], [3, -5, 3, -1], [3, -1, -1, -1], [5, 1, -3, -3]]

# print gfm.lineality_dim()
# 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 76 rays:
MaximalConesgfm = gfm.maximal_cones().values()[0]
# print MaximalConesgfm
# [[0, 1, 4], [0, 1, 5], [0, 2, 4, 9], [1, 3, 4], [1, 3, 6], [1, 5, 7], [1, 7, 8], [1, 6, 8, 10], [0, 5, 9], [2, 3, 4], [2, 3, 11], [3, 6, 15], [5, 7, 17], [7, 8, 17], [2, 9, 11], [5, 9, 12], [6, 10, 15], [8, 10, 16], [3, 11, 13], [5, 12, 17], [3, 13, 15], [8, 16, 17], [10, 15, 16], [12, 14, 17], [9, 11, 18], [9, 12, 14, 18], [11, 13, 19], [13, 15, 19], [14, 17, 19], [11, 18, 19], [14, 18, 19], [15, 16, 20], [16, 17, 20], [15, 19, 20], [17, 19, 20]]

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

PolyMaximalConesgfmDict = MaximalCones(gfm)
# save(PolyMaximalConesgfmDict, "PolyMaximalConesgfCommonRefinementDictyzxz.sobj")

# set([dim(x) for x in PolyMaximalConesgfmDict.values()])
# {4}