#This script computes the 270 equations that are satisfied by the Yoshidas-adapted anticanonical coordinates. These coodinates are obtained as  
# newCoordinate = oldCoordinate / quintic(oldCoordinate)


ds = [var("d%s"%i) for i in range(1,7)]
Sd =  PolynomialRing(QQ, ds)

#The following are the 45 new anticanonical coordinates.
# May be it would have been better to name the new variables something else.

# A list of variables X_12, ..., X_65
Xs = [var("X%s%s"%(i,j)) for i in range(1,7) for j in range(1,7) if i != j]

# A list of variables Y_123456, ..., Y_162534
Ys = [var("Y%s%s%s%s%s%s"%(i,j,k,l,m,n)) for i in range(1,7) for j in range(1,7) for k in range(1,7) for l in range(1,7) for m in range(1,7) for n in range(1,7) if i<j and k<l and m<n and i<k and k<m and len(uniq([i,j,k,l,m,n])) == 6]


S = PolynomialRing(QQ, ds+Xs+Ys)


############################################################
# Calculation of Linear equations in different coordinates #
############################################################

roots = load("Input/rootsE6.sobj")

quintics = load("Input/quintics.sobj")

##############################################
# Original equations from Ren-Shaw-Sturmfels #
##############################################
#This is the original linear equation (used in the file EquationsAnticanonicalEmbedding/Scripts/linear_equations.sage)

linear = set([expand((d3-d4)*(d1+d3+d4)*X21 - (d2-d4)*(d1+d2+d4)*X31 + (d2-d3)*(d1+d2+d3)*X41)])

# We multiply by 2 convenient roots and the quintic factors associated to each anticanonical triangle. The choice of these two roots is explained below.
# This is the resulting modified linear equation (that can be written in terms of Yoshida functions.

linear = set([expand((d5-d6)*(d1+d5+d6)*((d3-d4)*(d1+d3+d4)*quintics[S(X21)]*X21 - (d2-d4)*(d1+d2+d4)*quintics[S(X31)]*X31 + (d2-d3)*(d1+d2+d3)*quintics[S(X41)]*X41))])

# Unfortunately, it is too unwieldy for orbit calculations.
# So we will first convert it in terms of Yoshida functions.


##################################
# Equations in Yoshida functions #
##################################

# Denote by Yoshidai the Yoshida function corresponding to the ith row of the Yoshida matrix Yoshida_matrix.sobj where the roots are rootsE6.sobj, as we wrote before

# A list of variables Yoshida0, ..., Yoshida39
Yoshidas = [var("Yoshida%s"%i) for i in range(0,40)]

YM = load("../Output/Yoshida_matrix.sobj")
y_to_d = {Yoshidas[i]: prod([roots[j]^YM[i][j] for j in range(0,36)]) for i in range(0,40)}

T = PolynomialRing(QQ, Yoshidas + Xs + Ys)

# Then we have

# (d5-d6)*(d1+d5+d6)*(d3-d4)*(d1+d3+d4)*quintics[S(X21)] = Yoshida3-Yoshida37
# (d5-d6)*(d1+d5+d6)*(d2-d4)*(d1+d2+d4)*quintics[S(X31)] = Yoshida20-Yoshida8
# (d5-d6)*(d1+d5+d6)*(d2-d3)*(d1+d2+d3)*quintics[S(X41)] = Yoshida3-Yoshida5

# So the linear equation becomes:

linear = set([expand((Yoshida3-Yoshida37)*X21 - (Yoshida20-Yoshida8)*X31 + (Yoshida3-Yoshida5)*X41)])


##################################
# Equations in Crosses functions #
##################################

# Finally, we get the linear equation in terms of crosses

crosses = load("Input/binomial_yoshidas.sobj")
cross_variables = [var("Cross%s"%i) for i in range(0,135)]
cross_to_yoshidas = load("Input/cross_to_yoshidas.sobj")

# new dictionary:
linear = set([expand(Cross116*X21 - Cross2*X31 + Cross19*X41)])
# old dictionary:
#linear = set([expand(-Cross57*X21 + Cross17*X31 + Cross49*X41)])


######################################################################
# CALCULATION OF ORBITS of the linear equations in terms of Crosses. #
######################################################################


# Since the new X's = old X / quintic , the action of W(E6) on the Yoshida-adapted anticanonical coordinates changes: it acquires a sign.
# The action stored in "Input/Symbolic_action_on_everything.sobj" is on (old X / quintic) variables.

# We load the action of W(E6), previously computed on OperationsOnYoshidas.sage
allAction=load("Input/action_on_d_XY_yoshida_cross.sobj")

load("orbit_optimized.sage")

# linears = computeOrbit(linear, allAction)
# save(linears, "../Output/cross_linears.sobj")

linears = load("../Output/cross_linears.sobj")

# Remove the duplicates (up to sign) and store the result, both as a sage object and in plain text:
# load("../../General/RemoveDoubles.sage")
# only270 = removeDoubles(linears)
# save(only270, "../Output/270_cross_linears.sobj")

# f = file("../Output/270_cross_linears.txt", 'w')
# f.write(str(only270))
# f.close()



# We write the linear equations in terms of the Yoshida functions and store it. We only need to store 270 (the orbit calculation gave us 540 because of +/- sign). We use the removeDoubles function for this:

cross_to_y = {c : cross_to_yoshidas[c][0] for c in cross_to_yoshidas.keys()}

# modified_linears_Yoshidas = set([])
# i = 0
# for l in linears:
#     modified_linears_Yoshidas.add(l.substitute(cross_to_y))
#     i = i+1
#     print i
# save(modified_linears_Yoshidas, "../Output/yoshida_linears.sobj")

# We record the result as a plain text:
f = file("../Output/yoshida_linears.txt", 'w')
f.write(str(modified_linears_Yoshidas))
f.close()

# Remove the duplicates (up to sign) and store the result, both as a sage object and in plain text:
# # load("../../General/RemoveDoubles.sage")
# only270 = removeDoubles(modified_linears_Yoshidas)
# save(only270, "../Output/270_modified_linears_Yoshidas.sobj")

# f = file("../Output/270_modified_linears_Yoshidas.txt", 'w')
# f.write(str(only270))
# f.close()



###################################################
# Solution to the linear equations in terms of ds #
###################################################

# We find the solution to these equations in terms of the ds. For this, we first need to convert the linear equations from crosses to ds.


modified_linears_ds = set([])
i = 0
for l in linears:
    modified_linears_ds.add(l.substitute(cross_to_y).substitute(y_to_d))
    i = i+1
    print i

#save(modified_linears_ds, "../Output/modified_linear_ds.sobj")

# We run the script that solves this system of 270 linear equations in the variables Xs + Ys. The coefficients are rational functions in the ds. We remove all demonominators of the 4x45-solution matrix.

# load("../EquationsAnticanonicalEmbedding/Scripts/SolutionsInDsMatrix.sage")

# The calculation ran out of memory many times. We saved the output.
ModIntegerMM = load("ModIntegerMM.sobj")

# # Factor the entries and record the non-zero ones to generate a plain text file:
# ModIntegerMMFactored = {}
# for i in range(0, ModIntegerMM.nrows()):
#     for j in range(0, ModIntegerMM.ncols()):
#     	if ModIntegerMM[i][j] != 0:
# 	   ModIntegerMMFactored[(i,j)] = factor(ModIntegerMM[i][j])
#
# # Save the file as a plain text:
# f = file("../Output/SolutionMatrix_4x45_In_Ds.txt", 'w')
# f.write(str(ModIntegerMMFactored))
# f.close()


# The following checks that ModIntegerMM satisfies the linear equations.
# For faster computation, d's are evaluated at random values.
vald = {d: int(random()*100-50) for d in ds}
val_modified_linears_ds = [l.substitute(vald) for l in modified_linears_ds]

dring = valModIntegerMM[0][0].numerator().parent()
valdd = {dring(d): vald[d] for d in vald.keys()}
valModIntegerMM = ModIntegerMM.substitute(valdd)

def evaluate(l, v):
    dictXY = {(Xs+Ys)[i]:v[i] for i in range(0,45)}
    return expand(l.substitute(dictXY))

evals = [evaluate(l,v) for l in val_modified_linears_ds for v in valModIntegerMM]
all([e == 0 for e in evals])

#########################################################
# Solution to the linear equations in terms of Yoshidas #
#########################################################

# Convert ModIntegerMM to Yoshidas.
load("OperatorsOnYoshidas.sage")

# YoshidaIntegerMM =  {}
# def makeYoshidaIntegerMM():
#     for i in range(0,4):
#         for j in range(0,45):
#             if (i,j) in YoshidaIntegerMM.keys():
#                 continue
            
#             if ModIntegerMM[i][j] == 0:
#                 YoshidaIntegerMM[(i,j)] = 0
#             else:
#                 YoshidaIntegerMM[(i,j)] = SR(DtoY(ModIntegerMM[i][j].factor()))

# while len(YoshidaIntegerMM.keys()) < 180:                
#     try:
#         makeYoshidaIntegerMM()
#     except Exception:
#         print "retrying"

#save(YoshidaIntegerMM, "../Output/YoshidaIntegerMM.sobj")

YoshidaIntegerMM=load("../Output/YoshidaIntegerMM.sobj")

# # Save the Solution Matrix in terms of Yoshida functions as a plain text file:             
# f = file("../Output/YoshidaIntegerMM.txt", 'w')
# f.write(str(YoshidaIntegerMM))
# f.close()

                
# The following checks that YoshidaIntegerMM satisfies the linear equations.
# For faster computation, d's are evaluated at random values.
print 'Checking YoshidaIntegerMM'
vald = {d: int(random()*30-15) for d in ds}
valY = {Y:SR(y_to_d[Y]).substitute(vald) for Y in y_to_d.keys()}
YoshidaIntegerMM = load("../Output/YoshidaIntegerMM.sobj")
valYoshidaIntegerMM = {k:YoshidaIntegerMM[k].substitute(valY) for k in YoshidaIntegerMM.keys()}
valMatrixYoshidaIntegerMM = matrix(valYoshidaIntegerMM)

linears = list(load("../Output/cross_linears.sobj"))
# Convert from cross to Yoshidas and then evaluate
val_linears = [l.substitute(cross_to_y).substitute(valY) for l in linears]

evals = {(i,j):evaluate(val_linears[i], valMatrixYoshidaIntegerMM.rows()[j]) for i in range(0,540) for j in range(0,4)}
all([e==0 for e in evals])