| Download
Project: admcycles
Views: 723Visibility: Unlisted (only visible to those who know the link)
Image: ubuntu2004r"""1Code for computations involving Witten's r-spin class.2Authors: Felix Janda (main author), Aaron Pixton (code improvements), Johannes Schmitt (integration into admcycles)3"""45from __future__ import absolute_import, print_function67import itertools89from six.moves import range1011from admcycles.admcycles import Tautv_to_tautclass, tautclass12from admcycles.DR import *1314# from sage.combinat.subset import Subsets15from sage.arith.all import factorial, lcm16from sage.functions.other import floor, ceil17from sage.misc.misc_c import prod18from sage.rings.all import PolynomialRing, QQ, ZZ19from sage.modules.free_module_element import vector20# from sage.rings.polynomial.polynomial_ring import polygen21from sage.rings.polynomial.multi_polynomial_element import MPolynomial22from collections import Iterable23from sage.rings.integer import Integer24from sage.calculus.var import var25from sage.misc.functional import symbolic_sum262728# Computation of Witten's rspin class for large r2930def rspin_leg_factor(d,a):31if a < 0: a = X + a32return Pmpolynomial(d).subs(a=a)3334def rspin_edge_factor(w1,m,d1,d2):35R=PolynomialRing(QQ,1,'X')36X = R.gens()[0]37d = d1+d2+138S = 039for i in range(d+1):40S += R(rspin_leg_factor(i, (w1 + i) % (m-1))(X=m)*rspin_leg_factor(d-i, (m-2-i-w1) % (m-1))(X=m)*X**i)41S /= X+142assert(S.denominator() == 1)43S = S.numerator()44#print(-S[d1])45return -S[d1]4647def rspin_coeff_setup(num,g,r,n=0,dvector=(),moduli_type=MODULI_ST):48markings = tuple(range(1,n+1))49G = single_stratum(num,g,r,markings,moduli_type)50nr = G.M.nrows()51nc = G.M.ncols()52edge_list = []53exp_list = []54scalar_factor = 1/autom_count(num,g,r,markings,moduli_type)55given_weights = [1 - G.M[i+1,0][0] for i in range(nr-1)]56for i in range(1,nr):57for j in range(1,G.M[i,0].degree()+1):58scalar_factor /= factorial(G.M[i,0][j])59scalar_factor *= (-1)**G.M[i,0][j]60scalar_factor *= (rspin_leg_factor(j, 0))**(G.M[i,0][j])61given_weights[i-1] -= j*G.M[i,0][j]62for j in range(1,nc):63ilist = [i for i in range(1,nr) if G.M[i,j] != 0]64if G.M[0,j] == 0:65if len(ilist) == 1:66i1 = ilist[0]67i2 = ilist[0]68exp1 = G.M[i1,j][1]69exp2 = G.M[i1,j][2]70else:71i1 = ilist[0]72i2 = ilist[1]73exp1 = G.M[i1,j][1]74exp2 = G.M[i2,j][1]75edge_list.append([i1-1,i2-1])76exp_list.append(exp1)77exp_list.append(exp2)78else:79exp1 = G.M[ilist[0],j][1]80scalar_factor *= rspin_leg_factor(exp1, dvector[G.M[0,j][0]-1])81given_weights[ilist[0] - 1] += dvector[G.M[0,j][0] - 1] - exp182return edge_list,exp_list,given_weights,scalar_factor8384def rspin_coeff(num,g,r,n=0,dvector=(),r_coeff=None,step=1,m0given=-1,deggiven=-1,moduli_type=MODULI_ST):85markings = tuple(range(1,n+1))86G = single_stratum(num,g,r,markings,moduli_type)87nr = G.M.nrows()88nc = G.M.ncols()89edge_list,exp_list,given_weights,scalar_factor = rspin_coeff_setup(num,g,r,n,dvector,moduli_type)90if m0given == -1:91m0 = (ceil(sum([abs(i.subs(X=0)) for i in dvector])/2) + g*1+1)*step92else:93m0 = m0given94h0 = nc - nr - n + 195if deggiven == -1:96deg = 2*sum(exp_list) + 2*len(edge_list)97else:98deg = deggiven99if r_coeff is None:100mrange = list(range(m0 + step, m0 + step*deg + step + 1, step))101else:102mrange = [r_coeff+1] # just evaluate at a single value m = r_coeff103mvalues = []104for m in mrange:105given_weights_m = [ZZ(i.subs(X=m)) for i in given_weights]106total = 0107for weight_data in itertools.product(*[list(range(m-1)) for i in range(len(edge_list))]):108vertex_weights = copy(given_weights_m)109for i in range(len(edge_list)):110vertex_weights[edge_list[i][0]] += weight_data[i]111vertex_weights[edge_list[i][1]] -= weight_data[i]+2+exp_list[2*i]+exp_list[2*i+1]112if len([i for i in vertex_weights if i % (m-1) != 0]) > 0:113continue114term = 1115for i in range(len(edge_list)):116term *= rspin_edge_factor(weight_data[i], m, exp_list[2*i], exp_list[2*i+1])117total += term118if r_coeff is None:119mvalues.append(total*ZZ(m-1)**(-h0))120else:121#print(m-1)122mvalues.append(total*ZZ(m-1)**(-r-h0)) # undo the rescaling by r**degree123#print(mrange,mvalues, scalar_factor, r, h0)124mpoly = ZZ(-1)**(r-g)*(interpolate(mrange, mvalues)*scalar_factor).simplify_rational()125if r_coeff is not None:126mpoly = QQ(mpoly.subs(X=r_coeff))127return mpoly128129def rspin_compute(g,r,n=0,dvector=(),r_coeff=None,step=1,m0=-1,deg=-1,moduli_type=MODULI_ST):130answer = []131markings = tuple(range(1,n+1))132for i in range(num_strata(g,r,markings,moduli_type)):133answer.append(rspin_coeff(i,g,r,n,dvector,r_coeff,step,m0,deg,moduli_type))134return vector(answer)135136def rspin_degree_test(g,r,n=0,dvector=(),step=1,m0=-1,deg=-1,moduli_type=MODULI_ST):137answer = []138markings = tuple(range(1,n+1))139for i in range(num_strata(g,r,markings,moduli_type)):140answer.append(rspin_coeff(i,g,r,n,dvector,step,m0,deg,moduli_type).degree(X))141return answer142143def rspin_constant(g,r,n=0,dvector=(),step=1,m0=-1,deg=-1,moduli_type=MODULI_ST):144answer = []145markings = tuple(range(1,n+1))146for i in range(num_strata(g,r,markings,moduli_type)):147answer.append(rspin_coeff(i,g,r,n,dvector,step,m0,deg,moduli_type).subs(X=0))148return vector(answer)149150def Wittenrspin(g, Avector, r_coeff = None, d = None, rpoly = False):151r"""152Returns the polynomial limit of Witten's r-spin class in genus g with input Avector,153as discussed in the appendix of [Pandharipande-Pixton-Zvonkine '16].154155More precisely, Avector is expected to be a vector of linear polynomials in QQ[r],156with leading coefficients being rational numbers in the interval [0,1] and constant157coefficients being integers. Elements of Avector should sum to C * r + 2g-2 for some158nonnegative integer C. Then Wittenrspin returns the tautological class obtained as159the limit of160161r^(g-1+C) W_{g,n}^r(Avector)162163for r >> 0 sufficiently large and divisible, evaluated at r=0.164165INPUT::166167- ``g`` -- integer; underlying genus168169- ``Avector`` -- tuple ; a tuple of either integers, elements of a polynomial170ring QQ[r] in some variable r or tuples (C_i, D_i) which are interpreted as171polynomials C_i * r + D_i. For the entries with C_i = 1 we assume D_i <=-2.172173- ``r_coeff`` -- integer or None (default: `None`); if a particular integer174r_coeff = r is specified, the function will return the (unscaled) class W_{g,n}^r(Avector),175not taking a limit for large r.176177- ``d`` -- integer; desired degree in tautological ring; will be set to g-1+C by178default179180- ``rpoly`` -- bool (default: `False`); if True, return the limit of181r^(g-1+C) W_{g,n}^r(Avector) without evaluating at r=0, as a tautclass with182coefficients being polynomials in r183184EXAMPLES:185186We start by verifying the conjecture from the appendix of [Pandharipande-Pixton-Zvonkine '16]187for g = 2 and mu = (2)::188189sage: from admcycles import Wittenrspin, Strataclass190sage: H1 = Wittenrspin(2, (2,))191sage: H2 = Strataclass(2, 1, (2,))192sage: (H1-H2).is_zero()193True194195We can also verify a new conjecture for classes of strata of meromorphic differentials for196g = 1 and mu = (3,-1,-2). The argument (1/2,-1) stands for an insertion 1/2 * r -1::197198sage: H1 = Wittenrspin(1, (3, (1/2,-1), (1/2,-2)))199sage: H2 = Strataclass(1, 1, (3,-1,-2))200sage: (H1+H2).is_zero()201True202203As a variant of this, we also verify that insertions r-b stand for poles of order b with204vanishing residues::205206sage: R.<r> = PolynomialRing(QQ,1)207sage: H1 = Wittenrspin(1, (5,r-2, r-3))208sage: H2 = Strataclass(1, 1, (5,-2,-3), res_cond=(2,))209sage: (H1+H2).is_zero()210True211212We can also compute the (scaled) Witten's class without substituting r=0::213214sage: Wittenrspin(1,(2,r-2),rpoly=True)215Graph : [1] [[1, 2]] []216Polynomial : (1/12*r^2 - 5/24*r + 1/12)*(kappa_1^1 )_0217<BLANKLINE>218Graph : [1] [[1, 2]] []219Polynomial : (-1/12*r^2 + 29/24*r - 37/12)*psi_1^1220<BLANKLINE>221Graph : [1] [[1, 2]] []222Polynomial : (-1/12*r^2 + 17/24*r - 13/12)*psi_2^1223<BLANKLINE>224Graph : [0, 1] [[1, 2, 4], [5]] [(4, 5)]225Polynomial : (1/12*r^2 - 17/24*r + 13/12)*226<BLANKLINE>227Graph : [0] [[4, 5, 1, 2]] [(4, 5)]228Polynomial : (-1/48*r + 1/24)*229230Instead of calculating the asymptotic, polynomial behaviour of Witten's class,231we can also input a concrete value for r, using the option r_coeff. Below we232verify the CohFT property of Witten's class::233234sage: R.<r> = PolynomialRing(QQ,1)235sage: H1 = Wittenrspin(1, (5,r-2, r-3))236sage: H2 = Strataclass(1, 1, (5,-2,-3), res_cond=(2,))237sage: (H1+H2).is_zero()238True239"""240n = len(Avector)241242if r_coeff is None:243polyentries = [a for a in Avector if isinstance(a,MPolynomial)]244if len(polyentries) > 0:245R = polyentries[0].parent()246r = R.gens()[0]247else:248R = PolynomialRing(QQ,1,'r')249r = R.gens()[0]250251X = SR.var('X')252AvectorX = []253Cvector = []254Dvector = []255256# Extract entries of Avector into a unified format (AvectorX)257# Collect linear and constant coefficients of entries of AvectorX in Cvector, Dvector258for a in Avector:259if isinstance(a,Integer):260AvectorX.append(a)261Cvector.append(QQ(0))262Dvector.append(a)263elif isinstance(a,MPolynomial):264AvectorX.append(a[1]*X+a[0])265Cvector.append(QQ(a[1]))266Dvector.append(a[0])267elif isinstance(a,Iterable):268AvectorX.append(a[0]*X+a[1])269Cvector.append(QQ(a[0]))270Dvector.append(a[1])271else:272raise ValueError('Entries of Avector must be Integers, Polynomials or tuples (C_i, D_i)')273274# For C_i = 1, D_i = -1 we require a simple pole with vanishing residue => return zero class275if any((Cvector[i]==1) and (Dvector[i]==-1) for i in range(n)):276return tautclass([])277278step = lcm(c.denom() for c in Cvector)279C = sum(Cvector)280assert C in ZZ281assert sum(Dvector) == 2*g-2282else:283AvectorX=Avector284step=0285286if d is None:287if r_coeff is None:288d = ZZ(g-1+C)289else:290denom = ZZ((r_coeff-2)*(g-1) + sum(Avector))291if denom % r_coeff != 0:292return tautclass([]) # Witten's class vanishes unless this congruence is satisfied293else:294d = ZZ(denom/ZZ(r_coeff))295rvect = rspin_compute(g,d,n,tuple(AvectorX),r_coeff,step,m0=-1,deg=-1,moduli_type=MODULI_ST)296if not rpoly:297rvect = rvect.subs(X=0)298rvect = convert_vector_to_monomial_basis(rvect, g, d, tuple(range(1, n+1)), MODULI_ST)299300if rpoly:301rvect = vector((b.subs(X=r) for b in rvect))302303return Tautv_to_tautclass(rvect, g, n, d)304305@cached_function306def Pmpolynomial(m):307r"""308Returns the expression P_m(X,a) as defined in [Pandharipande-Pixton-Zvonkine '16, Section 4.5].309310TESTS::311312sage: from admcycles.witten import Pmpolynomial313sage: Pmpolynomial(0)3141315sage: Pmpolynomial(1)316-1/12*X^2 + 1/2*(X - 1)*a - 1/2*a^2 + 5/24*X - 1/12317sage: Pmpolynomial(2)3181/288*X^4 - 1/12*(5*X - 1)*a^3 + 1/8*a^4 + 7/288*X^3 + 1/48*(20*X^2 - 5*X - 4)*a^2 - 29/384*X^2 - 1/48*(6*X^3 + X^2 - 9*X + 2)*a + 7/288*X + 1/288319"""320(b,X,a) = var('b,X,a')321if m==0:322return 1+0*X323return ( QQ(1)/2 * symbolic_sum((2*m*X-X-2*b)*Pmpolynomial(m-1).subs(a=b-1),b,1,a) - QQ(1)/(4*m*X*(X-1)) * symbolic_sum((X-1-b)*(2*m*X-b)*(2*m*X-X-2*b)*Pmpolynomial(m-1).subs(a=b-1),b,1,X-2) ).simplify_rational()324325326327