| Download
Project: admcycles
Views: 724Visibility: Unlisted (only visible to those who know the link)
Image: ubuntu2004# -*- coding: utf-8 -*-1from __future__ import absolute_import, print_function2from six.moves import range34from copy import deepcopy5import itertools6import pickle78from sage.arith.all import factorial, bernoulli, gcd, lcm9from sage.combinat.all import Subsets10from sage.functions.other import floor, ceil, binomial11from sage.groups.perm_gps.permgroup import PermutationGroup12from sage.matrix.constructor import matrix13from sage.matrix.special import block_matrix14from sage.modules.free_module_element import vector15from sage.combinat.all import Permutations, Subsets, Partitions, Partition16from sage.functions.other import floor, ceil, binomial17from sage.arith.all import factorial, bernoulli, gcd, lcm, multinomial18from sage.misc.cachefunc import cached_function19from sage.misc.misc_c import prod20from sage.modules.free_module_element import vector21from sage.rings.all import Integer, Rational, QQ, ZZ2223from . import DR24#from . import double_ramification_cycle2526initialized = True27g = Integer(0)28n = Integer(3)2930def reset_g_n(gloc, nloc):31gloc = Integer(gloc)32nloc = Integer(nloc)33if gloc < 0 or nloc < 0:34raise ValueError35global g, n36g = gloc37n = nloc3839Hdatabase = {}4041# stgraph is a class modelling stable graphs G42# all vertices and all legs are uniquely globally numbered by integers, vertices numbered by 1,2...,m43# G.genera = list of genera of m vertices44# G.legs = list of length m, where ith entry is set of legs attached to vertex i45# CONVENTION: usually legs are labelled 1,2,3,4,....46# G.edges = edges of G, given as list of ordered 2-tuples, sorted in ascending order by smallest element47# TODO: do I really need sorted in ascending order ...48#49# G.g() = total genus of G5051# Examples52# --------53# Creating a stgraph with two vertices of genera 3,5 joined54# by an edge with a self-loop at the genus 3 vertex.55# sage: stgraph([3,5],[[1,3,5],[2]],[(1,2),(3,5)])56# [3, 5] [[1, 3, 5], [2]] [(1, 2), (3, 5)]5758from .stable_graph import StableGraph, GraphIsom59stgraph = StableGraph606162def trivgraph(g, n):63"""64Return the trivial graph in genus `g` with `n` markings.6566EXAMPLES::6768sage: from admcycles.admcycles import trivgraph69sage: trivgraph(1,1)70[1] [[1]] []71"""72return stgraph([g], [list(range(1, n + 1))], [])737475# gives a list of dual graphs for \bar M_{g,n} of codimension r, i.e. with r edges76#TODO: Huge potential for optimization, at the moment consecutively degenerate and check for isomorphisms77def list_strata(g,n,r):78if r==0:79return [stgraph([g],[list(range(1, n+1))],[])]80Lold=list_strata(g,n,r-1)81Lnew=[]82for G in Lold:83for v in range(G.num_verts()):84Lnew+=G.degenerations(v)85# now Lnew contains a list of all possible degenerations of graphs G with r-1 edges86# need to identify duplicates8788if not Lnew:89return []9091Ldupfree=[Lnew[0]]92count=193# duplicate-free list of strata9495for i in range(1, len(Lnew)):96newgraph=True97for j in range(count):98if Lnew[i].is_isomorphic(Ldupfree[j]):99newgraph=False100break101if newgraph:102Ldupfree+=[Lnew[i]]103count+=1104105return Ldupfree106107# degeneration_graph computes a data structure containing all isomorphism classes of stable graphs of genus g with n markings 1,...,n and the information about one-edge degenerations between them as well as a lookup table for graph-invariants, to quickly identify the index of a given graph108# optional argument: rmax = maximal number of edges that should be computed (so rmax < 3*g-3+n means only a partial degeneration graph is computed - if function had been called with higher rmax before, a larger degeneration_graph will be returned109# the function will look for previous calls in the dictionary of cached_function and restart from there110# output = (deglist,invlist)111# deglist is a list [L0,L1,L2,...] where Lk is itself a list [G0,G1,G2,...] with all isomorphism classes of stable graphs having k edges112# the entries Gj have the form [stgraph Gr, (tuple of the half-edges of Gr), set of morphisms to L(j-1), set of morphisms from L(j+1)]113# a morphism A->B is of the form (DT,halfedgeimages), where114# * DT is the number of B in L(j-1) or of A in L(j+1)115# * halfedgeimages is a tuple (of length #halfedges(B)) with entries 0,1,...,#halfedges(A)-1, where halfedgeimages[i]=j means halfedge number i of B corresponds to halfedge number j in A116# invlist is a list [I0,I1,I2,...] where Ik is a list [invariants,boundaries], where117# * invariants is a list of all Gr.invariant() for Gr in Lk118# * boundaries is a list [i0,i1,i2,...] of indices such that the Graphs in Lk with invariant number N (in invariants) are G(i(N)+1), ..., G(i(N+1))119120@cached_function121def degeneration_graph(g, n, rmax=None):122if rmax is None:123rmax = 3 *g-3 +n124if rmax>3 *g-3 +n:125rmax = 3 *g-3 +n126127# initialize deglist and invlist128try:129deglist, invlist = degeneration_graph.cached(g,n)130r0 = 3 *g-3 +n131except KeyError:132for r0 in range(3*g - 3 + n, -2, -1):133try:134deglist, invlist = degeneration_graph.cached(g,n,r0)135break136except KeyError:137pass138139if r0 == -1: # function has not been called before140deglist=[[[stgraph([g], [list(range(1, n+1))],[]), (), set(), set()]]]141invlist=[[[deglist[0][0][0].invariant()], [-1, 0]]]142143for r in range(r0+1, rmax+1):144deglist+=[[]] # add empty list L(r)145templist=[]146147# first we must degenerate all stgraphs in degree r-1 and record how the degeneration went148for i in range(len(deglist[r-1])):149Gr=deglist[r-1][i][0]150deg_Gr=Gr.degenerations()151for dGr in deg_Gr:152# since all degenerations deg_Gr have the same list of edges as Gr with the new edge appended and since stgraph.halfedges() returns the half-edges in this order, half-edge alpha of Gr corresponds to half-edge alpha in deg_Gr153lis=tuple(range(len(deglist[r-1][i][1])))154templist+=[[dGr, dGr.halfedges(), set([(i,lis)]), set([])]]155156# now templist contains a temporary list of all possible degenerations of deglist[r-1]157# we now proceed to consolidate this list to deglist[r], combining isomorphic entries158def inva(o):159return o[0].invariant()160templist.sort(key=inva)161162# first we reconstruct the blocks of entries having the same graph invariants163current_invariant=None164tempbdry=[]165tempinv=[]166count=0167for Gr in templist:168inv=Gr[0].invariant()169if inv != current_invariant:170tempinv+=[inv]171tempbdry+=[count-1]172current_invariant=inv173count+=1174tempbdry+=[len(templist)-1]175176# now go through blocks and check for actual isomorphisms177count=0 # counts the number of isomorphism classes already added to deglist[r]178179invlist+=[[[],[]]]180for i in range(len(tempinv)):181# go through templist[tempbdry[i]+1], ..., templist[tempbdry[i+1]]182isomcl=[]183for j in range(tempbdry[i]+1, tempbdry[i+1]+1):184Gr=templist[j][0] # must check if stgraph Gr is isomorphic to any of the graphs in isomcl185new_graph=True186for k in range(len(isomcl)):187Gr2=isomcl[k]188ans, Iso = Gr.is_isomorphic(Gr2[0], certificate=True)189if ans:190new_graph=False191dicl = Iso[1]192# modify isomcl[k] to add new graph upstairs, and record this isomorphism in the data of this upstairs graph193upstairs=list(templist[j][2])[0][0] # number of graph upstairs194# h runs through the half-edges of the graph upstairs, which are included via the same name in Gr, so dicl translates them to legs of Gr2195# and we record the index of this image leg in the list of half-edges of Gr2196lisnew=tuple((Gr2[1].index(dicl[h]) for h in deglist[r-1][upstairs][1]))197Gr2[2].add((upstairs, lisnew))198deglist[r-1][upstairs][3].add((count+k,lisnew))199break200if new_graph:201isomcl+=[templist[j]]202# now also add the info upstairs203upstairs=list(templist[j][2])[0][0]204lisnew=list(templist[j][2])[0][1]205deglist[r-1][upstairs][3].add((count+len(isomcl)-1, lisnew))206207deglist[r]+=isomcl208invlist[r][0]+=[isomcl[0][0].invariant()]209invlist[r][1]+=[count-1]210211count+=len(isomcl)212213invlist[r][1]+=[count-1]214215# remove old entry from cache_function dict if more has been computed now216if rmax > r0 and r0 != -1:217degeneration_graph.cache.pop(((g,n,r0),()))218219return (deglist,invlist)220221# finds Gr in the degeneration_graph of the corresponding (g,n), where markdic is a dictionary sending markings of Gr to 1,2,..,n222# returns a tuple (r,i,dicv,dicl), where223# * the graph isomorphic to Gr is in degeneration_graph(g,n)[0][r][i][0]224# * dicv, dicl are dictionaries giving the morphism of stgraphs from Gr to degeneration_graph(g,n)[0][r][i][0]225def deggrfind(Gr,markdic=None):226#global Graph1227#global Graph2228#global Grdeggrfind229#Grdeggrfind = deepcopy(Gr)230g = Gr.g()231n = len(Gr.list_markings())232r = Gr.num_edges()233234if markdic is not None:235Gr_rename = Gr.copy(mutable=True)236Gr_rename.rename_legs(markdic)237else:238markdic={i:i for i in range(1, n+1)}239Gr_rename=Gr240Gr_rename.set_immutable()241Inv=Gr_rename.invariant()242(deglist,invlist) = degeneration_graph(g,n,r)243244InvIndex=invlist[r][0].index(Inv)245for i in range(invlist[r][1][InvIndex]+1, invlist[r][1][InvIndex+1]+1):246#(Graph1,Graph2)= (Gr_rename,deglist[r][i][0])247ans, Isom = Gr_rename.is_isomorphic(deglist[r][i][0], certificate=True)248if ans:249dicl={l:Isom[1][markdic[l]] for l in markdic}250dicl.update({l:Isom[1][l] for l in Gr.halfedges()})251return (r,i,Isom[0],dicl)252print('Help, cannot find ' + repr(Gr))253print((g,n,r))254print((Inv,list_strata(2, 2, 0)[0].invariant()))255256257# returns the list of all A-structures on Gamma, for two stgraphs A,Gamma258# the output is a list [(dicv1,dicl1),(dicv2,dicl2),...], where259# * dicv are (surjective) dictionaries from the vertices of Gamma to the vertices of A260# * dicl are (injective) dictionaries from the legs of A to the legs of Gamma261# optional input: identGamma, identA of the form (rA,iA,dicvA,diclA) are result of deggrfind(A,markdic)262# TODO: for now assume that markings are 1,2,...,n263def Astructures(Gamma,A,identGamma=None,identA=None):264if A.num_edges() > Gamma.num_edges():265return []266g = Gamma.g()267mark = Gamma.list_markings()268n=len(mark)269270if g != A.g():271print('Error, Astructuring graphs of different genera')272print(Gamma)273print(A)274raise ValueError('A very specific bad thing happened')275return []276if set(mark) != set(A.list_markings()):277print('Error, Astructuring graphs with different markings')278return []279markdic={mark[i-1]:i for i in range(1, n+1)}280281282# first identify A and Gamma inside the degeneration graph283(deglist,invlist) = degeneration_graph(g, n, Gamma.num_edges())284285if identA is None:286(rA,iA,dicvA,diclA)=deggrfind(A,markdic)287else:288(rA,iA,dicvA,diclA)=identA289if identGamma is None:290(rG,iG,dicvG,diclG)=deggrfind(Gamma,markdic)291else:292(rG,iG,dicvG,diclG)=identGamma293294# go through the graph in deglist, starting from Gamma going upwards the necessary amount of steps and collect all morphisms in a set295# a morphism is encoded by a tuple (j,(i0,i1,...,im)), where j is the number of the current target of the morphism (in the list deglist[r])296# and i0,i1,..,im give the numbers of the half-edges in Gamma, to which the m+1 half-edges of the target correspond297298morphisms=set([(iG,tuple(range(len(deglist[rG][iG][1]))))]) # start with the identity on Gamma299300301for r in range(rG,rA,-1):302new_morphisms=set()303for mor in morphisms:304# add to new_morphisms all compositions of mor with a morphism starting at the target of mor305for (tar,psi) in deglist[r][mor[0]][2]:306composition=tuple((mor[1][j] for j in psi))307new_morphisms.add((tar,composition))308morphisms=new_morphisms309310# now pick out the morphisms actually landing in position (rA,iA)311new_morphisms = set([mor for mor in morphisms if mor[0]==iA])312morphisms=new_morphisms313314# at the end we have to precompose with the automorphisms of Gamma315Auto = GraphIsom(deglist[rG][iG][0],deglist[rG][iG][0])316new_morphisms=set()317for (Autov,Autol) in Auto:318for mor in morphisms:319# now composition is no longer of the format above320# instead composition[i] is the label of the half-edge in deglist[rG][iG][0] to which half-edge number i in the target corresponds321composition=tuple((Autol[deglist[rG][iG][1][j]] for j in mor[1]))322new_morphisms.add((mor[0],composition))323morphisms=new_morphisms324325326# now we must translate everything back to the original graphs Gamma,A and their labels,327# also reconstructing the vertex-action from the half-edge action328finalmorphisms=[]329dicmarkings={h:h for h in mark}330Ahalfedges=A.halfedges()331Ghalfedges=Gamma.halfedges()332dicAtild={deglist[rA][iA][1][i]:i for i in range(len(deglist[rA][iA][1]))}333diclGinverse={diclG[h]:h for h in Gamma.halfedges()}334335336for (targ,mor) in morphisms:337dicl={ h: diclGinverse[mor[dicAtild[diclA[h]]]] for h in Ahalfedges} # final dictionary for the half-edges338dicl.update(dicmarkings)339340#now reconstruct dicv from this341if A.num_verts() == 1:342dicv = {ver: 0 for ver in range(Gamma.num_verts())}343else:344# Now we pay the price for not recording the morphisms on the vertices before345# we need to reconstruct it by the known morphisms of half-edges by determining the connected components of the complement of beta(H_A)346# TODO: do we want to avoid this reconstruction by more bookkeeping?347348dicv = {Gamma.vertex(dicl[h]): A.vertex(h) for h in dicl}349remaining_vertices = set(range(Gamma.num_verts())) - set(dicv)350remaining_edges = set([e for e in Gamma.edges(copy=False) if e[0] not in dicl.values()])351current_vertices = set(dicv)352353while remaining_vertices:354newcurrent_vertices=set()355for e in deepcopy(remaining_edges):356if Gamma.vertex(e[0]) in current_vertices:357vnew=Gamma.vertex(e[1])358remaining_edges.remove(e)359# if vnew is in current vertices, we don't have to do anything (already know it)360# otherwise it must be a new vertex (cannot reach old vertex past the front of current vertices)361if vnew not in current_vertices:362dicv[vnew]=dicv[Gamma.vertex(e[0])]363remaining_vertices.discard(vnew)364newcurrent_vertices.add(vnew)365continue366if Gamma.vertex(e[1]) in current_vertices:367vnew=Gamma.vertex(e[0])368remaining_edges.remove(e)369# if vnew is in current vertices, we don't have to do anything (already know it)370# otherwise it must be a new vertex (cannot reach old vertex past the front of current vertices)371if vnew not in current_vertices:372dicv[vnew]=dicv[Gamma.vertex(e[1])]373remaining_vertices.discard(vnew)374newcurrent_vertices.add(vnew)375current_vertices=newcurrent_vertices376# now dicv and dicl are done, so they are added to the list of known A-structures377finalmorphisms+=[(dicv,dicl)]378return finalmorphisms379380381# find a list commdeg of generic (G1,G2)-stgraphs G3382# Elements of this list will have the form (G3,vdict1,ldict1,vdict2,ldict2), where383# * G3 is a stgraph384# * vdict1, vdict2 are the vertex-maps from vertices of G3 to G1 and G2385# * ldict1, ldict2 are the leg-maps from legs of G1 and G2 to G3386# modiso = True means we give a list of generic (G1,G2)-stgraphs G3 up to isomorphisms of G3387# if rename=True, the operation will also work if the markings of G1,G2 are not labelled 1,2,...,n; since this involves more bookkeeping, rename=False should run slightly faster388def common_degenerations(G1,G2,modiso=False,rename=False):389#TODO: fancy graph-theoretic implementation in case that G1 has a single edge?390391# almost brute force implementation below392g=G1.g()393mkings=G1.list_markings()394n=len(mkings)395396397# for convenience, we want r1 <= r2398switched = False399if G1.num_edges() > G2.num_edges():400temp = G1401G1 = G2402G2 = temp403switched = True404405if rename:406markdic={mkings[i]:i+1 for i in range(n)}407renamedicl1={l:l+n+1 for l in G1.leglist()}408renamedicl2={l:l+n+1 for l in G2.leglist()}409renamedicl1.update(markdic)410renamedicl2.update(markdic)411412G1 = G1.copy()413G2 = G2.copy()414G1.rename_legs(renamedicl1)415G2.rename_legs(renamedicl2)416G1.set_immutable()417G2.set_immutable()418419(r1,i1,dicv1,dicl1)=deggrfind(G1)420(r2,i2,dicv2,dicl2)=deggrfind(G2)421422(deglist,invlist)=degeneration_graph(g,n,r1+r2)423424commdeg=[]425426descendants1=set([i1])427for r in range(r1+1, r2+1):428descendants1=set([d[0] for i in descendants1 for d in deglist[r-1][i][3]])429430descendants2=set([i2])431432433for r in range(r2, r1+r2+1):434# look for intersection of descendants1, descendants2 and find generic (G1,G2)-struct435for i in descendants1.intersection(descendants2):436Gamma=deglist[r][i][0]437Gammafind=(r,i,{j:j for j in range(Gamma.num_verts())}, {l:l for l in Gamma.leglist()})438G1structures=Astructures(Gamma,G1,identGamma=Gammafind,identA=(r1,i1,dicv1,dicl1))439G2structures=Astructures(Gamma,G2,identGamma=Gammafind,identA=(r2,i2,dicv2,dicl2))440441if modiso is True:442AutGamma=GraphIsom(Gamma,Gamma)443AutGammatild=[(dicv,{dicl[l]:l for l in dicl}) for (dicv,dicl) in AutGamma]444tempdeg=[]445446# now need to identify the generic G1,G2-structures447numlegs=len(Gamma.leglist())448for (dv1,dl1) in G1structures:449for (dv2,dl2) in G2structures:450numcoveredlegs=len(set(list(dl1.values())+list(dl2.values())))451if numcoveredlegs==numlegs:452if switched:453if rename:454tempdeg.append((Gamma,dv2,{l:dl2[renamedicl2[l]] for l in renamedicl2},dv1,{l:dl1[renamedicl1[l]] for l in renamedicl1}))455else:456tempdeg.append((Gamma,dv2,dl2,dv1,dl1))457else:458if rename:459tempdeg.append((Gamma,dv1,{l:dl1[renamedicl1[l]] for l in renamedicl1},dv2,{l:dl2[renamedicl2[l]] for l in renamedicl2}))460else:461tempdeg.append((Gamma,dv1,dl1,dv2,dl2))462463464if modiso is True:465# eliminate superfluous isomorphic elements from tempdeg466while tempdeg:467(Gamma,dv1,dl1,dv2,dl2)=tempdeg.pop(0)468commdeg.append((Gamma,dv1,dl1,dv2,dl2))469for (dicv,dicltild) in AutGammatild:470try:471tempdeg.remove((Gamma,{v: dv1[dicv[v]] for v in dicv},{l: dicltild[dl1[l]] for l in dl1},{v: dv2[dicv[v]] for v in dicv},{l: dicltild[dl2[l]] for l in dl2}))472except:473pass474else:475commdeg+=tempdeg476477# (except in last step) compute new descendants478if r<r1+r2:479descendants1=set([d[0] for i in descendants1 for d in deglist[r][i][3]])480descendants2=set([d[0] for i in descendants2 for d in deglist[r][i][3]])481482return commdeg483484485# Gstgraph is a class modelling stable graphs Gamma with action of a group G and character data at legs486#487# Gamma.G = finite group acting on Gamma488# Gamma.gamma = underlying stable graph gamma of type stgraph489# Gamma.vertact and G.legact = dictionary, assigning to tuples (g,x) of g in G and x a vertex/leg a new vertex/leg490# Gamma.character = dictionary, assigning to leg l a tuple (h,e,k) of491# * a generator h in G of the stabilizer of l492# * the order e of the stabilizer493# * an integer 0<=k<e such that h acts on the tangent space of the curve at the leg l by exp(2 pi i/e * k)494495class Gstgraph(object):496def __init__(self,G,gamma,vertact,legact,character,hdata=None):497self.G=G498self.gamma=gamma499self.vertact=vertact500self.legact=legact501self.character=character502if hdata is None:503self.hdata=self.hurwitz_data()504else:505self.hdata=hdata506507def copy(self, mutable=True):508G = Gstgraph.__new__(Gstgraph)509G.G = self.G510G.gamma = self.gamma.copy(mutable=mutable)511G.vertact = self.vertact.copy()512G.legact = self.legact.copy()513G.character = self.character.copy()514return G515516def __repr__(self):517return repr(self.gamma)518519def __deepcopy__(self,memo):520raise RuntimeError("do not deepcopy")521return Gstgraph(self.G,deepcopy(self.gamma),deepcopy(self.vertact),deepcopy(self.legact),deepcopy(self.character),hdata=deepcopy(self.hdata))522523# returns dimension of moduli space at vertex v; if v is None, return dimension of entire stratum for self524def dim(self,v=None):525if v is None:526return self.quotient_graph().dim()527#TODO: optimization potential528return self.extract_vertex(v).dim()529530# renames legs according to dictionary di531# TODO: no check that this renaming is legal, i.e. produces well-defined graph532def rename_legs(self,di):533if not self.gamma.is_mutable():534self.gamma = self.gamma.copy()535self.gamma.rename_legs(di)536temp_legact={}537temp_character={}538for k in self.legact:539# the operation temp_legact[(k[0],k[1])]=self.legact[k] would produce copy of self.legact, di.get replaces legs by new name, if applicable, and leaves leg invariant otherwise540temp_legact[(k[0],di.get(k[1],k[1]))]=di.get(self.legact[k],self.legact[k])541for k in self.character:542temp_character[di.get(k,k)]=self.character[k]543544545# returns the stabilizer group of vertex i546def vstabilizer(self,i):547gen=[]548for g in self.G:549if self.vertact[(g,i)]==i:550gen+=[g]551try:552return self.G.ambient_group().subgroup(gen)553except AttributeError:554return self.G.subgroup(gen)555556# returns the stabilizer group of leg j557def lstabilizer(self,j):558try:559return self.G.ambient_group().subgroup([self.character[j][0]])560except AttributeError:561return self.G.subgroup([self.character[j][0]])562563564# converts self into a prodHclass (with gamma0 being the corresponding trivial graph)565def to_prodHclass(self):566return prodHclass(stgraph([self.gamma.g()],[list(self.gamma.list_markings())],[]),[self.to_decHstratum()])567568# converts self into a decHstratum (with gamma0 being the corresponding trivial graph)569def to_decHstratum(self):570dicv={}571dicvinv={}572dicl={}573quotgr=self.quotient_graph(dicv, dicvinv, dicl)574575spaces = {v: (self.gamma.genera(dicvinv[v]), HurwitzData(self.vstabilizer(dicvinv[v]),[self.character[l] for l in quotgr.legs(v, copy=False)])) for v in range(quotgr.num_verts())}576577masterdicl={}578579580581for v in range(quotgr.num_verts()):582n=1583gord=spaces[v][1].G.order()584tGraph=trivGgraph(*spaces[v])585vinv=dicvinv[v]586587# for l in quotgr.legs(v) find one leg l' in the orbit of l which belongs to vinv588legreps=[]589for l in quotgr.legs(v, copy=False):590for g in self.G:591if self.legact[(g,l)] in self.gamma.legs(vinv, copy=False):592legreps.append(self.legact[(g,l)])593break594595# make a list quotrep of representatives of self.G / spaces[v][1].G (quotient by Stab(vinv))596(quotrep,di)=leftcosetaction(self.G,spaces[v][1].G)597598for l in legreps:599for g in spaces[v][1].G:600for gprime in quotrep:601masterdicl[ self.legact[(gprime*g,l)]] = tGraph.legact[(g,n)]602n+=floor(QQ(gord)/self.character[l][1])603604605vertdata=[]606for w in range(self.gamma.num_verts()):607a=dicv[w]608609wdicl = {l:masterdicl[l] for l in self.gamma.legs(w, copy=False)}610611vertdata.append([a,wdicl])612613return decHstratum(deepcopy(self.gamma),spaces,vertdata)614615616# takes vertex v_i of the stable graph and returns a one-vertex Gstgraph with group G_{v_i} and all the legs and self-loops attached to v_i617def extract_vertex(self,i):618G_new=self.vstabilizer(i)619lgs = self.gamma.legs(i, copy=True) # takes list of legs attached to i620egs=self.gamma.edges_between(i,i)621gamma_new=stgraph([self.gamma.genera(i)],[lgs],egs)622623vertact_new={}624for g in G_new:625vertact_new[(g,0)]=0626627legact_new={}628for g in G_new:629for j in lgs:630legact_new[(g,j)]=self.legact[(g,j)]631632character_new={}633for j in lgs:634character_new[j]=deepcopy(self.character[j])635636return Gstgraph(G_new,gamma_new,vertact_new,legact_new,character_new)637638# glues the Gstgraph Gr (with group H) at the vertex i of a Gstgraph self639# and performs this operation equivariantly under the G-action on self640# all legs at the elements of the orbit of i are not renamed and the G-action and character on them is not changed641# optional arguments: if divGr/dil are given they are supposed to be a dictionary, which will be cleared and updated with the renaming-convention to pass from leg/vertex-names in Gr to leg/vertex-names in the glued graph (at former vertex i)642# similarly, divs will be a dictionary assigning vertex numbers in the old self (not in the orbit of i) the corresponding number in the new self643# necessary conditions:644# * the stabilizer of vertex i in self = H645# * every leg of i is also a leg in Gr646# note that legs of Gr which belong to an edge in Gr are not identified with legs of self, even if they have the same name647def equivariant_glue_vertex(self,i,Gr,divGr={},divs={},dil={}):648divGr.clear()649divs.clear()650dil.clear()651652for j in range(self.gamma.num_verts()):653divs[j]=j654655G=self.G656H=Gr.G657658# the orbit of vertex i is given by gH*i for g in L659(L,di)=leftcosetaction(G,H)660oldvertno = self.gamma.num_verts() - len(L) # number of old vertices that will remain in end661662# first we rename every leg in Gr that is not a leg of i such that the new leg-names do not appear in self663# TODO: to be changed (access to private _maxleg attribute)664Gr_mod = Gr.copy()665m = max(self.gamma._maxleg, Gr.gamma._maxleg)666667a = set().union(*Gr.gamma.legs(copy=False)) # legs of Gr668b = set(self.gamma.legs(i, copy=False)) # legs of self at vertex i669e = a-b # contains the set of legs of Gr that are not legs of vertex i670augment_dic={}671for l in e:672m+=1673augment_dic[l]=m674675Gr_mod.rename_legs(augment_dic)676677# collects the dictionaries translating legs in Gr_mod not belonging to vertex i to legs in the glued graph at elements g in L678legdiclist={}679680# records orbit of i under elements of L681orbi={}682oldlegsi = self.gamma.legs(i, copy=False)683684# go through the orbit of i685for g in L:686orbi[g]=self.vertact[(g,i)]687# create a translated copy of Gr using the information of the G-action on self688Gr_transl = Gr_mod.copy()689transl_dic={l:self.legact[(g,l)] for l in oldlegsi}690Gr_transl.rename_legs(transl_dic)691692# glue G_transl to vertex g*i of self (use divs to account for deletions in meantime) and remember how internal edges are relabelled693divGr_temp={}694divs_temp={}695dil_temp={}696self.gamma.glue_vertex(divs[self.vertact[(g,i)]], Gr_transl.gamma, divGr=divGr_temp, divs=divs_temp, dil=dil_temp)697698# update various dictionaries now699divs.pop(self.vertact[(g,i)]) # no longer needed700for j in divs:701divs[j]=divs_temp[divs[j]]702legdiclist.update({(g,j):dil_temp[j] for j in dil_temp})703704# adjust vertact705# ERROR here, from renaming vetices -> rather use vertact_reconstruct()706# for g in L:707# for h in G:708# self.vertact.pop((h,orbi[g])) # remove action of h on g*i709#710#711# for g1 in range(len(L)):712# for g2 in range(len(L)):713# for h in H:714# for j in range(len(Gr.gamma.genera)):715# # record where vertex number j in Gr(translated by g1) goes under g2*h716# self.vertact[(L[g2]*h,oldvertno+g1*len(Gr.gamma.genera)+j)]=oldvertno+di[(L[g2]*h,g1)]*len(Gr.gamma.genera)+j717718# adjust legact and character719for j in e:720for g1 in range(len(L)):721# we are adjusting here the data for the leg corresponding to j in Gr that was glued to vertex g1*i722for g2 in range(len(L)):723for h in H:724# the element g2*h*(g1)**{-1} acts on g1*i resulting in g2*h*i725# for g1 fixed, the expression g2*h*(g1)^{-1} runs through the elements of G726self.legact[(L[g2]*h*L[g1].inverse(),legdiclist[(L[g1],augment_dic[j])])]=legdiclist[(L[g2],augment_dic[Gr.legact[(h,j)]])]727# character given by conjugation728self.character[legdiclist[(L[g1],augment_dic[j])]]=(L[g1]*Gr.character[j][0]*L[g1].inverse(),Gr.character[j][1],Gr.character[j][2])729730self.vertact_reconstruct()731732dil.update({j:legdiclist[(L[0],augment_dic[j])] for j in e}) # here we assume that L[0]=id_G733divGr.update({j:oldvertno+j for j in range(Gr.gamma.num_verts())})734735def quotient_graph(self, dicv={}, dicvinv={}, dicl={}):736r"""737Computes the quotient graph of self under the group action.738739dicv, dicl are dictionaries that record the quotient map of vertices and legs,740dicvinv is a dictionary giving some section of dicv (going from vertices of741quotient to vertices of self).742743EXAMPLES::744745sage: from admcycles.admcycles import trivGgraph, HurData746sage: G = PermutationGroup([(1, 2)])747sage: H = HurData(G, [G[1], G[1], G[1], G[1]])748sage: tG = trivGgraph(1, H); tG.quotient_graph()749[0] [[1, 2, 3, 4]] []750sage: H = HurData(G, [G[1], G[1], G[1], G[1], G[0]])751sage: tG = trivGgraph(1, H); tG.quotient_graph()752[0] [[1, 2, 3, 4, 5]] []753sage: H = HurData(G, [G[1], G[1], G[1], G[1], G[1]])754sage: tG = trivGgraph(1, H); tG.quotient_graph() # no such cover exists (RH formula)755Traceback (most recent call last):756...757ValueError: Riemann-Hurwitz formula not satisfied in this Gstgraph.758759TESTS::760761sage: from admcycles.admcycles import trivGgraph, HurData762sage: G = PermutationGroup([(1, 2, 3)])763sage: H = HurData(G, [G[1]])764sage: tG = trivGgraph(2, H); tG.quotient_graph()765[1] [[1]] []766sage: H = HurData(G, [G[1] for i in range(6)])767sage: tG = trivGgraph(1, H); tG.quotient_graph() # quotient vertex would have genus -1768Traceback (most recent call last):769...770ValueError: Riemann-Hurwitz formula not satisfied in this Gstgraph.771"""772dicv.clear()773dicvinv.clear()774dicl.clear()775776vertlist = list(range(self.gamma.num_verts()))777leglist=[]778for v in vertlist:779leglist += self.gamma.legs(v, copy=False)780781countv=0782# compute orbits of G-action783for v in deepcopy(vertlist):784if v in vertlist: # has not yet been part of orbit of previous point785dicvinv[countv]=v786for g in self.G:787if self.vertact[(g,v)] in vertlist:788dicv[self.vertact[(g,v)]]=countv789vertlist.remove(self.vertact[(g,v)])790countv+=1791792for l in deepcopy(leglist):793if l in leglist: # has not yet been part of orbit of previous point794for g in self.G:795if self.legact[(g,l)] in leglist:796dicl[self.legact[(g,l)]]=l # note that leg-names in quotient graph equal the name of one preimage in self797leglist.remove(self.legact[(g,l)])798799# create new stgraph with vertices/legs given by values of dicv, dicl800quot_leg=[[] for j in range(countv)] # list of legs of quotient801legset=set(dicl.values())802for l in legset:803quot_leg[dicv[self.gamma.vertex(l)]]+=[l]804805quot_genera=[]806for v in range(countv):807Gv=self.vstabilizer(dicvinv[v]).order()808b=sum([1-QQ(1)/(self.character[l][1]) for l in quot_leg[v]]) # self.character[l][1] = e = order of Stab_l809# genus of quotient vertex by Riemann-Hurwitz formula810qgenus = ((2 *self.gamma.genera(dicvinv[v])-2)/QQ(Gv)+2 -b)/QQ(2)811if not (qgenus.is_integer() and qgenus>=0):812raise ValueError('Riemann-Hurwitz formula not satisfied in this Gstgraph.')813quot_genera+=[ZZ(qgenus)]814815quot_edges=[]816for e in self.gamma.edges(copy=False):817if e[0] in legset:818quot_edges+=[(e[0],dicl[e[1]])]819820quot_graph = stgraph(quot_genera,quot_leg, quot_edges, mutable=True)821quot_graph.tidy_up()822quot_graph.set_immutable()823824return quot_graph825826# returns Hurwitz data corresponding to the Hurwitz space in which self is a boundary stratum827# HOWEVER: since names of legs are not necessarily canonical, there is no guarantee that they will be compatible828# i.e. self is not necessarily a degeneration of trivGgraph(g,hurwitz_data(self))829def hurwitz_data(self):830quotgr=self.quotient_graph()831marks=list(quotgr.list_markings())832marks.sort()833return HurwitzData(self.G,[self.character[l] for l in marks])834835# TODO: currently ONLY works for cyclic groups G836def delta_degree(self,v):837r"""838Gives the degree of the delta-map from the Hurwitz space parametrizing the curve839placed on vertex v to the corresponding stable-maps space.840841Note: currently only works for cyclic groups G.842843EXAMPLES::844845sage: from admcycles.admcycles import trivGgraph, HurData846sage: G = PermutationGroup([(1, 2)])847sage: H = HurData(G, [G[1], G[1], G[1], G[1]])848sage: tG = trivGgraph(1, H); tG.delta_degree(0) # unique cover with 2 automorphisms8491/2850sage: H = HurData(G, [G[1], G[1], G[1], G[1], G[0]])851sage: tG = trivGgraph(1, H); tG.delta_degree(0) # unique cover with 1 automorphism8521853sage: H = HurData(G, [G[1], G[1], G[1], G[1], G[1]])854sage: tG = trivGgraph(1, H); tG.delta_degree(0) # no such cover exists (RH formula)8550856sage: G = PermutationGroup([(1, 2, 3)])857sage: H = HurData(G, [G[1]])858sage: tG = trivGgraph(2, H); tG.delta_degree(0) # no such cover exists (repr. theory)8590860"""861if self.gamma.num_verts() > 1:862return self.extract_vertex(v).delta_degree(0)863try:864quotgr=self.quotient_graph()865except ValueError: # Riemann-Hurwitz formula not satisfied, so space empty and degree 0866return 0867gprime = quotgr.genera(0)868n=self.G.order()869870# verify that Hurwitz space is nonempty #TODO: this assumes G cyclic (like the rest of the function)871if (n==2 and len([1 for l in quotgr.legs(0, copy=False) if self.character[l][1]==2])%2==1) or (n>2 and prod([self.character[l][0]**(ZZ(self.character[l][2]).inverse_mod(n)) for l in quotgr.legs(0, copy=False) if self.character[l][1]!=1]) != self.G.one()):872return 0873874# compute number of homomorphisms of pi_1(punctured base curve) to G giving correct Hurwitz datum by inclusion-exclusion type formula875sigmagcd = ceil(QQ(n)/lcm([self.character[l][1] for l in quotgr.legs(0, copy=False)]))876primfac=set(sigmagcd.prime_factors())877numhom=0878for S in Subsets(primfac):879numhom+=(-1)**len(S)*(QQ(n)/prod(S))**(2 *gprime)880881return numhom*prod([QQ(n)/self.character[l][1] for l in quotgr.legs(0, copy=False)])/QQ(n) #TODO: why not divide by gcd([self.character[l][1] for l in quotgr.legs(0)]) instead???882883# reconstructs action of G on vertices from action on legs. Always possible for connected stable graphs884def vertact_reconstruct(self):885if self.gamma.num_verts() == 1:886# only one vertex means trivial action887self.vertact={(g,0):0 for g in self.G}888else:889self.vertact={(g,v): self.gamma.vertex(self.legact[(g,self.gamma.legs(v, copy=False)[0])]) for g in self.G for v in range(self.gamma.num_verts())}890891892893# returns a list of all possible codimension 1 degenerations of the Gstgraph self occuring at the vertex v (and the elements of the orbit of v under G894# for this, v is extracted and seen as a Gstgraph with its own stabilizer group as symmetries and then split up according to the classification of boundary divisors in Hurwitz stacks. Then all possible degenerations of v are glued back in equivariantly to the original graph895# NOTE: we do not guarantee that (by global symmetries of the graph) all these degenerations are nonisomorphic896# TODO: Guarantee that for codim 0 to codim 1 all boundary components are only listed exactly once?897def degenerations(self,v):898# for more than one vertex, extract v and proceed as described above899if self.gamma.num_verts() > 1:900v0=self.extract_vertex(v)901L=v0.degenerations(0)902M = [self.copy() for j in range(len(L))]903for j in range(len(L)):904M[j].equivariant_glue_vertex(v,L[j])905return M906if self.gamma.num_verts() == 1:907dicv={}908dicvinv={}909dicl={}910quot_Graph=self.quotient_graph(dicv=dicv, dicvinv=dicvinv, dicl=dicl)911912quot_degen=quot_Graph.degenerations(0)913914degen=[] #list of degenerations915916if self.G.is_cyclic():917# first create dictionaries GtoNum, NumtoG translating elements of G to numbers 0,1, ..., G.order()-1918G=self.G919n=G.order()920921922# cumbersome method to obtain some generator. since G.gens() is not guaranteed to be minimal, go through and check order923Ggenerators=G.gens()924for ge in Ggenerators:925if ge.order()==n:926sigma=ge927break928929# now sigma is fixed generator of G, want to have dictionaries such that sigma <-> 1, sigma**2 <-> 2, ..., sigma**n <-> 0930GtoNum={}931NumtoG={}932mu=sigma933934for j in range(1, n+1):935GtoNum[mu]=j % n936NumtoG[j % n]=mu937mu=mu*sigma938939940# extract data of e_alpha, k_alpha, nu_alpha, where alpha runs through legs of dege941e={}942k={}943nu={}944945946for dege in quot_degen:947# all legs in dege come from legs of self, except the one added by the degeneration, which is last in the list dege.edges948last = dege.edges(copy=False)[dege.num_edges() - 1]949legs = set(dege.leglist())950legs.difference_update(last)951for alpha in legs:952e[alpha]=self.character[alpha][1]953r=floor(GtoNum[self.character[alpha][0]]*e[alpha]/QQ(n))954k[alpha]=(self.character[alpha][2]*r.inverse_mod(e[alpha]))% e[alpha]955nu[alpha]=k[alpha].inverse_mod(e[alpha]) # case e[alpha]=1 works, produces nu[alpha]=0956957if dege.num_verts() == 2:958I1 = dege.legs(0, copy=True)959last = dege.edges(copy=False)[dege.num_edges() - 1]960I1.remove(last[0])961I2 = dege.legs(1, copy=True)962I2.remove(last[1])963I = I1+I2964965n1_min=lcm([e[alpha] for alpha in I1])966n2_min=lcm([e[alpha] for alpha in I2])967968969#TODO: following for-loop has great optimization potential, but current form should be sufficient for now970for n1_add in range(1, floor(QQ(n)/n1_min + 1)):971for n2_add in range(1, floor(QQ(n)/n2_min + 1)):972n1=n1_add*n1_min # nj is the order of the stabilizer of the preimage of vertex j-1 in dege973n2=n2_add*n2_min974if n1.divides(n) and n2.divides(n) and lcm(n1,n2)==n: # lcm condition from connectedness of final graph975# we must make sure that the order e of the stabilizer of the preimage of the separating edge and the976# corresponding characters (described by k1,k2) at the legs produce nonempty Hurwitz spaces977# this amounts to a congruence relation mod n1, n2, which can be explicitly solved978A1=-floor(sum([QQ(n1)/e[alpha]*nu[alpha] for alpha in I1]))979e_new=floor(QQ(n1)/gcd(n1,A1))980nu1=floor(QQ(A1%n1)/(QQ(n1)/e_new))981nu2=e_new-nu1982k1=nu1.inverse_mod(e_new)983k2=nu2.inverse_mod(e_new)984985# now we list all ways to distribute the markings, avoiding duplicates by symmetries986I1_temp=deepcopy(I1)987I2_temp=deepcopy(I2)988possi={} #dictionary associating to a leg l in I the possible positions of l (from 0 to n/n1-1 if l in I1, ...)989990if I1:991possi[I1[0]]=[0] # use G-action to move marking I1[0] to zeroth vertex on the left992I1_temp.pop(0)993if I2: # if there is one more leg on the other vertex, use stabilizer of zeroth left vertex to move994possi[I2[0]]=list(range(gcd(n1,QQ(n)/n2)))995I2_temp.pop(0)996else:997if I2:998possi[I2[0]]=[0]999I2_temp.pop(0)1000for i in I1_temp:1001possi[i]=list(range(QQ(n)/n1))1002for i in I2_temp:1003possi[i]=list(range(QQ(n)/n2))100410051006# mark_dist is a list of the form [(v1,v2,...), ...] where leg I[0] goes to vertex number v1 in its orbit, ....1007mark_dist=itertools.product(*[possi[j] for j in I])10081009# now all choices are made, and we construct the degenerated graph, to add it to degen1010# first: use Riemann-Hurwitz to compute the genera of the vertices above 0,11011g1=floor((n1*(2 *dege.genera(0)-2)+n1*(sum([1-QQ(1)/e[alpha] for alpha in I1])+1-QQ(1)/e_new))/QQ(2)+1)1012g2=floor((n2*(2 *dege.genera(1)-2)+n2*(sum([1-QQ(1)/e[alpha] for alpha in I2])+1-QQ(1)/e_new))/QQ(2)+1)10131014# nonemptyness condition #TODO: is this the only thing you need to require?1015if g1<0 or g2<0:1016continue10171018new_genera=[g1 for j in range(QQ(n)/n1)]+[g2 for j in range(QQ(n)/n2)]1019new_legs=[[] for j in range(QQ(n)/n1+QQ(n)/n2)]1020new_legact=deepcopy(self.legact)1021new_edges = self.gamma.edges()1022new_character=deepcopy(self.character)10231024# now go through the orbit of the edge and adjust all necessary data1025# TODO: to be changed (access to private maxleg attribute)1026m0=self.gamma._maxleg+11027m=self.gamma._maxleg+110281029for j in range(QQ(n)/e_new):1030new_legs[j % (QQ(n)/n1)]+=[m]1031new_legs[(j % (QQ(n)/n2))+(QQ(n)/n1)]+=[m+1]1032new_edges+=[(m,m+1)]1033new_legact[(sigma,m)]=m0+((m+2 -m0)%(2 *QQ(n)/e_new))1034new_legact[(sigma,m+1)]=new_legact[(sigma,m)]+11035new_character[m]=(NumtoG[(QQ(n)/e_new)%n], e_new, k1)1036new_character[m+1]=(NumtoG[(QQ(n)/e_new)%n], e_new, k2)10371038m+=2103910401041# take care of remaining actions of g in G-{sigma} on legs belonging to orbit of new edge1042for jplus in range(2, n+1):1043for m in range(m0, m0 + 2*QQ(n)/e_new):1044new_legact[(NumtoG[jplus%n],m)]=new_legact[(sigma,new_legact[(NumtoG[(jplus-1)%n],m)])]10451046for md in mark_dist:1047new_legs_plus=deepcopy(new_legs)1048for alpha in range(len(I1)):1049for j in range(QQ(n)/e[I[alpha]]):1050new_legs_plus[(md[alpha]+j)%(QQ(n)/n1)]+=[self.legact[(NumtoG[j],I[alpha])]]1051for alpha in range(len(I1),len(I)):1052for j in range(QQ(n)/e[I[alpha]]):1053new_legs_plus[(md[alpha]+j)%(QQ(n)/n2)+(QQ(n)/n1)]+=[self.legact[(NumtoG[j],I[alpha])]]10541055# all the data is now ready to generate the new graph1056new_Ggraph=Gstgraph(G,stgraph(new_genera,new_legs_plus,new_edges),{},new_legact,new_character,hdata=self.hdata)1057new_Ggraph.vertact_reconstruct()1058degen+=[new_Ggraph]1059if dege.num_verts() == 1: # case of self-loop1060I = dege.legs(0, copy=True)1061last = dege.edges(copy=False)[dege.num_edges() - 1]1062I.remove(last[0]) # remove the legs corresponding to the degenerated edge1063I.remove(last[1])10641065n0_min=lcm([e[alpha] for alpha in I])10661067#TODO: following for-loop has great optimization potential, but current form should be sufficient for now1068for n0_add in range(1, floor(QQ(n)/n0_min + 1)):1069n0=n0_add*n0_min # n0 is the order of the stabilizer of any of the vertices1070if n0.divides(n):1071for e_new in n0.divisors():1072for g0 in [x for x in range(floor(QQ(n)/(2*n0)+1)) if gcd(x,QQ(n)/n0)==1]:1073for nu1 in [j for j in range(e_new) if gcd(j,e_new)==1]:1074if g0==0 and nu1>QQ(e_new)/2: # g0==0 means only one vertex and then there is a symmetry inverting edges1075continue1076nu2=e_new-nu11077k1=Integer(nu1).inverse_mod(e_new)1078k2=Integer(nu2).inverse_mod(e_new)10791080# now we list all ways to distribute the markings, avoiding duplicates by symmetries1081I_temp=deepcopy(I)10821083possi={} #dictionary associating to a leg l in I the possible positions of l (from 0 to n/n1-1 if l in I1, ...)10841085if I:1086possi[I[0]]=[0] # use G-action to move marking I1[0] to zeroth vertex1087I_temp.pop(0)1088for i in I_temp:1089possi[i]=list(range(QQ(n)/n0))109010911092# mark_dist is a list of the form [(v1,v2,...), ...] where leg I[0] goes to vertex number v1 in its orbit, ....1093mark_dist=itertools.product(*[possi[j] for j in I])109410951096# now all choices are made, and we construct the degenerated graph, to add it to degen1097# first: use Riemann-Hurwitz to compute the genera of the vertices above 0,11098gen0=floor((n0*(2 *dege.genera(0)-2)+n0*(sum([1-QQ(1)/e[alpha] for alpha in I])+2 -QQ(2)/e_new))/QQ(2)+1)10991100# nonemptyness condition #TODO: is this the only thing you need to require?1101if gen0<0:1102continue11031104new_genera=[gen0 for j in range(QQ(n)/n0)]1105new_legs=[[] for j in range(QQ(n)/n0)]1106new_legact=deepcopy(self.legact)1107new_edges = self.gamma.edges()1108new_character=deepcopy(self.character)11091110# now go through the orbit of the edge and adjust all necessary data1111# TODO: to be changed (access to private _maxleg attribute)1112m0=self.gamma._maxleg+11113m=self.gamma._maxleg+111141115for j in range(QQ(n)/e_new):1116new_legs[j % (QQ(n)/n0)]+=[m]1117new_legs[(j+g0) % (QQ(n)/n0)]+=[m+1]1118new_edges+=[(m,m+1)]1119new_legact[(sigma,m)]=m0+((m+2 -m0)%(2 *QQ(n)/e_new))1120new_legact[(sigma,m+1)]=new_legact[(sigma,m)]+11121new_character[m]=(NumtoG[(QQ(n)/e_new)%n],e_new, k1)1122new_character[m+1]=(NumtoG[(QQ(n)/e_new)%n],e_new, k2)11231124m+=211251126# take care of remaining actions of g in G-{sigma} on legs belonging to orbit of new edge1127for jplus in range(2, n+1):1128for m in range(m0, m0+2*QQ(n)/e_new):1129new_legact[(NumtoG[jplus%n],m)]=new_legact[(sigma,new_legact[(NumtoG[(jplus-1)%n],m)])]11301131for md in mark_dist:1132new_legs_plus=deepcopy(new_legs)1133for alpha in range(len(I)):1134for j in range(QQ(n)/e[I[alpha]]):1135new_legs_plus[(md[alpha]+j)%(QQ(n)/n0)]+=[self.legact[(NumtoG[j],I[alpha])]]113611371138# all the data is now ready to generate the new graph1139new_Ggraph=Gstgraph(G,stgraph(new_genera,new_legs_plus,new_edges),{},new_legact,new_character,hdata=self.hdata)1140new_Ggraph.vertact_reconstruct()1141degen+=[new_Ggraph]1142return degen1143else:1144print('Degenerations for non-cyclic groups not implemented!')1145return []11461147# returns the list of all G-isomorphisms of the Gstgraphs Gr1 and Gr21148# isomorphisms must respect markings (=legs that don't appear in edges)1149# they are given as [dictionary of images of vertices, dictionary of images of legs]1150#TODO:IDEA: first take quotients, compute isomorphisms between them, then lift??1151def equiGraphIsom(Gr1,Gr2):1152if Gr1.G != Gr2.G:1153return []1154Iso = GraphIsom(Gr1.gamma,Gr2.gamma) # list of all isomorphisms of underlying dual graphs1155GIso=[]11561157#TODO: once Gequiv=False, can abort other loops1158#TODO: check character data!!1159for I in Iso:1160Gequiv=True1161for g in Gr1.G.gens():1162# check if isomorphism I is equivariant with respect to g: I o g = g o I1163# action on vertices:1164for v in I[0]:1165if I[0][Gr1.vertact[(g,v)]] != Gr2.vertact[(g,I[0][v])]:1166Gequiv=False1167break1168# action on legs:1169for l in I[1]:1170if I[1][Gr1.legact[(g,l)]] != Gr2.legact[(g,I[1][l])]:1171Gequiv=False1172break1173# character:1174for l in I[1]:1175if not Gequiv:1176break1177for j in [t for t in range(Gr1.character[l][1]) if gcd(t,Gr1.character[l][1])==1]:1178if Gr1.character[l][0]**j == Gr2.character[I[1][l]][0] and not ((j*Gr1.character[l][2]-Gr2.character[I[1][l]][2])%Gr1.character[l][1])==0:1179Gequiv=False1180break1181if Gequiv:1182GIso+=[I]11831184return GIso11851186# HurwitzData models the ramification data D of a Galois cover of curves in the following way:1187#1188# D.G = finite group1189# D.l = list of elements of the form (h,e,k) corresponding to orbits of markings with1190# * a generator h in G of the stabilizer of the element p of an orbit1191# * the order e of the stabilizer1192# * an integer 0<=k<e such that h acts on the tangent space of the curve at p by exp(2 pi i/e * k)1193class HurwitzData:1194def __init__(self,G,l):1195self.G=G1196self.l=l1197def __eq__(self,other):1198if not isinstance(other,HurwitzData):1199return False1200return (self.G==other.G) and (self.l==other.l)1201def __hash__(self):1202return hash((self.G,tuple(self.l)))1203def __repr__(self):1204return repr((self.G,self.l))1205def __deepcopy__(self,memo):1206return HurwitzData(self.G,deepcopy(self.l))1207def nummarks(self):1208g=self.G.order()1209N=sum([QQ(g)/r[1] for r in self.l])1210return floor(N)12111212# HurData takes a group G and a list l of group elements (giving monodromy around the corresponding markings) and gives back the corresponding Hurwitz data1213# This ist just a more convenient way to enter HurwitzData1214def HurData(G,l):1215return HurwitzData(G,[(h,h.order(),min(h.order(),1)) for h in l])12161217# returns the trivial stable graph of genus gen with group action corresponding to the HurwitzData D1218def trivGgraph(gen,D):1219G=D.G12201221#graph has only one vertex with trivial G-action1222vertact={}1223for g in D.G:1224vertact[(g,0)]=012251226n=0 #counts number of legs as we go through Hurwitz data1227legact={}1228character={}12291230for e in D.l:1231# h=e[0] generates the stabilizer of (one element of) the orbit corresponding to e1232H=G.subgroup([e[0]])1233# markings corresponding to entry e are labelled by cosets gH and G-action is given by left-action of G1234(L,dic)=leftcosetaction(G,H)1235for g in G:1236for j in range(len(L)):1237legact[(g,j+1+n)] = dic[(g,j)]+1+n #shift comes from fact that the n legs 1,2,...,n have already been assigned12381239for j in range(len(L)):1240# Stabilizer at point gH generated by ghg**{-1}, of same order e, acts still by exp(2 pi i/e * k)1241character[j+1+n]=(L[j]*e[0]*L[j].inverse(), e[1],e[2])1242n+=len(L)124312441245gamma=stgraph([gen],[list(range(1, n+1))],[])12461247return Gstgraph(G,gamma,vertact,legact,character,hdata=D)124812491250# gives a (duplicate-free) list of Gstgraphs for the Hurwitz space of genus g with HurwitzData H in codimension r1251#TODO: Huge potential for optimization, at the moment consecutively degenerate and check for isomorphisms12521253@cached_function1254def list_Hstrata(g,H,r):1255if r==0:1256return [trivGgraph(g,H)]1257Lold=list_Hstrata(g,H,r-1)1258Lnew=[]1259for Gr in Lold:1260for v in range(Gr.gamma.num_verts()):1261Lnew+=Gr.degenerations(v)1262# now Lnew contains a list of all possible degenerations of graphs G with codimension r-11263# need to identify duplicates12641265if not Lnew:1266return []12671268Ldupfree=[Lnew[0]]1269count=11270# duplicate-free list of strata12711272for i in range(1, len(Lnew)):1273newgraph=True1274for j in range(count):1275if equiGraphIsom(Lnew[i],Ldupfree[j]):1276newgraph=False1277break1278if newgraph:1279Ldupfree+=[Lnew[i]]1280count+=112811282return Ldupfree12831284# gives a list of the quotgraphs of the entries of list_Hstrata(g,H,r) and the corresponding dictionaries1285# list elements are of the form [quotient graph,deggrfind(quotient graph),dicv,dicvinv,dicl]1286# if localize=True, finds the quotient graphs in the corresponding degeneration_graph and records their index as well as an isomorphism to the standard-graph1287# then the result has the form: [quotient graph,index in degeneration_graph,dicv,dicvinv,dicl,diclinv], where1288# * dicv is the quotient map sending vertices of Gr in list_Hstrata(g,H,r) to vertices IN THE STANDARD GRAPH inside degeneration_graph1289# * dicvinv gives a section of dicv1290# * dicl sends leg names in Gr to leg names IN THE STANDARD GRAPH1291# * diclinv gives a section of dicl12921293@cached_function1294def list_quotgraphs(g,H,r,localize=True):1295Hgr=list_Hstrata(g,H,r)1296result=[]1297for gr in Hgr:1298dicv={}1299dicvinv={}1300dicl={}1301qgr=gr.quotient_graph(dicv,dicvinv,dicl)1302if localize:1303dgfind=deggrfind(qgr)1304defaultdicv=dgfind[2]1305defaultdicl=dgfind[3]1306defaultdicvinv={defaultdicv[v]:v for v in defaultdicv}1307defaultdiclinv={defaultdicl[l]:l for l in defaultdicl}13081309result.append([qgr,dgfind[1],{v: defaultdicv[dicv[v]] for v in dicv},{w:dicvinv[defaultdicvinv[w]] for w in defaultdicvinv},{l:defaultdicl[dicl[l]] for l in dicl},defaultdiclinv])1310else:1311result.append([qgr,dicv,dicvinv,dicl])1312return result13131314# cyclicGstgraph is a convenient method to create a Gstgraph for a cyclic group action. It takes as input1315# * a stgraph Gr1316# * a number n giving the order of the cyclic group G1317# * a tuple perm of the form ((1,2),(3,),(4,5,6)) describing the action of a generator sigma of G on the legs; all legs need to appear to fix their order for interpreting cha1318# * a tuple cha=(2,3,1) giving for each of the cycles of legs l in perm the number k such that sigma**{n/|G_l|} acts on the tangent space at l by exp(2 pi i/|G_l| *k)1319# and returns the corresponding Gstgraph1320# optionally, the generator sigma can be given as input directly and G is computed as its parent1321# Here we use that the G-action on a connected graph Gr is uniquely determined by the G-action on the legs1322# TODO: user-friendliness: allow (3) instead of (3,)1323def cyclicGstgraph(Gr,n,perm,cha,sigma=None):1324if sigma is None:1325G=PermutationGroup(tuple(range(1, n+1)))1326sigma=G.gens()[0]1327else:1328G=sigma.parent()132913301331# first legaction and character1332perm_dic={} # converts action of sigma into dictionary1333character={}13341335for cycleno in range(len(perm)):1336e=floor(QQ(n)/len(perm[cycleno])) # order of stabilizer1337for j in range(len(perm[cycleno])):1338perm_dic[perm[cycleno][j]] = perm[cycleno][(j+1)%len(perm[cycleno])]1339character[perm[cycleno][j]]=(sigma**(len(perm[cycleno])), e, cha[cycleno])134013411342legact={(sigma,k):perm_dic[k] for k in perm_dic}13431344mu=sigma1345for j in range(2, n+1):1346mu2=mu*sigma1347legact.update({(mu2,k):perm_dic[legact[(mu,k)]] for k in perm_dic})1348mu=mu213491350#TODO: insert appropriate hdata1351Gr_new=Gstgraph(G,Gr,{},legact,character) # we don't bother to fill in vertact, since it is reconstructable anyway1352Gr_new.vertact_reconstruct()13531354return Gr_new13551356#returns a tuple (L,d) of1357# * a list L =[g_0, g_1, ...] of representatives g_i of the cosets g_i H of H in G1358# * a dictionary d associating to tuples (g,n) of g in G and integers n the integer m, such that g*(g_n H) = g_m H1359#TODO: assumes first element of L is id_G1360def leftcosetaction(G,H):1361C=G.cosets(H, side='left')1362L=[C[i][0] for i in range(len(C))]1363d={}13641365for g in G:1366for i in range(len(C)):1367gtild=g*L[i]1368for j in range(len(C)):1369if gtild in C[j]:1370d[(g,i)]=j1371return (L,d)137213731374#stores list of decstratums, together with methods for scalar multiplication from left and sum1375class tautclass(object):1376r"""1377A tautological class on the moduli space \Mbar_{g,n} of stable curves.13781379Internally, it is represented by a list ``terms`` of objects of type1380``decstratum``. In most cases, a ``tautclass`` should be created using1381functions like ``kappaclass``, ``psiclass`` and1382``StableGraph.boundary_pushforward``, but it is possible to create it1383from a list of elements of type ``decstratum``.13841385EXAMPLES::13861387sage: from admcycles.admcycles import *1388sage: gamma = StableGraph([1,2],[[1,2],[3]],[(2,3)])1389sage: ds1 = decstratum(gamma, kappa=[[1],[]]); ds11390Graph : [1, 2] [[1, 2], [3]] [(2, 3)]1391Polynomial : 1*(kappa_1^1 )_01392sage: ds2 = decstratum(gamma, kappa=[[],[1]]); ds21393Graph : [1, 2] [[1, 2], [3]] [(2, 3)]1394Polynomial : 1*(kappa_1^1 )_11395sage: t = tautclass([ds1, ds2])1396sage: (t - gamma.to_tautclass() * kappaclass(1,3,1)).is_zero()1397True1398"""1399def __init__(self, terms):1400self.terms = terms14011402# returns the pullback of self under the forgetful map which forgot the markings in the list legs1403# rename: if True, checks if legs contains leg names already taken for edges1404def forgetful_pullback(self,legs,rename=True):1405r"""1406Returns the pullback of a given tautological class under the map pi : \bar M_{g,A cup B} --> \bar M_{g,A}.14071408INPUT:14091410legs : list1411List B of legs that are forgotten by the map pi.141214131414EXAMPLES::14151416sage: from admcycles import *14171418sage: psiclass(2,1,2).forgetful_pullback([3])1419Graph : [1] [[1, 2, 3]] []1420Polynomial : 1*psi_2^11421<BLANKLINE>1422Graph : [1, 0] [[1, 4], [2, 3, 5]] [(4, 5)]1423Polynomial : (-1)*1424"""1425return sum([tautclass([])]+[c.forgetful_pullback(legs,rename) for c in self.terms])14261427# computes the pushforward under the map forgetting the given legs1428def forgetful_pushforward(self,legs):1429r"""1430Returns the pushforward of a given tautological class under the map pi : \bar M_{g,n} --> \bar M_{g,{1,...,n} \ A}.14311432INPUT:14331434legs : list1435List A of legs that are forgotten by the map pi.14361437EXAMPLES::14381439sage: from admcycles import *14401441sage: s1=psiclass(3,1,3)^2;s1.forgetful_pushforward([2,3])1442Graph : [1] [[1]] []1443Polynomial : 1*14441445::14461447sage: t=tautgens(2,2,1)[1]+2*tautgens(2,2,1)[3]1448sage: t.forgetful_pushforward([1])1449Graph : [2] [[2]] []1450Polynomial : 3*1451<BLANKLINE>1452Graph : [2] [[2]] []1453Polynomial : 2*1454"""1455return tautclass([t.forgetful_pushforward(legs) for t in self.terms]).consolidate()14561457def rename_legs(self,dic,rename=False):1458for t in self.terms:1459t.rename_legs(dic,rename)1460return self14611462def __add__(self,other):1463if other==0:1464return deepcopy(self)1465new=deepcopy(self)1466new.terms+=deepcopy(other.terms)1467return new.consolidate()14681469def __neg__(self):1470return (-1)*self14711472def __sub__(self,other):1473return self + (-1)*other14741475def __iadd__(self,other):1476if other==0:1477return self1478if isinstance(other,tautclass):1479self.terms+=other.terms1480return self1481if isinstance(other,decstratum):1482self.terms.append(other)1483return self14841485def __radd__(self,other):1486if other==0:1487return deepcopy(self) #changed this adding deepcopy1488else:1489return self+other14901491def __mul__(self,other):1492return self.__rmul__(other)14931494def __rmul__(self,other):1495if isinstance(other,tautclass):1496result=tautclass([])1497for t1 in self.terms:1498for t2 in other.terms:1499result+=t1*t21500return result1501#if isinstance(other,sage.rings.integer.Integer) or isinstance(other,sage.rings.rational.Rational) or isinstance(other,int):1502else:1503new=deepcopy(self)1504for i in range(len(new.terms)):1505new.terms[i]=other*new.terms[i]1506return new15071508def __pow__(self, exponent):1509if isinstance(exponent, (Integer, Rational, int)):1510if exponent == 0:1511L = self.gnr_list()1512gnset = set((g, n) for g, n,_ in L)1513if len(gnset) == 1:1514return fundclass(*gnset.pop())1515else:1516return 11517else:1518return prod(exponent*[self])15191520def coeff_subs(self,dic):1521r"""1522If coefficients of self are polynomials, it tries to substitute variable assignments1523given by dictionary ``dic``. This is done inplace, so the class is changed by this1524operation.15251526EXAMPLES::15271528sage: from admcycles import psiclass1529sage: R.<a1, a2> = PolynomialRing(QQ,2)1530sage: t = a1 * psiclass(1,2,2) + a2 * psiclass(2,2,2); t1531Graph : [2] [[1, 2]] []1532Polynomial : a1*psi_1^11533<BLANKLINE>1534Graph : [2] [[1, 2]] []1535Polynomial : a2*psi_2^11536sage: t.coeff_subs({a2:1-a1})1537Graph : [2] [[1, 2]] []1538Polynomial : a1*psi_1^11539<BLANKLINE>1540Graph : [2] [[1, 2]] []1541Polynomial : (-a1 + 1)*psi_2^11542"""1543safetycopy=deepcopy(self)1544try:1545for t in self.terms:1546coelist=t.poly.coeff1547for i in range(len(coelist)):1548coelist[i]=coelist[i].subs(dic)1549return self1550except:1551print('Substitution failed!')1552return safetycopy155315541555def __repr__(self):1556s=''1557for i in range(len(self.terms)):1558s+=repr(self.terms[i])+'\n\n'1559return s.rstrip('\n\n')15601561def toTautvect(self,g=None,n=None,r=None):1562if g is None or n is None or r is None:1563l = self.gnr_list()1564if len(l) == 1:1565g, n, r = l[0]1566return converttoTautvect(self,g,n,r)15671568def toTautbasis(self,g=None,n=None,r=None,moduli='st'):1569r"""1570Computes vector expressing the class in a basis of the tautological ring.15711572If moduli is given, computes expression for the tautological ring of an open subset1573of \Mbar_{g,n}.15741575Options:15761577- 'st' : all stable curves1578- 'tl' : treelike curves (all cycles in the stable graph have length 1)1579- 'ct' : compact type (stable graph is a tree)1580- 'rt' : rational tails (there exists vertex of genus g)1581- 'sm' : smooth curves15821583TESTS::15841585sage: from admcycles import StableGraph, psiclass1586sage: b = StableGraph([1],[[1,2,3]],[(2,3)]).to_tautclass()1587sage: b.toTautbasis()1588(10, -10, -14)1589sage: b.toTautbasis(moduli='ct')1590(0, 0)1591sage: c = psiclass(1,2,2)**21592sage: for mod in ('st', 'tl', 'ct', 'rt', 'sm'):1593....: print(c.toTautbasis(moduli = mod))1594(0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0)1595(0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0)1596(5/6, -1/6, 1/3, 0, 0)1597(1)1598()1599"""1600if g is None or n is None or r is None:1601l = self.gnr_list()1602if len(l) == 1:1603g, n, r = l[0]1604return Tautvecttobasis(converttoTautvect(self,g,n,r),g,n,r,moduli)16051606def toprodtautclass(self,g,n):1607return prodtautclass(stgraph([g],[list(range(1, n+1))],[]), [[t] for t in deepcopy(self.terms)])16081609def dimension_filter(self):1610for i in range(len(self.terms)):1611self.terms[i].dimension_filter()1612return self.consolidate()16131614#remove all terms of degree greater than dmax1615def degree_cap(self,dmax):1616for i in range(len(self.terms)):1617self.terms[i].degree_cap(dmax)1618return self.consolidate()16191620# returns degree d part of self1621def degree_part(self,d):1622result=tautclass([])1623for i in range(len(self.terms)):1624result+=self.terms[i].degree_part(d)1625return result.consolidate()16261627def simplify(self,g=None,n=None,r=None):1628r"""1629Simplifies self by combining terms with same tautological generator, returns self.16301631EXAMPLES::16321633sage: from admcycles import psiclass1634sage: t = psiclass(1,2,1) + 11*psiclass(1,2,1); t1635Graph : [2] [[1]] []1636Polynomial : 1*psi_1^11637<BLANKLINE>1638Graph : [2] [[1]] []1639Polynomial : 11*psi_1^11640sage: t.simplify()1641Graph : [2] [[1]] []1642Polynomial : 12*psi_1^11643sage: t1644Graph : [2] [[1]] []1645Polynomial : 12*psi_1^11646"""1647L = self.gnr_list()1648if not L:1649return self1650if r is not None:1651# Assume that self is not combination from different g,n1652g, n, _ = L[0]1653self.terms= Tautv_to_tautclass(self.toTautvect(g,n,r),g,n,r).terms1654return self1655else:1656result=tautclass([])1657for gnr in L:1658result+=Tautv_to_tautclass(self.toTautvect(*gnr),*gnr)1659self.terms=result.terms1660return self16611662def simplified(self,g=None,n=None,r=None):1663r"""1664Returns a simplified version of self by combining terms with same tautological generator,1665leaving self invariant. If r is specified, only return the corresponding degree r part of self.1666"""1667L = self.gnr_list()1668if not L:1669return deepcopy(self)16701671if r is not None:1672# Assume that self is not combination from different g,n1673g, n, _ = L[0]1674return Tautv_to_tautclass(self.toTautvect(g,n,r),g,n,r)1675else:1676result=tautclass([])1677for gnr in L:1678result+=Tautv_to_tautclass(self.toTautvect(*gnr),*gnr)1679return result16801681#newself=Tautv_to_tautclass(self.toTautvect(g,n,r),g,n,r)1682#self.terms=newself.terms1683#return self16841685# uses FZ relations to convert self in a sum of the basis tautclasses1686def FZsimplify(self,r=None):1687r"""1688Returns representation of self as a tautclass formed by a linear combination of1689the preferred tautological basis.1690If r is given, only take degree r part.16911692EXAMPLES::16931694sage: from admcycles import fundclass, psiclass1695sage: t = 7*fundclass(0,4) + psiclass(1,0,4) + 3 * psiclass(2,0,4) - psiclass(3,0,4)1696sage: s = t.FZsimplify(); s1697Graph : [0] [[1, 2, 3, 4]] []1698Polynomial : 3*(kappa_1^1 )_01699<BLANKLINE>1700Graph : [0] [[1, 2, 3, 4]] []1701Polynomial : 7*1702sage: u = t.FZsimplify(r=0); u1703Graph : [0] [[1, 2, 3, 4]] []1704Polynomial : 7*1705"""1706L = self.gnr_list()1707if not L:1708return deepcopy(self)17091710if r is not None:1711# Assume that self is not combination from different g,n1712g, n, _ = L[0]1713return Tautvb_to_tautclass(self.toTautbasis(g,n,r),g,n,r)1714else:1715result=tautclass([])1716for gnr in L:1717result+=Tautvb_to_tautclass(self.toTautbasis(*gnr),*gnr)1718return result17191720def consolidate(self):1721#TODO: Check for isomorphisms of stgraphs for terms, add corresponding kppolys1722for i in range(len(self.terms) - 1, -1, -1):1723self.terms[i].consolidate()1724if not self.terms[i].poly.monom:1725self.terms.pop(i)1726return self1727# computes integral against the fundamental class of the corresponding moduli space1728# will not complain if terms are mixed degree or if some of them do not have the right codimension1729def evaluate(self):1730r"""1731Computes integral against the fundamental class of the corresponding moduli space,1732e.g. the degree of the zero-cycle part of the tautological class.17331734EXAMPLES::17351736sage: from admcycles import kappaclass1737sage: t = kappaclass(1,1,1)1738sage: t.evaluate()17391/241740"""1741return sum([t.evaluate() for t in self.terms])17421743def fund_evaluate(self, g=None, n=None):1744r"""1745Computes degree zero part of the class as multiple of the fundamental class.17461747EXAMPLES::17481749sage: from admcycles import psiclass1750sage: t = psiclass(1,2,1)1751sage: s = t.forgetful_pushforward([1]); s1752Graph : [2] [[]] []1753Polynomial : 2*1754sage: s.fund_evaluate()175521756"""1757if g is None or n is None:1758l = self.gnr_list()1759if not l:1760return 01761g, n, _ = l[0]17621763return self.toTautvect(g,n,0)[0]17641765def gnr_list(self):1766r"""1767Returns a list [(g,n,r), ...] of all genera g, number n of markings and degrees r1768for the terms appearing in the class.17691770EXAMPLES::17711772sage: from admcycles import psiclass1773sage: a = psiclass(1,2,3)1774sage: t = a + a*a1775sage: t.gnr_list()1776[(2, 3, 2), (2, 3, 1)]1777"""1778# why not a set ?1779return list(set(s for t in self.terms for s in t.gnr_list()))17801781def is_zero(self, moduli='st'):1782r"""1783Return whether this class is a known tautological relation (using1784Pixton's implementation of the generalized Faber-Zagier relations).17851786If optional argument `moduli` is given, it checks the vanishing on1787an open subset of \Mbar_{g,n}.17881789Options:17901791- 'st' : all stable curves1792- 'tl' : treelike curves (all cycles in the stable graph have length 1)1793- 'ct' : compact type (stable graph is a tree)1794- 'rt' : rational tails (there exists vertex of genus g)1795- 'sm' : smooth curves17961797EXAMPLES::17981799sage: from admcycles import kappaclass, lambdaclass1800sage: diff = kappaclass(1,3,0) - 12*lambdaclass(1,3,0)1801sage: diff.is_zero()1802False1803sage: diff.is_zero(moduli='sm')1804True1805"""1806return all(self.toTautbasis(*gnr, moduli=moduli).is_zero() for gnr in self.gnr_list())180718081809# TODO: storing the data as (unordered) lists is very inefficient!1810class KappaPsiPolynomial(object):1811r"""1812Polynomial in kappa and psi-classes on a common stable graph.18131814The data is stored as a list monomials of entries (kappa,psi) and a list1815coeff of their coefficients. Here (kappa,psi) is a monomial in kappa, psi1816classes, represented as18171818- ``kappa``: list of length self.gamma.num_verts() of lists of the form [3,0,2]1819meaning that this vertex carries kappa_1**3*kappa_3**218201821- ``psi``: dictionary, associating nonnegative integers to some legs, where1822psi[l]=3 means that there is a psi**3 at this leg for a kppoly p. The values1823can not be zero.18241825If ``p`` is such a polynomial, ``p[i]`` is of the form ``(kappa,psi,coeff)``.18261827EXAMPLES::18281829sage: from admcycles.admcycles import kppoly1830sage: p1 = kppoly([([[0, 1], [0, 0, 2]], {1:2, 2:3}), ([[], [0, 1]], {})], [-3, 5])1831sage: p11832(-3)*(kappa_2^1 )_0 (kappa_3^2 )_1 psi_1^2 psi_2^3 +5*(kappa_2^1 )_11833"""1834def __init__(self, monom, coeff):1835self.monom = monom1836self.coeff = coeff1837assert len(self.monom) == len(self.coeff)18381839def copy(self):1840r"""1841TESTS::18421843sage: from admcycles.admcycles import kppoly1844sage: p = kppoly([([[], [0, 1]], {0:3, 1:2})], [5])1845sage: p.copy()18465*(kappa_2^1 )_1 psi_0^3 psi_1^21847"""1848res = KappaPsiPolynomial.__new__(KappaPsiPolynomial)1849res.monom = [([x[:] for x in k], p.copy()) for k,p in self.monom]1850res.coeff = self.coeff[:]1851return res18521853def __iter__(self): #TODO fix1854return iter([self[i] for i in range(len(self))])18551856def __getitem__(self,i):1857return self.monom[i]+(self.coeff[i],)18581859def __len__(self):1860return len(self.monom)18611862def __neg__(self):1863return (-1)*self18641865def __add__(self,other): #TODO: Add __radd__1866r"""1867TESTS:18681869Check that mutable data are not shared between the result and the terms::18701871sage: from admcycles.admcycles import kappacl, psicl1872sage: A = kappacl(0,2,2)1873sage: B = psicl(3,2)1874sage: C = A + B1875sage: C.monom[1][1].update({2:3})1876sage: A18771*(kappa_2^1 )_01878sage: B18791*psi_3^11880"""1881if other==0:1882return self.copy()1883new = self.copy()1884for (k,p,c) in other:1885try:1886ind = new.monom.index((k,p))1887new.coeff[ind] += c1888if new.coeff[ind] == 0:1889new.monom.pop(ind)1890new.coeff.pop(ind)1891except ValueError:1892if c != 0:1893new.monom.append((k[:], p.copy()))1894new.coeff.append(c)1895return new18961897# returns the degree of the ith term of self (default: i=0). If self is empty, give back None1898def deg(self,i=0):1899if not self.monom:1900return None1901else:1902(kappa,psi)=self.monom[i]1903return sum([sum((j+1)*kvec[j] for j in range(len(kvec))) for kvec in kappa])+sum(psi.values())19041905# computes the pullback of self under a graph morphism described by dicv, dicl1906# returns a kppoly on this graph1907def graphpullback(self,dicv,dicl):1908if len(self.monom)==0:1909return kppoly([],[])19101911numvert_self=len(self.monom[0][0])1912numvert = len(dicv)19131914preim = [[] for i in range(numvert_self)]1915for v in dicv:1916preim[dicv[v]].append(v)19171918resultpoly = kppoly([], [])1919for (kappa,psi,coeff) in self:1920# pullback of psi-classes1921psipolydict = {dicl[l]: psi[l] for l in psi}1922psipoly = kppoly([([[] for i in range(numvert)],psipolydict)],[1])19231924# pullback of kappa-classes1925kappapoly = prod([prod([sum([kappacl(w, k+1, numvert) for w in preim[v]])**kappa[v][k] for k in range(len(kappa[v]))]) for v in range(numvert_self)])19261927resultpoly += coeff * psipoly * kappapoly1928return resultpoly19291930def rename_legs(self,dic):1931for count in range(len(self.monom)):1932self.monom[count]=(self.monom[count][0],{dic.get(l,l): self.monom[count][1][l] for l in self.monom[count][1]})1933return self19341935# takes old kappa-data and embeds self in larger graph with numvert vertices on vertices start, start+1, ...1936def expand_vertices(self,start,numvert):1937for count in range(len(self.monom)):1938self.monom[count] = ([[] for i in range(start)]+self.monom[count][0] + [[] for i in range(numvert-start-len(self.monom[count][0]))], self.monom[count][1])1939return self19401941def __radd__(self, other):1942if other == 0:1943return self.copy()19441945def __mul__(self,other):1946if isinstance(other,kppoly):1947return kppoly([([kappaadd(kappa1[i],kappa2[i]) for i in range(len(kappa1))],{l:psi1.get(l,0)+psi2.get(l,0) for l in set(list(psi1)+list(psi2))}) for (kappa1,psi1) in self.monom for (kappa2,psi2) in other.monom],[a*b for a in self.coeff for b in other.coeff]).consolidate()1948#if isinstance(other,sage.rings.integer.Integer) or isinstance(other,sage.rings.rational.Rational):1949else:1950return self.__rmul__(other)19511952def __rmul__(self,other):1953if isinstance(other,kppoly):1954return self.__mul__(other)1955#if isinstance(other,sage.rings.integer.Integer) or isinstance(other,sage.rings.rational.Rational) or isinstance(other,int):1956else:1957new = self.copy()1958for i in range(len(self)):1959new.coeff[i] *= other1960return new19611962def __pow__(self,exponent):1963if isinstance(exponent,Integer) or isinstance(exponent,Rational) or isinstance(exponent,int):1964return prod(exponent*[self])19651966def __repr__(self):1967s=''1968for (kappa,psi,coeff) in self:1969coestring=repr(coeff)1970if coestring.find('+')!=-1 or coestring.find('-')!=-1: # coefficient contains sum or minus (sign)1971coestring='('+coestring+')'1972s+=coestring+'*'1973count=-11974for l in kappa:1975count+=11976if len(l)==0:1977continue1978s+='('1979for j in range(len(l)):1980if l[j]==0:1981continue1982s+='kappa_'+ repr(j+1)+'^'+repr(l[j])+' '1983s+=')_' + repr(count)+' '1984for l in psi:1985if psi[l]==0:1986continue1987s+='psi_'+repr(l)+'^'+repr(psi[l])+' '1988s+='+'1989return s.rstrip('+')19901991# TODO: this should rather be called "normalize"1992def consolidate(self):1993r"""1994Remove trailing zeroes in kappa and l with psi[l]=0 and things with coeff=01995and sum up again.19961997TESTS::19981999sage: from admcycles.admcycles import kppoly2000sage: kppoly([([[], [0, 1]], {})], [5]).consolidate()20015*(kappa_2^1 )_12002sage: kppoly([([[0, 0], [0,1,0]], {})], [3]).consolidate()20033*(kappa_2^1 )_12004sage: kppoly([([[], [0,1]], {})], [0]).consolidate() # known bug200502006"""2007for (kappa,psi) in self.monom:2008for kap in kappa:2009remove_trailing_zeros(kap)2010for l in list(psi):2011if psi[l] == 0:2012psi.pop(l)2013# now we must check for newly identified duplicates2014newself=kppoly([],[])+self2015self.monom=newself.monom2016self.coeff=newself.coeff2017return self20182019kppoly = KappaPsiPolynomial20202021# some functions for entering kappa and psi classes20222023# returns (kappa_index)_vertex - must give total number of vertices2024def kappacl(vertex,index,numvert,g=None,n=None):2025if index==0:2026return (2 *g-2 +n)*onekppoly(numvert)2027if index<0:2028return kppoly([],[])2029li=[[] for i in range(numvert)]2030li[vertex]=(index-1)*[0]+[1]2031return kppoly([(li,{})],[1])20322033# returns (psi)_leg - must give total number of vertices2034def psicl(leg,numvert):2035li=[[] for i in range(numvert)]2036return kppoly([(li,{leg:1})],[1])20372038def onekppoly(numvert):2039return kppoly([([[] for cnt in range(numvert)],{})],[1])20402041# def clean(self):2042# for k in self.coeff.keys():2043# if self.coeff[k]==0:2044# self.coeff.pop(k)2045# return self2046#TODO: Multiplication with other kppoly204720482049class decstratum(object):2050r"""2051A tautological class given by a boundary stratum, decorated by a polynomial in2052kappa-classes on the vertices and psi-classes on the legs (i.e. half-edges2053and markings)20542055The internal structure is as follows20562057- ``gamma``: underlying stgraph for the boundary stratum2058- ``kappa``: list of length gamma.num_verts() of lists of the form [3,0,2]2059meaning that this vertex carries kappa_1^3*kappa_3^22060- ``psi``: dictionary, associating nonnegative integers to some legs, where2061psi[l]=3 means that there is a psi^3 at this leg20622063We adopt the convention that: kappa_a = pi_*(psi_{n+1}^{a+1}), where2064pi is the universal curve over the moduli space, psi_{n+1} is the psi-class2065at the marking n+1 that is forgotten by pi.2066"""2067def __init__(self,gamma,kappa=None,psi=None,poly=None):2068self.gamma = gamma.copy(mutable=False) # pick an immutable copy (possibly avoids copy)2069if kappa is None:2070kappa=[[] for i in range(gamma.num_verts())]2071if psi is None:2072psi={}2073if poly is None:2074self.poly=kppoly([(kappa,psi)],[1])2075else:2076if isinstance(poly,kppoly):2077self.poly=poly2078#if isinstance(poly,sage.rings.integer.Integer) or isinstance(poly,sage.rings.rational.Rational):2079else:2080self.poly=poly*onekppoly(gamma.num_verts())20812082#def deg(self):2083#return len(self.gamma.edges)+sum([sum([k[i]*(i+1) for i in range(len(k))]) for k in kappa])+sum(psi.values())20842085def copy(self, mutable=True):2086S = decstratum.__new__(decstratum)2087S.gamma = self.gamma.copy(mutable)2088S.poly = self.poly.copy()2089return S20902091def automorphism_number(self):2092r"""Returns number of automorphisms of underlying graph fixing decorations by kappa and psi-classes.2093Currently assumes that self has exaclty one nonzero term."""2094kappa,psi=self.poly.monom[0]2095g,n,r = self.gnr_list()[0]2096markings=range(1,n+1)2097num=DR.num_of_stratum(Pixtongraph(self.gamma,kappa,psi),g,r,markings)2098return DR.autom_count(num,g,r,markings,DR.MODULI_ST)20992100def rename_legs(self,dic,rename=False):2101if rename:2102legl=self.gamma.leglist()2103shift=max(legl+list(dic.values())+[0])+12104dicenh={l:dic.get(l,l+shift) for l in legl}2105else:2106dicenh=dic2107if not self.gamma.is_mutable():2108self.gamma = self.gamma.copy()2109self.gamma.rename_legs(dicenh)2110self.poly.rename_legs(dicenh)2111return self21122113def __repr__(self):2114return 'Graph : ' + repr(self.gamma) +'\n'+ 'Polynomial : ' + repr(self.poly)21152116# returns a list of lists of the form [kp1,kp2,...,kpr] where r is the number of vertices in self.gamma and kpj is a kppoly being the part of self.poly living on the jth vertex (we put coeff in the first vertex)2117def split(self):2118result=[]2119numvert = self.gamma.num_verts()2120lvdict={l:v for v in range(numvert) for l in self.gamma.legs(v, copy=False)}21212122for (kappa,psi,coeff) in self.poly:2123psidiclist=[{} for v in range(numvert)]2124for l in psi:2125psidiclist[lvdict[l]][l]=psi[l]2126newentry=[kppoly([([kappa[v]],psidiclist[v])],[1]) for v in range(numvert)]2127newentry[0]*=coeff2128result.append(newentry)2129return result21302131def __neg__(self):2132return (-1)*self21332134def __mul__(self,other):2135return self.__rmul__(other)21362137def __rmul__(self,other):2138if isinstance(other,Hdecstratum):2139return other.__mul__(self)21402141if isinstance(other,decstratum):2142# TODO: to be changed (direct access to _edges attribute)2143if other.gamma.num_edges() > self.gamma.num_edges():2144return other.__rmul__(self)2145if self.gamma.num_edges() == other.gamma.num_edges() == 0:2146# both are just vertices with some markings and some kappa-psi-classes2147# multiply those kappa-psi-classes and return the result2148return tautclass([decstratum(self.gamma, poly=self.poly*other.poly)])2149p1=self.convert_to_prodtautclass()2150p2=self.gamma.boundary_pullback(other)2151p=p1*p22152return p.pushforward()2153#if isinstance(other,sage.rings.integer.Integer) or isinstance(other,sage.rings.rational.Rational):2154if isinstance(other,tautclass):2155return tautclass([self*t for t in other.terms])2156else:2157new=self.copy()2158new.poly*=other2159return new21602161def toTautvect(self,g=None,n=None,r=None):2162return converttoTautvect(self,g,n,r)2163def toTautbasis(self,g=None,n=None,r=None):2164return Tautvecttobasis(converttoTautvect(self,g,n,r),g,n,r)21652166# remove all terms that must vanish for dimension reasons2167def dimension_filter(self):2168for i in range(len(self.poly.monom)):2169(kappa,psi)=self.poly.monom[i]2170for v in range(self.gamma.num_verts()):2171# local check: we are not exceeding dimension at any of the vertices2172if sum([(k+1)*kappa[v][k] for k in range(len(kappa[v]))]) + sum([psi[l] for l in self.gamma.legs(v, copy=False) if l in psi]) > 3 *self.gamma.genera(v) - 3 + self.gamma.num_legs(v):2173self.poly.coeff[i]=02174break2175# global check: the total codimension of the term of the decstratum is not higher than the total dimension of the ambient space2176#if self.poly.deg(i)+len(self.gamma.edges)>3*self.gamma.g()-3+len(self.gamma.list_markings()):2177# self.poly.coeff[i]=02178#break2179return self.consolidate()21802181# remove all terms of degree higher than dmax2182def degree_cap(self,dmax):2183for i in range(len(self.poly.monom)):2184if self.poly.deg(i) + self.gamma.num_edges() > dmax:2185self.poly.coeff[i]=02186return self.consolidate()21872188# returns degree d part of self2189def degree_part(self,d):2190result=self.copy()2191for i in range(len(result.poly.monom)):2192if result.poly.deg(i) + result.gamma.num_edges() != d:2193result.poly.coeff[i]=02194return result.consolidate()21952196def consolidate(self):2197self.poly.consolidate()2198return self21992200# computes integral against the fundamental class of the corresponding moduli space2201# will not complain if terms are mixed degree or if some of them do not have the right codimension2202def evaluate(self):2203answer=02204for (kappa,psi,coeff) in self.poly:2205temp=12206for v in range(self.gamma.num_verts()):2207psilist=[psi.get(l,0) for l in self.gamma.legs(v, copy=False)]2208kappalist=[]2209for j in range(len(kappa[v])):2210kappalist+=[j+1 for k in range(kappa[v][j])]2211if sum(psilist+kappalist) != 3 *self.gamma.genera(v)-3 +len(psilist):2212temp = 02213break2214temp*=DR.socle_formula(self.gamma.genera(v),psilist,kappalist)2215answer+=coeff*temp2216return answer22172218# computes the pushforward under the map forgetting the given markings2219# this returns a decstratum on the stgraph obtained by forgetting and stabilizing2220# currently works by forgetting one marking at a time2221def forgetful_pushforward(self,markings,dicv=False):2222result = self.copy(mutable=True)2223result.dimension_filter()2224vertim = list(range(self.gamma.num_verts()))2225for m in markings:2226# take current result and forget m2227v = result.gamma.vertex(m)2228if result.gamma.genera(v) == 0 and len(result.gamma.legs(v, copy=False)) == 3:2229# this vertex will become unstable, but the previous dimension filter took care that no term in result.poly has kappa or psi on vertex v2230# thus we can simply forget this vertex and go on, but we must adjust the kappa entries in result.poly22312232# If the vertex v carries another marking m' and a leg l, then while m' and l cannot have psi-data, the psi-data from the complementary leg of l must be transferred to m'2233# TODO: use more elegant (dicv,dicl)-return of .stabilize()2234vlegs = result.gamma.legs(v, copy=True)2235vlegs.remove(m)2236vmarks=set(vlegs).intersection(set(result.gamma.list_markings()))22372238if vmarks:2239mprime=list(vmarks)[0]2240vlegs.remove(mprime)2241leg=vlegs[0]2242legprime=result.gamma.leginversion(leg)2243for (kappa,psi,coeff) in result.poly:2244if legprime in psi:2245temp=psi.pop(legprime)2246psi[mprime]=temp22472248result.gamma.forget_markings([m])2249result.gamma.stabilize()2250vertim.pop(v)2251for (kappa,psi,coeff) in result.poly:2252kappa.pop(v)2253else:2254# no vertex becomes unstable, hence result.gamma is unchanged2255# we proceed to push forward on the space corresponding to marking v2256g=result.gamma.genera(v)2257n = result.gamma.num_legs(v) - 1 # n of the target space2258numvert = result.gamma.num_verts()22592260newpoly=kppoly([],[])2261for (kappa,psi,coeff) in result.poly:2262trunckappa=deepcopy(kappa)2263trunckappa[v]=[]2264truncpsi=deepcopy(psi)2265for l in result.gamma.legs(v, copy=False):2266truncpsi.pop(l,0)2267trunckppoly=kppoly([(trunckappa,truncpsi)],[coeff]) # kppoly-term supported away from vertex v22682269kappav=kappa[v]2270psiv={l:psi.get(l,0) for l in result.gamma.legs(v, copy=False) if l !=m}2271psivpoly=kppoly([([[] for i in range(numvert)],psiv)],[1])2272a_m = psi.get(m,0)22732274cposs=[]2275for b in kappav:2276cposs.append(list(range(b+1)))2277# if kappav=[3,0,1] then cposs = [[0,1,2,3],[0],[0,1]]22782279currpoly=kppoly([],[])22802281for cvect in itertools.product(*cposs):2282currpoly+=prod([binomial(kappav[i],cvect[i]) for i in range(len(cvect))])*prod([kappacl(v,i+1, numvert)**cvect[i] for i in range(len(cvect))])*psivpoly*kappacl(v,a_m+sum([(kappav[i]-cvect[i])*(i+1) for i in range(len(cvect))])-1, numvert, g=g, n=n)2283if a_m==0:2284kappapoly=prod([kappacl(v, i+1, numvert)**kappav[i] for i in range(len(kappav))])22852286def psiminuspoly(l):2287psivminus=deepcopy(psiv)2288psivminus[l]-=12289return kppoly([([[] for i in range(numvert)],psivminus)],[1])22902291currpoly+=sum([psiminuspoly(l) for l in psiv if psiv[l]>0])*kappapoly22922293newpoly+=currpoly*trunckppoly22942295result.poly=newpoly2296result.gamma.forget_markings([m])229722982299result.dimension_filter()23002301if dicv:2302return (result,{v:vertim[v] for v in range(len(vertim))})2303else:2304return result230523062307# returns a tautclass giving the pullback under the forgetful map which forgot the markings in the list newmark2308def forgetful_pullback(self,newmark,rename=True):2309if newmark==[]:2310return tautclass([self])2311if rename: #TODO: at the moment rename is not really functional, since total overhaul would be necessary at the first calling2312# make sure that markings in newmark do not already appear as leg-labels in some edge2313# TODO: to be changed (access to private attribute _maxleg)2314mleg=max(self.gamma._maxleg+1, max(newmark)+1)2315rndic={}2316rnself = self.copy()23172318for e in self.gamma.edges(copy=False):2319if e[0] in newmark:2320rndic[e[0]]=mleg2321mleg+=12322if e[1] in newmark:2323rndic[e[1]]=mleg2324mleg+=12325rnself.gamma.rename_legs(rndic)23262327#now must fix psi-data2328for l in rndic:2329for (kappa,psi,coeff) in rnself.poly:2330if l in psi:2331psi[rndic[l]]=psi.pop(l)2332else:2333rnself=self233423352336# TODO: change recursive definition to explicit, direct definition?2337# mdist = itertools.product(*[list(range(len(rnself.gamma.genera))) for j in range(len(newmark))])23382339# Strategy: first compute pullback under forgetting newmark[0], then by recursion pull this back under forgetting other elements of newmark2340a=newmark[0]2341nextnewmark=deepcopy(newmark)2342nextnewmark.pop(0)2343partpullback=[] # list of decstratums to construct tautological class obtained by forgetful-pulling-back a23442345for v in range(self.gamma.num_verts()):2346# v gives the vertex where the marking a goes23472348# first the purely kappa-psi-term2349# grv is the underlying stgraph2350# TODO: to be changed (direct access to private attributes _legs)2351grv = rnself.gamma.copy()2352grv._legs[v] += [a]2353grv.tidy_up()2354grv.set_immutable()23552356# polyv is the underlying kppoly2357polyv=kppoly([],[])2358for (kappa,psi,coeff) in rnself.poly:2359binomdist=[list(range(0, m+1)) for m in kappa[v]] # for kappa[v]=(3,0,2) gives [[0,1,2,3],[0],[0,1,2]]23602361for bds in itertools.product(*binomdist):2362kappa_new=deepcopy(kappa)2363kappa_new[v]=list(bds)23642365psi_new=deepcopy(psi)2366psi_new[a]=sum([(i+1)*(kappa[v][i]-bds[i]) for i in range(len(bds))])2367polyv=polyv+kppoly([(kappa_new,psi_new)],[((-1)**(sum([(kappa[v][i]-bds[i]) for i in range(len(bds))])))*coeff*prod([binomial(kappa[v][i],bds[i]) for i in range(len(bds))])])23682369# now the terms with boundary divisors2370# TODO: to be changed (access to private attribute _maxleg)2371rnself.gamma.tidy_up()2372newleg = max(rnself.gamma._maxleg+1, a+1) #this is the leg landing on vertex v, connecting it to the rational component below2373grl={} # dictionary: leg l -> graph with rational component attached where l was2374for l in rnself.gamma.legs(v, copy=False):2375# create a copy of rnself.gamma with leg l replaced by a rational component attached to v, containing a and l2376# TODO: use functions!!2377new_gr = rnself.gamma.copy()2378new_gr._legs[v].remove(l)2379new_gr._legs[v] += [newleg]2380new_gr._genera += [0]2381new_gr._legs += [[newleg+1, a, l]]2382new_gr._edges += [(newleg, newleg+1)]2383new_gr.tidy_up()2384new_gr.set_immutable()23852386grl[l] = new_gr23872388# contains kappa, psi polynomials attached to the graphs above2389polyl = {l:kppoly([],[]) for l in rnself.gamma.legs(v, copy=False)}23902391for (kappa,psi,coeff) in rnself.poly:2392for j in set(psi).intersection(set(rnself.gamma.legs(v, copy=False))):2393if psi[j]==0: # no real psi class here2394continue2395kappa_new=deepcopy(kappa)+[[]] # rational component, which is placed last in the graphs above, doesn't get kappa classes23962397psi_new=deepcopy(psi)2398psi_new[newleg]=psi_new.pop(j)-1 # psi**b at leg l is replaced by -psi^(b-1) at the leg connecting v to the new rational component23992400polyl[j]=polyl[j]+kppoly([(kappa_new,psi_new)],[-coeff])2401partpullback+=[decstratum(grv,poly=polyv)]+[decstratum(grl[l], poly=polyl[l]) for l in rnself.gamma.legs(v, copy=False) if len(polyl[l])>0]2402# now partpullback is list of decstratums; create tautclass from this and return pullback under nextnewmark2403T=tautclass(partpullback)2404return T.forgetful_pullback(nextnewmark)24052406# converts the decstratum self to a prodtautclass on self.gamma2407def convert_to_prodtautclass(self):2408# TODO: to be changed (direct access to private attributes _genera, _legs)2409terms = []2410for (kappa,psi,coeff) in self.poly:2411currterm=[]2412for v in range(self.gamma.num_verts()):2413psiv={(lnum+1):psi[self.gamma._legs[v][lnum]] for lnum in range(len(self.gamma._legs[v])) if self.gamma._legs[v][lnum] in psi}2414currterm.append(decstratum(stgraph([self.gamma.genera(v)], [list(range(1, len(self.gamma._legs[v])+1))], []), kappa=[kappa[v]],psi=psiv))2415currterm[0].poly*=coeff2416terms.append(currterm)2417return prodtautclass(self.gamma,terms)24182419# returns a list [(g,n,r), ...] of all genera g, number n of markings and degrees r for the terms appearing in self2420def gnr_list(self):2421g = self.gamma.g()2422n = self.gamma.n()2423e = self.gamma.num_edges()2424L = ((g, n, e + self.poly.deg(i)) for i in range(len(self.poly.monom)))2425return list(set(L))24262427#stores list of Hdecstratums, together with methods for scalar multiplication from left and sum2428class Htautclass(object):2429def __init__(self, terms):2430self.terms = terms24312432def __neg__(self):2433return (-1)*self24342435def __add__(self,other): # TODO: add += operation2436if other==0:2437return deepcopy(self)2438new=deepcopy(self)2439new.terms+=deepcopy(other.terms)2440return new.consolidate()2441def __radd__(self,other):2442if other==0:2443return self2444else:2445return self+other2446def __rmul__(self,other):2447new=deepcopy(self)2448for i in range(len(new.terms)):2449new.terms[i]=other*new.terms[i]2450return new2451def __mul__(self,other):2452if isinstance(other,tautclass):2453return sum([a*b for a in self.terms for b in other.terms])24542455def __repr__(self):2456s=''2457for i in range(len(self.terms)):2458s+=repr(self.terms[i])+'\n\n'2459return s.rstrip('\n\n')2460def consolidate(self): #TODO: Check for isomorphisms of Gstgraphs for terms, add corresponding kppolys2461for i in range(len(self.terms)):2462self.terms[i].consolidate()2463return self24642465# returns self as a prodHclass on the correct trivial graph2466# TODO: if self.terms==[], this does not produce a meaningful result2467def to_prodHclass(self):2468if self.terms==[]:2469raise ValueError('Htautclass does not know its graph, failed to convert to prodHclass')2470else:2471gamma0=stgraph([self.terms[0].Gr.gamma.g()],[self.terms[0].Gr.gamma.list_markings()],[])2472terms=[t.to_decHstratum() for t in self.terms]2473return prodHclass(gamma0,terms)24742475# computes the pushforward under the map \bar H -> M_{g',b}, sending C to C/G2476# the result is a tautclass2477def quotient_pushforward(self):2478res=tautclass([])2479for c in self.terms:2480res+=c.quotient_pushforward()2481return res24822483# pull back tautclass other on \bar M_{g',b} under the quotient map delta and return multiplication with self2484# effectively: pushforward self (divide by deg(delta)), multiply with other, then pull back entire thing2485# TODO: verify that this is allowed2486def quotient_pullback(self, other):2487if self.terms==[]:2488return deepcopy(self) #self=0, so pullback still zero2489g=self.terms[0].Gr.gamma.g()2490hdata=self.terms[0].Gr.hdata2491trivGg=trivGgraph(g,hdata)2492Gord=hdata.G.order()2493deltadeg=trivGg.delta_degree(0)24942495push=(QQ(1)/deltadeg)*self.quotient_pushforward()2496pro=push*other2497pro.dimension_filter()24982499# result = quotient pullback of pro to trivGg2500result=Htautclass([])2501edgenums=set(t.gamma.num_edges() for t in pro.terms)2502Hstrata={i:list_Hstrata(g,hdata,i) for i in edgenums}2503Hquots={i:list_quotgraphs(g,hdata,i,True) for i in edgenums}250425052506for t in pro.terms:2507# t is a decstratum on \bar M_{g',b} living on a graph t.gamma2508# list all Hstrata having quotient graph t.gamma, then pull back t.poly and scale by right multiplicity2509(r,ti,tdicv,tdicl)=deggrfind(t.gamma)2510# tdicv and tdicl send vertex/leg names in t.gamma to the names of the vertices/legs in the standard graph25112512preims=[j for j in range(len(Hstrata[r])) if Hquots[r][j][1]==ti]25132514for j in preims:2515# this corresponds to a Gstgraph, which needs to be decorated by pullbacks of t.poly and added to result2516newGr = Hstrata[r][j].copy()2517numvert = newGr.gamma.num_verts()2518multiplicity=prod([newGr.lstabilizer(e0).order() for (e0,e1) in Hquots[r][j][0].edges(copy=False)]) * t.gamma.automorphism_number() / QQ(len(equiGraphIsom(newGr,newGr)))2519(quotientgraph,inde,dicv,dicvinv,dicl,diclinv) = Hquots[r][j]25202521# dv is a dictionary associating to vertex numbers in t.gamma a preimage vertex in newGr2522dv={v: tdicv[v] for v in range(t.gamma.num_verts())}25232524newpoly=kppoly([],[])2525for (kappa,psi,coeff) in t.poly:2526newcoeff=coeff2527newkappa=[[] for u in range(numvert)]2528for v in range(len(kappa)):2529newkappa[dicvinv[tdicv[v]]]=kappa[v]2530newcoeff*=(QQ(1)/newGr.vstabilizer(dicvinv[tdicv[v]]).order())**(sum(kappa[v]))2531newpsi={diclinv[tdicl[l]]:psi[l] for l in psi}2532newcoeff*=prod([(newGr.lstabilizer(diclinv[tdicl[l]]).order())**(psi[l]) for l in psi])25332534newpoly+=kppoly([(newkappa,newpsi)],[newcoeff])253525362537newpoly*=multiplicity253825392540result.terms.append(Hdecstratum(newGr,poly=newpoly))2541return result2542254325442545# H-tautological class given by boundary stratum of Hurwitz space, i.e. a Gstgraph, decorated by a polynomial in kappa-classes on the vertices and psi-classes on the legs (i.e. half-edges and markings)2546# self.Gr = underlying Gstgraph for the boundary stratum2547# self.kappa = list of length len(self.gamma.genera) of lists of the form (3,0,2) meaning that this vertex carries kappa_1^3*kappa_3^22548# self.psi = dictionary, associating nonnegative integers to some legs, where psi[l]=3 means that there is a psi^3 at this leg2549# Convention: kappa_a = pi_*(c_1(omega_pi(sum p_i))^{a+1}), where pi is the universal curve over the moduli space, omega_pi the relative dualizing sheaf and p_i are the (images of the) sections of pi given by the marked points2550class Hdecstratum(object):2551def __init__(self,Gr,kappa=None,psi=None,poly=None):2552self.Gr=Gr2553if kappa is None:2554kappa = [[] for i in range(Gr.gamma.num_verts())]2555if psi is None:2556psi={}2557if poly is None:2558self.poly=kppoly([(kappa,psi)],[1])2559else:2560self.poly=poly2561#def deg(self):2562#return len(self.gamma.edges)+sum([sum([k[i]*(i+1) for i in range(len(k))]) for k in kappa])+sum(psi.values())25632564def __neg__(self):2565return (-1)*self25662567def __rmul__(self,other):2568if isinstance(other,decstratum):2569return self.__mul__(other)2570#if isinstance(other,sage.rings.integer.Integer) or isinstance(other,sage.rings.rational.Rational):2571else:2572new=deepcopy(self)2573new.poly*=other2574return new257525762577def to_decHstratum(self):2578result=self.Gr.to_decHstratum()2579result.poly*=self.poly2580return result25812582def __mul__(self,other):2583if isinstance(other,decstratum):2584# Step 1: find a list commdeg of generic (self,other)-Gstgraphs Gamma32585# Elements of this list will have the form (Gamma3,currentdeg) with currentdeg a list of elements (vdict1,ldict1,vdict2,ldict2), where2586# * Gamma3 is a Gstgraph with same group G as self.Gr2587# * vdict1, vdict2 are the vertex-maps from vertices of Gamma3 to self.Gr.gamma and other.gamma2588# * ldict1, ldict2 are the leg-maps from legs of self.Gr.gamma, other.gamma to Gamma32589if self.Gr.gamma.num_edges() == 0: # TODO: implement intersection with non-trivial Gstgraph2590mindeg = ceil(QQ(other.gamma.num_edges())/self.Gr.G.order()) # minimal number of edge-orbits to hope for generic graph2591maxdeg = other.gamma.num_edges()2592numotheraut = other.gamma.automorphism_number()25932594Ggrdegenerations=[list_Hstrata(self.Gr.gamma.g(),self.Gr.hdata, r) for r in range(maxdeg+1)]25952596commdeg=[]2597multiplicityvector=[] # gives for every entry in commdeg the number of degenerations isomorphic to this entry via Aut_G(Gamma3)2598for deg in range(mindeg, maxdeg+1):2599for Gamma3 in Ggrdegenerations[deg]:2600currentdeg=[]2601currentmultiplicity=[]2602otstructures=Astructures(Gamma3.gamma,other.gamma)2603if otstructures==[]:2604continue2605AutGr=equiGraphIsom(Gamma3,Gamma3)26062607# Now single out representatives of the AutGr-action on otstructures2608while otstructures:2609(dicvcurr,diclcurr)=otstructures[0]2610multiplicity=02611for (dicv,dicl) in AutGr:2612dicvnew={v:dicvcurr[dicv[v]] for v in dicv}2613diclnew={l:dicl[diclcurr[l]] for l in diclcurr}2614try:2615otstructures.remove((dicvnew,diclnew))2616multiplicity+=12617except:2618pass2619# now must test if we are dealing with a generic structure26202621coveredlegs=set()2622for g in Gamma3.G:2623for l in diclcurr.values():2624coveredlegs.add(Gamma3.legact[(g,l)])2625# in general case would need to include half-edges from self.Gr.gamma2626if set(Gamma3.gamma.halfedges()).issubset(coveredlegs):2627currentdeg += [({v:0 for v in range(Gamma3.gamma.num_verts())}, {l:l for l in self.Gr.gamma.leglist()}, dicvcurr, diclcurr)]2628multiplicity*=QQ(1)/len(AutGr) #1*numotheraut/len(AutGr)2629currentmultiplicity+=[multiplicity]2630commdeg+=[(Gamma3,currentdeg)]2631multiplicityvector+=[currentmultiplicity]26322633#Step 1 finished263426352636# Step2: go through commdeg, use vdicts to transfer kappa-classes and ldicts to transfer psi-classes2637# from self and other to Gamma3; also look for multiple edges in G-orbits on Gamma3 and put2638# the contribution (-psi-psi') coming from excess intersection there2639# TODO: figure out multiplicities from orders of orbits, automorphisms, etc.2640# TODO: eliminate early any classes that must vanish for dimension reasons2641# sum everything up and return it2642result = Htautclass([])2643for count in range(len(commdeg)):2644if commdeg[count][1]==[]:2645continue2646Gamma3 = commdeg[count][0]2647gammaresultpoly=kppoly([],[])2648for count2 in range(len(commdeg[count][1])):2649(vdict1,ldict1,vdict2,ldict2) = commdeg[count][1][count2]2650# compute preimages of vertices in self, other under the chosen morphisms to assign kappa-classes2651preim1 = [[] for i in range(self.Gr.gamma.num_verts())]2652for v in vdict1:2653preim1[vdict1[v]]+=[v]2654preim2 = [[] for i in range(other.gamma.num_verts())]2655for v in vdict2:2656preim2[vdict2[v]]+=[v]2657numvert = Gamma3.gamma.num_verts()26582659# excess intersection classes2660quotdicl={}2661quotgraph=Gamma3.quotient_graph(dicl=quotdicl)2662numpreim={l:0 for l in quotgraph.halfedges()}26632664for l in self.Gr.gamma.halfedges():2665numpreim[quotdicl[ldict1[l]]]+=12666for l in other.gamma.halfedges():2667numpreim[quotdicl[ldict2[l]]]+=126682669excesspoly=onekppoly(numvert)2670for (e0,e1) in quotgraph.edges(copy=False):2671if numpreim[e0]>1:2672excesspoly*=( (-1)*(psicl(e0,numvert)+psicl(e1,numvert)))**(numpreim[e0]-1)267326742675resultpoly=kppoly([],[])2676for (kappa1,psi1,coeff1) in self.poly:2677for (kappa2,psi2,coeff2) in other.poly:2678# pullback of psi-classes2679psipolydict={ldic1t[l]: psi1[l] for l in psi1}2680for l in psi2:2681psipolydict[ldict2[l]]=psipolydict.get(ldict2[l],0)+psi2[l]2682psipoly=kppoly([([[] for i in range(numvert)],psipolydict)],[1])26832684# pullback of kappa-classes2685kappapoly1=prod([prod([sum([kappacl(w, k+1, numvert) for w in preim1[v]])**kappa1[v][k] for k in range(len(kappa1[v]))]) for v in range(self.Gr.gamma.num_verts())])2686kappapoly2=prod([prod([sum([kappacl(w, k+1, numvert) for w in preim2[v]])**kappa2[v][k] for k in range(len(kappa2[v]))]) for v in range(other.gamma.num_verts())])2687resultpoly+=coeff1*coeff2*psipoly*kappapoly1*kappapoly226882689resultpoly*=multiplicityvector[count][count2]*excesspoly #TODO: divide by |Aut_G(Gamma3)| ??? if so, do in definition of multiplicityvector2690gammaresultpoly+=resultpoly26912692result.terms+=[Htautclass([Hdecstratum(Gamma3,poly=gammaresultpoly)])]2693return result26942695def __repr__(self):2696return 'Graph : ' + repr(self.Gr) +'\n'+ 'Polynomial : ' + repr(self.poly)26972698def consolidate(self):2699self.poly.consolidate()2700return self27012702# computes the pushforward under the map \bar H -> M_{g',b}, sending C to C/G2703# the result is a tautclass2704def quotient_pushforward(self):2705Gord=self.Gr.G.order()27062707qdicv={}2708qdicvinv={}2709qdicl={}2710quotgr=self.Gr.quotient_graph(dicv=qdicv, dicvinv=qdicvinv, dicl=qdicl)27112712preimlist = [[] for i in range(quotgr.num_verts())]2713for v in range(self.Gr.gamma.num_verts()):2714preimlist[qdicv[v]]+=[v]27152716# first, for each vertex w in quotgr compute the degree of the delta-map \bar H(preim(v)) -> \bar M_{g(v),n(v)}2717deltadeg={w: self.Gr.delta_degree(qdicvinv[w]) for w in qdicvinv}27182719# result will be decstratum on the graph quotgr2720# go through all terms of self.poly and add their pushforwards2721# for each term, go through vertices of quotient graph and collect all classes above it,2722# i.e. kappa-classes from vertex-preimages and psi-classes from leg-preimages27232724resultpoly=kppoly([],[])2725for (kappa,psi,coeff) in self.poly:2726kappanew=[]2727coeffnew=coeff2728psinew={}2729for w in range(quotgr.num_verts()):2730kappatemp=[]2731for v in preimlist[w]:2732kappatemp=kappaadd(kappatemp,kappa[v])2733kappanew+=[kappatemp]2734coeffnew*=(QQ(Gord)/len(preimlist[w]))**sum(kappatemp) # Gord/len(preimlist[w]) is the order of the stabilizer of any preimage v of w2735for l in psi:2736psinew[qdicl[l]]=psinew.get(qdicl[l],0)+psi[l]2737coeffnew*=(QQ(1)/self.Gr.character[l][1])**psi[l]2738coeffnew*=prod(deltadeg.values())2739resultpoly+=kppoly([(kappanew,psinew)],[coeffnew])2740return tautclass([decstratum(quotgr,poly=resultpoly)])27412742def remove_trailing_zeros(l):2743r"""2744Remove the trailing zeroes t the end of l.27452746EXAMPLES::27472748sage: from admcycles.admcycles import remove_trailing_zeros2749sage: l = [0, 1, 0, 0]2750sage: remove_trailing_zeros(l)2751sage: l2752[0, 1]2753sage: remove_trailing_zeros(l)2754sage: l2755[0, 1]27562757sage: l = [0, 0]2758sage: remove_trailing_zeros(l)2759sage: l2760[]2761sage: remove_trailing_zeros(l)2762sage: l2763[]2764"""2765while l and l[-1] == 0:2766l.pop()27672768# returns sum of kappa-vectors for same vertex, so [0,2]+[1,2,3]=[1,4,3]2769def kappaadd(a,b):2770if len(a)<len(b):2771aprime=a+(len(b)-len(a))*[0]2772else:2773aprime=a2774if len(b)<len(a):2775bprime=b+(len(a)-len(b))*[0]2776else:2777bprime=b2778return [aprime[i]+bprime[i] for i in range(len(aprime))]27792780# converts a Pixton-style Graph (Matrix with univariate-Polynomial entries) into a decstratum2781# assume that markings on Graph are called 1,2,3,...,n2782def Graphtodecstratum(G):2783# first care about vertices, i.e. genera and kappa-classes2784genera=[]2785kappa=[]2786for i in range(1, G.M.nrows()):2787genera.append(G.M[i,0][0])2788kappa.append([G.M[i,0][j] for j in range(1, G.M[i,0].degree()+1)])2789# now care about markings as well as edges together with the corresponding psi-classes2790legs=[[] for i in genera]2791edges=[]2792psi={}2793legname=G.M.ncols() # legs corresponding to edges are named legname, legname+1, ...; cannot collide with marking-names2794for j in range(1, G.M.ncols()):2795for i in range(1, G.M.nrows()):2796if G.M[i,j] != 0:2797break2798# now i is the index of the first nonzero entry in column j2799if G.M[0,j] != 0: # this is a marking2800legs[i-1].append(ZZ(G.M[0,j]))2801if G.M[i,j].degree()>=1:2802psi[G.M[0,j]]=G.M[i,j][1]28032804if G.M[0,j] == 0: # this is a leg, possibly self-node2805if G.M[i,j][0] == 2: #self-node2806legs[i-1]+=[legname, legname+1]2807edges.append((legname,legname+1))2808if G.M[i,j].degree()>=1:2809psi[legname]=G.M[i,j][1]2810if G.M[i,j].degree()>=2:2811psi[legname+1]=G.M[i,j][2]2812legname+=22813if G.M[i,j][0] == 1: #edge to different vertex2814if G.M.nrows()<=i+1:2815print('Attention')2816print(G.M)2817for k in range(i+1, G.M.nrows()):2818if G.M[k,j] != 0:2819break2820# now k is the index of the second nonzero entry in column j2821legs[i-1]+=[legname]2822legs[k-1]+=[legname+1]2823edges.append((legname,legname+1))2824if G.M[i,j].degree()>=1:2825psi[legname]=G.M[i,j][1]2826if G.M[k,j].degree()>=1:2827psi[legname+1]=G.M[k,j][1]2828legname+=22829return decstratum(stgraph(genera,legs,edges),kappa=kappa,psi=psi)28302831# returns the indices of a generating set of R^r \bar M_{g,n} in Pixton's all_strata(g,r,(1,..,n),MODULI_ST)2832# usually this is computed by using the FZ-relations2833# however, if r is greater than half the dimension of \bar M_{g,n} and all cohomology is tautological, we rather compute a generating set using Poincare duality and the intersection pairing28342835@cached_function2836def generating_indices(g,n,r,FZ=False):2837if FZ or r<=(3 *g-3 +n)/QQ(2) or not (g <=1 or (g==2 and n <= 19)):2838rel=matrix(DR.list_all_FZ(g,r,tuple(range(1, n+1))))2839relconv=matrix([DR.convert_vector_to_monomial_basis(rel.row(i),g,r,tuple(range(1, n+1))) for i in range(rel.nrows())])2840del(rel)2841reverseindices=list(range(relconv.ncols()))2842reverseindices.reverse()2843reordrelconv=relconv.matrix_from_columns(reverseindices)2844del(relconv)2845reordrelconv.echelonize() #(algorithm='classical')28462847reordnongens=reordrelconv.pivots()2848nongens=[reordrelconv.ncols()-i-1 for i in reordnongens]2849gens = [i for i in range(reordrelconv.ncols()) if not i in nongens]28502851# compute result of genstobasis here, for efficiency reasons2852gtob={gens[j]: vector(QQ,len(gens),{j:1}) for j in range(len(gens))}2853for i in range(len(reordnongens)):2854gtob[nongens[i]]=vector([-reordrelconv[i,reordrelconv.ncols()-j-1] for j in gens])2855genstobasis.set_cache([gtob[i] for i in range(reordrelconv.ncols())],2856*(g,n,r))2857return gens2858else:2859gencompdeg=generating_indices(g,n,3 *g-3 +n-r,True)2860maxrank=len(gencompdeg)2861currrank = 02862M=matrix(maxrank,0)2863gens=[]2864count=02865while currrank<maxrank:2866Mnew=block_matrix([[M,matrix(DR.pairing_submatrix(tuple(gencompdeg),(count,),g,3 *g-3 +n-r,tuple(range(1, n+1))))]],subdivide=False)2867if Mnew.rank()>currrank:2868gens.append(count)2869M=Mnew2870currrank=Mnew.rank()2871count+=12872return gens28732874def Hintnumbers(g,dat,indices=None, redundancy=False):2875Hbar = Htautclass([Hdecstratum(trivGgraph(g,dat),poly=onekppoly(1))])2876dimens = Hbar.terms[0].Gr.dim()2877markins=tuple(range(1, 1 + len(Hbar.terms[0].Gr.gamma.list_markings())))28782879strata=DR.all_strata(g,dimens,markins)2880decst=[Graphtodecstratum(Grr) for Grr in strata]2881intnumbers=[]28822883if indices is None or redundancy:2884effindices=list(range(len(strata)))2885else:2886effindices=indices28872888for i in effindices:2889intnumbers+=[(Hbar*(tautclass([ Graphtodecstratum(strata[i])]))).quotient_pushforward().evaluate()]28902891if not redundancy:2892return intnumbers28932894M=DR.FZ_matrix(g,dimens,markins)2895Mconv=matrix([DR.convert_vector_to_monomial_basis(M.row(i),g,dimens,markins) for i in range(M.nrows())])2896relcheck= Mconv*vector(intnumbers)2897if relcheck==2 *relcheck: # lazy way to check relcheck==(0,0,...,0)2898print('Intersection numbers are consistent, number of checks: '+repr(len(relcheck)))2899else:2900print('Intersection numbers not consistent for '+repr((g,dat.G,dat.l)))2901print('relcheck = '+repr(relcheck))29022903if indices is None:2904return intnumbers2905else:2906return [intnumbers[j] for j in indices]2907290829092910def Hidentify(g,dat,method='pull',vecout=False,redundancy=False,markings=None):2911r"""2912Identifies the (pushforward of the) fundamental class of \bar H_{g,dat} in2913terms of tautological classes.29142915INPUT:29162917- ``g`` -- integer; genus of the curve C of the cover C -> D29182919- ``dat`` -- HurwitzData; ramification data of the cover C -> D29202921- ``method`` -- string (default: `'pull'`); method of computation29222923method='pull' means we pull back to the boundary strata and recursively2924identify the result, then use linear algebra to restrict the class of \bar H2925to an affine subvectorspace (the corresponding vector space is the2926intersection of the kernels of the above pullback maps), then we use2927intersections with kappa,psi-classes to get additional information.2928If this is not sufficient, return instead information about this affine2929subvectorspace.29302931method='pair' means we compute all intersection pairings of \bar H with a2932generating set of the complementary degree and use this to reconstruct the2933class.29342935- ``vecout`` -- bool (default: `False`); return a coefficient-list with respect2936to the corresponding list of generators tautgens(g,n,r).2937NOTE: if vecout=True and markings is not the full list of markings 1,..,n, it2938will return a vector for the space with smaller number of markings, where the2939markings are renamed in an order-preserving way.29402941- ``redundancy`` -- bool (default: `True`); compute pairings with all possible2942generators of complementary dimension (not only a generating set) and use the2943redundant intersections to check consistency of the result; only valid in case2944method='pair'.29452946- ``markings`` -- list (default: `None`); return the class obtained by the2947fundamental class of the Hurwitz space by forgetting all points not in2948markings; for markings = None, return class without forgetting any points.29492950EXAMPLES::29512952sage: from admcycles.admcycles import trivGgraph, HurData, Hidentify2953sage: G = PermutationGroup([(1, 2)])2954sage: H = HurData(G, [G[1], G[1]])2955sage: t = Hidentify(2, H, markings=[]); t # Bielliptic locus in \Mbar_22956Graph : [2] [[]] []2957Polynomial : 30*(kappa_1^1 )_02958<BLANKLINE>2959Graph : [1, 1] [[6], [7]] [(6, 7)]2960Polynomial : (-9)*2961sage: G = PermutationGroup([(1, 2, 3)])2962sage: H = HurData(G, [G[1]])2963sage: t = Hidentify(2, H, markings=[]); t.is_zero() # corresponding space is empty2964True2965"""2966Hbar = trivGgraph(g,dat)2967r = Hbar.dim()2968n = len(Hbar.gamma.list_markings())29692970#global Hexam2971#Hexam=(g,dat,method,vecout,redundancy,markings)29722973if markings is None:2974markings=list(range(1, n+1))29752976N=len(markings)2977msorted=sorted(markings)2978markingdictionary={i+1:msorted[i] for i in range(N)}2979invmarkingdictionary={markingdictionary[j]:j for j in markingdictionary}29802981# Section to deal with the case that the group G is trivial quickly2982if dat.G.order()==1:2983if len(markings)<n:2984return tautclass([])2985else:2986return tautclass([decstratum(stgraph([g],[list(range(1, n+1))],[]))])298729882989if not cohom_is_taut(g,len(markings),2 *(3 *g-3 +len(markings)-r)):2990print('Careful, Hurwitz cycle '+repr((g,dat))+' might not be tautological!\n')29912992# first: try to look up if we already computed this class (or a class with even more markings) before2993found=False2994try:2995Hdb=Hdatabase[(g,dat)]2996for marks in Hdb:2997if set(markings).issubset(set(marks)):2998forgottenmarks=set(marks)-set(markings)2999result=deepcopy(Hdb[marks])3000result=result.forgetful_pushforward(list(forgottenmarks))3001found=True3002#TODO: potentially include this result in Hdb30033004except KeyError:3005pass300630073008# if not, proceed according to method given3009if method=='pair' and not found:3010if not (g <= 1 or (g==2 and n <= 19)):3011print('Careful, potentially outside of Gorenstein range!\n')3012genscompl= generating_indices(g,n,r) # indices of generators in rank with which \bar H will be paired3013gens = generating_indices(g,n,3 *g-3 +n-r) # indices of generators in rank in which \bar H sits30143015intnumbers = Hintnumbers(g,dat,indices=genscompl,redundancy=redundancy)3016Mpairing = matrix(DR.pairing_submatrix(tuple(genscompl),tuple(gens),g,r,tuple(range(1, n+1))))30173018# in Mpairing, the rows correspond to generators in genscompl, the columns to generators in gens3019# we want to solve the equation Mpairing * x = intnumbers, which should have a unique solution3020# then \bar H = sum_i x_i * generator(gens[i]) in R^(3g-3+n-r) \bar M_{g,n}3021x = Mpairing.solve_right(vector(intnumbers))30223023strata = DR.all_strata(g,3 *g-3 +n-r, tuple(range(1, n+1)))3024stratmod=[]3025for i in range(len(gens)):3026if x[i]==0:3027continue3028currstrat=Graphtodecstratum(strata[gens[i]])3029currstrat.poly*=x[i]3030stratmod.append(currstrat)30313032result= tautclass(stratmod)3033if (g,dat) not in Hdatabase:3034Hdatabase[(g,dat)]={}3035Hdatabase[(g,dat)][tuple(range(1, n+1))]=result30363037if len(markings)<n:3038result=Hidentify(g,dat,method,False,redundancy,markings)3039if method=='pull' and not found:3040gens = generating_indices(g,N,3 *g-3 +N-r) # indices of generators in rank in which forgetful_* \bar H sits3041ngens=len(gens)3042bdrydiv=list_strata(g,N,1)30433044# first compute matrix describing the pullbacks of the generators in gens to the boundary3045# NEW: omit pullback to irreducible boundary, since this leads to computing FZ-relations on spaces with many markings3046A=pullback_matrix(g,N,3 *g-3 +N-r,irrbdry=False)30473048irrpullback = False3049if g >= 1 and (A.rank() < ngens or redundancy):3050irrpullback = True3051irrbd = stgraph([g-1], [list(range(1, N+3))], [(N+1, N+2)])3052A = block_matrix([[A], [pullback_matrix(g, N, 3*g - 3 + N - r, bdry=irrbd)]], subdivide=False)30533054kpinters=False3055if A.rank() < ngens or redundancy:3056# pullback to boundary not injective, try to throw in intersection numbers with kappa-psi-polynomials of complementary degree as well3057kpinters=True3058(B,kppolynomials)=kpintersection_matrix(g,N,3 *g-3 +N-r)3059A=block_matrix([[A],[B]],subdivide=False)3060306130623063if A.rank()<ngens:3064print('Matrix rank in computation of '+repr((g,dat))+' not full!\n')306530663067# now we must compute the pullbacks of ( forgetful_* \bar H) to the boundary, collect in vector b3068# if kpinters, we also must compute the intersection numbers with polynomials in kppolynomials and concatenate to b3069b=[]3070forgetmarks=list(set(range(1, n+1))-set(markings))3071bar_H=prodHclass(stgraph([g],[list(range(1, N+1))],[]),[decHstratum(stgraph([g],[list(range(1, N+1))],[]), {0: (g,dat)}, [[0,markingdictionary]])])30723073for bdry in bdrydiv:3074if bdry.num_verts()!=1:3075b+=list(bar_H.gamma0_pullback(bdry).toprodtautclass().totensorTautbasis(3 *g-3 +N-r,vecout=True))30763077if irrpullback:3078b+=list(bar_H.gamma0_pullback(irrbd).toprodtautclass().totensorTautbasis(3 *g-3 +N-r,vecout=True))307930803081if kpinters:3082b+=[(bar_H*tautclass([kppo])).evaluate() for kppo in kppolynomials]308330843085# finally, the coefficients x of ( forgetful_* \bar H) in basis gens are the solutions of A*x=b3086# the matrix A likely has many more rows than ngens, providing some redundancy3087# if its rank is smaller than ngens, we cannot obtain ( forgetful_* \bar H) like this (at most restrict it to affine subvectorspace)3088# this CANNOT happen if both degree r and degree 3*g-3+N-r are completely tautological!3089try:3090x = A.solve_right(vector(b))3091except ValueError:3092global examA3093examA=A3094global examb3095examb=b3096global Hexam3097Hexam=(g,dat,method,vecout,redundancy,markings)3098raise ValueError('In Hidentify, matrix equation has no solution')30993100result=Tautvb_to_tautclass(x,g,N,3 *g-3 +N-r)3101result.rename_legs(markingdictionary,rename=True)31023103# return section3104if not found:3105if (g,dat) not in Hdatabase:3106Hdatabase[(g,dat)]={}3107Hdatabase[(g,dat)][tuple(markings)]=deepcopy(result)31083109if vecout:3110result.rename_legs(invmarkingdictionary,rename=True)3111return result.toTautvect(g,len(markings),3 *g-3 +len(markings)-r)3112else:3113return result31143115# if bdry is given (as a stgraph with exactly one edge) it computes the matrix whose columns are the pullback of our preferred generating_indices(g,n,d) to the divisor bdry, identified in the (tensor product) Tautbasis3116# if bdry=None, concatenate all matrices for all possible boundary divisors (in their order inside list_strata(g,n,1))3117# if additionally irrbdry=False, omit the pullback to the irreducible boundary, since there computations get more complicated3118def pullback_matrix(g,n,d,bdry=None,irrbdry=True):3119genindices = generating_indices(g,n,d) # indices of generators corresponding to columns of A3120strata = DR.all_strata(g, d, tuple(range(1, n+1)))3121gens=[Graphtodecstratum(strata[i]) for i in genindices]3122ngens=len(genindices)31233124if bdry is None:3125A=matrix(0,ngens)3126bdrydiv=list_strata(g,n,1)31273128for bdry in bdrydiv:3129if irrbdry or bdry.num_verts() != 1:3130A=block_matrix([[A],[pullback_matrix(g,n,d,bdry)]],subdivide=False)3131return A3132else:3133pullbacks=[bdry.boundary_pullback(cl).totensorTautbasis(d,True) for cl in gens]3134return matrix(pullbacks).transpose()313531363137def kpintersection_matrix(g,n,d):3138r"""3139Computes the matrix B whose columns are the intersection numbers of our preferred3140generating_indices(g,n,d) with the list kppolynomials of kappa-psi-polynomials of3141opposite degree. Returns the tuple (B,kppolynomials), where the latter are decstrata.31423143TESTS::31443145sage: from admcycles.admcycles import kpintersection_matrix3146sage: kpintersection_matrix(0,5,2)[0]3147[1]3148"""3149genindices = generating_indices(g,n,d) # indices of generators corresponding to columns of A3150strata = DR.all_strata(g,d,tuple(range(1, n+1)))3151opstrata = DR.all_strata(g,3 *g-3 +n-d,tuple(range(1, n+1)))3152opgenindices = list(range(len(opstrata))) #generating_indices(g,n,3*g-3+n-d) # indices of generators in opposite degree31533154gens=[Graphtodecstratum(strata[i]) for i in genindices]3155ngens=len(genindices)3156kppolynomials=[Graphtodecstratum(opstrata[i]) for i in opgenindices]3157kppolynomials=[s for s in kppolynomials if s.gamma.num_edges()==0]31583159B=[[(gen*s).evaluate() for gen in gens] for s in kppolynomials]3160return (matrix(B),kppolynomials)31613162class prodtautclass(object):3163r"""3164A tautological class on the product of several spaces \bar M_{g_i, n_i},3165which correspond to the vertices of a stable graph.31663167One way to construct such a tautological class is to pull back an other3168tautological class under a boundary gluing map.31693170Internally, the product class is stored as a list ``self.terms``, where3171each entry is of the form ``[ds(0),ds(1), ..., ds(m)]``, where ``ds(j)``3172are decstratums.31733174Be careful that the marking-names on the ``ds(j)`` are 1,2,3, ...3175corresponding to legs 0,1,2, ... in gamma.legs(j).3176If argument prodtaut is given, it is a list (of length gamma.num_verts())3177of tautclasses with correct leg names and this should produce the3178prodtautclass given as pushforward of the product of these classes that is,3179``self.terms`` is the list of all possible choices of a decstratum from the3180various factors.3181"""3182def __init__(self, gamma, terms=None, protaut=None):3183self.gamma = gamma.copy(mutable=False)3184if terms is not None:3185self.terms = terms3186elif protaut is not None:3187declists=[t.terms for t in protaut]3188self.terms=[deepcopy(list(c)) for c in itertools.product(*declists)]3189else:3190self.terms = [[decstratum(StableGraph([self.gamma.genera(i)], [list(range(1, self.gamma.num_legs(i)+1))], [])) for i in range(self.gamma.num_verts())]]31913192# TODO: implement a copy function instead of relying on deepcopy!3193# (e.g. there is no need to copy the stable graph)31943195# returns the pullback of self under the forgetful map which forgot the markings in the list newmark3196# rename: if True, checks if newmark contains leg names already taken for edges3197#def forgetful_pullback(self,newmark,rename=True):3198# return sum([c.forgetful_pullback(newmark,rename) for c in self.terms])31993200# returns the tautclass obtained by pushing forward under the gluing map associated to gamma3201def pushforward(self):3202result=tautclass([])3203for t in self.terms:3204# we have to combine the graphs from all terms such that there are no collisions of half-edge names3205# * genera will be the concatenation of all t[i].gamma.genera3206# * legs will be the concatenation of renamed versions of the t[i].gamma.legs3207# * edges will be the union of the renamed edge-lists from the elements t[i] and the edge-list of self.gamma3208maxleg=max(self.gamma.leglist()+[0])+13209genera=[]3210legs=[]3211# TODO: to be changed (direct access to "hidden" attributes)3212edges = self.gamma.edges()3213numvert = sum(s.gamma.num_verts() for s in t) # total number of vertices in new graph3214poly = onekppoly(numvert)32153216for i in range(len(t)):3217# s is a decstratum3218s = t[i]3219gr = s.gamma3220start = len(genera)3221genera.extend(gr.genera())32223223# create the renaming-dictionary3224# TODO: to be changed (direct access to _legs and _edges attribute)3225renamedic = {j+1: self.gamma._legs[i][j] for j in range(len(self.gamma._legs[i]))}3226for (e0,e1) in gr._edges:3227renamedic[e0] = maxleg3228renamedic[e1] = maxleg+13229maxleg += 232303231legs += [[renamedic[l] for l in ll] for ll in gr._legs]3232edges += [(renamedic[e0], renamedic[e1]) for (e0,e1) in gr._edges]32333234grpoly = deepcopy(s.poly)3235grpoly.rename_legs(renamedic)3236grpoly.expand_vertices(start,numvert)3237poly*=grpoly3238result+=tautclass([decstratum(stgraph(genera,legs,edges),poly=poly)])3239return result32403241# given a stgraph gamma0 and (dicv,dicl) a morphism from self.gamma to gamma0 (i.e. self.gamma is a specialization of gamma0)3242# returns the prodtautclass on gamma0 obtained by the partial pushforward along corresponding gluing maps3243def partial_pushforward(self,gamma0,dicv,dicl):3244gamma=self.gamma3245preimv=[[] for v in gamma0.genera(copy=False)]3246for w in dicv:3247preimv[dicv[w]].append(w)3248usedlegs=dicl.values()32493250# now create list of preimage-graphs, where leg-names and orders are preserved so that decstratums can be used unchanged3251preimgr=[]3252for v in range(gamma0.num_verts()):3253preimgr.append(stgraph([gamma.genera(w) for w in preimv[v]], [gamma.legs(w) for w in preimv[v]], [e for e in gamma.edges(copy=False) if (gamma.vertex(e[0]) in preimv[v] and gamma.vertex(e[1]) in preimv[v]) and not (e[0] in usedlegs or e[1] in usedlegs)]))32543255# now go through self.terms, split decstrata onto preimgr (creating temporary prodtautclasses)3256# push those forward obtaining tautclasses on the moduli spaces corresponding to the vertices of gamma0 (but wrong leg labels)3257# rename them according to the inverse of dicl, then multiply all factors together using distributivity3258result=prodtautclass(gamma0,[])3259rndics=[{dicl[gamma0.legs(v, copy=False)[j]]:j+1 for j in range(gamma0.num_legs(v))} for v in range(gamma0.num_verts())]3260for t in self.terms:3261tempclasses=[prodtautclass(preimgr[v],[[t[i] for i in preimv[v]]]) for v in range(gamma0.num_verts())]3262temppush=[s.pushforward() for s in tempclasses]3263for v in range(gamma0.num_verts()):3264temppush[v].rename_legs(rndics[v],rename=True)3265for comb in itertools.product(*[tp.terms for tp in temppush]):3266result.terms.append(list(comb))3267return result32683269# converts self into a prodHclass on gamma0=self.gamma, all spaces being distinct with trivial Hurwitz datum (for the trivial group) and not having any identifications among themselves3270def toprodHclass(self):3271gamma0 = self.gamma.copy(mutable=False)3272terms=[]3273trivgp=PermutationGroup([()])3274trivgpel=trivgp[0]3275result=prodHclass(gamma0,[])32763277for t in self.terms:3278# we have to combine the graphs from all terms such that there are no collisions of half-edge names3279# * genera will be the concatenation of all t[i].gamma.genera3280# * legs will be the concatenation of renamed versions of the t[i].gamma.legs3281# * edges will be the union of the renamed edge-lists from the elements t[i] and the edge-list of self.gamma3282dicv0={}3283dicl0={l:l for l in self.gamma.leglist()}3284spaces={}3285vertdata=[]32863287maxleg=max(self.gamma.leglist()+[0])+13288genera=[]3289legs=[]3290edges = self.gamma.edges()3291numvert = sum(s.gamma.num_verts() for s in t) # total number of vertices in new graph3292poly = onekppoly(numvert)32933294for i in range(len(t)):3295# s is a decstratum3296s = t[i]3297gr = s.gamma3298start = len(genera)3299genera += gr.genera(copy=False)3300dicv0.update({j:i for j in range(start, len(genera))})33013302# create the renaming-dictionary3303renamedic = {j+1: self.gamma.legs(i, copy=False)[j] for j in range(self.gamma.num_legs(i))}3304for (e0,e1) in gr.edges():3305renamedic[e0]=maxleg3306renamedic[e1]=maxleg+13307maxleg+=233083309legs+=[[renamedic[l] for l in ll] for ll in gr._legs]3310edges+=[(renamedic[e0],renamedic[e1]) for (e0,e1) in gr._edges]3311spaces.update({j:[genera[j],HurwitzData(trivgp,[(trivgpel,1,0) for counter in range(len(legs[j]))])] for j in range(start, len(genera))})3312vertdata+=[[j,{legs[j][lcount]:lcount+1 for lcount in range(len(legs[j]))}] for j in range(start,len(genera))]33133314grpoly=deepcopy(s.poly)3315grpoly.rename_legs(renamedic)3316grpoly.expand_vertices(start,numvert)3317poly*=grpoly33183319gamma=stgraph(genera,legs,edges)33203321result.terms.append(decHstratum(gamma,spaces,vertdata,dicv0,dicl0,poly))3322return result332333243325# takes a set of vertices of self.gamma and a prodtautclass prodcl on a graph with the same number of vertices as the length of this list - returns the multiplication of self with a pullback (under a suitable projection) of prodcl3326# more precisely: let vertices=[v_1, ..., v_r], then for j=1,..,r the vertex v_j in self.gamma has same genus and number of legs as vertex j in prodcl.gamma3327# multiply the pullback of these classes accordingly3328def factor_pullback(self,vertices,prodcl):3329other=prodtautclass(self.gamma) # initialize with 1-class3330defaultterms = [other.terms[0][i] for i in range(self.gamma.num_verts())] #collect 1-terms on various factors here3331other.terms=[]3332for t in prodcl.terms:3333newterm=deepcopy(defaultterms)3334for j in range(len(vertices)):3335newterm[vertices[j]]=t[j]3336other.terms.append(newterm)3337return self*other33383339def __neg__(self):3340return (-1)*self33413342def __add__(self,other):3343if other==0:3344return deepcopy(self)3345new=deepcopy(self)3346new.terms+=deepcopy(other.terms)3347return new.consolidate()33483349def __radd__(self,other):3350if other==0:3351return self3352else:3353return self+other33543355def __iadd__(self,other):3356if other==0:3357return self3358if not isinstance(other,prodtautclass) or other.gamma != self.gamma:3359raise ValueError("summing two prodtautclasses on different graphs!")3360self.terms+=other.terms3361return self33623363def __mul__(self,other):3364return self.__rmul__(other)33653366def __imul__(self,other):3367if isinstance(other,prodtautclass):3368self.terms=(self.__rmul__(other)).terms3369return self3370#if isinstance(other,sage.rings.integer.Integer) or isinstance(other,sage.rings.rational.Rational) or isinstance(other,int):3371else:3372for t in self.terms:3373t[0]*=other3374return self33753376def __rmul__(self,other):3377if isinstance(other,prodtautclass):3378# multiply the decstratums on each vertex together separately3379if other.gamma != self.gamma:3380raise ValueError("product of two prodtautclass on different graphs!")33813382result = prodtautclass(self.gamma,[])3383for T1 in self.terms:3384for T2 in other.terms:3385# go through vertices, multiply decstratums -> obtain tautclasses3386# in the end: must output all possible combinations of those tautclasses3387vertextautcl=[]3388for v in range(self.gamma.num_verts()):3389vertextautcl.append((T1[v]*T2[v]).terms)3390for choice in itertools.product(*vertextautcl):3391result+=prodtautclass(self.gamma,[list(choice)])3392return result3393#if isinstance(other,sage.rings.integer.Integer) or isinstance(other,sage.rings.rational.Rational) or isinstance(other,int):3394else:3395# new = self.copy()3396new = deepcopy(self)3397for t in new.terms:3398t[0]=other*t[0]3399return new34003401def dimension_filter(self):3402for t in self.terms:3403for s in t:3404s.dimension_filter()3405return self.consolidate()34063407def consolidate(self):3408revrange=list(range(len(self.terms)))3409revrange.reverse()3410for i in revrange:3411for s in self.terms[i]:3412if len(s.poly.monom)==0:3413self.terms.pop(i)3414break3415return self34163417def __repr__(self):3418s='Outer graph : ' + repr(self.gamma) +'\n'3419for i in range(len(self.terms)):3420for j in range(len(self.terms[i])):3421s+='Vertex '+repr(j)+' :\n'+repr(self.terms[i][j])+'\n'3422s+='\n\n'3423return s.rstrip('\n\n')34243425def factor_reconstruct(self,m, otherfactors):3426r"""3427Assuming that the prodtautclass is a product v_0 x ... x v_(s-1) of tautclasses v_i on the3428vertices, this returns the mth factor v_m assuming that otherfactors=[v_1, ..., v_(m-1), v_(m+1), ... v_(s-1)]3429is the list of factors on the vertices not equal to m.34303431EXAMPLES::34323433sage: from admcycles.admcycles import StableGraph, prodtautclass, psiclass, kappaclass, fundclass3434sage: Gamma = StableGraph([1,2,1],[[1,2],[3,4],[5]],[(2,3),(4,5)])3435sage: pt = prodtautclass(Gamma,protaut=[psiclass(1,1,2),kappaclass(2,2,2),fundclass(1,1)])3436sage: res = pt.factor_reconstruct(1,[psiclass(1,1,2),fundclass(1,1)])3437sage: res3438Graph : [2] [[1, 2]] []3439Polynomial : 1*(kappa_2^1 )_034403441NOTE::34423443This method works by computing intersection numbers in the other factors. This could be replaced3444with computing a basis in the tautological ring of these factors.3445In principle, it is also possible to reconstruct all factors v_i (up to scaling) just assuming3446that the prodtautclass is a pure tensor (product of pullbacks from factors).3447"""3448s = self.gamma.num_verts()3449if not len(otherfactors) == s-1:3450raise ValueError('otherfactors must have exactly one entry for each vertex not equal to m')3451otherindices = list(range(s))3452otherindices.remove(m)34533454othertestclass={}3455otherintersection=[]3456for i in range(s-1):3457gv, nv, rv = otherfactors[i].gnr_list()[0]3458rvcomplement = 3*gv-3+nv-rv3459for tc in tautgens(gv,nv,rvcomplement):3460intnum = (tc*otherfactors[i]).evaluate()3461if intnum != 0:3462othertestclass[otherindices[i]]=tc3463otherintersection.append(intnum)3464break3465result = tautclass([])3466for t in self.terms:3467result+=prod([(t[oi]*othertestclass[oi]).evaluate() for oi in otherindices])*tautclass([t[m]])3468result.simplify()3469return 1/prod(otherintersection,QQ(1)) * result347034713472# returns the cohomological degree 2r part of self, where we see self as a cohomology vector in3473# H*(\bar M_{g1,n1} x \bar M_{g2,n2} x ... x \bar M_{gm,nm}) with m the number of vertices of self.gamma3474# currently only implemented for m=1 (returns vector) or m=2 (returns list of matrices of length r+1)3475# if vecout=True, in the case m=2 a vector is returned that is the concatenation of all row vectors of the matrices3476# TODO: make r optional argument?3477def totensorTautbasis(self,r,vecout=False):3478if self.gamma.num_verts() == 1:3479g = self.gamma.genera(0)3480n = self.gamma.num_legs(0)3481result = vector(QQ,len(generating_indices(g,n,r)))3482for t in self.terms:3483# t is a one-element list containing a decstratum3484result+=t[0].toTautbasis(g,n,r)3485return result3486if self.gamma.num_verts() == 2:3487g1 = self.gamma.genera(0)3488n1 = self.gamma.num_legs(0)3489g2 = self.gamma.genera(1)3490n2 = self.gamma.num_legs(1)3491rmax1 = 3*g1 - 3 + n13492rmax2 = 3*g2 - 3 + n234933494rmin=max(0,r-rmax2) # the degree r1 in the first factor can only vary between rmin and rmax3495rmax=min(r,rmax1)34963497result=[matrix(QQ,0,0) for i in range(rmin)]34983499for r1 in range(rmin, rmax+1):3500M=matrix(QQ,len(generating_indices(g1,n1,r1)),len(generating_indices(g2,n2,r-r1)))3501for t in self.terms:3502# t is a list [d1,d2] of two decstratums3503vec1=t[0].toTautbasis(g1,n1,r1)3504vec2=t[1].toTautbasis(g2,n2,r-r1)3505M+=matrix(vec1).transpose()*matrix(vec2)3506result.append(M)35073508result+=[matrix(QQ,0,0) for i in range(rmax+1, r+1)]35093510if vecout:3511return vector([M[i,j] for M in result for i in range(M.nrows()) for j in range(M.ncols())])3512else:3513return result35143515if self.gamma.num_verts() > 2:3516print('totensorTautbasis not yet implemented on graphs with more than two vertices')3517return 035183519# converts a decstratum or tautclass to a vector in Pixton's basis3520# tries to reconstruct g,n,r from first nonzero term on the first decstratum3521# if they are explicitly given, only extract the part of D that lies in right degree, ignore others3522def converttoTautvect(D,g=None,n=None,r=None):3523#global copyD3524#copyD=(D,g,n,r)3525if isinstance(D,tautclass):3526if len(D.terms)==0:3527if (g is None or n is None or r is None):3528print('Unable to identify g, n, r for empty tautclass')3529return 03530else:3531return vector(QQ,DR.num_strata(g,r,tuple(range(1, n+1))))3532if g is None:3533g=D.terms[0].gamma.g()3534if n is None:3535n=len(D.terms[0].gamma.list_markings())3536if r is None:3537polydeg=D.terms[0].poly.deg()3538if polydeg is not None:3539r = D.terms[0].gamma.num_edges() + polydeg3540result=converttoTautvect(D.terms[0],g,n,r)3541for i in range(1, len(D.terms)):3542result+=converttoTautvect(D.terms[i],g,n,r)3543return result3544if isinstance(D,decstratum):3545D.dimension_filter()3546if g is None:3547g=D.gamma.g()3548if n is None:3549n=len(D.gamma.list_markings())3550if len(D.poly.monom)==0:3551if (r is None):3552print('Unable to identify r for empty decstratum')3553return 03554else:3555return vector(QQ,DR.num_strata(g,r,tuple(range(1, n+1))))3556if r is None:3557polydeg=D.poly.deg()3558if polydeg is not None:3559r = D.gamma.num_edges() + polydeg3560# just assume now that g,n,r are given3561length=DR.num_strata(g,r,tuple(range(1, n+1)))3562markings=tuple(range(1, n+1))3563result = vector(QQ,length)3564graphdegree = D.gamma.num_edges()3565for (kappa,psi,coeff) in D.poly:3566if graphdegree + sum([sum((j+1)*kvec[j] for j in range(len(kvec))) for kvec in kappa])+sum(psi.values()) == r:3567try:3568pare=coeff.parent()3569except:3570pare=QQ3571result+=vector(pare,length,{DR.num_of_stratum(Pixtongraph(D.gamma,kappa,psi),g,r,markings):coeff})3572return result35733574# returns a Pixton-style graph given a stgraph G and data kappa and psi as in kppoly3575def Pixtongraph(G,kappa,psi):3576firstcol=[-1]3577for v in range(G.num_verts()):3578entry = G.genera(v)3579for k in range(len(kappa[v])):3580entry+=kappa[v][k]*(DR.X)**(k+1)3581firstcol.append(entry)3582columns=[firstcol]35833584for l in G.list_markings():3585vert=G.vertex(l)3586entry = 13587if l in psi:3588entry+=psi[l]*DR.X3589columns.append([l]+[0 for i in range(vert)]+[entry]+[0 for i in range(G.num_verts()-vert-1)])3590for (e0,e1) in G.edges(copy=False):3591v0=G.vertex(e0)3592v1=G.vertex(e1)3593if v0==v1:3594entry=23595if (psi.get(e0,0)!=0) and (psi.get(e1,0)==0):3596entry+=psi[e0]*DR.X3597if (psi.get(e1,0)!=0) and (psi.get(e0,0)==0):3598entry+=psi[e1]*DR.X3599if (psi.get(e0,0)!=0) and (psi.get(e1,0)!=0):3600if psi[e0]>=psi[e1]:3601entry+=psi[e0]*DR.X + psi[e1]*(DR.X)**23602else:3603entry+=psi[e1]*DR.X + psi[e0]*(DR.X)**23604columns.append([0]+[0 for i in range(v0)]+[entry]+[0 for i in range(G.num_verts()-v0-1)])3605else:3606col = [0 for i in range(G.num_verts()+1)]3607entry = 13608if e0 in psi:3609entry+=psi[e0]*DR.X3610col[v0+1]=entry36113612entry = 13613if e1 in psi:3614entry+=psi[e1]*DR.X3615col[v1+1]=entry3616columns.append(col)36173618M=matrix(columns).transpose()3619return DR.Graph(M)36203621# converts vector v in generating set all_strata(g,r,(1,..,n)) to preferred basis computed by generating_indices(g,n,r)3622def Tautvecttobasis(v,g,n,r,moduli='st'):3623if (2 *g-2 +n<=0) or (r<0) or (r>3 *g-3 +n):3624return vector([])3625vecs=genstobasis(g,n,r)3626res=vector(QQ,len(vecs[0]))3627for i in range(len(v)):3628if v[i]!=0:3629res+=v[i]*vecs[i]3630if moduli == 'st':3631return res3632else:3633W = op_subset_space(g,n,r,moduli)3634if not res:3635return W.zero()3636reslen = len(res)3637return sum(res[i] * W(vector(QQ,reslen,{i:1})) for i in range(reslen) if not res[i].is_zero())36383639# gives a list whose ith entry is the representation of the ith generator in all_strata(g,r,(1,..,n)) in the preferred basis computed by generating_indices(g,n,r)36403641@cached_function3642def genstobasis(g, n, r):3643generating_indices(g, n, r, True)3644return genstobasis(g, n, r)364536463647# converts Tautvector to tautclass3648def Tautv_to_tautclass(v,g,n,r):3649strata = DR.all_strata(g,r, tuple(range(1, n+1)))3650stratmod=[]3651for i in range(len(v)):3652if v[i]==0:3653continue3654currstrat=Graphtodecstratum(strata[i])3655currstrat=v[i]*currstrat3656stratmod.append(currstrat)36573658return tautclass(stratmod)36593660# converts vector in generating_indices-basis to tautclass3661def Tautvb_to_tautclass(v,g,n,r):3662genind=generating_indices(g,n,r)3663strata = DR.all_strata(g, r, tuple(range(1, n+1)))3664stratmod=[]3665for i in range(len(v)):3666if v[i]==0:3667continue3668currstrat=Graphtodecstratum(strata[genind[i]])3669currstrat.poly*=v[i]3670stratmod.append(currstrat)36713672return tautclass(stratmod)367336743675class prodHclass(object):3676r"""3677A sum of gluing pushforwards of pushforwards of fundamental classes of3678Hurwitz spaces under products of forgetful morphisms.36793680This is all relative to a fixed stable graph gamma0.3681"""3682def __init__(self, gamma0, terms):3683self.gamma0 = gamma0.copy(mutable=False)3684self.terms = terms36853686def __neg__(self):3687return (-1)*self36883689def __add__(self,other):3690if other==0:3691return deepcopy(self)3692if isinstance(other,prodHclass) and other.gamma0==self.gamma0:3693new=deepcopy(self)3694new.terms+=deepcopy(other.terms)3695return new.consolidate()36963697def __radd__(self,other):3698if other==0:3699return self3700else:3701return self+other37023703def __iadd__(self,other):3704if other==0:3705return self3706if isinstance(other,prodHclass) and other.gamma0==self.gamma0:3707self.terms+=other.terms #TODO: should we use a deepcopy here?3708else:3709raise ValueError("sum of prodHclasses on different graphs")3710return self37113712def __mul__(self,other):3713if isinstance(other,tautclass) or isinstance(other,decstratum):3714return self.__rmul__(other)3715#if isinstance(other,sage.rings.integer.Integer) or isinstance(other,sage.rings.rational.Rational) or isinstance(other,int):3716else:3717new=deepcopy(self)3718for t in new.terms:3719t[0]=other*t[0]3720return new372137223723def __rmul__(self,other):3724if isinstance(other,tautclass) or isinstance(other,decstratum):3725pbother = self.gamma0.boundary_pullback(other) # this is a prodtautclass on self.gamma0 now3726pbother = pbother.toprodHclass() # this is a prodHclass with gamma0 = self.gamma0 now, but we know that all spaces are \bar M_{g_i,n_i}, which we use below!3727result = prodHclass(deepcopy(self.gamma0),[])3728for t in pbother.terms:3729# t is a decHstratum and (t.dicv0,t.dicl0) is a morphism from t.gamma to self.gamma03730# we pull back self along this morphism to get a prodHclass on t.gamma3731temppullback = self.gamma0_pullback(t.gamma,t.dicv0,t.dicl0)37323733# now it remains to distribute the kappa and psi-classes from t.poly (living on t.gamma) to the terms of temppullback3734for s in temppullback.terms:3735# s is now a decHstratum and (s.dicv0,s.dicl0) describes a map s.gamma -> t.gamma3736s.poly*=t.poly.graphpullback(s.dicv0,s.dicl0)3737# rearrange s to be again a decHstratum with respect to gamma03738s.dicv0={v: t.dicv0[s.dicv0[v]] for v in s.dicv0}3739s.dicl0={l: s.dicl0[t.dicl0[l]] for l in t.dicl0}3740result.terms+=temppullback.terms3741return result3742#if isinstance(other,sage.rings.integer.Integer) or isinstance(other,sage.rings.rational.Rational) or isinstance(other,int):3743else:3744new=deepcopy(self)3745for t in new.terms:3746t[0]=other*t[0]3747return new374837493750# evaluates self against the fundamental class of the ambient space (of self.gamma0), returns a rational number3751def evaluate(self):3752return sum([t.evaluate() for t in self.terms])3753375437553756# tries to convert self into a prodtautclass on the graph gamma03757# Note: since not all Hurwitz spaces have tautological fundamental class and since not all diagonals have tautological Kunneth decompositions, this will not always work3758def toprodtautclass(self):3759result = prodtautclass(self.gamma0, [])3760for t in self.terms:3761# create a prodtautclass on t.gamma, then do partial pushforward under (t.dicv0,t.dicl0)3762tempres=decstratum(deepcopy(t.gamma),poly=deepcopy(t.poly))3763tempres=tempres.convert_to_prodtautclass()37643765# now group vertices of t.gamma according to the spaces they access3766# create diagonals and insert the identification of the Hurwitz stack fundamental class3767# make sure that for this you use the version forgetting the most markings possible3768spacelist={a:[] for a in t.spaces}3769for v in range(t.gamma.num_verts()):3770spacelist[t.vertdata[v][0]].append(v)3771for a in spacelist:3772if spacelist[a]==[]:3773spacelist.pop(a)37743775for a in spacelist:3776# find out which of the markings 1, .., N of the space containing the Hurwitz cycle are used3777usedmarks=set()3778for v in spacelist[a]:3779usedmarks.update(set(t.vertdata[v][1].values()))3780usedmarks=list(usedmarks)3781usedmarks.sort()3782rndic={usedmarks[i]:i+1 for i in range(len(usedmarks))}37833784Hclass=Hidentify(t.spaces[a][0],t.spaces[a][1],markings=usedmarks)3785Hclass.rename_legs(rndic,rename=True)37863787legdics = [{l: rndic[t.vertdata[v][1][l]] for l in t.gamma.legs(v, copy=False)} for v in spacelist[a]]3788diagclass = forgetful_diagonal(t.spaces[a][0], len(usedmarks), [t.gamma.legs(v, copy=False) for v in spacelist[a]], legdics, T=Hclass)37893790tempres=tempres.factor_pullback(spacelist[a],diagclass)379137923793tempres=tempres.partial_pushforward(self.gamma0,t.dicv0,t.dicl0)3794result+=tempres3795return result37963797# given a different graph gamma1 carrying a gamma0-structure encoded by dicv, dicl, pull back self under this structure and return a prodHclass with underlying graph gamma13798def gamma0_pullback(self,gamma1,dicv=None,dicl=None):3799gamma0=self.gamma038003801if dicv is None:3802if gamma0.num_verts() == 1:3803dicv={v:0 for v in range(gamma1.num_verts())}3804else:3805raise RuntimeError('dicv not uniquely determined')3806if dicl is None:3807if gamma0.num_edges() == 0:3808dicl={l:l for l in gamma0.leglist()}3809else:3810raise RuntimeError('dicl not uniquely determined')38113812# Step 1: factor the map gamma1 -> gamma0 in a series of 1-edge degenerations3813# TODO: choose a smart order to minimize computations later (i.e. start w/ separating edges with roughly equal genus on both sides38143815imdicl=dicl.values()3816extraedges=[e for e in gamma1.edges(copy=False) if e[0] not in imdicl]3817delta_e = gamma1.num_edges() - gamma0.num_edges()3818if delta_e != len(extraedges):3819print('Warning: edge numbers')3820global exampullb3821exampullb=(deepcopy(self),deepcopy(gamma1),deepcopy(dicv),deepcopy(dicl))3822raise ValueError('Edge numbers')3823edgeorder=list(range(delta_e))3824contredges=[extraedges[i] for i in edgeorder]38253826contrgraph=[gamma1]3827actvert=[]3828edgegraphs=[]3829vnumbers=[]3830contrdicts=[]3831vimages=[dicv[v] for v in range(gamma1.num_verts())]3832count=03833for e in contredges:3834gammatild = contrgraph[count].copy()3835(av,edgegraph,vnum,diccv) = gammatild.contract_edge(e,adddata=True)3836gammatild.set_immutable()3837contrgraph.append(gammatild)3838actvert.append(av)3839edgegraphs.append(edgegraph)3840vnumbers.append(vnum)3841contrdicts.append(diccv)3842if len(vnum)==2:3843vimages.pop(vnum[1])3844count+=13845contrgraph.reverse()3846actvert.reverse()3847edgegraphs.reverse()3848vnumbers.reverse()3849contrdicts.reverse()385038513852# Step 2: pull back the decHstratum t one edge-degeneration at a time3853# Note that this will convert the decHstratum into a prodHclass, so we need to take care of all the different terms in each step38543855# First we need to adapt (a copy of) self so that we take into account that gamma0 is isomorphic to the contracted version of gamma13856# The gamma0-structure on all terms of self must be modified to be a contrgraph[0]-structure3857result=deepcopy(self)3858result.gamma0=contrgraph[0]3859dicvisom={dicv[vimages[v]]:v for v in range(gamma0.num_verts())} # gives isomorphism V(gamma0) -> V(contrgraph[0])3860diclisom={dicl[l]:l for l in dicl} # gives isomorphism L(contrgraph[0]) -> L(gamma0)3861for t in result.terms:3862t.dicv0={v:dicvisom[t.dicv0[v]] for v in t.dicv0}3863t.dicl0={l:t.dicl0[diclisom[l]] for l in diclisom}3864# dicv, dicl should have served their purpose and should no longer appear below38653866for edgenumber in range(delta_e):3867newresult=prodHclass(contrgraph[edgenumber+1],[])3868gam0=contrgraph[edgenumber] # current gamma0, we now pull back under contrgraph[edgenumber+1] -> gam03869av=actvert[edgenumber] # vertex in gam0 where action happens3870diccv=contrdicts[edgenumber] # vertex dictionary for contrgraph[edgenumber+1] -> gam03871diccvinv={diccv[v]:v for v in diccv} # section of diccv: bijective if loop is added, otherwise assigns one of the preimages to the active vertex3872vextractdic={v:vnumbers[edgenumber][v] for v in range(len(vnumbers[edgenumber]))}387338743875for t in result.terms:3876gamma=t.gamma38773878# Step 2.1: extract the preimage of the active vertex inside t.gamma3879actvertices=[v for v in range(gamma.num_verts()) if t.dicv0[v]==av]3880#global exsubvar3881#exsubvar=(t,gamma,actvertices,gam0, av)3882#avloops=[e for e in gam0.edges if (e[0] in gam0.legs(av, copy=False) and e[1] in gam0.legs(av, copy=False))]3883#outlegs=set(gam0.legs(av, copy=False))-set([e[0] for e in avloops] + [e[1] for e in avloops])3884(gammaprime,dicvextract,diclextract) =gamma.extract_subgraph(actvertices,outgoing_legs=[t.dicl0[l] for l in gam0.legs(av, copy=False)], rename=False)38853886dicvextractinv={dicvextract[v]:v for v in dicvextract}3887# dicvextract maps vertex numbers in gamma to vertex-numbers in gammaprime3888# diclextract not needed since rename=False38893890# Step 2.2: find the common degenerations of this preimage and the one-edge graph relevant for this step3891# We need to rename edgegraphs[edgenumber] to match the leg-names in gammaprime3892egraph = edgegraphs[edgenumber].copy()3893eshift=max(list(t.dicl0.values())+[0])+13894egraph.rename_legs(t.dicl0,shift=eshift)3895egraph.set_immutable()38963897#(gammaprime,dicvextract,diclextract) =gamma.extract_subgraph(actvertices,outgoing_legs=[l for l in egraph.list_markings()],rename=False)38983899#dicvextractinv={dicvextract[v]:v for v in dicvextract}3900# dicvextract maps vertex numbers in gamma to vertex-numbers in gammaprime3901# diclextract not needed since rename=False39023903commdeg=common_degenerations(gammaprime,egraph,modiso=True,rename=True)390439053906(e0,e1) = egraph.edges(copy=False)[0]3907for (gammadeg,dv1,dl1,dv2,dl2) in commdeg:3908if gammadeg.num_edges() == gammaprime.num_edges():3909# Step 2.3a: if gammadeg isomorphic to gammaprime, we have excess intersection, introducing -psi - psi'3910term=deepcopy(t)3911numvert = gamma.num_verts()3912dl1inverse={dl1[l]:l for l in dl1}3913term.poly*=-(psicl(dl1inverse[dl2[e0]],numvert)+psicl(dl1inverse[dl2[e1]],numvert))3914# adapt term.dicv0 and term.dicl039153916# term.dicv0 must be lifted along the contraction map using diccvinv outside of the active area gammaprime; there we use dv23917dv1inverse={dv1[v]:v for v in dv1} #TODO: is this what we want?3918term.dicv0={v:diccvinv[term.dicv0[v]] for v in term.dicv0}3919term.dicv0.update({v:vextractdic[dv2[dv1inverse[dicvextract[v]]]] for v in actvertices})39203921# term.dicl0 should now assign to the two half-edges that are new the corresponding edge in term.gamma3922term.dicl0[e0-eshift]= dl1inverse[dl2[e0]]3923term.dicl0[e1-eshift]= dl1inverse[dl2[e1]]39243925newresult.terms.append(term)3926else:3927# Step 2.3b: else there is a real degeneration going on3928# This means we need to find the vertex of gamma where this new edge would appear,3929# then check if the Hurwitz space glued to this spot admits a compatible degeneration.3930# Note that here we must take note that we possibly forgot some markings of this Hurwitz space.393139323933vactiveprime=dv1[gammadeg.vertex(dl2[e0])] # vertex in gammaprime where degeneration occurs3934vactive=dicvextractinv[vactiveprime] # vertex in gamma where degeneration occurs393539363937#TODO: it is likely much faster to go through elements Gr of Hbdry, go through edges e of Gr, contract everything except e and check if3938#compatible with the boundary partition3939(a,spacedicl)=t.vertdata[vactive]3940(gsp,Hsp)=t.spaces[a]3941numHmarks = trivGgraph(gsp,Hsp).gamma.num_legs(0)39423943#Hbdry=list_Hstrata(gsp,Hsp,1) # these are the Gstgraphs which might get glued in at vertex vactive (+ other vertices accessing space a)3944#if Hbdry==[]:3945# continue # we would require a degeneration, but none can occur3946#numHmarks=len(Hbdry[0].gamma.list_markings())3947#presentmarks=spacedicl.values()3948#forgetmarks=[i for i in range(1,numHmarks) if i not in presentmarks]39493950# now create list divlist of divisors D living in the space of Hbdry such that we can start enumerating D-structures on elements of Hbdry3951# first we must extract the original boundary graph, which is a subgraph of gammadeg on the preimages of vactiveprime3952vactiveprimepreim=[gammadeg.vertex(dl2[e0]),gammadeg.vertex(dl2[e1])]3953if vactiveprimepreim[0]==vactiveprimepreim[1]:3954vactiveprimepreim.pop(1)3955# old, possibly faulty: (egr,dvex,dlex) =gammadeg.extract_subgraph(vactiveprimepreim,outgoing_legs=gammaprime.legs(vactiveprime),rename=False)3956(egr,dvex,dlex) = gammadeg.extract_subgraph(vactiveprimepreim,outgoing_legs=[dl1[le] for le in gammaprime.legs(vactiveprime, copy=False)], rename=False, mutable=True)39573958# rename egr to fit in the space of Hbdry3959rndic = {dl1[l]: spacedicl[l] for l in gammaprime.legs(vactiveprime, copy=False)}3960egr.rename_legs(rndic,shift=numHmarks+1)3961egr.set_immutable()39623963degenlist = Hbdrystructures(gsp,Hsp,egr)39643965for (Gr,mult,degendicv,degendicl) in degenlist:3966term=deepcopy(t)3967e = egr.edges(copy=False)[0]3968speciale_im = term.replace_space_by_Gstgraph(a,Gr,specialv=vactive,speciale=(degendicl[e[0]], degendicl[e[1]]))3969term.poly*=mult3970# term is still a decHstratum with correct gamma, spaces, vertdata and poly, BUT still referencing to ga03971# now repair dicv0, dicl039723973# term.dicl0 should now assign to the two half-edges that are new the corresponding edge in term.gamma3974dl1inverse={dl1[l]:l for l in dl1}3975eindex = e.index(dl2[e0]+numHmarks+1)3976term.dicl0[e0-eshift]= speciale_im[eindex]3977term.dicl0[e1-eshift]= speciale_im[1 -eindex]39783979# term.dicv0 is now reconstructed from dicl03980# TODO: add more bookkeeping to do this directly???3981term.dicv0=dicv_reconstruct(term.gamma,contrgraph[edgenumber+1],term.dicl0)39823983newresult.terms.append(term)39843985result=newresult398639873988return result39893990# def dimension_filter(self):3991# for t in self.terms:3992# for s in t:3993# s.dimension_filter()3994# return self.consolidate()39953996def consolidate(self):3997# revrange=range(len(self.terms))3998# revrange.reverse()3999# for i in revrange:4000# for s in self.terms[i]:4001# if len(s.poly.monom)==0:4002# self.terms.pop(i)4003# break4004return self40054006def __repr__(self):4007s='Outer graph : ' + repr(self.gamma0) +'\n'4008for t in self.terms:4009s+=repr(t) + '\n\n'4010# for i in range(len(self.terms)):4011# for j in range(len(self.terms[i])):4012# s+='Vertex '+repr(j)+' :\n'+repr(self.terms[i][j])+'\n'4013# s+='\n\n'4014return s.rstrip('\n\n')40154016# takes a genus g, a HurwitzData H and a stgraph bdry with exactly one edge, whose markings are a subset of the markings 1,2,..., N of trivGgraph(g,H)4017# returns a list with entries (Gr,mult,dicv,dicl) giving all boundary Gstgraphs derived from (g,H) which are specializations of (the forgetful pullback of) Gr, where4018# * Gr is a Gstgraph4019# * mult is a rational number giving an appropriate multiplicity for the intersection4020# * dicv is a surjective dictionary of the vertices of Gr.gamma to the vertices of bdry4021# * dicl is an injection from the legs of bdry to the legs of Gr.gamma4022def Hbdrystructures(g,H,bdry):4023preHbdry, N = preHbdrystructures(g, H)40244025if bdry.num_verts() == 1:4026try:4027tempresult=preHbdry[(g)]4028except:4029return []4030(e0,e1) = bdry.edges(copy=False)[0]4031result=[]4032for (Gr,mult,dicv,dicl) in tempresult:4033result.append([Gr,mult,dicv,{e0:dicl[N+1], e1:dicl[N+2]}])4034return result403540364037if bdry.num_verts() == 2:4038presentmarks=bdry.list_markings()4039forgetmarks=[i for i in range(1, N+1) if i not in presentmarks]4040possi=[[0,1] for j in forgetmarks]404140424043(e0,e1) = bdry.edges(copy=False)[0]4044lgs = bdry.legs(copy=True)4045ve0 = bdry.vertex(e0)4046lgs[ve0].remove(e0) #leglists stripped of the legs forming the edge of bdry4047lgs[1 -ve0].remove(e1)40484049result=[]40504051for distri in itertools.product(*possi):4052# distri = (0,0,1,0,1,1) means that forgetmarks number 0,1,3 go to erg.legs(0) and 2,4,5 go to erg.legs(1)4053lgscpy=deepcopy(lgs)4054lgscpy[0]+=[forgetmarks[j] for j in range(len(forgetmarks)) if distri[j]==0]4055lgscpy[1]+=[forgetmarks[j] for j in range(len(forgetmarks)) if distri[j]==1]40564057if bdry.genera(0) > bdry.genera(1) or (bdry.genera(0) == bdry.genera(1) and 1 in lgscpy[1]):4058vert=14059else:4060vert=04061ed=(ve0+vert)%2406240634064try:4065tempresult = preHbdry[(bdry.genera(vert), tuple(sorted(lgscpy[vert])))]4066for (Gr,mult,dicv,dicl) in tempresult:4067result.append((Gr,mult,{j:(dicv[j]+vert)%2 for j in dicv},{e0:dicl[N+1 +ed], e1:dicl[N+2 -ed]}))4068except KeyError:4069pass4070return result40714072# returns (result,N), where N is the number of markings in Gstgraphs with (g,H) and result is a dictionary sending some tuples (g',(m_1, ..., m_l)), which correspond to boundary divisors, to the list of all entries of the form [Gr,mult,dicv,dicl] like in the definition of Hbdrystructures. Here dicv and dicl assume that the boundary divisor bdry is of the form stgraph([g',g-g'],[[m_1, ..., m_l,N+1],[complementary legs, N+2],[(N+1,N+2)]]4073# Convention: we choose the key above such that g'<=g-g'. If there is equality, either there are no markings OR we ask that m_1=1. Note that m_1, ... are assumed to be ordered4074# For the unique boundary divisor parametrizing irreducible curves (having a self-node), we reserve the key (g)40754076@cached_function4077def preHbdrystructures(g, H):4078Hbdry=list_Hstrata(g,H,1)4079if len(Hbdry)==0:4080return ({},len(trivGgraph(g,H).gamma.list_markings()))4081N=len(Hbdry[0].gamma.list_markings())4082result={}4083for Gr in Hbdry:4084AutGr=equiGraphIsom(Gr,Gr)40854086edgs = Gr.gamma.edges()40874088while edgs:4089econtr=edgs[0] # this edge will remain after contraction40904091# now go through AutGr and eliminate all other edges which are in orbit of econtr, record multiplicity4092multiplicity=04093for (dicv,dicl) in AutGr:4094try:4095edgs.remove((dicl[econtr[0]],dicl[econtr[1]]))4096multiplicity+=14097except:4098pass409941004101# now contract all except econtr and add the necessary data to the dictionary result4102contrGr = Gr.gamma.copy(mutable=True)4103diccv = {j:j for j in range(contrGr.num_verts())}4104for e in Gr.gamma.edges(copy=False):4105if e != econtr:4106(v0, d1, d2,diccv_new) = contrGr.contract_edge(e,True)4107diccv = {j:diccv_new[diccv[j]] for j in diccv}41084109# OLD: multiplicity*=len(GraphIsom(contrGr,contrGr))4110# NOTE: multiplicity must be multiplied by the number of automorphisms of the boundary stratum, 2 in this case4111# this compensates for the fact that these automorphisms give rise to different bdry-structures on Gr41124113if contrGr.num_verts() == 1:4114# it remains a self-loop4115if g not in result:4116result[(g)]=[]4117result[(g)].append([Gr,QQ(multiplicity)/len(AutGr),{j:0 for j in range(Gr.gamma.num_verts())},{N+1: econtr[0], N+2: econtr[1]}])4118result[(g)].append([Gr,QQ(multiplicity)/len(AutGr),{j:0 for j in range(Gr.gamma.num_verts())},{N+2: econtr[0], N+1: econtr[1]}])41194120else:4121# we obtain a separating boundary; now we must distinguish cases to ensure that g' <= g-g' and for equality m_1=14122if contrGr.genera(0) > contrGr.genera(1) or (contrGr.genera(0) == contrGr.genera(1) and 1 in contrGr.legs(1, copy=False)):4123switched=14124else:4125switched=041264127eswitched=(contrGr.vertex(econtr[0])+switched)%2 # eswitched=0 means econtr[0] is in the vertex with g',[m_1,..]; eswitched=1 means it is at the other vertex41284129gprime = contrGr.genera(switched)4130markns = contrGr.legs(switched, copy=True)4131markns.remove(econtr[eswitched])4132markns.sort()4133markns=tuple(markns)413441354136try:4137result[(gprime,markns)].append([Gr,QQ(multiplicity)/len(AutGr),{v:(diccv[v]+switched)%2 for v in diccv},{N+1: econtr[eswitched], N+2: econtr[1 -eswitched]}])4138except KeyError:4139result[(gprime,markns)]=[[Gr,QQ(multiplicity)/len(AutGr),{v:(diccv[v]+switched)%2 for v in diccv},{N+1: econtr[eswitched], N+2: econtr[1 -eswitched]}]]41404141if contrGr.genera(0) == contrGr.genera(1) and N == 0:4142# in this case the resulting boundary divisor has an automorphism, so we need to record also the switched version of the above bdry-structure4143result[(gprime,markns)].append([Gr,QQ(multiplicity)/len(AutGr),{v:(diccv[v]+1 -switched)%2 for v in diccv},{N+2: econtr[eswitched], N+1: econtr[1 -eswitched]}])41444145return (result,N)414641474148414941504151415241534154# summands in prodHclass4155# * gamma is a stgraph on which this summand lives4156# * spaces is a dictionary, assigning entries of the form [g,Hdata] of a genus and HurwitzData Hdata to integer labels a - since these change a lot and are replaced frequently by multiple new labels, this data structure hopefully allows for flexibility4157# * vertdata gives a list of entries of the form [a, dicl] for each vertex v of gamma, where a is a label (sent to (g,Hdata) by spaces) and dicl is an injective dictionary mapping leg-names of v to the corresponding (standard) numbers of the Hurwitz space of (g,Hdata)4158# * poly is a kppoly on gamma4159# * dicv0 is a surjective dictionary mapping vertices of gamma to vertices of the including prodHclass4160# * dicl0 is an injective dictionary mapping legs of gamma0 to legs of gamma4161class decHstratum(object):4162def __init__(self,gamma,spaces,vertdata,dicv0=None,dicl0=None,poly=None):4163self.gamma=gamma4164self.spaces=spaces4165self.vertdata=vertdata4166if dicv0 is None:4167self.dicv0={v:0 for v in range(gamma.num_verts())}4168else:4169self.dicv0=dicv041704171if dicl0 is None:4172self.dicl0={l:l for l in gamma.list_markings()}4173else:4174self.dicl0=dicl041754176if poly is None:4177self.poly = onekppoly(gamma.num_verts())4178else:4179self.poly=poly41804181def __repr__(self):4182return repr((self.gamma,self.spaces,self.vertdata,self.dicv0,self.dicl0,self.poly))41834184def __neg__(self):4185return (-1)*self41864187def __rmul__(self,other):4188#if isinstance(other,sage.rings.integer.Integer) or isinstance(other,sage.rings.rational.Rational) or isinstance(other,int):4189new=deepcopy(self)4190new.poly*=other4191return new41924193# self represents a morphism H_1 x ... x H_r -i-> M_1 x ... x M_r -pf-> M_1'x ... x M_g', where the first map i is the inclusion of Hurwitz stacks and the second map pf is a product of forgetful maps of some of the M_i. The spaces M_j' are those belonging to the vertices of self.gamma4194# prodforgetpullback takes a dictionary spacelist mapping space labels a_1, ..., a_r to the lists of vertices j using these spaces (i.e. the M_j'-component of pf depends on the M_a_i argument if j in spacelist[a_i])4195# it returns the pullback of self.poly under pf, which is a prodtautclass on a stgraph being the disjoint union of M_1, ..., M_r (in the order of spacelist.keys())4196# this function is meant as an auxilliary function for self.evaluate()4197def prodforgetpullback(self,spacelist):4198alist=list(spacelist)4199decs=decstratum(self.gamma,poly=self.poly)4200splitdecs=decs.split()42014202# find number of markings of spaces M_i4203N={a: self.spaces[a][1].nummarks() for a in spacelist}4204forgetlegs = [list(set(range(1, N[self.vertdata[v][0]]+1))-set(self.vertdata[v][1].values())) for v in range(self.gamma.num_verts())]4205resultgraph=stgraph([self.spaces[a][0] for a in spacelist],[list(range(1, N[a]+1)) for a in N],[])42064207# rename splitted kppolys to match the standard names of markings on the M_i and wrap them in decstrata4208trivgraphs = [stgraph([self.gamma.genera(i)],[list(self.vertdata[i][1].values())],[]) for i in range(self.gamma.num_verts())]4209splitdecs=[[s[i].rename_legs(self.vertdata[i][1]) for i in range(len(s))] for s in splitdecs]4210splitdecs=[[decstratum(trivgraphs[i],poly=s[i]) for i in range(len(s))] for s in splitdecs]42114212result=prodtautclass(resultgraph,[])4213for s in splitdecs:4214term=prodtautclass(resultgraph)4215for j in range(len(alist)):4216t=prod([s[i].forgetful_pullback(forgetlegs[i]) for i in spacelist[alist[j]]])4217t.dimension_filter()4218t=t.toprodtautclass(self.spaces[alist[j]][0],N[alist[j]])4219term=term.factor_pullback([j],t)4220result+=term4221return result422242234224# evaluates self against the fundamental class of the ambient space (of self.gamma), returns a rational number4225def evaluate(self):4226spacelist={a:[] for a in self.spaces}4227for v in range(self.gamma.num_verts()):4228spacelist[self.vertdata[v][0]].append(v)4229for a in spacelist:4230if spacelist[a]==[]:4231spacelist.pop(a)4232alist=list(spacelist)42334234pfp=self.prodforgetpullback(spacelist)4235Gr=[trivGgraph(*self.spaces[a]) for a in spacelist]4236prodGr=[prodHclass(gr.gamma,[gr.to_decHstratum()]) for gr in Gr]4237result=04238for t in pfp.terms:4239tempresult=14240for i in range(len(alist)):4241if t[i].gamma.edges()==[]: # decstratum t[i] is a pure kppoly4242tempresult*=(Hdecstratum(Gr[i],poly=t[i].poly)).quotient_pushforward().evaluate()4243else:4244tempresult*=(prodGr[i]*t[i]).evaluate()4245result+=tempresult4246return result42474248# replaces the fundamental class of the Hurwitz space with label a with the (fundamental class of the stratum corresponding to) Gstgraph Gr4249# in particular, self.gamma is changed by replacing all vertices with vertdata a with the graph Gr (more precisely: with the graph obtained by forgetting the corresponding markings in Gr and then stabilizing4250# if the vertex specialv in self and the edge speciale in Gr are given, also return the edge corresponding to the version of speciale that was glued to the vertex specialv (which by assumption has space labelled by a)4251# IMPORTANT: here we assume that speciale is not contracted in the forgetful-stabilization process!4252def replace_space_by_Gstgraph(self,a,Gr,specialv=None,speciale=None):4253#global examinrep4254#examinrep=(deepcopy(self),a,deepcopy(Gr),specialv,speciale)42554256avertices = [v for v in range(self.gamma.num_verts()) if self.vertdata[v][0] == a]4257avertices.reverse()425842594260gluein=Gr.to_decHstratum()426142624263maxspacelabel=max(self.spaces)4264self.spaces.pop(a)4265self.spaces.update({l+maxspacelabel+1: gluein.spaces[l] for l in gluein.spaces})42664267for v in avertices:4268# First we must create the graph to be glued in at v. For this we take the graph from gluein and forget unnecessary markings, stabilize and then rename the remaining markings as dictated by self.vertdata[v][1]. During all steps we must record which vertices/legs go where42694270gluegraph = gluein.gamma.copy()4271dicl_v=self.vertdata[v][1]42724273usedlegs=dicl_v.values()4274unusedlegs=[l for l in gluegraph.list_markings() if l not in usedlegs]42754276gluegraph.forget_markings(unusedlegs)427742784279(forgetdicv,forgetdicl,forgetdich)=gluegraph.stabilize()42804281forgetdicl.update({l:l for l in gluegraph.leglist() if l not in forgetdicl})428242834284# now we must rename the legs in gluegraph so they fit with the leg-names in self.gamma, using the inverse of dicl_v4285dicl_v_inverse={dicl_v[l]:l for l in dicl_v}4286shift=max(list(dicl_v)+[0])+14287gluegraph.rename_legs(dicl_v_inverse,shift)428842894290divGr={}4291divs={}4292dil={}4293numvert_old = self.gamma.num_verts()4294num_newvert = gluegraph.num_verts()4295dic_voutgoing = {l:l for l in self.gamma.legs(v, copy=False)}4296self.gamma = self.gamma.copy()4297self.gamma.glue_vertex(v,gluegraph,divGr,divs,dil)4298self.gamma.set_immutable()4299dil.update(dic_voutgoing)4300dil_inverse={dil[l]:l for l in dil}430143024303if specialv==v:4304speciale_im=(dil[forgetdich.get(speciale[0],speciale[0])+shift],dil[forgetdich.get(speciale[1],speciale[1])+shift])430543064307# now repair vertdata: note that new vertex list is obtained by deleting v and appending vertex list of gluegraph4308# also note that in the gluing-procedure ALL leg-labels at vertices different from v have NOT changed, so those vertdatas do not have to be adapted4309self.vertdata.pop(v)4310for w in range(num_newvert):4311# forgetdicv[w] is the vertex number that w corresponded to before stabilizing, inside gluein4312(pre_b,pre_diclb)=gluein.vertdata[forgetdicv[w]]431343144315b=pre_b+maxspacelabel+14316# pre_diclb is injective dictionary from legs of vertex forgetdicv[w] in gluein.gamma to the standard-leg-numbers in self.spaces[b]4317# thus we need to reverse the shift/rename and the glue_vertex-rename and compose with diclb4318diclb={l:pre_diclb[forgetdicl[dicl_v.get(dil_inverse[l],dil_inverse[l]-shift)]] for l in self.gamma.legs(w+numvert_old-1, copy=False)}43194320self.vertdata.append([b,diclb])432143224323# now repair dicv0, the surjective map from vertices of self.gamma to the gamma0 of the containing prodHclass4324# vertex v went to some vertex w, now all vertices that have replaced v must go to w4325# unfortunately, all vertices v+1, v+2, ... have shifted down by 1, so this must be taken care of4326w=self.dicv0.pop(v)4327for vnumb in range(v, numvert_old-1):4328self.dicv0[vnumb]=self.dicv0[vnumb+1]4329for vnumb in range(numvert_old-1, numvert_old - 1 + num_newvert):4330self.dicv0[vnumb]=w43314332# now repair dicl0, the injective map from leg-names of gamma0 to leg-names of self.gamma4333# ACTUALLY: since gluing does not change any of the old leg-names, dicl0 is already up to date!433443354336# now repair self.poly, term by term, collecting the results in newpoly, which replaces self.poly at the end4337newpoly=kppoly([],[])4338for (kappa,psi,coeff) in self.poly:4339trunckappa=deepcopy(kappa)4340kappav=trunckappa.pop(v)4341trunckappa+=[[] for i in range(num_newvert)]4342trunckppoly=kppoly([(trunckappa,psi)],[coeff])4343trunckppoly*=prod([(sum([kappacl(numvert_old-1+k, j+1, numvert_old+num_newvert-1) for k in range(num_newvert)]))**kappav[j] for j in range(len(kappav))])4344newpoly+=trunckppoly43454346self.poly=newpoly43474348spacesused = {a: False for a in self.spaces}4349for a, _ in self.vertdata:4350spacesused[a] = True4351unusedspaces = [a for a in spacesused if not spacesused[a]]43524353degree=14354for a in unusedspaces:4355trivGraph=trivGgraph(*self.spaces.pop(a))4356if trivGraph.dim()==0:4357degree*=trivGraph.delta_degree(0)4358else:4359degree=04360break43614362self.poly*=degree43634364if specialv is not None:4365return speciale_im43664367# takes genus g, number n of markings and lists leglists, legdics of length r (and possibly a tautclass T on \bar M_{g,n})4368# for i=0,..,r-1 the dictionary legdics[i] is an injection of the list leglists[i] in {1,..,n}4369# computes the class DELTA of the small diagonal in (\bar M_{g,n})^r, intersects with T on the first factor and then pushes down under a forgetful map4370# this map exactly forgets in the ith factor all markings not in the image of legdics[i]4371# afterwards it renames the markings on the ith factor such that marking legdics[i][j] is called j4372# it returns a prodtautclass on a graph that is the disjoint union of genus g vertices with leglists given by leglists (assumed to be disjoint)4373def forgetful_diagonal(g,n,leglists,legdics,T=None):4374#global inpu4375#inpu = (g,n,leglists,legdics,T)4376r=len(leglists)4377if T is None:4378T=tautclass([decstratum(stgraph([g],[list(range(1, n+1))],[]))])43794380gamma=stgraph([g for i in range(r)],[list(range(1, n+1)) for i in range(r)],[])4381result=prodtautclass(gamma,[[decstratum(stgraph([g],[list(range(1, n+1))],[])) for i in range(r)]])43824383# We want to obtain small diagonal by intersecting Delta_12, Delta_13, ..., Delta_1r4384# since we will forget markings in the factors 2,..,r, only components of Delta_1i are relevant which have sufficient cohomological degree in the second factor (all others vanish in the pushforward having positive-dimensional fibres)4385forgetlegs=[[j for j in range(1, n+1) if j not in legdics[i].values()] for i in range(r)]4386fl=[len(forg) for forg in forgetlegs]4387rndics=[{legdics[i][leglists[i][j]]: j+1 for j in range(len(leglists[i]))} for i in range(r)]438843894390if r >1:4391#TODO: check that everything is tautological more carefully, only checking those degrees that are really crucial4392for rdeg in range(0, 2 *(3 *g-3 +n)+1):4393if not cohom_is_taut(g,n,rdeg):4394print("In computation of diagonal : H^%s(bar M_%s,%n) not necessarily tautological" % (rdeg,g,n))43954396# now compute diagonal4397Dgamma=stgraph([g,g],[list(range(1, n+1)),list(range(1, n+1))],[])4398Delta=prodtautclass(Dgamma,[])4399for deg in range(0, (3 *g-3 +n)+1):4400gi1=generating_indices(g,n,deg)4401gi2=generating_indices(g,n,3 *g-3 +n-deg)4402strata1 = DR.all_strata(g,deg,tuple(range(1, n+1)))4403gens1=[Graphtodecstratum(strata1[j]) for j in gi1]4404strata2 = DR.all_strata(g,3 *g-3 +n-deg,tuple(range(1, n+1)))4405gens2=[Graphtodecstratum(strata2[j]) for j in gi2]44064407# compute inverse of intersection matrix4408invintmat=inverseintmat(g,n,tuple(gi1),tuple(gi2),deg)4409for a in range(len(gi1)):4410for b in range(len(gi2)):4411if invintmat[a,b] != 0:4412Delta.terms.append([invintmat[a,b]*gens1[a],gens2[b]])44134414for i in range(1, r):4415result=result.factor_pullback([0,i],Delta)4416result=result.factor_pullback([0],T.toprodtautclass(g,n))44174418# now we take care of forgetful pushforwards and renaming4419for t in result.terms:4420for i in range(r):4421t[i]=t[i].forgetful_pushforward(forgetlegs[i])4422t[i].rename_legs(rndics[i]) #note: there can be no half-edge names in the range 1,..,n so simply renaming legs is possible44234424result.gamma=stgraph([g for i in range(r)],deepcopy(leglists),[])44254426return result.consolidate()442744284429@cached_function4430def inverseintmat(g,n,gi1,gi2,deg):4431if deg <= QQ(3 *g-3 +n)/2:4432intmat=matrix(DR.pairing_submatrix(gi1,gi2,g,deg,tuple(range(1, n+1))))4433else:4434intmat=matrix(DR.pairing_submatrix(gi2,gi1,g,3 *g-3 +n-deg,tuple(range(1, n+1)))).transpose()4435return intmat.transpose().inverse()44364437def inverseintmat2(g,n,gi1,gi2,deg):4438tg1=tautgens(g,n,deg)4439tg2=tautgens(g,n,3 *g-3 +n-deg)4440intmat=matrix([[(tg1[i]*tg2[j]).evaluate() for i in gi1] for j in gi2])4441return intmat.transpose().inverse()44424443""" # remove all terms that must vanish for dimension reasons4444def dimension_filter(self):4445for i in range(len(self.poly.monom)):4446(kappa,psi)=self.poly.monom[i]4447for v in range(len(self.gamma.genera)):4448# local check: we are not exceeding dimension at any of the vertices4449if sum([(k+1)*kappa[v][k] for k in range(len(kappa[v]))])+sum([psi[l] for l in self.gamma.legs(v) if l in psi])>3*self.gamma.genera(v)-3+len(self.gamma.legs(v)):4450self.poly.coeff[i]=04451break4452# global check: the total codimension of the term of the decstratum is not higher than the total dimension of the ambient space4453#if self.poly.deg(i)+len(self.gamma.edges)>3*self.gamma.g()-3+len(self.gamma.list_markings()):4454# self.poly.coeff[i]=04455#break4456return self.consolidate()44574458def consolidate(self):4459self.poly.consolidate()4460return self """44614462""" # computes integral against the fundamental class of the corresponding moduli space4463# will not complain if terms are mixed degree or if some of them do not have the right codimension4464def evaluate(self):4465answer=04466for (kappa,psi,coeff) in self.poly:4467temp=14468for v in range(len(self.gamma.genera)):4469psilist=[psi.get(l,0) for l in self.gamma.legs(v)]4470kappalist=[]4471for j in range(len(kappa[v])):4472kappalist+=[j+1 for k in range(kappa[v][j])]4473if sum(psilist+kappalist) != 3*self.gamma.genera(v)-3+len(psilist):4474temp = 04475break4476temp*=DR.socle_formula(self.gamma.genera(v),psilist,kappalist)4477answer+=coeff*temp4478return answer """4479""" def rename_legs(self,dic):4480self.gamma.rename_legs(dic)4481self.poly.rename_legs(dic)4482return self4483def __repr__(self):4484return 'Graph : ' + repr(self.gamma) +'\n'+ 'Polynomial : ' + repr(self.poly) """448544864487# old code snippet for boundary pullback of decHstrata4488"""4489if len(erg.genera)==1:4490# loop graph, add all missing markings to it4491egr.legs(0)+=forgetmarks4492divlist=[egr]4493if len(erg.genera)==2:4494possi=[[0,1] for j in forgetmarks]4495divlist=[]4496for distri in itertools.product(*possi):4497# distri = (0,0,1,0,1,1) means that forgetmarks number 0,1,3 go to erg.legs(0) and 2,4,5 go to erg.legs(1)4498ergcpy=deepcopy(erg)4499ergcpy.legs(0)+=[forgetmarks[j] for j in range(len(forgetmarks)) if distri[j]==0]4500ergcpy.legs(1)+=[forgetmarks[j] for j in range(len(forgetmarks)) if distri[j]==1]4501divlist.append(ergcpy)4502"""45034504# Knowing that there exists a morphism Gamma -> A, we want to reconstruct the corresponding dicv (from vertices of Gamma to vertices of A) from the known dicl4505#TODO: we could give dicv as additional argument if we already know part of the final dicv4506def dicv_reconstruct(Gamma, A, dicl):4507if A.num_verts() == 1:4508return {ver:0 for ver in range(Gamma.num_verts())}4509else:4510dicv={Gamma.vertex(dicl[h]): A.vertex(h) for h in dicl}4511remaining_vertices=set(range(Gamma.num_verts()))-set(dicv)4512remaining_edges=set([e for e in Gamma.edges(copy=False) if e[0] not in dicl.values()])4513current_vertices=set(dicv)45144515while remaining_vertices:4516newcurrent_vertices=set()4517for e in deepcopy(remaining_edges):4518if Gamma.vertex(e[0]) in current_vertices:4519vnew=Gamma.vertex(e[1])4520remaining_edges.remove(e)4521# if vnew is in current vertices, we don't have to do anything (already know it)4522# otherwise it must be a new vertex (cannot reach old vertex past the front of current vertices)4523if vnew not in current_vertices:4524dicv[vnew]=dicv[Gamma.vertex(e[0])]4525remaining_vertices.discard(vnew)4526newcurrent_vertices.add(vnew)4527continue4528if Gamma.vertex(e[1]) in current_vertices:4529vnew=Gamma.vertex(e[0])4530remaining_edges.remove(e)4531# if vnew is in current vertices, we don't have to do anything (already know it)4532# otherwise it must be a new vertex (cannot reach old vertex past the front of current vertices)4533if vnew not in current_vertices:4534dicv[vnew]=dicv[Gamma.vertex(e[1])]4535remaining_vertices.discard(vnew)4536newcurrent_vertices.add(vnew)4537current_vertices=newcurrent_vertices4538return dicv453945404541# approximation to the question if H^{r}(\overline M_{g,n},QQ) is generated by tautological classes4542# if True is returned, the answer is yes, if False is returned, the answer is no; if None is returned, we make no claim4543# we give a reference for each of the claims, though we do not claim that this is the historically first such reference4544def cohom_is_taut(g, n, r):4545if g == 0:4546return True # [Keel - Intersection theory of moduli space of stable N-pointed curves of genus zero]4547if g == 1 and r % 2 == 0:4548return True # [Petersen - The structure of the tautological ring in genus one]4549if g == 1 and n < 10:4550return True # [Graber, Pandharipande - Constructions of nontautological classes on moduli spaces of curves] TODO: also for n=10 probably, typo in GP ?4551if g == 1 and n >= 11 and r == 11:4552return False # [Graber, Pandharipande]4553if g == 2 and r % 2 == 0 and n < 20:4554return True # [Petersen - Tautological rings of spaces of pointed genus two curves of compact type]4555if g == 2 and n <= 3:4556return True # [Getzler - Topological recursion relations in genus 2 (above Prop. 16+e_11=-1,e_2=0)] + [Yang - Calculating intersection numbers on moduli spaces of curves]4557if g == 2 and n == 4:4558return True # [Bini,Gaiffi,Polito - A formula for the Euler characteristic of \bar M_{2,4}] + [Yang] + [Petersen - Tautological ...]: all even cohomology is tautological and dimensions from Yang sum up to Euler characteristic4559# WARNING: [BGP] contains error, later found in [Bini,Harer - Euler characteristics of moduli spaces of curves], does not affect the numbers above4560if g == 3 and n == 0:4561return True # [Getzler - Topological recursion relations in genus 2 (Proof of Prop. 16)] + [Yang]4562if g == 3 and n == 1:4563return True # [Getzler,Looijenga - The hodge polynomial of \bar M_3,1] + [Yang]4564if g == 3 and n == 2:4565return True # [Bergström - Cohomology of moduli spaces of curves of genus three via point counts]4566if g == 4 and n == 0:4567return True # [Bergström, Tommasi - The rational cohomology of \bar M_4] + [Yang]45684569#TODO: Getzler computes Serre characteristic in genus 1 ('96)4570#TODO: [Bini,Harer - Euler characteristics of moduli spaces of curves]: recursive computation of chi(\bar M_{g,n})4571# This suggests that g=3,n=2 should also work (Yang's numbers add up to Euler char.)45724573if r in [0, 2]:4574return True # r=0: fundamental class is tautological, r=1: [Arbarello, Cornalba - The Picard groups of the moduli spaces of curves]4575if r in [1, 3, 5]:4576return True # [Arbarello, Cornalba - Calculating cohomology groups of moduli spaces of curves via algebraic geometry]4577if 2 * (3 * g - 3 + n) - r in [0, 1, 2, 3, 5]:4578return True4579# dim=(3*g-3+n), then by Hard Lefschetz, there is isomorphism H^(dim-k) -> H^(dim+k) by multiplication with ample divisor4580# this divisor is tautological, so if all of H^(dim-k) is tautological, also H^(dim+k) is tautological45814582return None45834584# approximation to the question if the Faber-Zagier relations generate all tautological relations in A^d(\bar M_{g,n})4585# if True is returned, the answer is yes, if False is returned, it means that no proof has been registered yet4586# we give a reference for each of the claims, though we do not claim that this is the historically first such reference4587def FZ_conjecture_holds(g, n, d):4588if g == 0:4589return True # [Keel - Intersection theory of moduli space of stable N-pointed curves of genus zero]4590if g == 1:4591return True # [Petersen - The structure of the tautological ring in genus one] %TODO: verify FZ => Getzler's rel. + WDVV45924593############4594#4595# Lookup and saving functions for Hdatabase4596#4597############45984599# To ensure future compatibility: only append functions to this list! Do not change order of existing functions.4600PICKLED_FUNCTIONS = [generating_indices, genstobasis]460146024603def save_FZrels():4604"""4605Saving previously computed Faber-Zagier relations to file new_geninddb.pkl.46064607TESTS::46084609sage: from admcycles import *4610sage: generating_indices(0,3,0)4611[0]4612sage: admcycles.genstobasis(0,3,0)4613[(1)]4614sage: cache1=deepcopy(dict(generating_indices.cache))4615sage: cache2=deepcopy(dict(admcycles.genstobasis.cache))4616sage: save_FZrels()4617sage: generating_indices.clear_cache()4618sage: admcycles.genstobasis.clear_cache()4619sage: load_FZrels()4620sage: cache1==dict(generating_indices.cache)4621True4622sage: cache2==dict(admcycles.genstobasis.cache)4623True4624"""4625with open('new_geninddb.pkl', 'wb') as file:4626dumped = {fun[0]: dict(fun[1].cache)4627for fun in enumerate(PICKLED_FUNCTIONS)}4628pickle.dump(dumped, file)462946304631def load_FZrels():4632"""4633Loading values4634"""4635try:4636with open('new_geninddb.pkl', 'rb') as file:4637dumped = pickle.load(file)4638for k, cache in dumped.items():4639PICKLED_FUNCTIONS[k].cache.update(cache)4640except IOError:4641raise IOError('Could not find file new_geninddb.pkl\nIf you previously saved FZ relations in geninddb.pkl, try old_load_FZrels() followed by save_FZrels() to go to new format.\nType old_load_FZrels? for more info.')46424643def old_load_FZrels():4644"""4645Loading values from old format in 'geninddb.pkl'.4646If you previously saved FZ relations to this file, you should first call old_load_FZrels(), followed by save_FZrels() to save them to the new location 'new_geninddb.pkl'. From now on, you can load the relations from this new file using load_FZrels().4647Once you are certain that the transfer worked, you can delete the file 'geninddb.pkl'.4648"""4649pkl_file=open('geninddb.pkl','rb')4650dumpdic=pickle.load(pkl_file)4651# Dictionary of the form4652# {('generating_indices', (0, 3, 0)): [0],4653# ('genstobasis', (2, 1, 0)): [(1)], ... }4654for k, cache in dumpdic.items():4655if k[0]=='generating_indices':4656generating_indices.set_cache(cache, *(k[1]))4657if k[0]=='genstobasis':4658genstobasis.set_cache(cache, *(k[1]))4659pkl_file.close()466046614662# looks up if the Hurwitz cycle for (g,dat) has already been computed. If yes, returns a fresh tautclass representing this4663# it takes care of the necessary reordering and possible pushforwards by forgetting markings4664# if this cycle has not been computed before, it raises a KeyError4665def Hdb_lookup(g,dat,markings):4666# first get a sorted version of dat, remembering the reordering4667li=dat.l4668sort_li=sorted(enumerate(li),key=lambda x:x[1])4669ind_reord=[i[0] for i in sort_li] # [2,0,1] means that li looking like [a,b,c] now is reordered as li_reord=[c,a,b]4670li_reord=[i[1] for i in sort_li]4671dat_reord=HurwitzData(dat.G,li_reord)46724673# now create a dictionary legpermute sending the marking-names in the space you want to the corresponding names in the standard-ordered space4674mn_list=[]4675num_marks=04676Gord=dat.G.order()4677for count in range(len(li)):4678new_marks=QQ(Gord)/li[count][1]4679mn_list.append(list(range(num_marks+1, num_marks+new_marks+1)))46804681mn_reord=[]4682for i in ind_reord:4683mn_reord.append(mn_list[i])46844685legpermute={mn_reord[j]:j+1 for j in range(len(mn_reord))}4686#markings_reord=46874688Hdb=Hdatabase[(g,dat_reord)]4689for marks in Hdb:4690if set(markings).issubset(set(marks)):4691forgottenmarks=set(marks)-set(markings)4692result=deepcopy(Hdb[marks])4693result=result.forgetful_pushforward(list(forgottenmarks))4694found=True4695#TODO: potentially include this result in Hdb46964697469846994700############4701#4702# Library of interesting cycle classes4703#4704############47054706def fundclass(g,n):4707r"""4708Return the fundamental class of \Mbar_g,n.4709"""4710if not all(isinstance(i, (int, Integer)) and i >= 0 for i in (g,n)):4711raise ValueError("g,n must be non-negative integers")4712if 2*g-2+n <= 0:4713raise ValueError("g,n must satisfy 2g-2+n>0")4714return trivgraph(g,n).to_tautclass()47154716# returns a list of generators of R^r(\bar M_{g,n}) as tautclasses in the order of Pixton's program4717# if decst=True, returns generators as decstrata, else as tautclasses4718def tautgens(g,n,r,decst=False):4719r"""4720Returns a lists of all tautological classes of degree r on \bar M_{g,n}.47214722INPUT:47234724g : integer4725Genus g of curves in \bar M_{g,n}.4726n : integer4727Number of markings n of curves in \bar M_{g,n}.4728r : integer4729The degree r of of the classes.4730decst : boolean4731If set to True returns generators as decorated strata, else as tautological classes.47324733EXAMPLES::47344735sage: from admcycles import *47364737sage: tautgens(2,0,2)[1]4738Graph : [2] [[]] []4739Polynomial : 1*(kappa_1^2 )_047404741::47424743sage: L=tautgens(2,0,2);2*L[3]+L[4]4744Graph : [1, 1] [[2], [3]] [(2, 3)]4745Polynomial : 2*psi_2^14746<BLANKLINE>4747Graph : [1] [[2, 3]] [(2, 3)]4748Polynomial : 1*(kappa_1^1 )_04749"""4750L = DR.all_strata(g,r,tuple(range(1, n+1)))4751if decst:4752return [Graphtodecstratum(l) for l in L]4753else:4754return [tautclass([Graphtodecstratum(l)]) for l in L]47554756# prints a list of generators of R^r(\bar M_{g,n}) as tautclasses in the order of Pixton's program4757def list_tautgens(g,n,r):4758r"""4759Lists all tautological classes of degree r on \bar M_{g,n}.47604761INPUT:47624763g : integer4764Genus g of curves in \bar M_{g,n}.4765n : integer4766Number of markings n of curves in \bar M_{g,n}.4767r : integer4768The degree r of of the classes.47694770EXAMPLES::47714772sage: from admcycles import *47734774sage: list_tautgens(2,0,2)4775[0] : Graph : [2] [[]] []4776Polynomial : 1*(kappa_2^1 )_04777[1] : Graph : [2] [[]] []4778Polynomial : 1*(kappa_1^2 )_04779[2] : Graph : [1, 1] [[2], [3]] [(2, 3)]4780Polynomial : 1*(kappa_1^1 )_04781[3] : Graph : [1, 1] [[2], [3]] [(2, 3)]4782Polynomial : 1*psi_2^14783[4] : Graph : [1] [[2, 3]] [(2, 3)]4784Polynomial : 1*(kappa_1^1 )_04785[5] : Graph : [1] [[2, 3]] [(2, 3)]4786Polynomial : 1*psi_2^14787[6] : Graph : [0, 1] [[3, 4, 5], [6]] [(3, 4), (5, 6)]4788Polynomial : 1*4789[7] : Graph : [0] [[3, 4, 5, 6]] [(3, 4), (5, 6)]4790Polynomial : 1*4791"""4792L=tautgens(g,n,r)4793for i in range(len(L)):4794print('['+repr(i)+'] : '+repr(L[i]))47954796def kappaclass(a, g=None, n=None):4797r"""4798Returns the (Arbarello-Cornalba) kappa-class kappa_a on \bar M_{g,n} defined by47994800kappa_a= pi_*(psi_{n+1}^{a+1})48014802where pi is the morphism \bar M_{g,n+1} --> \bar M_{g,n}.48034804INPUT:48054806a : integer4807The degree a of the kappa class.4808g : integer4809Genus g of curves in \bar M_{g,n}.4810n : integer4811Number of markings n of curves in \bar M_{g,n}.48124813EXAMPLES::48144815sage: from admcycles import *48164817sage: kappaclass(2,3,1)4818Graph : [3] [[1]] []4819Polynomial : 1*(kappa_2^1 )_048204821When working with fixed g and n for the moduli space \bar M_{g,n},4822it is possible to specify the desired value of the global variables g4823and n using ``reset_g_n`` to avoid giving them as an argument each time::48244825sage: reset_g_n(3,2)4826sage: kappaclass(1)4827Graph : [3] [[1, 2]] []4828Polynomial : 1*(kappa_1^1 )_04829"""4830if g is None:4831g=globals()['g']4832if n is None:4833n=globals()['n']4834return tautclass([decstratum(trivgraph(g,n), poly=kappacl(0,a,1,g,n))])48354836# returns psi_i on \bar M_{gloc,nloc}4837def psiclass(i,g=None,n=None):4838r"""4839Returns the class psi_i on \bar M_{g,n}.48404841INPUT:48424843i : integer4844The leg i associated to the psi class.4845g : integer4846Genus g of curves in \bar M_{g,n}.4847n : integer4848Number of markings n of curves in \bar M_{g,n}.48494850EXAMPLES::48514852sage: from admcycles import *48534854sage: psiclass(2,2,3)4855Graph : [2] [[1, 2, 3]] []4856Polynomial : 1*psi_2^148574858When working with fixed g and n for the moduli space \bar M_{g,n},4859it is possible to specify the desired value of the global variables g4860and n using ``reset_g_n`` to avoid giving them as an argument each time::48614862sage: reset_g_n(3,2)4863sage: psiclass(1)4864Graph : [3] [[1, 2]] []4865Polynomial : 1*psi_1^148664867TESTS::48684869sage: psiclass(3,2,1)4870Traceback (most recent call last):4871...4872ValueError: Index of marking for psiclass smaller than number of markings4873"""4874if g is None:4875g=globals()['g']4876if n is None:4877n=globals()['n']4878if i > n:4879raise ValueError('Index of marking for psiclass smaller than number of markings')4880return tautclass([decstratum(trivgraph(g,n), poly=psicl(i,1))])488148824883# returns the pushforward of the fundamental class under the boundary gluing map to \bar M_{gloc,nloc} associated to the partition (g1,A|g-g1, {1,...,n} \ A) of (g,n)4884def sepbdiv(g1,A,g=None,n=None):4885r"""4886Returns the pushforward of the fundamental class under the boundary gluing map \bar M_{g1,A} X \bar M_{g-g1,{1,...,n} \ A} --> \bar M_{g,n}.48874888INPUT:48894890g1 : integer4891The genus g1 of the first vertex.4892A: list4893The list A of markings on the first vertex.4894g : integer4895The total genus g of the graph.4896n : integer4897The total number of markings n of the graph.48984899EXAMPLES::49004901sage: from admcycles import *49024903sage: sepbdiv(1,(1,3,4),2,5)4904Graph : [1, 1] [[1, 3, 4, 6], [2, 5, 7]] [(6, 7)]4905Polynomial : 1*49064907When working with fixed g and n for the moduli space \bar M_{g,n},4908it is possible to specify the desired value of the global variables g4909and n using ``reset_g_n`` to avoid giving them as an argument each time::49104911sage: reset_g_n(3,3)4912sage: sepbdiv(1,(2,))4913Graph : [1, 2] [[2, 4], [1, 3, 5]] [(4, 5)]4914Polynomial : 1*4915"""4916if g is None:4917g=globals()['g']4918if n is None:4919n=globals()['n']4920return tautclass([decstratum(stgraph([g1,g-g1],[list(A)+[n+1], list(set(range(1, n+1)) - set(A))+[n+2]],[(n+1,n+2)]))])49214922# returns the pushforward of the fundamental class under the irreducible boundary gluing map \bar M_{g-1,n+2} -> \bar M_{g,n}4923def irrbdiv(g=None, n=None):4924r"""4925Returns the pushforward of the fundamental class under the irreducible boundary gluing map \bar M_{g-1,n+2} -> \bar M_{g,n}.49264927INPUT:49284929g : integer4930The total genus g of the graph.4931n : integer4932The total number of markings n of the graph.49334934EXAMPLES::49354936sage: from admcycles import *49374938sage: irrbdiv(2,5)4939Graph : [1] [[1, 2, 3, 4, 5, 6, 7]] [(6, 7)]4940Polynomial : 1*49414942When working with fixed g and n for the moduli space \bar M_{g,n},4943it is possible to specify the desired value of the global variables g4944and n using ``reset_g_n`` to avoid giving them as an argument each time::49454946sage: reset_g_n(3,0)4947sage: irrbdiv()4948Graph : [2] [[1, 2]] [(1, 2)]4949Polynomial : 1*4950"""4951if g is None:4952g=globals()['g']4953if n is None:4954n=globals()['n']4955return tautclass([decstratum(stgraph([g-1], [list(range(1, n+3))], [(n+1, n+2)]))])495649574958495949604961496249634964# tries to return the class lambda_d on \bar M_{g,n}. We hope that this is a variant of the DR-cycle4965# def lambdaish(d,g,n):4966# dvector=vector([0 for i in range(n)])4967# return (-2)**(-g) * DR_cycle(g,dvector,d)49684969# Returns the chern character ch_d(EE) of the Hodge bundle EE on \bar M_{g,n} as a mixed degree tautological class (up to maxim. degree dmax)4970# implements the formula from [Mumford - Towards an enumerative geometry ...]49714972@cached_function4973def hodge_chern_char(g, n, d):4974bdries=list_strata(g,n,1)4975irrbdry=bdries.pop(0)4976autos = [bd.automorphism_number() for bd in bdries]4977result = tautgens(g,n,0)[0] #=149784979if d==0 :4980return deepcopy(tautgens(g,n,0)[0])4981elif (d%2 == 0) or (d<0):4982return tautclass([])4983else:4984psipsisum_onevert=sum([((-1)**i)*(psicl(n+1, 1)**i)*(psicl(n+2,1)**(d-1-i)) for i in range(d)])4985psipsisum_twovert=sum([((-1)**i)*(psicl(n+1, 2)**i)*(psicl(n+2, 2)**(d-1-i)) for i in range(d)])49864987contrib=kappaclass(d,g,n)-sum([psiclass(i,g,n)**d for i in range(1, n+1)])498849894990# old: contrib=kappaclass(d,g,n)-sum([psiclass(i,g,n) for i in range(1,n+1)])4991contrib+=(QQ(1)/2)*tautclass([decstratum(irrbdry,poly=psipsisum_onevert)])4992contrib+=sum([(QQ(1)/autos[ind])*tautclass([decstratum(bdries[ind], poly=psipsisum_twovert)])for ind in range(len(bdries))])49934994contrib.dimension_filter()49954996return (bernoulli(d+1)/factorial(d+1)) * contrib49974998# converts the function chclass (sending m to ch_m) to the corresponding chern polynomial, up to degree dmax4999def chern_char_to_poly(chclass,dmax,g,n):5000result=deepcopy(tautgens(g,n,0)[0])50015002argum=sum([factorial(m-1)*chclass(m) for m in range(1, dmax+1)])5003expo=deepcopy(argum)5004result=result+argum50055006for m in range(2, dmax+1):5007expo*=argum5008expo.degree_cap(dmax)5009result+=(QQ(1)/factorial(m))*expo50105011return result5012501350145015def chern_char_to_class(t,char,g=None,n=None):5016r"""5017Turns a Chern character into a Chern class in the tautological ring of \Mbar_{g,n}.50185019INPUT:50205021t : integer5022degree of the output Chern class5023char : tautclass or function5024Chern character, either represented as a (mixed-degree) tautological class or as5025a function m -> ch_m50265027EXAMPLES:50285029Note that the output of generalized_hodge_chern takes the form of a chern character::50305031sage: from admcycles import *5032sage: from admcycles.GRRcomp import *5033sage: g=2;n=2;l=0;d=[1,-1];a=[[1,[1],-1]]5034sage: chern_char_to_class(1,generalized_hodge_chern(l,d,a,1,g,n))5035Graph : [2] [[1, 2]] []5036Polynomial : 1/12*(kappa_1^1 )_05037<BLANKLINE>5038Graph : [2] [[1, 2]] []5039Polynomial : (-13/12)*psi_1^15040<BLANKLINE>5041Graph : [2] [[1, 2]] []5042Polynomial : (-1/12)*psi_2^15043<BLANKLINE>5044Graph : [0, 2] [[1, 2, 4], [5]] [(4, 5)]5045Polynomial : 1/12*5046<BLANKLINE>5047Graph : [1, 1] [[4], [1, 2, 5]] [(4, 5)]5048Polynomial : 1/12*5049<BLANKLINE>5050Graph : [1, 1] [[2, 4], [1, 5]] [(4, 5)]5051Polynomial : 13/12*5052<BLANKLINE>5053Graph : [1] [[4, 5, 1, 2]] [(4, 5)]5054Polynomial : 1/24*5055"""5056if g is None:5057g=globals()['g']5058if n is None:5059n=globals()['n']50605061if isinstance(char,tautclass):5062arg = [(-1)**(s-1)*factorial(s-1)*char.simplified(r=s) for s in range(1,t+1)]5063else:5064arg = [(-1)**(s-1)*factorial(s-1)*char(s) for s in range(1,t+1)]5065exp = sum(multinomial(s.to_exp())/factorial(len(s))*prod(arg[k-1] for k in s) for s in Partitions(t))5066if t==0:5067return exp * fundclass(g, n)5068return exp.simplify(r=t)5069## the degree t part of the exponential function5070## we sum over the partitions of degree element to give a degree t element5071## The length of a partition is the power in the exponent we are looking at5072## Then we need to see the number of ways to devide the degrees over the x's to get the right thing.5073507450755076# returns lambda-class lambda_d on \bar M_{g,n}5077# r"""Returns lambda-class lambda_d on \\bar M_{g,n}"""5078def lambdaclass_old(d,g=None,n=None):5079r"""5080Old implementation of lambda_d on \bar M_{g,n} defined as the d-th Chern class50815082lambda_d = c_d(E)50835084of the Hodge bundle E. The result is represented as a sum of stable graphs with kappa and psi classes.50855086INPUT:50875088d : integer5089The degree d.5090g : integer5091Genus g of curves in \bar M_{g,n}.5092n : integer5093Number of markings n of curves in \bar M_{g,n}.50945095EXAMPLES::50965097sage: from admcycles.admcycles import lambdaclass_old50985099sage: lambdaclass_old(1,2,1)5100Graph : [2] [[1]] []5101Polynomial : 1/12*(kappa_1^1 )_05102<BLANKLINE>5103Graph : [2] [[1]] []5104Polynomial : (-1/12)*psi_1^15105<BLANKLINE>5106Graph : [1, 1] [[3], [1, 4]] [(3, 4)]5107Polynomial : 1/12*5108<BLANKLINE>5109Graph : [1] [[3, 4, 1]] [(3, 4)]5110Polynomial : 1/24*51115112When working with fixed g and n for the moduli space \bar M_{g,n},5113it is possible to specify the desired value of the global variables g5114and n using ``reset_g_n`` to avoid giving them as an argument each time::51155116sage: from admcycles import reset_g_n5117sage: reset_g_n(1,1)5118sage: lambdaclass_old(1)5119Graph : [1] [[1]] []5120Polynomial : 1/12*(kappa_1^1 )_05121<BLANKLINE>5122Graph : [1] [[1]] []5123Polynomial : (-1/12)*psi_1^15124<BLANKLINE>5125Graph : [0] [[3, 4, 1]] [(3, 4)]5126Polynomial : 1/24*51275128TESTS::51295130sage: from admcycles.admcycles import lambdaclass, lambdaclass_old5131sage: L1 = lambdaclass_old(3, 3, 0)5132sage: L2 = lambdaclass(3, 3, 0)5133sage: (L1 - L2).toTautvect().is_zero()5134True5135"""5136if g is None:5137g=globals()['g']5138if n is None:5139n=globals()['n']51405141if d>g:5142return tautclass([])51435144def chcl(m):5145return hodge_chern_char(g,n,m)51465147cherpoly=chern_char_to_poly(chcl,d,g,n)51485149result=cherpoly.degree_part(d)5150result.simplify(g,n,d)5151return result515251535154# returns lambda-class lambda_d on \bar M_{g,n}5155# r"""Returns lambda-class lambda_d on \\bar M_{g,n}"""5156def lambdaclass(d, g=None, n=None, pull=True):5157r"""5158Returns the tautological class lambda_d on \bar M_{g,n} defined as the d-th Chern class51595160lambda_d = c_d(E)51615162of the Hodge bundle E. The result is represented as a sum of stable graphs with kappa and psi classes.51635164INPUT:51655166d : integer5167The degree d.5168g : integer5169Genus g of curves in \bar M_{g,n}.5170n : integer5171Number of markings n of curves in \bar M_{g,n}.5172pull : boolean, default = True5173Whether lambda_d on \bar M_{g,n} is computed as pullback from \bar M_{g}51745175EXAMPLES::51765177sage: from admcycles import *51785179sage: lambdaclass(1,2,0)5180Graph : [2] [[]] []5181Polynomial : 1/12*(kappa_1^1 )_05182<BLANKLINE>5183Graph : [1, 1] [[2], [3]] [(2, 3)]5184Polynomial : 1/24*5185<BLANKLINE>5186Graph : [1] [[2, 3]] [(2, 3)]5187Polynomial : 1/24*51885189When working with fixed g and n for the moduli space \bar M_{g,n},5190it is possible to specify the desired value of the global variables g5191and n using ``reset_g_n`` to avoid giving them as an argument each time::51925193sage: reset_g_n(1,1)5194sage: lambdaclass(1)5195Graph : [1] [[1]] []5196Polynomial : 1/12*(kappa_1^1 )_05197<BLANKLINE>5198Graph : [1] [[1]] []5199Polynomial : (-1/12)*psi_1^15200<BLANKLINE>5201Graph : [0] [[3, 4, 1]] [(3, 4)]5202Polynomial : 1/24*52035204TESTS::52055206sage: from admcycles import lambdaclass5207sage: inputs = [(0,0,4), (1,1,3), (1,2,1), (2,2,1), (3,2,1), (-1,2,1), (2,3,2)]5208sage: for d,g,n in inputs:5209....: assert (lambdaclass(d, g, n)-lambdaclass(d, g, n, pull=False)).is_zero()52105211"""5212if g is None:5213g=globals()['g']5214if n is None:5215n=globals()['n']52165217if d > g or d < 0:5218return tautclass([])52195220if n > 0 and pull:5221if g == 0:5222return fundclass(g,n) # Note: previous check on d forces d == 0 in this case5223if g == 1:5224newmarks = list(range(2, n+1))5225return lambdaclass(d, 1, 1, pull=False).forgetful_pullback(newmarks)5226newmarks = list(range(1, n+1))5227return lambdaclass(d, g, 0, pull=False).forgetful_pullback(newmarks)52285229def chcl(m):5230return hodge_chern_char(g,n,m)52315232return chern_char_to_class(d,chcl,g,n)5233523452355236def barH(g,dat,markings=None):5237"""Returns \\bar H on genus g with Hurwitz datum dat as a prodHclass on the trivial graph, remembering only the marked points given in markings"""5238Hbar = trivGgraph(g,dat)5239n = len(Hbar.gamma.list_markings())52405241#global Hexam5242#Hexam=(g,dat,method,vecout,redundancy,markings)52435244if markings is None:5245markings = list(range(1, n + 1))52465247N=len(markings)5248msorted=sorted(markings)5249markingdictionary={i+1: msorted[i] for i in range(N)}5250return prodHclass(stgraph([g],[list(range(1, N+1))],[]),[decHstratum(stgraph([g],[list(range(1 ,N+1))],[]), {0: (g,dat)}, [[0 ,markingdictionary]])])525152525253def Hyperell(g,n=0 ,m=0):5254"""Returns the cycle class of the hyperelliptic locus of genus g curves with n marked fixed points and m pairs of conjugate points in \\barM_{g,n+2m}.52555256TESTS::52575258sage: from admcycles import *5259sage: H=Hyperell(3) # long time5260sage: H.toTautbasis() # long time5261(3/4, -9/4, -1/8)5262"""5263if n>2 *g+2:5264print('A hyperelliptic curve of genus '+repr(g)+' can only have '+repr(2 *g+2)+' fixed points!')5265return 05266if 2 *g-2 +n+2 *m<=0:5267print('The moduli space \\barM_{'+repr(g)+','+repr(n+2 *m)+'} does not exist!')5268return 052695270G=PermutationGroup([(1, 2)])5271H=HurData(G,[G[1] for i in range(2 *g+2)]+[G[0] for i in range(m)])5272marks=list(range(1, n+1))+list(range(2 *g+2 + 1, 2 *g+2 + 1 +2 *m))5273factor=QQ(1)/factorial(2 *g+2 -n)5274result= factor*Hidentify(g,H,markings=marks)52755276# currently, the markings on result are named 1,2,..,n, 2*g+3, ..., 2*g+2+2*m5277# shift the last set of markings down by 2*g+3-(n+1)5278rndict={i:i for i in range(1, n+1)}5279rndict.update({j:j-(2 *g+2 -n) for j in range(2*g+2+1, 2*g+2+1+2*m)})5280result.rename_legs(rndict,True)52815282return result52835284def Biell(g,n=0 ,m=0):5285r"""5286Returns the cycle class of the bielliptic locus of genus ``g`` curves with ``n`` marked fixed points and ``m`` pairs of conjugate points in `\bar M_{g,n+2m}`.52875288TESTS::52895290sage: from admcycles import *5291sage: B=Biell(2) # long time5292sage: B.toTautbasis() # long time5293(15/2, -9/4)5294"""5295if g==0:5296print('There are no bielliptic curves of genus 0!')5297return 05298if n>2 *g-2:5299print('A bielliptic curve of genus '+repr(g)+' can only have '+repr(2 *g-2)+' fixed points!')5300return 05301if 2 *g-2 +n+2 *m<=0:5302print('The moduli space \\barM_{'+repr(g)+','+repr(n+2 *m)+'} does not exist!')5303return 053045305G=PermutationGroup([(1, 2)])5306H=HurData(G,[G[1] for i in range(2 * g - 2)]+[G[0] for i in range(m)])5307marks=list(range(1, n+1))+list(range(2*g-2+1, 2*g-2+1+2*m))5308factor= QQ((1,factorial(2 *g-2 -n)))5309if g==2 and n==0 and m==0:5310factor/=25311result= factor*Hidentify(g,H,markings=marks)53125313# currently, the markings on result are named 1,2,..,n, 2*g-1, ..., 2*g-2+2*m5314# shift the last set of markings down by 2*g-1-(n+1)5315rndict={i:i for i in range(1, n+1)}5316rndict.update({j:j-(2 *g-2 -n) for j in range(2*g-2+1, 2*g-2+1+2*m)})5317result.rename_legs(rndict,True)53185319return result5320############5321#5322# Transfer maps5323#5324############53255326def Hpullpush(g,dat,alpha):5327"""Pulls the class alpha to the space \\bar H_{g,dat} via map i forgetting the action, then pushes forward under delta."""5328Hb = Htautclass([Hdecstratum(trivGgraph(g,dat),poly=onekppoly(1))])5329Hb = Hb * alpha5330return Hb.quotient_pushforward()53315332def FZreconstruct(g,n,r):5333genind=generating_indices(g,n,r)5334gentob=genstobasis(g,n,r)5335numgens=len(gentob)5336M=[]53375338for i in range(numgens):5339v=insertvec(gentob[i],numgens,genind)5340v[i]-=15341M.append(v)53425343return matrix(M)534453455346def insertvec(w,length,positions):5347v=vector(QQ,length)5348for j in range(len(positions)):5349v[positions[j]]=w[j]5350return v53515352############5353#5354# Test routines5355#5356############535753585359def pullpushtest(g,dat,r):5360r"""Test if for Hurwitz space specified by (g,dat), pulling back codimension r relations under the source map and pushing forward under the target map gives relations.5361"""5362Gr=trivGgraph(g,dat)5363n=Gr.gamma.n()53645365Grquot=Gr.quotient_graph()5366gquot=Grquot.g()5367nquot=Grquot.n()53685369M=FZreconstruct(g,n,r)53705371N=matrix([Hpullpush(g,dat,alpha).toTautbasis(gquot,nquot,r) for alpha in tautgens(g,n,r)])53725373return (N.transpose()*M.transpose()).is_zero()53745375for v in M.rows():5376if not v.is_zero():5377alpha=Tautv_to_tautclass(v,g,n,r)5378pb=Hpullpush(g,dat,alpha)5379print(pb.is_zero())538053815382# compares intersection numbers of Pixton generators computed by Pixton's program and by own program5383# computes all numbers between generators for \bar M_{g,n} in degrees r, 3g-3+n-r, respectively5384def checkintnum(g,n,r):5385markings=tuple(range(1, n+1))5386Mpixton=DR.pairing_matrix(g,r,markings)53875388strata1=[Graphtodecstratum(Grr) for Grr in DR.all_strata(g,r,markings)]5389strata2=[Graphtodecstratum(Grr) for Grr in DR.all_strata(g,3 *g-3 +n-r,markings)]5390Mself=[[(s1*s2).evaluate() for s2 in strata2] for s1 in strata1]53915392return Mpixton==Mself539353945395# pulls back all generators of degree r to the boundary divisors and identifies them, checks if results are compatible with FZ-relations between generators5396def pullandidentify(g,n,r):5397markings=tuple(range(1, n+1))5398M=DR.FZ_matrix(g,r,markings)5399Mconv=matrix([DR.convert_vector_to_monomial_basis(M.row(i),g,r,markings) for i in range(M.nrows())])54005401L=list_strata(g,n,1)5402strata = DR.all_strata(g, r, markings)5403decst=[Graphtodecstratum(Grr) for Grr in strata]54045405for l in L:5406pullbacks=[((l.boundary_pullback(st)).dimension_filter()).totensorTautbasis(r,True) for st in decst]5407for i in range(Mconv.nrows()):5408vec=0 *pullbacks[0]5409for j in range(Mconv.ncols()):5410if Mconv[i,j]!=0:5411vec+=Mconv[i,j]*pullbacks[j]5412if not vec.is_zero():5413return False5414return True541554165417# computes all intersection numbers of generators of R^(r1)(\bar M_{g,n1}) with pushforwards of generators of R^(r2)(\bar M_{g,n2}), where n1<n25418# compares with result when pulling back and intersecting5419# condition: r1 + r2- (n2-n1) = 3*g-3+n15420def pushpullcompat(g,n1,n2,r1):5421#global cu5422#global cd54235424r2=3 *g-3 +n2-r154255426strata_down=[tautclass([Graphtodecstratum(Grr)]) for Grr in DR.all_strata(g,r1,tuple(range(1, n1+1)))]5427strata_up=[tautclass([Graphtodecstratum(Grr)]) for Grr in DR.all_strata(g,r2,tuple(range(1, n2+1)))]54285429fmarkings=list(range(n1+1, n2+1))54305431for class_down in strata_down:5432for class_up in strata_up:5433intersection_down = (class_down*class_up.forgetful_pushforward(fmarkings)).evaluate()5434intersection_up = (class_down.forgetful_pullback(fmarkings) * class_up).evaluate()54355436if not intersection_down == intersection_up:5437return False5438return True543954405441# test delta-pullback composed with delta-pushforward is multiplication with delta-degree for classes in R^r(\bar M_{g',b}) on space of quotient curves5442def deltapullpush(g,H,r):5443trivGg=trivGgraph(g,H)5444trivGclass=Htautclass([Hdecstratum(trivGg)])5445ddeg=trivGg.delta_degree(0)5446quotgr=trivGg.quotient_graph()5447gprime = quotgr.genera(0)5448bprime = quotgr.num_legs(0)54495450gens=tautgens(gprime,bprime,r)54515452for cl in gens:5453clpp=trivGclass.quotient_pullback(cl).quotient_pushforward()5454if not (clpp.toTautvect(gprime,bprime,r)==ddeg*cl.toTautvect(gprime,bprime,r)):5455return False54565457#print 'The class \n'+repr(cl)+'\npulls back to\n'+repr(trivGclass.quotient_pullback(cl))+'\nunder delta*.\n'5458return True545954605461## Lambda-classes must satisfy this:5462def lambdaintnumcheck(g):5463intnum=((lambdaclass(g,g,0)*lambdaclass(g-1, g, 0)).simplify()*lambdaclass(g-2 ,g, 0)).evaluate()5464result= QQ((-1,2)) / factorial(2 *(g-1)) * bernoulli(2 *(g-1)) / (2*(g-1))*bernoulli(2 *g) / (2 *g)5465return intnum==result54665467@cached_function5468def op_subset_space(g,n,r,moduli):5469r"""5470Returns a vector space W over QQ which is the quotient W = V / U, where5471V = space of vectors representing classes in R^r(\Mbar_{g,n}) in a basis5472U = subspace of vectors representing classes supported outside the open subset5473of \Mbar_{g,n} specified by moduli (= 'sm', 'rt', 'ct' or 'tl')5474Thus for v in V, the vector W(v) is the representation of the class represented by v5475in a basis of R^r(M^{moduli}_{g,n}).5476"""5477L = tautgens(g,n,r)5478Msm=[(0*L[0]).toTautbasis(g,n,r)] #start with zero vector to get dimensions right later54795480if moduli == 'st':5481def modtest(gr):5482return False5483if moduli == 'sm':5484def modtest(gr):5485return bool(gr.num_edges())5486if moduli == 'rt':5487def modtest(gr):5488return g not in gr.genera(copy = False)5489if moduli == 'ct':5490def modtest(gr):5491return gr.num_edges() - gr.num_verts() + 1 != 05492if moduli == 'tl':5493def modtest(gr):5494return gr.num_edges() - gr.num_loops() - gr.num_verts() + 1 != 054955496for t in L:5497if modtest(t.terms[0].gamma):5498# t is supported on a graph outside the locus specified by moduli5499Msm.append(t.toTautbasis())55005501MsmM=matrix(QQ,Msm)5502U=MsmM.row_space()5503W=U.ambient_vector_space()/U5504return W55055506############5507#5508# Further doctests5509#5510############5511r"""5512EXAMPLES::55135514sage: from admcycles import *5515sage: from admcycles.admcycles import pullpushtest, deltapullpush, checkintnum, pullandidentify, pushpullcompat, lambdaintnumcheck5516sage: G=PermutationGroup([(1,2)])5517sage: H=HurData(G,[G[1],G[1]])5518sage: pullpushtest(2,H,1)5519True5520sage: pullpushtest(2,H,2) # long time5521True5522sage: deltapullpush(2,H,1)5523True5524sage: deltapullpush(2,H,2) # long time5525True5526sage: G=PermutationGroup([(1,2,3)])5527sage: H=HurData(G,[G[0]])5528sage: c=Hidentify(1,H)5529sage: c.toTautbasis()5530(5, -3, 2, 2, 2)5531sage: checkintnum(2,0,1)5532True5533sage: checkintnum(2,1,2)5534True5535sage: pullandidentify(2,1,2)5536True5537sage: pullandidentify(3,0,2)5538True5539sage: pushpullcompat(2,0,1,2) # long time5540True5541sage: lambdaintnumcheck(2)5542True5543sage: lambdaintnumcheck(3) # long time5544True55455546"""554755485549