| Download
Project: admcycles
Views: 724Visibility: Unlisted (only visible to those who know the link)
Image: ubuntu2004r"""1Double ramification cycle2"""3from __future__ import absolute_import, print_function45import itertools67from six.moves import range89from admcycles.admcycles import Tautv_to_tautclass, Tautvb_to_tautclass, Tautvecttobasis, tautgens, psiclass, tautclass10from admcycles.stable_graph import StableGraph11import admcycles.DR as DR1213from sage.combinat.subset import Subsets14from sage.arith.all import factorial15from sage.functions.other import floor, ceil16from sage.misc.misc_c import prod17from sage.rings.all import PolynomialRing, QQ, ZZ18from sage.modules.free_module_element import vector19from sage.rings.polynomial.polynomial_ring import polygen20from sage.rings.polynomial.multi_polynomial_element import MPolynomial2122############23#24# Old DR-cycle implementation25#26############272829def DR_cycle_old(g, Avector, d=None, k=None, tautout=True, basis=False):30r"""Returns the k-twisted Double ramification cycle in genus g and codimension d31for the partition Avector of k*(2g-2+n).3233In the notation of the [JPPZ]-paper, the output is 2^(-d) * P_g^{d,k}(Avector).3435Note: This is based on the old implementation DR_compute by Pixton. A new one,36which can be faster in many examples is DR_cycle.3738INPUT:3940- ``tautout` -- bool (default: `True`); if False, returns a vector41(in all generators for basis=false, in the basis of the ring for basis=true)4243- ``basis` -- bool (default: `False`); if True, use FZ relations to give out the44DR cycle in a basis of the tautological ring45"""46if d is None:47d = g48n = len(Avector)49if k is None:50k = floor(sum(Avector)/(2*g-2+n))51if sum(Avector) != k*(2*g-2+n):52raise ValueError('2g-2+n must divide the sum of entries of Avector.')5354v = DR.DR_compute(g, d, n, Avector, k)55v1 = vector([QQ(a) for a in DR.convert_vector_to_monomial_basis(v, g, d, tuple(range(1, n+1)), DR.MODULI_ST)])5657if basis:58v2 = Tautvecttobasis(v1, g, n, d)59if tautout:60return Tautvb_to_tautclass(v2, g, n, d)61return v262else:63if tautout:64return Tautv_to_tautclass(v1, g, n, d)65return v1666768############69#70# New DR-cycle implementation71#72############7374def DR_cycle(g, Avector, d=None, k=None, rpoly=False, tautout=True, basis=False):75r"""Returns the k-twisted Double ramification cycle in genus g and codimension d76for the partition Avector of k*(2g-2+n). If some elements of Avector are elements of a77polynomial ring, compute DR-polynomial and substitute the entries of Avector - the result78then has coefficients in a polynomial ring.7980In the notation of the [JPPZ]-paper, the output is 2^(-d) * P_g^{d,k}(Avector).8182Note: This is a reimplementation of Pixton's DR_compute which can be faster in many examples.83To access the old version, use DR_cycle_old - it might be faster on older SageMath-versions.8485INPUT:8687- ``rpoly`` -- bool (default: `False`); if True, return tautclass 2^(-d) * P_g^{d,r,k}(Avector)88whose coefficients are polynomials in the variable r.8990- ``tautout` -- bool (default: `True`); if False, returns a vector91(in all generators for basis=false, in the basis of the ring for basis=true)9293- ``basis` -- bool (default: `False`); if True, use FZ relations to give out the94DR cycle in a basis of the tautological ring9596EXAMPLES::9798sage: from admcycles import DR_cycle, DR_cycle_old99sage: DR_cycle(1, (2, 3, -5))100Graph : [1] [[1, 2, 3]] []101Polynomial : 2*psi_1^1102<BLANKLINE>103Graph : [1] [[1, 2, 3]] []104Polynomial : 9/2*psi_2^1105<BLANKLINE>106Graph : [1] [[1, 2, 3]] []107Polynomial : 25/2*psi_3^1108<BLANKLINE>109Graph : [0, 1] [[2, 3, 5], [1, 6]] [(5, 6)]110Polynomial : (-2)*111<BLANKLINE>112Graph : [0, 1] [[1, 3, 5], [2, 6]] [(5, 6)]113Polynomial : (-9/2)*114<BLANKLINE>115Graph : [0, 1] [[1, 2, 5], [3, 6]] [(5, 6)]116Polynomial : (-25/2)*117<BLANKLINE>118Graph : [0] [[5, 6, 1, 2, 3]] [(5, 6)]119Polynomial : (-1/24)*120sage: D = DR_cycle(1, (1, 3, -4), d=2)121sage: D2 = DR_cycle_old(1, (1, 3, -4), d=2) # long time122sage: D.toTautvect() == D2.toTautvect() # long time123True124sage: D.is_zero() # Clader-Janda: P_g^{d,k}(Avector)=0 for d>g125True126sage: v = DR_cycle(2, (3,), d=4, rpoly=true).evaluate()127sage: v # Conjecture by Longting Wu predicts that r^1-term vanishes128-1/1728*r^6 + 1/576*r^5 + 37/34560*r^4 - 13/6912*r^3 - 1/1152*r^2129130We can create a polynomial ring and use an Avector with entries in this ring.131The result is a tautological class with polynomial coefficients::132133sage: R.<a1,a2>=PolynomialRing(QQ,2)134sage: Q=DR_cycle(1,(a1,a2,-a1-a2))135sage: Q136Graph : [1] [[1, 2, 3]] []137Polynomial : 1/2*a1^2*psi_1^1138<BLANKLINE>139Graph : [1] [[1, 2, 3]] []140Polynomial : 1/2*a2^2*psi_2^1141<BLANKLINE>142Graph : [1] [[1, 2, 3]] []143Polynomial : (1/2*a1^2 + a1*a2 + 1/2*a2^2)*psi_3^1144<BLANKLINE>145Graph : [0, 1] [[2, 3, 5], [1, 6]] [(5, 6)]146Polynomial : (-1/2*a1^2)*147<BLANKLINE>148Graph : [0, 1] [[1, 3, 5], [2, 6]] [(5, 6)]149Polynomial : (-1/2*a2^2)*150<BLANKLINE>151Graph : [0, 1] [[1, 2, 5], [3, 6]] [(5, 6)]152Polynomial : (-1/2*a1^2 - a1*a2 - 1/2*a2^2)*153<BLANKLINE>154Graph : [0] [[5, 6, 1, 2, 3]] [(5, 6)]155Polynomial : (-1/24)*156157158TESTS::159160sage: from admcycles import DR_cycle, DR_cycle_old161sage: D = DR_cycle(2, (3, 4, -2), k=1) # long time162sage: D2 = DR_cycle_old(2, (3, 4, -2), k=1) # long time163sage: D.toTautvect() == D2.toTautvect() # long time164True165sage: D = DR_cycle(2, (1, 4, 5), k=2) # long time166sage: D2 = DR_cycle_old(2, (1, 4, 5), k=2) # long time167sage: D.toTautvect() == D2.toTautvect() # long time168True169sage: D = DR_cycle(3, (2,4), k=1) # long time170sage: D2 = DR_cycle_old(3, (2,4), k=1) # long time171sage: D.toTautvect() == D2.toTautvect() # long time172True173174sage: D = DR_cycle(2, (3,),tautout=False); D175(0, 1/8, -9/4, 81/8, 1/4, -9/4, -1/8, 1/8, 1/4, -1/8, 1/48, 1/240, -3/16, -41/240, 1/48, 1/48, 1/1152)176sage: D2=DR_cycle_old(2, (3,), tautout=False); D2177(0, 1/8, -9/4, 81/8, 1/4, -9/4, -1/8, 1/8, 1/4, -1/8, 1/48, 1/240, -3/16, -41/240, 1/48, 1/48, 1/1152)178sage: D3 = DR_cycle(2, (3,), tautout=False, basis=True); D3179(-5, 2, -8, 18, 3/2)180sage: D4 = DR_cycle_old(2, (3,), tautout=False, basis=True); D4181(-5, 2, -8, 18, 3/2)182"""183if d is None:184d = g185n = len(Avector)186187if any(isinstance(ai, MPolynomial) for ai in Avector):188# compute DRpoly and substitute the ai in there189pol, gens = DRpoly(g, d, n, tautout=tautout, basis=basis, gensout=True)190subsdict = {gens[i]: Avector[i] for i in range(n)}191192if tautout:193pol.coeff_subs(subsdict)194pol.simplify()195return pol196return pol.apply_map(lambda x: x.subs(subsdict))197198if k is None:199k = floor(sum(Avector)/(2*g-2+n))200if sum(Avector) != k*(2*g-2+n):201raise ValueError('2g-2+n must divide the sum of entries of Avector.')202203gens = tautgens(g, n, d, decst=True)204v = vector([DR_coeff_new(decs, g, n, d, Avector, k, rpoly) for decs in gens]) / ZZ(2)**d205#v1=vector([QQ(a) for a in DR.convert_vector_to_monomial_basis(v,g,d,tuple(range(1, n+1)),DR.MODULI_ST)])206207if basis:208v1 = Tautvecttobasis(v, g, n, d)209if tautout:210return Tautvb_to_tautclass(v1, g, n, d)211return v1212else:213if tautout:214return Tautv_to_tautclass(v, g, n, d)215return v216217218def interpolate(A, B):219r"""220Univariate Lagrange interpolation over the rationals.221222EXAMPLES::223224sage: from admcycles.double_ramification_cycle import interpolate225sage: p = interpolate([1/2, -2, 3], [4/5, 2/3, -7/6])226sage: p(1/2)2274/5228sage: p(-2)2292/3230sage: p(3)231-7/6232233TESTS::234235sage: from admcycles.double_ramification_cycle import interpolate236sage: parent(interpolate([], []))237Univariate Polynomial Ring in r over Rational Field238"""239R = polygen(QQ, 'r').parent()240return R.lagrange_polynomial(zip(A, B))241242243def DR_coeff_setup(G, g, n, d, Avector, k):244gamma = G.gamma245kappa, psi = G.poly.monom[0]246exp_list = []247scalar_factor = QQ((1, G.automorphism_number()))248given_weights = [-k * (2 * gv - 2 + len(gamma._legs[i]))249for i, gv in enumerate(gamma._genera)]250251# contributions to scalar_factor from kappa-classes252for v in range(gamma.num_verts()):253if len(kappa[v]) == 1: # some kappa_1^e at this vertex254scalar_factor *= (-k**2)**(kappa[v][0]) / factorial(kappa[v][0])255256# contributions to scalar_factor and given_weights from markings257for i in range(1, n + 1):258v = gamma.vertex(i)259psipow = psi.get(i, 0) # if i in dictionary, have psi_i^(psi[i]), otherwise have psi_i^0260given_weights[v] += Avector[i - 1]261scalar_factor *= (Avector[i - 1])**(2 * psipow) / factorial(psipow)262263# contributions to scalar_factor and explist from edges264for (lu, lv) in gamma._edges:265psipowu = psi.get(lu, 0)266psipowv = psi.get(lv, 0)267exp_list.append(psipowu + psipowv + 1)268scalar_factor *= ((-1)**(psipowu+psipowv)) / factorial(psipowv) / factorial(psipowu) / (psipowu+psipowv+1)269270return exp_list, given_weights, scalar_factor271272273def DR_coeff_new(G, g, n, d, Avector, k, rpoly):274gamma = G.gamma # underlying stable graph of decstratum G275kappa, _ = G.poly.monom[0]276# kappa is a list of lenght = # vertices, of entries like [3,0,2] meaning that at this vertex there is a kappa_1^3 * kappa_3^2277# _ = psi is a dictionary and psi[3]=4 means there is a decoration psi^4 at marking/half-edge 3278279if any(len(kalist) > 1 for kalist in kappa):280return 0 # vertices can only carry powers of kappa_1, no higher kappas allowed281282m0 = ceil(sum([abs(i) for i in Avector]) / ZZ(2)) + g*abs(k) # value from which on the Pixton-sum below will be polynomial in m283h0 = gamma.num_edges() - gamma.num_verts() + 1 # first Betti number of the graph Gamma284exp_list, given_weights, scalar_factor = DR_coeff_setup(G, g, n, d, Avector, k)285286deg = 2 * sum(exp_list) # degree of polynomial in m287# R = PolynomialRing(QQ, len(gamma.edges) + 1, 'd')288# P=prod([(di*(d0-di))**exp_list[i] for i,di in enumerate(R.gens()[1:])])/(d0**h0)289u = gamma.flow_solve(given_weights)290Vbasis = gamma.cycle_basis()291292#mpoly = mpoly_affine_sum(P, u, Vbasis,m0+1,deg)293mpoly = mpoly_special(exp_list, h0, u, Vbasis, m0+1, deg)294295# mpoly = interpolate(mrange, mvalues)296if rpoly:297return mpoly * scalar_factor298return mpoly.subs(r=0) * scalar_factor299300301def mpoly_special(exp_list, h0, u, Vbasis, start, degr):302r"""303Return the sum of the rational function in DR_coeff_new over the304affine space u + V modulo r in ZZ^n as a univariate polynomial in r.305V is the lattice with basis Vbasis.306Use that this sum is polynomial in r starting from r=start and the polynomial has307degree degr (for now).308"""309mrange = list(range(start, start + degr + 2))310mvalues = []311rank = len(Vbasis)312ulen = len(u)313314for m in mrange:315total = 0316for coeff in itertools.product(*[list(range(m)) for i in range(rank)]):317v = u + sum([coeff[i]*Vbasis[i] for i in range(rank)])318v = [vi % m for vi in v]319total += prod([(v[i]*(m-v[i]))**exp_list[i] for i in range(ulen)])320mvalues.append(total/QQ(m**h0))321322return interpolate(mrange, mvalues)323324325def mpoly_affine_sum(P, u, Vbasis, start, degr):326r"""327Return the sum of the polynomial P in variables r, d1, ..., dn over the328affine space u + V modulo r in ZZ^n as a univariate polynomial in r.329V is the lattice with basis Vbasis.330Use that this sum is polynomial in r starting from r=start and the polynomial has331degree degr (for now).332"""333mrange = list(range(start, start + degr + 2))334mvalues = []335rank = len(Vbasis)336337for m in mrange:338total = 0339for coeff in itertools.product(*[list(range(m)) for i in range(rank)]):340v = u + sum([coeff[i] * Vbasis[i] for i in range(rank)])341v = [vi % m for vi in v]342total += P(m, *v)343mvalues.append(total)344345return interpolate(mrange, mvalues)346347############348#349# DR polynomial350#351############352353354def multivariate_interpolate(f, d, n, gridwidth=1, R=None, generator=None):355r"""Takes a vector/number-valued function f on n variables and interpolates it on a grid around 0.356Returns a vector with entries in a polynomial ring.357358INPUT:359360- ``d`` -- integer; maximal degree of output-polynomial in any of the variables361362- ``gridwidth` -- integer (default: `1`); width of the grid used for interpolation363364- ``R` -- polynomial ring (default: `None`); expected to be polynomial ring in365at least n variables; if None, use ring over QQ in variables z0,...,z(n-1)366367- ``generator` -- list (default: `None`); list of n variables in the above polynomial368ring to be used in output; if None, use first n generators of ring369370EXAMPLES::371372sage: from admcycles.double_ramification_cycle import multivariate_interpolate373sage: from sage.modules.free_module_element import vector374sage: def foo(x,y):375....: return vector((x**2+y,2*x*y))376....:377sage: multivariate_interpolate(foo,2,2)378(z0^2 + z1, 2*z0*z1)379"""380if R is None:381R = PolynomialRing(QQ, 'z', n)382if generator is None:383generator = R.gens()384385if n == 0: # return constant386return R.one() * f()387388cube = [list(range(d + 1))] * n389points = itertools.product(*cube)390# shift cube containing evaluation points in negative direction to reduce abs value of evaluation points391shift = floor((d + 1) / QQ(2)) * gridwidth392result = 0393394for p in points:395# compute Lagrange polynomial not vanishing exactly at gridwidth*p396lagr = prod([prod([generator[i]-(j*gridwidth-shift) for j in range(d+1) if j != p[i]]) for i in range(n)])397pval = [gridwidth*coord-shift for coord in p]398#global fex,lagrex, pvalex399#fex=f400#lagrex=lagr401#pvalex=pval402value = lagr.subs({generator[i]: pval[i] for i in range(n)})403if n == 1: # avoid problems with multiplication of polynomial with vector404mult = lagr/value405result += vector((e*mult for e in f(*pval)))406else:407result += f(*pval) / value * lagr408return result409410411def DRpoly(g, d, n, dplus=0, tautout=True, basis=False, ring=None, gens=None, gensout=False):412r"""413Return the Double ramification cycle in genus g with n markings in degree d as a414tautclass with polynomial coefficients. Evaluated at a point (a_1, ..., a_n) with415sum a_i = k(2g-2+n) it equals DR_cycle(g,(a_1, ..., a_n),d).416417The computation uses interpolation and the fact that DR is a polynomial in the a_i418of degree 2*d.419420INPUT:421422- ``dplus`` -- integer (default: `0`); if dplus>0, the interpolation is performed423on a larger grid as a consistency check424425- ``tautout` -- bool (default: `True`); if False, returns a polynomial-valued vector426(in all generators for basis=false, in the basis of the ring for basis=true)427428- ``basis` -- bool (default: `False`); if True, use FZ relations to give out the429DR cycle in a basis of the tautological ring430431- ``ring` -- polynomial ring (default: `None`); expected to be polynomial ring in432at least n variables; if None, use ring over QQ in variables a1,...,an433434- ``gens` -- list (default: `None`); list of n variables in the above polynomial435ring to be used in output; if None, use first n generators of ring436437- ``gensout` -- bool (default: `False`); if True, return (DR polynomial, list of generators438of coefficient ring)439440EXAMPLES::441442sage: from admcycles import DRpoly, DR_cycle, psiclass443sage: D, (a1, a2) = DRpoly(1, 1, 2, gensout=True)444sage: D.coeff_subs({a2:-a1}).simplify()445Graph : [1] [[1, 2]] []446Polynomial : 1/2*a1^2*psi_1^1447<BLANKLINE>448Graph : [1] [[1, 2]] []449Polynomial : 1/2*a1^2*psi_2^1450<BLANKLINE>451Graph : [0] [[4, 5, 1, 2]] [(4, 5)]452Polynomial : (-1/24)*453sage: (D*psiclass(1,1,2)).evaluate() # intersection number with psi_14541/24*a1^2 - 1/24455sage: (D.coeff_subs({a1:4})-DR_cycle(1,(4,-4))).is_zero() # polynomial agrees with DR_cycle at (4,-4)456True457sage: DRpoly(1,2,2).is_zero() # Clader-Janda: DR vanishes in degree d>g458True459460TESTS::461462sage: R.<a1,a2,a3,b1,b2,b3> = PolynomialRing(QQ, 6)463sage: D = DRpoly(1, 1, 3, ring=R, gens=[a1, a2, a3])464sage: Da = deepcopy(D)465sage: Db = deepcopy(D).coeff_subs({a1:b1, a2:b2, a3:b3})466sage: Dab = deepcopy(D).coeff_subs({a1:a1+b1, a2:a2+b2, a3:a3+b3})467sage: diff = Da*Db - Da*Dab # this should vanish in compact type by [Holmes-Pixton-Schmitt]468sage: diff.toTautbasis(moduli='ct')469(0)470"""471def f(*Avector):472return DR_cycle(g, Avector, d, tautout=False, basis=basis)473#k=floor(QQ(sum(Avector))/(2 * g - 2 + n))474#return vector([QQ(e) for e in DR_red(g,d,n,Avector, k, basis)])475#return vector(DR_compute(g,d,n,Avector, k))476477if ring is None:478ring = PolynomialRing(QQ, ['a%s' % i for i in range(1, n + 1)])479if gens is None:480gens = ring.gens()[0:n]481482gridwidth = 2*g - 2 + n483interp = multivariate_interpolate(f, 2*d + dplus, n, gridwidth, R=ring, generator=gens)484485if not tautout:486ans = interp487else:488if not basis:489ans = Tautv_to_tautclass(interp, g, n, d)490else:491ans = Tautvb_to_tautclass(interp, g, n, d)492if gensout:493return (ans, gens)494return ans495496497def degree_filter(polyvec, d):498r"""Takes a vector of polynomials in several variables and returns a vector containing the (total)499degree d part of these polynomials.500501INPUT:502503- vector of polynomials504505- integer degree506507EXAMPLES::508509sage: from admcycles.double_ramification_cycle import degree_filter510sage: R.<x,y> = PolynomialRing(QQ, 2)511sage: v = vector((x+y**2+2*x*y,1+x**3+x*y))512sage: degree_filter(v, 2)513(2*x*y + y^2, x*y)514"""515resultvec = []516for pi in polyvec:517s = 0518for c, m in pi:519if m.degree() == d:520s += c * m521resultvec.append(s)522return vector(resultvec)523524525def Hain_divisor(g, A):526r"""527Returns a divisor class D extending the pullback of the theta-divisor under the Abel-Jacobi map (on compact type) given by partition A of zero.528529Note: D^g/g! agrees with the Double-Ramification cycle in compact type.530531EXAMPLES::532533sage: from admcycles import *534535sage: R = PolynomialRing(QQ, 'z', 3)536sage: z0, z1, z2 = R.gens()537sage: u = Hain_divisor(2, (z0, z1, z2))538sage: g = DRpoly(2, 1, 3, ring=R, gens=[z0, z1, z2]) #u,g should agree inside compact type539sage: (u.toTautvect() - g.toTautvect()).subs({z0: -z1-z2})540(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1/24)541"""542n = len(A)543544def delt(l, g, n, I):545# takes genus l and subset I of markings 1,...,n and returns generalized boundary divisor Delta_{l,I}546I = set(I)547if l == 0 and len(I) == 1:548return -psiclass(list(I)[0], g, n)549if l == g and len(I) == n-1:550i_set = set(range(1, n+1)) - I551return -psiclass(list(i_set)[0], g, n)552if (l == 0 and len(I) == 0) or (l == g and len(I) == n):553return tautclass([])554555Icomp = set(range(1, n+1)) - I556gra = StableGraph([l, g-l], [list(I) + [n+1], list(Icomp) + [n+2]],557[(n+1, n+2)])558return QQ((1, gra.automorphism_number())) * gra.to_tautclass()559560result = sum([sum([(sum([A[i - 1] for i in I])**2) * delt(l, g, n, I)561for I in Subsets(list(range(1, n + 1)))])562for l in range(g + 1)])563# note index i-1 since I subset {1,...,n} but array indices subset {0,...,n-1}564565return QQ((-1, 4)) * result566567568