# This is a sage script for computing the 25 fundamental Hodge classes
# for the semiample classes \bE_i

WE6 = WeylGroup('E6').as_permutation_group()
oWE6 = WE6.order(); # 51840

tmp_chartable = WE6.character_table() 
tmp_ccs = WE6.conjugacy_classes()

# We reorder the rows in the character table and then the conjugacy classes
# so that they match the character table obtained by
# the command 'Display(CharacterTable("W(E6)")' in GAP

# Characters of WE6 are uniquely determined by (dim, sum(values on conj classes)
dim_sum_to_gap_row = { 
    (1, 25): 0,
    (1, 5): 1,
    (10, 8): 2,
    (6, 10): 3,
    (6, 4): 4,
    (20, 20): 5,
    (15, 31): 6,
    (15, 21): 7,
    (15, 33): 8,
    (15, 15): 9,
    (20, 50): 10,
    (20, 18): 11,
    (24, 50): 12,
    (24, 36): 13,
    (30, 32): 14,
    (30, 20): 15,
    (60, 67): 16,
    (80, 57): 17,
    (90, 87): 18,
    (60, 71): 19,
    (60, 51): 20,
    (64, 71): 21,
    (64, 45): 22,
    (81, 92): 23,
    (81, 76): 24,
}

L = [[] for _ in range(25)]

for row in tmp_chartable:
    d, s = row[0], sum(row)
    L[ dim_sum_to_gap_row[(d, s)] ] = row

tmp_chartable2 = matrix(L)    

# These are the columns in GAP, for each conjugacy class
columns_to_gap_columns = { 
    (1,1,10,6,6,20,15,15,15,15,20,20,24,24,30,30,60,80,90,60,60,64,64,81,81): 0,
    (1,1,-6,-2,-2,4,-1,-1,7,7,4,4,8,8,-10,-10,12,-16,-6,-4,-4,0,0,9,9): 1,
    (1,1,2,2,2,-4,-1,-1,3,3,4,4,0,0,2,2,4,0,-6,4,4,0,0,-3,-3): 2,
    (1,1,1,-3,-3,-7,6,6,-3,-3,2,2,6,6,3,3,-3,-10,9,6,6,-8,-8,0,0): 3,
    (1,1,-2,3,3,2,3,3,0,0,5,5,0,0,3,3,-6,-4,0,-3,-3,4,4,0,0): 4,
    (1,1,4,0,0,2,0,0,3,3,-1,-1,3,3,3,3,0,2,0,-3,-3,-2,-2,0,0): 5,
    (1,1,2,2,2,4,3,3,-1,-1,0,0,0,0,-2,-2,4,0,2,0,0,0,0,-3,-3): 6,
    (1,1,-2,0,0,0,-1,-1,1,1,0,0,0,0,0,0,0,0,2,0,0,0,0,-1,-1): 7,
    (1,1,0,1,1,0,0,0,0,0,0,0,-1,-1,0,0,0,0,0,0,0,-1,-1,1,1): 8,
    (1,1,-3,1,1,1,2,2,1,1,-2,-2,2,2,-1,-1,-3,2,-3,2,2,0,0,0,0): 9,
    (1,1,0,1,1,-2,-1,-1,-2,-2,1,1,2,2,-1,-1,0,2,0,-1,-1,0,0,0,0): 10,
    (1,1,0,-2,-2,-2,2,2,1,1,1,1,-1,-1,-1,-1,0,2,0,-1,-1,0,0,0,0): 11,
    (1,1,2,-1,-1,2,-1,-1,0,0,1,1,0,0,-1,-1,-2,0,0,1,1,0,0,0,0): 12,
    (1,1,1,0,0,-1,0,0,0,0,-1,-1,0,0,0,0,0,-1,0,0,0,1,1,0,0): 13,
    (1,1,-1,-1,-1,1,0,0,-1,-1,0,0,0,0,1,1,1,0,-1,0,0,0,0,0,0): 14,
    (1,-1,0,4,-4,0,5,-5,5,-5,10,-10,4,-4,10,-10,0,0,0,10,-10,16,-16,9,-9): 15,
    (1,-1,0,0,0,0,-3,3,1,-1,2,-2,4,-4,-2,2,0,0,0,2,-2,0,0,-3,3): 16,
    (1,-1,0,-2,2,0,1,-1,3,-3,2,-2,0,0,-4,4,0,0,0,-2,2,0,0,3,-3): 17,
    (1,-1,0,2,-2,0,1,-1,-1,1,2,-2,0,0,0,0,0,0,0,-2,2,0,0,-1,1): 18,
    (1,-1,0,1,-1,0,-1,1,2,-2,1,-1,-2,2,1,-1,0,0,0,1,-1,-2,2,0,0): 19,
    (1,-1,0,-2,2,0,2,-2,-1,1,1,-1,1,-1,1,-1,0,0,0,1,-1,-2,2,0,0): 20,
    (1,-1,0,0,0,0,0,0,1,-1,-1,1,1,-1,1,-1,0,0,0,-1,1,0,0,0,0): 21,
    (1,-1,0,0,0,0,-1,1,-1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,1,-1): 22,
    (1,-1,0,-1,1,0,0,0,0,0,0,0,-1,1,0,0,0,0,0,0,0,1,-1,-1,1): 23,
    (1,-1,0,1,-1,0,1,-1,0,0,-1,1,0,0,-1,1,0,0,0,1,-1,0,0,0,0): 24,
}

ccs = ['' for _ in range(25)]
L = [[] for _ in range(25)]

for k, cc in enumerate(tmp_ccs):
    column = tuple( tmp_chartable2.column(k) )
    gap_k = columns_to_gap_columns[column]
    L[gap_k] = column
    ccs[gap_k] = cc
    
chartable = matrix.column(L)
ccreps = [cc[0] for cc in ccs]

# Labels used by GAP for the characters and conjugacy classes
chi_lbls = [str(k+1) for k in range(25)]
chi_dim_lbls = ['1', '1-', '10', '6', '6-', '20a',
              '15a', '15a-', '15b', '15b-', '20b', '20b-',
              '24', '24-', '30', '30-', '60a', '80', '90',
              '60b', '60b-', '64', '64-', '81', '81-']
def chi_dim_lbl_to_tex(x):
    if x[-1] == '-':
        return '\\tov{' + x[:-1] + '}'
    else:
        return x
tex_chi_dim_lbls = [chi_dim_lbl_to_tex(x) for x in chi_dim_lbls]    

cc_lbls = ['1a', '2a', '2b', '3a', '3b', '3c', '4a', '4b', '5a', '6a', '6b',
           '6c', '6d', '9a', '12a', '2c', '2d', '4c', '4d', '6e', '6f', '6g',
           '8a', '10a', '12b']

# Create dictionaries for easier access to chars, ccs, ccreps

Dchi = {x: chartable[k] for k, x in enumerate(chi_dim_lbls)}
Dcc = {x: ccs[k] for k, x in enumerate(cc_lbls)}
Dccrep = {x: ccreps[k] for k, x in enumerate(cc_lbls)}

def print_chartable():
    """Prints the character table of W(E_6), produced by GAP using Display(CharacterTable("W(E6)"))"""
    print('\n\\begin{table}[htp!]\\centering')
    print('\\setlength{\\tabcolsep}{1.5pt}')
    print('\\begin{tabular}{|rl|rrrrrrrrrrrrrrr|rrrrrrrrrr|} \\hline')
    print('$\\chi$ & name &',  ' & '.join(cc_lbls), '\\\\')
    print('\\hline')
    
    def view(x):
        if x == 0:
            return '.'
        else:
            return str(x)
    
    for k, row in enumerate(chartable.rows()):
        str_row = [view(x) for x in row]
        print(chi_lbls[k], '&', tex_chi_dim_lbls[k], r'\ {} &', ' & '.join(str_row), '\\\\')
        
    print('\\hline \\end{tabular}')
    print('\\medskip\\caption{The character table of $W(E_6)$}')
    print('\\label{tab:E6-chartable}')
    print('\\end{table}')

# Three special conjugacy classes
cc2c = ccs[cc_lbls.index('2c')]
cc2b = ccs[cc_lbls.index('2b')]
cc3b = ccs[cc_lbls.index('3b')]
cc2crep, cc2brep, cc3brep = cc2c[0], cc2b[0], cc3b[0] # representatives

def irrep_mults(G):
    r"""Let G be a subgroup of WE6. Returns a list of length 25. Entries are dim( (\rho_\alpha)^G ) = \sum_{g\in G} Tr_{\rho_\alpha} (g) / |G| for the 25 irreps \alpha of WE6"""
    answer = vector([0 for _ in range(25)])
    for g in G:
        for k in range(25):
            if g in ccs[k]: 
                answer += chartable.column(k)
    return answer / G.order()

WE6_cyclic_subgroups = [WE6.subgroup([g]) for g in ccreps]

def irrep_mults_B():
    r"""Return the matrix of multiplicities of \bE_i in H^0(K_{B_\alpha})
for the 25 curves B_\alpha = \wC / W_\alpha, where W_\alpha is the cyclic group generated by a representative w_\alpha of a conjugacy class. 
Notation: rows for characters, columns for conjugacy classes."""
    return matrix.column(
        [irrep_mults(G) for G in WE6_cyclic_subgroups]
    )

def irrep_mults_A():
    r"""Return the matrix of multiplicities of \bE_i in H^0(K_{A_\alpha})
for the 25 curves A_\alpha = \wC / Z_\alpha, where Z_\alpha is the centralizer 
of w_\alpha. Notation: rows for characters, columns for conjugacy classes."""
    return matrix.column(
        [ [1] + [0 for _ in range(24)] ] +  # the first column is easy
        [ irrep_mults(WE6.centralizer(g)) for g in ccreps[1:] ]
    )

def degree(G):
    r"""Degree of \wC / G over P^1"""
    return oWE6 / G.order()

def genus(G):
    r"""Genus of \wC/ G"""
    [a2c, b2c, a2b, b2b, a3b, b3b] = ramif_profile(G)
    return 10 * a2c - b2c + 1

def ramif_profile(G):
    r"""Return [a_2c, b_2c, a_2b, b_2b, a_3b, b_3b], the ramification profile for a curve C = \wC/G"""
    answer = []
    deg = oWE6 / G.order()
    SG = set(G)
    
    for c in [cc2crep, cc2brep, cc3brep]:
        cc = WE6.conjugacy_class(c)
        Scc = set(cc)
        Sintersection = Scc & SG
        Zc = WE6.centralizer(c)
        b = len(Sintersection) * Zc.order() / G.order()
        a = (deg - b) / c.order()
        answer += [a, b]
    return answer
 
def divisor(G):
    r"""For the curve \wC/G, return a list [c0, c_syz, c_azy] for the correponding lambda-divisor c0 D0 + c_syz D_syz + c_azy D_azy on \overline{Hur}"""
    [a_2c, b_2c, a_2b, b_2b, a_3b, b_3b] = ramif_profile(G)
    return divisor_from_as(a_2c, a_2b, a_3b)

def divisor_from_as(a_2c, a_2b, a_3b):
    r"""Starting with a_2c, a_2b, a_3b, return a list [c0, c_syz, c_azy] for the correponding lambda-divisor c0 D0 + c_syz D_syz + c_azy D_azy"""
    c0 = 11 * a_2c / 92
    c_syz = (66 * a_2c /23 - 3 * a_2b / 2) / 6
    c_azy = (66 * a_2c /23 - 8 * a_3b / 3) / 8
    return [c0, c_syz, c_azy]

def a2c2b3b():
    mults = irrep_mults_B()
    inv = mults.transpose().inverse()
    L = []
    for G in WE6_cyclic_subgroups:
        [a2c, b2c, a2b, b2b, a3b, b3b] = ramif_profile(G)
        L.append([a2c, a2b, a3b])
    return inv * matrix(L)

def degs_ranks_lambdas_a2c2b3b():
    r"""Return 25 x 7 matrix of ranks, lambda classes, and a2c, a2b, a3b of the fundamental Hodge bundles \bE_i"""
    mults = irrep_mults_B()
    inv = mults.transpose().inverse()

    L = []
    for G in WE6_cyclic_subgroups:
        [a2c, b2c, a2b, b2b, a3b, b3b] = ramif_profile(G)
        L.append(
            [genus(G)] +
            divisor(G) +
            [a2c, a2b, a3b]
        )
    return inv * matrix(L)

def ranks_and_lambdas():
    r"""Return 25 x 7 matrix of ranks and lambda classes (in two bases) of the fundamental Hodge bundles \bE_i"""
    mults = irrep_mults_B()
    inv = mults.transpose().inverse()
    g = inv * matrix.column([genus(G) for G in WE6_cyclic_subgroups])
    div = inv * matrix([ divisor(G) for G in WE6_cyclic_subgroups ])

    # Now recompute in a new basis
    m = matrix([ div[3], [0, 1/2, 0], div[1] ])
    div2 = div * m.inverse()

    answer = matrix.column(g.columns() + div.columns() + div2.columns())
    return answer

def print_matrix_latex(m, header, footer):
    print(header)
    for k, lbl in enumerate(tex_chi_dim_lbls):
        print(str(k+1), '&', lbl, '&', ' & '.join([str(x) for x in m[k]]), '\\\\')
    print(footer)

def print_degs_ranks_lambdas_a2c2b3b():
    r"""Print 25 x 7 matrix of ranks and lambda classes (in two bases) of the fundamental Hodge bundles \bE_i"""
    m = degs_ranks_lambdas_a2c2b3b()
    header = (
        '\n' + r'\newcommand\tov[1]{$\overline{\text{#1}}$}' + 
        '\n\\begin{table}[htp!]\\centering' + 
        '\n\\begin{tabular}{|rl|c|ccc|ccc|} \\hline' +
        '\n$\\chi$ & name & $\\rk\\bE_i$ & $D_0$ & $D\\syz$ & $D\\azy$ & $a_{2c}$ & $a_{2b}$ & $a_{3b}$  \\\\' +
        '\n\\hline\n'
        )
    footer = (
        '\\hline \\end{tabular}' +
        '\n\\medskip\\caption{$\\chi_i$, $\\rk\\bE_i$, $\\lambda_i$, and $(a_{2c},a_{2b},a_{3b})(\\chi_i)$}' +
        '\n\\label{tab:rk-lambda}' +
        '\n\\end{table}'
    )
    print_matrix_latex(m, header, footer)

def define_divisors(coef_of_dsyz=1, coef_from_jets=1):
    r"""Define all the main divisors. Assume multiplying by k in the formula for kappa"""
    global lp, div
    global d0,    ds,    da,    lam,   la5,   fn,    fq,   ka6, kap, gam, fd
    global nd0,   nds,   nda,   nlam,  nla5,  nfn,   nfq,  nka6, nkap, ngam, nfd
    
    ral = ranks_and_lambdas()
    lp = { k+1: ral[k][1:4] for k in range(25) }

    d0, ds, da = vector([1,0,0]), vector([0,1,0]), vector([0,0,1])
    la5 = lp[4]
    lam = lp[1] + lp[4] + lp[11]
    ka6 = 7 * la5 - d0

    kap = 12 * lam - 6 * d0 - 2 * ds * coef_of_dsyz  # coef of ds is crucial!
    gam = da - kap / 3 + d0 + ds * coef_of_dsyz * coef_from_jets

    fn = - lam - la5 + gam
    #fn = -5 * lam - la5 + 3 * d0 + 5/3 * ds + da

    fd = - lam - kap + 15/2 * gam
    
    fq = 39/11 * la5 + 12/11 * d0 + 39/11 * fn - lam
    
    div = {
        'd0': d0, 'ds': ds, 'da': da,
        'lam': lam, 'la5': la5,
        'kap': kap, 'gam': gam, 'fn': fn, 'fd': fd,
        'fq': fq, 'ka6': ka6,
        }

    print('Assuming that the coef of Dsyz in kapps is', coef_of_dsyz)

    print('\nIn the basis d0, ds, da:')
    for x in [
            'd0', 'ds', 'da', 'lam', 'la5',
            'kap', 'gam', 'fn', 'fd', 'ka6',
            #'fq', 
    ]:
        y = div[x]
        #print('debug y=', y, 'type(y)=', type(y))
        #print(x, '\t', "{:<27}".format(y), normalize_vec(y, 2).n(digits=3))
        print(x, '\t', y, normalize_vec(y, 2).n(digits=3))

    print('\nAlso available: lp[k] for k = 1...25')

    mi = matrix([la5, d0, fn]).transpose().inverse()
    print('\nIn the basis la5, d0, fn:')
    for x in [
            'd0', 'ds', 'da', 'lam', 'la5',
            'kap', 'gam', 'fn', 'fd', 'ka6',
            #'fq', 
    ]:
        y = mi * div[x]
        #print(x, '\t', "{:<27}".format(y), normalize_vec(y, 0).n(digits=3))
        print(x, '\t', y, normalize_vec(y, 0).n(digits=3))

def div_ifnzero(v):
    """Assuming coef_of_dsyz=2 and n=0, return divisor in basis ds, da"""
    return v + v[1] * vector([39/28, -1, -59/140])
    
def normalize_vec(v, place):
    """Starting with (c0, c1, c2), return (c0/c2, c1/c2, 1)"""
    if v[place] == 0:
        return v
    else:
        return v / abs(v[place])

print('Run print_degs_ranks_lambdas_a2c2b3b()')
#print('Run find_a6_slope() or find_a6_slope(2)')

WE6_maximal_subgroups = [
    WE6.subgroup([ '(1,16,18,26,25,22,17,2,3)(4,7,14,13,23,10,24,11,19)(5,9,8,6,12,20,27,21,15)', '(4,6,8)(5,7,9)(10,13,11)(12,16,14)(15,19,17)(25,26,27)' ]),
    WE6.subgroup([ '(6,8)(7,9)(10,11)(12,14)(15,17)(26,27)', '(1,18,13,10,6)(2,4,8,11,21)(3,14,27,22,15)(5,23,9,25,19)(7,12,16,20,26)' ]),
    WE6.subgroup([ '(1,22)(2,20)(3,25,6,24,4,26)(5,10,7)(8,27)(9,11,13)(12,15)(14,21,16,17,18,19)', '(1,25)(2,22)(3,24)(4,20)(6,27)(7,17)(8,26)(9,15)(10,11)(12,14)(13,21)(16,23)' ]), 
  WE6.subgroup([ '(1,6,15,25,27,9)(2,10,26,24,14,8)(3,4,21,22,20,16)(5,18,19)(7,12,23,17,11,13)', '(1,9,19,23,25,10)(2,16,17,26,4,20)(3,24,15,12,5,7)(6,14,22)(8,27,21,18,11,13)' ]),
    WE6.subgroup([ '(1,12,7,10,15,2,4,3)(5,6)(8,14,16,20,26,22,21,17)(9,18,24,13,25,19,11,23)', '(1,9,24,26,21,4)(2,3,16,27,25,10)(5,7,19,23,18,11)(6,8,14,20,22,15)(12,13,17)' ]),
    WE6.subgroup([ '(1,22)(3,15,6,19)(4,21,8,17)(5,10,13,9)(12,26,16,24)(14,25,18,27)', '(2,13,9,7)(3,6,8,16)(4,18,14,12)(5,10,11,20)(15,21,27,24)(17,26)', '(1,23)(2,17,9,24,20,26,10,21)(3,14,16,12,18,6,4,8)(5,27,7,25,13,15,11,19)' ])
]

# sage: [oWE6/G.order() for G in WE6_maximal_subgroups]
# [2, 27, 36, 40, 40, 45]

# sage: [irrep_mults(G) for G in WE6_maximal_subgroups]
# [(1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),
#  (1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),
#  (1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),
#  (1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),
#  (1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),
#  (1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)]

def find_a6_slope(coef_of_dsyz=1, coef_from_jets=1):
    define_divisors(coef_of_dsyz, coef_from_jets)
    p = Polyhedron(base_ring=QQ, rays=[d0,ds,da,lam,la5,fn,fd])

    var('s')
    fns = [vector(x[1:]) for x in p.inequalities_list()]
    ineqs = [fns[k] * (s * la5 - d0) for k in range(len(fns))]
    answer = solve([f >= 0 for f in ineqs], s)
    print('\nConclusion: s * la5 - d0 is effective if:', answer)