# This script computes relevant functions on the Yoshida coordinates:
# 1) Functions involving the action of W(E6) (computed by WE6_action.sage) and compositions of actions.
# 2) the intersection points of suitable pairs of extremal lines. We use this in "CalculationOfAllPds.sage"
# 3) Functions to convert from crosses to Yoshidas and from Yoshidas to polynomials in ds. Also, whenever possible, a function that converts a rational function in the ds to a rational function in the Yoshidas.

# The calculations of all cross functions and triples of four-factors associated to each Eckardt quintic is done in "four_factors.sage".

###########################################################################
# Setting up the ambient rings with the W(E6) action.
###########################################################################

# A list of variables d_1, ..., d_6 for the coefficients
ds = [var("d%s"%i) for i in range(1,7)]

# A list of variables E_1, ..., E_6
Es = [var("E%s"%i) for i in range(1,7)]

# A list of variables F_12, ..., F_56
Fs = [var("F%s%s"%(i,j)) for i in range(1,7) for j in range(1,7) if i < j]

# A list of variables G_1, ..., G_6
Gs = [var("G%s"%i) for i in range(1,7)]

# 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]

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

# A polynomial ring on all these generators
basering = FractionField(PolynomialRing(QQ, ds))
S = PolynomialRing(basering, Xs+Ys)
T = PolynomialRing(QQ, ds+Es+Fs+Gs)

YM = load("../Output/Yoshida_matrix.sobj")
YRing = PolynomialRing(QQ, Yoshidas)
YField = FractionField(YRing)

load("Yoshida.sage")


# Correspondence between X_{ij}, Y_{ijklmn} and anti-canonical triangles.
# X_{ij} = E_i F_{ij} G_j
# Y_{ijklmn} = F_{ij} F_{kl} F_{mn}
# correspondence = dict([])
# for i in range(1,7):
#     for j in range(1,7):
#         if i < j:
#             correspondence[S("X%s%s"%(i,j))] = T("E%s"%i)*T("F%s%s"%(i,j))*T("G%s"%j)
#         if j < i:
#             correspondence[S("X%s%s"%(i,j))] = T("E%s"%i)*T("F%s%s"%(j,i))*T("G%s"%j)
    
# 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:
#                             correspondence[S("Y%s%s%s%s%s%s"%(i,j,k,l,m,n))] =  T("F%s%s"%(i,j)) *  T("F%s%s"%(k,l)) * T("F%s%s"%(m,n))

#Note: the variables are missing the quintic denominators.

#save(correspondence, "Input/XYEFG_correspondence.sobj")
correspondence = load("Input/XYEFG_correspondence.sobj")



################################################
# 1) Functions involved in the action of W(E6) #
################################################

# allAction is a dictionary indexed by the 6 generators of W(E6). Working with dictionaries is faster than working with homomorphisms.  The commands SR() converts all the variables and expressions to Symbolic Ring. It is easier to manipulate symbolic expressions than keeping track of the various parental rings.
# The action was computed with the Sage script "WE6_action.sage"

allAction = load("Input/Symbolic_action_on_everything.sobj")

# We record the 6 entries of allAction (at the level of the ds action):
#
# allAction[1] = (12)
# allAction[2] = (23)
# allAction[3] = (34)
# allAction[4] = (45)
# allAction[5] = (56)
# allAction[6] = Exceptional switch
# [d1 + 1/3*d4 + 1/3*d5 + 1/3*d6,
#  d2 + 1/3*d4 + 1/3*d5 + 1/3*d6,
#  d3 + 1/3*d4 + 1/3*d5 + 1/3*d6,
#  1/3*d4 - 2/3*d5 - 2/3*d6,
#  -2/3*d4 + 1/3*d5 - 2/3*d6,
#  -2/3*d4 - 2/3*d5 + 1/3*d6]   


# In order to act by other elements with W(E6), we must be able to compose. The following functions allows us to compose pairs dictionaries:

def compose(dict1, dict2):
    answer = {}
    for k in dict1.keys():
        answer[k] = k.substitute(dict1).substitute(dict2)
    return answer

# The next function allows for inductive composition of multiple elements. We will use it to compose multiple dictionaries.
# By convention:
# Composition of [d1, d2, ..., dn] is dn o d(n-1) o ... o d2 o d1

def composeAll(l):
    answer = {}
    for d in l:
        answer = compose(d, answer)
    return answer


# From these, we can generate all actions involved in the script "CalculationOfAllPds.sage" 
a12 = allAction[1]
a13 = composeAll([allAction[1], allAction[2], allAction[1]])
a14 = composeAll([a13, allAction[3], a13])
a15 = composeAll([a14, allAction[4], a14])
a16 = composeAll([a15, allAction[5], a15])
a23 = allAction[2]
a24 = composeAll([allAction[2], allAction[3], allAction[2]])
a25 = composeAll([a24, allAction[4], a24])
a26 = composeAll([a25, allAction[5], a25])
a34 = allAction[3]
a35 = composeAll([allAction[3], allAction[4], allAction[3]])
a36 = composeAll([a35, allAction[5], a35])
a45 = allAction[4]
a46 = composeAll([allAction[4], allAction[5], allAction[4]])
a56 = allAction[5]
exc = allAction[6]



############################################
# 3) Functions used for Calculation of Pds #
############################################
 
# We use these functions as inputs for "CalculationsOfAllPds.sage"
# Pij = Intersection of Ei and Fij
# Gives a list of 0/Infinity entries where Infinity corresponds to anticanonical coordinates divisible by Ei or Fij. 
def p(i,j):
    line1 = T("E%s"%i)
    if i < j:
        line2 = T("F%s%s"%(i,j))
    if j < i:
        line2 = T("F%s%s"%(j,i))

    answer = dict()
    for a in correspondence.keys():
        if line1.divides(correspondence[a]) or line2.divides(correspondence[a]):
            answer[a] = Infinity
        else:
            answer[a] = 0
    return answer

# We now use the above to get actual coordinates for the Pij

# We load the 4x45 matrix whose rows span the P^3 spanned by the del Pezzo cubic X inside P^44. 
ModIntegerMM = load("Input/ModIntegerMM.sobj")

def pd(i,j):
    pij = p(i,j)
    # The pij above will give us which entries should be zero (indicated by Infinity in pij)

    Mij = ModIntegerMM.matrix_from_columns(sorted([S.gens().index(x) for x in pij.keys() if pij[x] == Infinity]))
    # Pick the submatrix of ModIntegerMM by selecting columns that should vanish.
        
    b = Mij.left_kernel().basis()[0]
    # Find the linear combination of the rows that makes those corresponding column entry vanish.
    
    v = ModIntegerMM.transpose() * vector(b)
    # Take the linear combination of the original matrix.

    # Return the result as a family (indexed by anti-canonical triangles X_kl, Y_klmnpq).
    return {S.gens()[i]:v[i] for i in range(0,len(S.gens()))}


# # We check that the matrix M12 has rank 3, hence the left-kernel is 1-dimensional:
# p12 = p(1,2)
# M12 = ModIntegerMM.matrix_from_columns(sorted([S.gens().index(x) for x in p12.keys() if p12[x] == Infinity]))
# print rank(M12)
# # 3
# M12.left_kernel().basis()
# [
# (0, 0, 0, 1)
# ]



################################################################################################################################
# 4) Functions to convert Crosses to Yoshidas an back, and expressions in roots and quintics to rational functions in Yoshidas #
################################################################################################################################

# We load the function that gives the possible 4 roots that we need to multiply each quintic to obtain a Cross function. There are precisely 3 choices.

load("four_factors.sage")

# input c = expression in d's
# output in terms of yoshida functions
roots=load("Input/rootsE6.sobj")

Y = list(load("../Output/Yoshida_functions.sobj"))

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)}

# Input: exp = A monomial expression that contains linears and quintics
# Output: an expression in terms of Yoshidas (Yoshidas[0] to Yoshidas[39]) correct up to a constant factor.

def DtoYHelper(linear_factors, quintic_factors):
    #if there are no quintics, we are done
    if len(quintic_factors) == 0:
        linear_vector = sum([e * vector(exponent_vector(expand(f),roots)) for (f,e) in linear_factors])
#       linear_vector = sum([e * vector(YoshidaToVector(f.factor())) for (f,e) in linear_factors])
        if linear_vector in span(ZZ, YM.rows()):
            yoshida_vector = YM.solve_left(linear_vector)
            return prod([Yoshidas[i]^yoshida_vector[i] for i in range(0, 40)])
        else:
            return None
    else:
        # Convert the first quintic factors into yoshidas in all possible ways.
        (quintic, e) = quintic_factors[0]
        ffs = four_factors(quintic)
        for f in ffs:
            # Express f * quintic as a difference
            product = f * quintic
            def expand_product(product):
                for i in range(0,40):
                    for j in range(0,i):
                        vi = VectorToYoshida(YM.rows()[i], roots)
                        vj = VectorToYoshida(YM.rows()[j], roots)
                        
                        work = [signi * Yoshidas[i] + signj*Yoshidas[j] for signi in [-1,1] for signj in [-1,1] if basering(signi * vi + signj*vj) == basering(product)]

                        if len(work) > 0:
                            return work[0]
                print "Nothing worked! Wrong four factors?"
                return None
                            
            product_yoshida = expand_product(basering(product))

            # Now express the "rest". "rest" is obtained by removing (quintic, e) from the quintics and adding the factors of f with exponent -e to the linears.
            rest_linear_factors = linear_factors + [(new_linear, -e * new_exp) for (new_linear, new_exp) in f]
            rest_quintic_factors = quintic_factors[1:]
            do_the_rest = DtoYHelper(rest_linear_factors, rest_quintic_factors)

            # If the rest succeeds, then return the answer combined with our quintic factor
            if do_the_rest != None:
                return product_yoshida^e*do_the_rest

        #If we are here, then it means that none of the ffs worked. This is bad
        print "Nothing worked! Was the expression invariant?"
        return None

# We test that the conversion is correct by evaluating at random integer values of ds. Since we are dividing by the quintics, we have to make sure that the quintic doesn't vanish for this particular choice. Whenever we get an error from this vanishing we have to repeat, so we add the while loop that exits whenever the evaluation was possible (success == True).

def DtoY(exp):
    linear_factors = []
    quintic_factors = []

    for (f,e) in exp:
        if f.total_degree() == 1:
            linear_factors.append((f,e))
            
        elif f.total_degree() == 5:
            quintic_factors.append((f,e))

    answer = DtoYHelper(linear_factors, quintic_factors)

    success = False
    while success == False:
        try:
            vald = {SR(d): int(random()*100-50) for d in ds}
            answer_val = SR(answer.substitute(y_to_d)).substitute(vald)
            exp_val = SR(expand(exp)).substitute(vald)
            success = True
        except Exception:
            print "Division by zero. Retrying."
            
    if answer_val == exp_val:
        return answer
    elif answer_val == -exp_val:
        return -1*answer
    else:
        print "DtoY helper is wrong! "
        print "Answer_val = " + str(answer_val)
        print "exp_val = " + str(exp_val)
        return None


#Tests
        
#DtoY(VectorToYoshida(YM[0], roots).factor())
#DtoY(-1*VectorToYoshida(YM[0], roots).factor() * VectorToYoshida(YM[17], roots).factor())
#DtoY(quintic12.factor() * four_factors(quintic12)[2])


# troubleMakers = []
# for k in range(0,40):
#     print str(k)
#     if DtoY(Y[k].factor()) == Yoshidas[k]:
#        print 'True with plus'
#     if DtoY(Y[k].factor()) == -1*Yoshidas[k]:
#        print 'True with minus'
#     else:
# 	print 'error'
# 	troubleMakers = troubleMakers + [k]


# DtoY(Y[40].factor())
# -Yoshida11^2*Yoshida12^2*Yoshida13*Yoshida3*Yoshida7/(Yoshida0*Yoshida10*Yoshida15*Yoshida17*Yoshida4*Yoshida8)
# Y[40].factor()
# (d5 - d6) * (d3 + d5 + d6) * (-d2 + d4) * (d2 + d4 + d6) * (d2 + d4 + d5) * (-d1 + d3) * (d1 + d5 + d6) * (d1 + d3 + d4) * (d1 + d2 + d3)




##########################################################
# 5) Correspondence Yoshida matrix and Yoshida functions #
##########################################################

# # As given, the list of 80 Yoshida functions (40 up to sign) does not correspond to the rows of the Yoshida matrix. We therefore need to reorder them in a consisitent way, namely [Y0,..., Y40, -Y0,..., -Y40].

# # The following functions create the approapriate order.

# # Y is the ordered list of Yoshida functions we have
# # YM is the Yoshida matrix

# factoredYoshidas = dict()
# for i in range(0,80):
#     print str(i)
#     factoredYoshidas[i] = YField(DtoY(basering(Y[i]).factor()))
    
# orderedYoshidas = dict()
# for i in range(0,40):
#      print '\n'
#      print str(i)
#      a = DtoY(VectorToYoshida(YM[i], roots).factor())
#      b = [x for x in factoredYoshidas.keys() if factoredYoshidas[x] == YField(a)][0]
#      print str(b)
#      c = [x for x in factoredYoshidas.keys() if factoredYoshidas[x] == YField(-1*a)][0]
#      print str(c)
#      orderedYoshidas[i] = (b, Y[b], +1)
#      orderedYoshidas[40+i] = (c, Y[c], -1)

# orderedYoshidaFunctions = {i: orderedYoshidas[i][1] for i in sorted(orderedYoshidas.keys())}

# save(orderedYoshidaFunctions, "Input/orderedYoshidaFunctionsDictionary.sobj")