| Download
Project: admcycles
Views: 725Visibility: Unlisted (only visible to those who know the link)
Image: ubuntu2004r"""1Recursive computation of strata of k-differentials23This is following [Farkas-Pandharipande; Schmitt]4"""5from __future__ import absolute_import, print_function6from six.moves import range7from admcycles.admcycles import (Tautv_to_tautclass, fundclass, tautclass)8from admcycles.stable_graph import StableGraph9from admcycles.double_ramification_cycle import DR_cycle10from admcycles.diffstrata.sig import Signature11from admcycles.diffstrata.levelstratum import GeneralisedStratum1213import itertools14from copy import deepcopy15from collections import Iterable1617from sage.misc.misc import subsets18from sage.combinat.integer_vector import IntegerVectors19from sage.combinat.partition import Partitions20from sage.misc.misc_c import prod21from sage.rings.all import ZZ22from sage.misc.cachefunc import cached_function23from sage.combinat.words.word import Word24from sage.functions.other import ceil2526# Given a genus g, an order k and a partition mu of k(2g-2), twistenum27# returns a list of graphs with twists2829# Elements of this list are of the following form30# [31# [g_0, [markings of v_0] ],32# [g_1, [markings of v_1], #edges of v_1, [twists of edges to v_1] ],33# ...,34# [g_r, [markings of v_r], #edges of v_r, [twists of edges to v_r] ]35# ]36# NOTE: Twists are counted in multiples of k, so a twist of 3 actually means 3k3738# ordering of graphs (always descending)39# genus g0 of central vertex > outlying vertex data40# ordering of outlying vertices41# total genus > genus > total edge twist424344def classes(l):45"""46INPUT:4748- l -- a list4950EXAMPLES::5152sage: from admcycles.stratarecursion import classes53sage: classes([])54[]55sage: classes([4,4,3,3,3,1,1])56[[0, 1], [2, 3, 4], [5, 6]]57"""58if not l:59return []60indices = []61current = [0]62currentitem = l[0]63for i in range(1, len(l)):64if l[i] == currentitem:65current += [i]66else:67indices += [current]68current = [i]69currentitem = l[i]70indices += [current]71return indices727374def SetPart(S, s):75"""76Return a list of partitions of the set S into the disjoint union of s sets.7778A partition is a list of sets.7980EXAMPLES::8182sage: from admcycles.stratarecursion import SetPart83sage: S = set(['b','o','f'])84sage: SetPart(S,2)85[[set(), {'b', 'f', 'o'}],86...87[{'b', 'f', 'o'}, set()]]88sage: SetPart(S,1)89[[{'b', 'f', 'o'}]]90sage: SetPart(S,-1)91[]92sage: SetPart(S,0)93[]94"""95if s < 0:96return []97if s == 0:98if S:99return []100else:101return [[]]102if s == 1:103return [[S]]104105resu = []106for T in subsets(S):107for r in SetPart(S - set(T), s - 1):108resu.append([set(T)] + r)109return resu110111112def twistenum(g, k, mu):113k = ZZ(k)114lis = [] # list of all graphs, returned in the end115116# collect all negative or non-divisible-by-k markings in A, all others in B117A = []118B = []119for i in range(len(mu)):120if mu[i] < 0 or (mu[i] % k != 0):121A.append(i)122else:123B.append(i)124Atot = sum([mu[i] for i in A])125126for g0 in range(g+1): # g0 = genus of center vertex127# iterate over possible partitions of the rest of the genera to vertices128for totgenlist in Partitions(g - g0):129# list containing lists of indices with equal total genus130cl = classes(totgenlist)131132# we can have 1 up to total_genus many edges, distribute the len(c) many possibilities to those133# each of the len(c) vertices can have 1, ..., total_genus many edges; vector (1,0,2) means134# 1 vertex has 1 edge, no vertex has 2 edges, 2 vertices have 3 edges135par = [IntegerVectors(len(c), totgenlist[c[0]]) for c in cl]136137for p in itertools.product(*par):138# gra contains the information that will eventually be added to lis139gra = [[g0, A]]140numed = 0 # total number of edges141for j in range(len(cl)):142# p[j] is now an integer vector, encoding edge numbers of vertices c in cl[j] with totgenlist[c[0]] total genus143for z in range(len(p[j])):144gra += [[totgenlist[cl[j][0]] - z, [], z + 1, []]145for _ in range(p[j][z])]146numed += p[j][z]*(z+1)147148# iterate over additional markings sent to v_0149for markings0 in subsets(B):150# vertex v0 unstable151if g0 == 0 and numed+len(A)+len(markings0) <= 2:152continue153Btot = sum([mu[i] for i in markings0])154if not k.divides(Atot+Btot):155continue156157Brest = [b for b in B if b not in markings0]158159Itotal = (Atot + Btot) // k - 2*numed - (2*g0-2)160161# Itotal is the total number of ADDITIONAL twists k that can be distributed to the various edges162# (all have k as default anyway, so each edge has at least a pole of order 2*k on the central vertex)163164# TODO: break if Itotal<0165166vertcl = classes(gra)167168for coarsesttwistdist in IntegerVectors(Itotal, len(vertcl)-1):169# distribute the additional twists on the classes of outlying vertices170twist = []171# twist will be a list for every class of outlying vertex recording a list of all possible twist->edge distr.172# element of the list for the class {3,4,5} have the form173# [[a,b],[c],[d,e,f]] if v_3 has two edges, v_4 one and v_5 three174175for i in range(len(vertcl)-1):176# for each of the classes collect the possibilities to distribute inside this class177vertcltwistlist = []178# will be added as an element to twist179for coarsetwistdist in Partitions(coarsesttwistdist[i]+len(vertcl[i+1]), length=len(vertcl[i+1])):180# coarsetwistdist describes how twists are distributed to the vertices in vertcl[i+1]181# artificially added ones to be able to use Partitions182183inditwist = [Partitions(184coarsetwistdist[v]+gra[vertcl[i+1][v]][2]-1, length=gra[vertcl[i+1][v]][2]) for v in range(len(vertcl[i+1]))]185# list for every vertex in the class giving all possible twists of the edges of this vertex186187vertcltwistlist += itertools.product(188*inditwist)189190twist.append(vertcltwistlist)191for inde in itertools.product(*twist):192# inde is now of the form [ [[a,b], [c], [d,e,f]] , [[g], [h,i]], ... ]193grap = deepcopy(gra)194# this will be updated now with the edge twists determined above and the markings sent to v_0195grap[0][1] += markings0196count = 1197for i in inde:198for j in i:199grap[count][3] = j200count += 1201202twicla = classes(grap)203# print((twicla,len(twicla)-1,Brest))204205for pa in SetPart(set(Brest), len(twicla)-1):206mpar = [SetPart(pa[c], len(twicla[c+1]))207for c in range(len(twicla)-1)]208if not mpar and Brest:209continue210for part in itertools.product(*mpar):211# part is now of the form [ [ set([0,2]), set([1]) ] , [ set([3]) ] , ... ]212graph = deepcopy(grap)213count = 1214adm = True # check if graph satisfies twist conditions at all vertices215for i in part:216for j in i:217# we are now at the vertex v_count and need to check if the twist-condition is satisfied218graph[count][1] = list(j)219if (2*graph[count][0]-2)*k + k*sum(-l+1 for l in graph[count][3]) != sum(mu[l] for l in graph[count][1]):220adm = False221count += 1222if adm:223sgraph = (224(graph[0][0], tuple(sorted(m+1 for m in graph[0][1]))),)225sgraph += tuple(sorted(((gv, tuple(sorted(m+1 for m in marki)), tuple(226sorted(k*t for t in etwist))) for gv, marki, enu, etwist in graph[1:])))227lis.append(sgraph)228229return list(set(lis))230231232def Strataclass(g, k, mu, virt=False, res_cond = (), xi_power = 0, method = 'pull'):233r"""234Returns the fundamental class of the closure of the stratum of k-differentials in genus g in Mbar_{g,n}235with vanishing and pole orders mu.236237INPUT:238239- ``g`` -- integer ; genus of the curves240241- ``k`` -- integer ; power of the canonical line bundle in definition of stratum242243- ``mu`` -- tuple ; tuple of integers of length n giving zero and pole multiplicities244of k-differential, required to sum up to k*(2g-2)245246- ``virt`` -- bool (default: `False`); if True, k=1 and all entries of mu nonnegative, this247computes the virtual codimension g class supported on the codimension g-1 stratum of248holomorphic 1-differentials.249250- ``res_cond`` -- tuple (default: `()`); tuple of residue conditions. Each entry of251res_cond can be of one of the following two types:252253- an integer i from 1, ..., n indicating a marking pi with mu[i]<0, such that the254differential eta is required to have a pole with vanishing residue at the marking.255256- a tuple (c1, ..., cn) of rational numbers indicating that a condition257c1 * Res_{p1}(eta) + ... + cn * Res_{pn}(eta) = 0258is imposed. Currently only implemented for ci in {0,1}.259260The function then computes the class of the closure of the locus of smooth curves261having such a differential eta. Currently only implemented for k=1.262263- ``xi_power`` -- integer (default: `0`); if positive, returns the pushforward of264the corresponding power of the first Chern class xi = c_1(O(-1)) of the tautological265bundle on the moduli space of multi-scale differentials from [BCGGM3].266Currently only implemented for k=1 and with method = 'diffstrata'.267268- ``method`` -- string (default: `'pull'`); when computing a stratum of 1-differentials269with residue conditions, there are two choices here: 'pull' will compute it via boundary270pullbacks of higher genus strata, 'diffstrata' will use the package `diffstrata`, which271iteratively replaces residue conditions with equivalent divisorial conditions. The two272results should be equal.273274NOTE::275276The class is computed using a formula from papers by Farkas-Pandharipande and Schmitt,277proven by [Holmes-Schmitt 19] and [Bae-Holmes-Pandharipande-Schmitt-Schwarz 20].278The formula for differentials with residue conditions is based on unpublished work279relying on [Bainbridge-Chen-Gendron-Grushevsky-Moeller 16].280281WARNING::282283Imposing residue conditions at poles of high order leads to very long computations,284since the method works by computing strata of differentials on high-genus moduli285spaces.286287TESTS::288289sage: from admcycles import Hyperell, Biell290sage: from admcycles.stratarecursion import Strataclass291sage: L=Strataclass(2,1,(3,-1)); L.is_zero()292True293sage: L=Strataclass(3,1,(5,-1)); L.is_zero() # doctest: +SKIP294True295sage: L=Strataclass(2,1,(2,)); (L-Hyperell(2,1)).is_zero()296True297298In g=2, the locus Hbar_2(2,2) decomposes into the codimension 1 set of Hbar_2(1,1) and the299codimension 2 set of curves (C,p,q) with p,q Weierstrass points. The latter is equal to the cycle300Hyperell(2,2). We can obtain it by subtracting the virtual cycle for the partition (1,1) from the301virtual cycle for the partition (2,2)::302303sage: H1 = Strataclass(2, 2, (2, 2), virt = True)304sage: H2 = Strataclass(2, 1, (1, 1), virt = True)305sage: T = H1 - H2 - Hyperell(2, 2)306sage: T.is_zero()307True308309In g=1, the locus Hbar_1(2,-2) is the locus of curves (E,p,q) with p-q being 2-torsion in E.310Equivalently, this is the locus of bielliptic curves with a pair of bielliptic conjugate points::311312sage: (Strataclass(1,1,(2,-2)) - Biell(1,0,1)).is_zero()313True314315Some tests of computations involving residue conditions::316317sage: from admcycles import Strataclass318sage: OmegaR = Strataclass(1,1,(6,-4,-2),res_cond=(3,))319sage: OmegaRalt = Strataclass(1,1,(6,-4,-2),res_cond=(2,)) # long time320sage: (OmegaR - OmegaRalt).is_zero() # long time321True322sage: (OmegaR.forgetful_pushforward([2,3])).fund_evaluate()32342324sage: a=4; (a+2)**2 + a**2 - 10 # formula from [Castorena-Gendron, Cor. 5.5]32542326sage: OmegaR2 = Strataclass(1,1,(4,-2,-2),res_cond=(3,))327sage: (OmegaR2.forgetful_pushforward([2,3])).fund_evaluate()32810329sage: OmegaR3 = Strataclass(1,1,(5,-3,-2),res_cond=(2,)) # not tested330sage: (OmegaR3.forgetful_pushforward([2,3])).fund_evaluate() # not tested33124332sage: a=3; (a+2)**2 + a**2 - 10 # formula from [Castorena-Gendron, Cor. 5.5] # not tested33324334sage: OmegaR5 = Strataclass(2,1,(5,-1,-2),res_cond=(3,)) # not tested335sage: OmegaR5.is_zero() # vanishes by residue theorem # not tested336True337338We can also check that the two ways of computing residue conditions (via pullbacks and339via the package diffstrata) coincide::340341sage: a = Strataclass(1,1,(4,-2,-2), res_cond=(2,))342sage: b = Strataclass(1,1,(4,-2,-2), res_cond=(2,), method='diffstrata')343sage: (a-b).is_zero()344True345346The following computes the locus of genus 1 curves admitting a differential with multiplicity347vector (8,-2,-2,-2,-2) at the markings such that the sum of residues at p2 and p3 equals zero::348349sage: c = Strataclass(1,1,(8,-2,-2,-2,-2), res_cond=((0,1,1,0,0),))350"""351n = len(mu)352k = ZZ(k)353if sum(mu) != k*(2*g-2):354raise ValueError('mu must be a partition of k*(2g-2).')355356357358if method == 'diffstrata' or xi_power > 0 or any(isinstance(rc,Iterable) for rc in res_cond):359# preprocessing residue conditions360fancy_res_cond = []361for rc in res_cond:362if isinstance(rc,Iterable):363if not all(a == 0 or a == 1 for a in rc):364raise NotImplementedError('Only residue conditions with coefficients 0,1 implemented.')365fancy_res_cond.append([(0,i) for i, a in enumerate(rc) if a == 1])366else:367fancy_res_cond.append([(0,rc-1)])368369X = GeneralisedStratum(sig_list = [Signature(tuple(mu))], res_cond = fancy_res_cond)370return (X.xi**xi_power).to_prodtautclass().pushforward()371372373if len(res_cond)>0:374# res_cond = tuple(res_cond) # make sure it is of type tuple375if not k==1:376raise NotImplementedError('Residue conditions only implemented for k=1')377if virt:378raise ValueError('Residue conditions not compatible with virt=True')379380res_cond = sorted(res_cond)381poles = [i for (i,mui) in enumerate(mu) if mui<0] # indices of poles382383if len(poles)==1: # residue conditions are automatic384return Strataclass(g, k, mu, virt=virt)385386if any(mu[i-1]>=0 for i in res_cond): # try imposing residue cond. at non-pole387raise ValueError('Residue conditions can only be imposed at poles')388if any(mu[i-1]==-1 for i in res_cond): # try imposing residue cond. at simple pole389return tautclass([])390391num_add_marks = len([1 for i in res_cond if mu[i-1]%2]) # additional markings392maxleg = n + num_add_marks+2393394genera = [g] + [ceil(-mu[i-1]/2) for i in res_cond]395munew = []396outermus = []397markcounter = 1398legs = [[]]399edges = [(maxleg+2*j, maxleg+2*j+1) for j in range(len(res_cond))]400401for i in range(1,n+1):402if i not in res_cond:403legs[0].append(markcounter)404munew.append(mu[i-1])405markcounter+=1406else:407legs[0].append(maxleg)408if mu[i-1]%2 == 0:409legs.append([maxleg+1])410outermus.append([-mu[i-1]-2])411else:412legs.append([maxleg+1,markcounter])413outermus.append([-mu[i-1]-2,1])414munew.append(1)415markcounter+=1416maxleg+=2417gamma = StableGraph(genera, legs, edges)418419outerg = sum(genera)420outerclass = Strataclass(outerg, k, munew, virt=False)421422pullb = gamma.boundary_pullback(outerclass)423return pullb.factor_reconstruct(0,[Strataclass(genera[j+1],1,outermus[j]) for j in range(len(res_cond))])424425if (g == 0) or all(m == 0 for m in mu):426return fundclass(g, n)427428# Hbar has codimension g429meromorphic = any(not k.divides(m) or m < 0 for m in mu)430431if k > 1 and not meromorphic and not virt:432# all entries divisible by k and nonnegative AND user wants a codim g-1 class433# return the corresponding stratum of abelian differentials434return Strataclass(g, 1, tuple(m // k for m in mu))435436ordering_permutation = Word(mu).standard_permutation().inverse()437ordering_dict = {i+1: j for i, j in enumerate(ordering_permutation)}438# this is the dictionary such that renaming the answer with standard-ordering according to ordering_dict gives correct answer439sortmu = tuple(sorted(mu))440441try:442v = StrataDB.cached(g, k, sortmu, virt)443if meromorphic or virt: # return codim g class444ans_preord = Tautv_to_tautclass(v, g, n, g)445else:446ans_preord = Tautv_to_tautclass(v, g, n, g-1)447ans_preord.rename_legs(ordering_dict, rename=True)448return ans_preord449except KeyError:450pass451452# at this point, we actually have to compute453# we will compute the answer with standard ordering, store its vector in StrataDB and then return the renamed answer454if meromorphic or virt:455# return codimension g class456bdry_list = twistenum(g, k, sortmu)457indfind1 = tuple((i for i, l in enumerate(bdry_list) if l[0][0] == g))458assert len(indfind1) == 1, (g, k, mu, tuple(indfind1))459bdry_list.pop(indfind1[0]) # remove the trivial graph460461# right hand side of recursion462result = DR_cycle(g, tuple((m+k for m in sortmu)))463result -= Strataboundarysum(g, k, sortmu, bdry_list)464465v = result.toTautvect(g, n, g)466StrataDB.set_cache(v, *(g, k, sortmu, virt))467result = Tautv_to_tautclass(v, g, n, g) # gives simplified version of result468result.rename_legs(ordering_dict, rename=True)469return result470else:471# return codimension g-1 class472# by previous check can assume that k=1473assert k == 1, (g, k, mu)474475sortmuprime = list(sortmu) + [-1]476sortmuprime[n-1] += 1477sortmuprime = tuple(sortmuprime)478479bdry_list = twistenum(g, k, sortmuprime)480481indfind1 = tuple((i for i, l in enumerate(bdry_list) if l[0][0] == g))482assert len(indfind1) == 1, (g, k, mu, tuple(indfind1))483bdry_list.pop(indfind1[0]) # remove the trivial graph484485indfind2 = tuple((i for i, l in enumerate(bdry_list) if len(486l) == 2 and l[0][0] == 0 and l[0][1] == (n, n+1) and l[1][0] == g))487assert len(indfind2) == 1, (g, k, mu)488bdry_list.pop(indfind2[0])489490preresult = DR_cycle(g, tuple(m + k for m in sortmuprime))491preresult -= Strataboundarysum(g, k, sortmuprime, bdry_list)492result = preresult.forgetful_pushforward([n+1])493result *= (1/sortmuprime[n-1])494495v = result.toTautvect(g, n, g-1)496StrataDB.set_cache(v, *(g, k, sortmu, virt))497result = Tautv_to_tautclass(v, g, n, g-1) # gives simplified version of result498result.rename_legs(ordering_dict, rename=True)499return result500501502def Strataboundarysum(g, k, mu, bdry_list, termsout=False):503r""" Returns sum of boundary terms in Conjecture A for entries of bdry_list. These entries have the format from the output of twistenum, but might be a subset of these.504"""505resultterms = []506n = len(mu)507for vdata in bdry_list:508genera = [b[0] for b in vdata]509legs = [list(b[1]) for b in vdata] # just markings at the moment510edges = []511maxleg = n+1512# stores twist at legs and half-edges513twist = {i: mu[i-1] for i in range(1, n+1)}514for v, (_, _, twistvect) in list(enumerate(vdata))[1:]:515for I in twistvect:516legs[0].append(maxleg)517legs[v].append(maxleg+1)518edges.append((maxleg, maxleg+1))519twist[maxleg] = -I-k520twist[maxleg+1] = I-k521maxleg += 2522bdry_graph = StableGraph(genera, legs, edges)523524vertterms = [Strataclass(genera[0], k, [twist[l] for l in legs[0]])]525vertterms += [Strataclass(genera[v], 1, [twist[l] // k for l in legs[v]])526for v in range(1, len(genera))]527528bdry_term = bdry_graph.boundary_pushforward(vertterms)529coeff = prod([I for (_, _, twistvect) in vdata[1:] for I in twistvect])530coeff /= k**(len(genera)-1) * automorphism_number_fixing_twist(bdry_graph, twist)531resultterms += [coeff * bdry_term]532if termsout:533return resultterms534else:535result = sum(resultterms)536result.simplify()537return result538539540@cached_function541def StrataDB(g, k, mu):542raise NotImplementedError('StrataDB is just an internal database '543'for strata classes, use Strataclass instead')544545546def automorphism_number_fixing_twist(gr, I):547r"""548Return number of automorphisms of gr leaving the twist I on gr invariant.549550EXAMPLES::551552sage: from admcycles import StableGraph553sage: from admcycles.stratarecursion import automorphism_number_fixing_twist554sage: gr = StableGraph([0,0],[[1,2,3],[4,5,6]],[(1,2),(3,4),(5,6)])555sage: twist = {i:0 for i in range(7)}556sage: automorphism_number_fixing_twist(gr,twist)5578558sage: twist[1] = 1559sage: automorphism_number_fixing_twist(gr,twist)5602561sage: twist[2] = 1562sage: automorphism_number_fixing_twist(gr,twist)5634564"""565halfedges = gr.halfedges()566G = gr.leg_automorphism_group()567return len([1 for g in G if all(I[g(h)] == I[h] for h in halfedges)])568569570