Source code for admcycles.DR

# This file was *autogenerated* from the file DR.sage
from __future__ import absolute_import

from copy import copy, deepcopy

from sage.rings.all import Integer, RealNumber, PolynomialRing, ZZ
from sage.arith.all import factorial, binomial, bernoulli
from sage.matrix.constructor import matrix
from sage.combinat.all import Combinations, IntegerVectors, Partitions, Permutations, Subsets
from sage.functions.other import floor

_sage_const_3 = Integer(3); _sage_const_2 = Integer(2); _sage_const_1 = Integer(1); _sage_const_0 = Integer(0); _sage_const_7 = Integer(7); _sage_const_6 = Integer(6); _sage_const_5 = Integer(5); _sage_const_4 = Integer(4); _sage_const_8 = Integer(8); _sage_const_28 = Integer(28); _sage_const_40 = Integer(40); _sage_const_2520 = Integer(2520); _sage_const_84 = Integer(84); _sage_const_112 = Integer(112); _sage_const_420 = Integer(420); _sage_const_90 = Integer(90); _sage_const_840 = Integer(840); _sage_const_12 = Integer(12); _sage_const_10 = Integer(10); _sage_const_1000 = Integer(1000); _sage_const_15 = Integer(15); _sage_const_30 = Integer(30); _sage_const_70 = Integer(70); _sage_const_168 = Integer(168); _sage_const_100 = Integer(100); _sage_const_280 = Integer(280)# Sage code by Aaron Pixton (apixton@mit.edu)
#
# This file contains code to compute the double ramification cycle and simplify
# it using the 3-spin tautological relations.
#
# The functions are nearly all undocumented. The main function to use to
# compute DR cycles is the following:
# - DR_compute(g,r,n,dvector,kval=0,moduli_type=MODULI_ST):
#     g = genus
#     r = cohomological degree (set to g for the actual DR cycle, should give
#                               tautological relations for r > g)
#     n = number of marked points
#     dvector = vector of weights to place on the marked points, should have
#               sum zero in the case of the DR cycle
#     kval: 0 by default (DR cycles), can be set to other integers to twist
#           by copies of the log canonical bundle (e.g. 1 for differentials)
#     moduli_type: MODULI_ST by default (moduli of stable curves), can be
#                  set to moduli spaces containing less of the boundary (e.g.
#                  MODULI_CT for compact type) to compute there instead
#
#   Example: to compute the genus 1 DR cycle with weight vector (2,-2):
#
#   sage: DR_compute(1,1,2,(2,-2))
#   (0, 2, 2, 0, -1/24)
#
#   The function returns a vector of rational numbers giving the coefficients
#   of the standard tautological generators (xi_Gamma)_*(kappa-psi monomial).
#   To see how the generators are ordered, you can view a list of them by the
#   command:
#
#   sage: list_strata(1,1,2)
#   generator 0
#   [   -1     1     2]
#   [X + 1     1     1]
#   -------------------------
#   generator 1
#   [   -1     1     2]
#   [    1 X + 1     1]
#   -------------------------
#   generator 2
#   [   -1     1     2]
#   [    1     1 X + 1]
#   -------------------------
#   generator 3
#   [-1  1  2  0]
#   [ 0  1  1  1]
#   [ 1  0  0  1]
#   -------------------------
#   generator 4
#   [-1  0  1  2]
#   [ 0  2  1  1]
#   -------------------------
#
#   Here each generator is represented by an augmented vertex-edge incidence
#   matrix:
#   - each row after the first one corresponds to a vertex
#   - each column after the first one corresponds to an edge or leg
#   - the first column gives the genera of the vertices
#   - the first row gives the markings on the legs
#   - the other cells are 0 if the vertex and edge/leg are not incident
#                         1 if the vertex and edge/leg are incident once
#                         2 if the edge is a loop at the vertex
#   - entries with polynomials in X describe kappa/psi decorations:
#     - in the first column, each X^n term corresponds to a kappa_n at that
#       vertex
#     - in other locations, each X term corresponds to a psi at that half-edge
#     - at loops, 2 + aX + bX^2 corresponds to having psi^a on one side of the
#       loop and psi^b on the other
#
#   In the example above, the five classes described on R^1(Mbar_{1,2}) are:
#     kappa_1, psi_1, psi_2, delta_{12}, delta_{irr}.
#   (Here by delta_{irr} we mean the pushforward of 1 under the gluing map
#   Mbar_{0,4} -> Mbar_{1,2}, so twice the class of the physical locus.)
#
#
#   There are a couple other functions that might be convenient when computing
#   DR cycles:
#   - DR_sparse: same as DR_compute except that it returns the answer as a
#                sparse vector, e.g.
#                sage: DR_sparse(1,1,2,(2,-2))
#                [[1, 2], [2, 2], [4, -1/24]]
#   - DR_reduced: same as DR_compute except that it only requires two arguments
#                 (g and dvector) and simplifies the answer using the 3-spin
#                 tautological relations, e.g.
#                 sage: DR_reduced(1,(2,-2))
#                 (0, 0, 0, 4, 1/8)
import itertools

from .cache import capply

R = PolynomialRing(ZZ,_sage_const_1 ,order='lex', names=('X',)); (X,) = R._first_ngens(1)
A_list = [factorial(_sage_const_6 *n)/(factorial(_sage_const_3 *n)*factorial(_sage_const_2 *n)) for n in range(_sage_const_100 )]
B_list = [factorial(_sage_const_6 *n+_sage_const_1 )/((_sage_const_6 *n-_sage_const_1 )*factorial(_sage_const_3 *n)*factorial(_sage_const_2 *n)) for n in range(_sage_const_100 )]
MODULI_SMALL = -_sage_const_1
MODULI_SM = _sage_const_0
MODULI_RT = _sage_const_1
MODULI_CT = _sage_const_2
MODULI_ST = _sage_const_3
ENABLE_DPRINT = True
ENABLE_DSAVE = False

[docs]def dprint(str,*args): if ENABLE_DPRINT: print str % args
[docs]def dsave(str,*args): if ENABLE_DSAVE: save(_sage_const_0 ,str % args)
[docs]class Graph: def __init__(self,M=None,genus_list=None): if M: self.M = copy(M) elif genus_list: self.M = matrix(R,len(genus_list)+_sage_const_1 ,_sage_const_1 ,[-_sage_const_1 ]+genus_list) else: self.M = matrix(R,_sage_const_1 ,_sage_const_1 ,-_sage_const_1 )
[docs] def num_vertices(self): return self.M.nrows() - _sage_const_1
[docs] def num_edges(self): return self.M.ncols() - _sage_const_1
[docs] def h1(self): return self.M.ncols()-self.M.nrows()+_sage_const_1
[docs] def add_vertex(self,g): self.M = self.M.stack(matrix(_sage_const_1 ,self.M.ncols())) self.M[-_sage_const_1 ,_sage_const_0 ] = g
[docs] def add_edge(self,i1,i2,marking=_sage_const_0 ): self.M = self.M.augment(matrix(self.M.nrows(),_sage_const_1 )) self.M[_sage_const_0 ,-_sage_const_1 ] = marking if i1 > _sage_const_0 : self.M[i1,-_sage_const_1 ] += _sage_const_1 if i2 > _sage_const_0 : self.M[i2,-_sage_const_1 ] += _sage_const_1
[docs] def del_vertex(self,i): self.M = self.M[:i].stack(self.M[(i+_sage_const_1 ):])
[docs] def del_edge(self,i): if i == self.num_edges(): self.M = self.M[_sage_const_0 :,:i] else: self.M = self.M[_sage_const_0 :,:i].augment(self.M[_sage_const_0 :,(i+_sage_const_1 ):])
[docs] def compute_degree_vec(self): self.degree_vec = [_sage_const_0 for i in range(_sage_const_1 ,self.M.nrows())] for i in range(_sage_const_1 ,self.M.nrows()): for j in range(_sage_const_1 ,self.M.ncols()): self.degree_vec[i-_sage_const_1 ] += self.M[i,j][_sage_const_0 ]
[docs] def degree(self,i): return self.degree_vec[i-_sage_const_1 ]
[docs] def split_vertex(self,i,row1,row2): self.M = self.M.stack(matrix(_sage_const_2 ,self.M.ncols(),row1+row2)) self.add_edge(self.M.nrows()-_sage_const_2 , self.M.nrows()-_sage_const_1 ) self.del_vertex(i)
[docs] def set_target_parity(self): self.target_parity = _sage_const_0 for i in range(_sage_const_1 ,self.M.nrows()): local_parity = _sage_const_1 + self.M[i,_sage_const_0 ][_sage_const_0 ] for j in range(_sage_const_1 ,self.M.ncols()): local_parity += self.M[i,j][_sage_const_1 ] + self.M[i,j][_sage_const_2 ] for j in range(_sage_const_1 ,self.M[i,_sage_const_0 ].degree()+_sage_const_1 ): local_parity += j*self.M[i,_sage_const_0 ][j] local_parity %= _sage_const_2 self.target_parity += (local_parity << (i-_sage_const_1 ))
[docs] def replace_vertex_with_graph(self,i,G): nv = self.num_vertices() ne = self.num_edges() #i should have degree d, there should be no classes near i, and G should have markings 1,...,d and genus equal to the genus of i hedge_list = [] for k in range(_sage_const_1 ,self.M.ncols()): for j in range(self.M[i,k]): hedge_list.append(k) self.del_vertex(i) for j in range(G.num_edges() - len(hedge_list)): self.add_edge(_sage_const_0 ,_sage_const_0 ) for j in range(G.num_vertices()): self.add_vertex(G.M[j+_sage_const_1 ,_sage_const_0 ]) col = ne+_sage_const_1 for k in range(_sage_const_1 ,G.M.ncols()): if G.M[_sage_const_0 ,k] > _sage_const_0 : mark = ZZ(G.M[_sage_const_0 ,k]) for j in range(G.num_vertices()): if self.M[nv+j,hedge_list[mark-_sage_const_1 ]] == _sage_const_0 : self.M[nv+j,hedge_list[mark-_sage_const_1 ]] = G.M[j+_sage_const_1 ,k] elif G.M[j+_sage_const_1 ,k] != _sage_const_0 : a = self.M[nv+j,hedge_list[mark-_sage_const_1 ]][_sage_const_1 ] b = G.M[j+_sage_const_1 ,k][_sage_const_1 ] self.M[nv+j,hedge_list[mark-_sage_const_1 ]] = _sage_const_2 + max(a,b)*X + min(a,b)*X**_sage_const_2 else: for j in range(G.num_vertices()): self.M[nv+j,col] = G.M[j+_sage_const_1 ,k] col += _sage_const_1
[docs] def compute_invariant(self): nr,nc = self.M.nrows(),self.M.ncols() self.invariant = [[self.M[i,_sage_const_0 ], [], [], [[] for j in range(_sage_const_1 ,nr)]] for i in range(_sage_const_1 ,nr)] for k in range(_sage_const_1 ,nc): L = [i for i in range(_sage_const_1 ,nr) if self.M[i,k] != _sage_const_0 ] if len(L) == _sage_const_1 : if self.M[_sage_const_0 ,k] != _sage_const_0 : self.invariant[L[_sage_const_0 ]-_sage_const_1 ][_sage_const_2 ].append((self.M[_sage_const_0 ,k],self.M[L[_sage_const_0 ],k])) else: self.invariant[L[_sage_const_0 ]-_sage_const_1 ][_sage_const_1 ].append(self.M[L[_sage_const_0 ],k]) else: self.invariant[L[_sage_const_0 ]-_sage_const_1 ][_sage_const_3 ][L[_sage_const_1 ]-_sage_const_1 ].append((self.M[L[_sage_const_0 ],k],self.M[L[_sage_const_1 ],k])) self.invariant[L[_sage_const_1 ]-_sage_const_1 ][_sage_const_3 ][L[_sage_const_0 ]-_sage_const_1 ].append((self.M[L[_sage_const_1 ],k],self.M[L[_sage_const_0 ],k])) for i in range(_sage_const_1 ,nr): self.invariant[i-_sage_const_1 ][_sage_const_3 ] = [term for term in self.invariant[i-_sage_const_1 ][_sage_const_3 ] if len(term) > _sage_const_0 ] for k in range(len(self.invariant[i-_sage_const_1 ][_sage_const_3 ])): self.invariant[i-_sage_const_1 ][_sage_const_3 ][k].sort() self.invariant[i-_sage_const_1 ][_sage_const_3 ][k] = tuple(self.invariant[i-_sage_const_1 ][_sage_const_3 ][k]) self.invariant[i-_sage_const_1 ][_sage_const_3 ].sort() self.invariant[i-_sage_const_1 ][_sage_const_3 ] = tuple(self.invariant[i-_sage_const_1 ][_sage_const_3 ]) self.invariant[i-_sage_const_1 ][_sage_const_2 ].sort() self.invariant[i-_sage_const_1 ][_sage_const_2 ] = tuple(self.invariant[i-_sage_const_1 ][_sage_const_2 ]) self.invariant[i-_sage_const_1 ][_sage_const_1 ].sort() self.invariant[i-_sage_const_1 ][_sage_const_1 ] = tuple(self.invariant[i-_sage_const_1 ][_sage_const_1 ]) self.invariant[i-_sage_const_1 ] = tuple(self.invariant[i-_sage_const_1 ]) vertex_invariants = [[i,self.invariant[i-_sage_const_1 ]] for i in range(_sage_const_1 ,nr)] self.invariant.sort() self.invariant = tuple(self.invariant) vertex_invariants.sort(key=lambda x: x[_sage_const_1 ]) self.vertex_groupings = [] for i in range(nr-_sage_const_1 ): if i == _sage_const_0 or vertex_invariants[i][_sage_const_1 ] != vertex_invariants[i-_sage_const_1 ][_sage_const_1 ]: self.vertex_groupings.append([]) self.vertex_groupings[-_sage_const_1 ].append(vertex_invariants[i][_sage_const_0 ])
[docs] def purify(self): for i in range(self.M.nrows()): for j in range(self.M.ncols()): self.M[i,j] = R(self.M[i,j][_sage_const_0 ])
[docs] def contract(self,i,vlist,elist): # assumes graph is undecorated if self.M[_sage_const_0 ,i] != _sage_const_0 : print "ERROR: cannot contract a marking" return S = [row for row in range(_sage_const_1 ,self.M.nrows()) if self.M[row,i] != _sage_const_0 ] if len(S) == _sage_const_1 : self.M[S[_sage_const_0 ],_sage_const_0 ] += _sage_const_1 self.del_edge(i) elist = elist[:(i-_sage_const_1 )] + elist[i:] else: self.del_edge(i) elist = elist[:(i-_sage_const_1 )] + elist[i:] self.add_vertex(_sage_const_0 ) self.M[-_sage_const_1 ] += self.M[S[_sage_const_0 ]] self.M[-_sage_const_1 ] += self.M[S[_sage_const_1 ]] self.del_vertex(S[_sage_const_1 ]) self.del_vertex(S[_sage_const_0 ]) vlist = vlist[:(S[_sage_const_0 ]-_sage_const_1 )] + vlist[S[_sage_const_0 ]:(S[_sage_const_1 ]-_sage_const_1 )] + vlist[S[_sage_const_1 ]:] + [vlist[S[_sage_const_0 ]-_sage_const_1 ] + vlist[S[_sage_const_1 ]-_sage_const_1 ]] return vlist,elist
[docs]def graph_isomorphic(G1,G2): if G1.invariant != G2.invariant: return False else: return isomorphic(G1.M,G2.M,G1.vertex_groupings,G2.vertex_groupings)
[docs]def isomorphic(M1,M2,group1,group2): nr,nc = M1.nrows(),M1.ncols() PermList = [Permutations(range(len(group))) for group in group1] for sigma_data in itertools.product(*PermList): sigma = [_sage_const_0 for i in range(nr-_sage_const_1 )] for i in range(len(group1)): for j in range(len(group1[i])): sigma[group1[i][j]-_sage_const_1 ] = group2[i][sigma_data[i][j]] good = True for i in range(_sage_const_1 ,nr): ii = sigma[i-_sage_const_1 ] for j in range(_sage_const_1 ,i): jj = sigma[j-_sage_const_1 ] L1 = [] for k in range(_sage_const_1 ,nc): if M1[i,k] != _sage_const_0 and M1[j,k] != _sage_const_0 : L1.append([M1[i,k],M1[j,k]]) L1.sort() L2 = [] for k in range(_sage_const_1 ,nc): if M2[ii,k] != _sage_const_0 and M2[jj,k] != _sage_const_0 : L2.append([M2[ii,k],M2[jj,k]]) L2.sort() if L1 != L2: good = False break if good == False: break if good: return True return False
[docs]def graph_count_automorphisms(G,vertex_orbits=False): return count_automorphisms(G.M,G.vertex_groupings,vertex_orbits)
[docs]def count_automorphisms(M,grouping,vertex_orbits=False): nr,nc = M.nrows(),M.ncols() count = _sage_const_0 PermList = [Permutations(range(len(group))) for group in grouping] if vertex_orbits: isom_list = [] for sigma_data in itertools.product(*PermList): sigma = [_sage_const_0 for i in range(nr-_sage_const_1 )] for i in range(len(grouping)): for j in range(len(grouping[i])): sigma[grouping[i][j]-_sage_const_1 ] = grouping[i][sigma_data[i][j]] good = True for i in range(_sage_const_1 ,nr): ii = sigma[i-_sage_const_1 ] for j in range(_sage_const_1 ,i): jj = sigma[j-_sage_const_1 ] L1 = [] for k in range(_sage_const_1 ,nc): if M[i,k] != _sage_const_0 and M[j,k] != _sage_const_0 : L1.append([M[i,k],M[j,k]]) L1.sort() L2 = [] for k in range(_sage_const_1 ,nc): if M[ii,k] != _sage_const_0 and M[jj,k] != _sage_const_0 : L2.append([M[ii,k],M[jj,k]]) L2.sort() if L1 != L2: good = False break if good == False: break if good: count += _sage_const_1 if vertex_orbits: isom_list.append(sigma) if vertex_orbits: orbit_list = [] vertices_used = [] while len(vertices_used) < nr-_sage_const_1 : i = [ii for ii in range(_sage_const_1 ,nr) if ii not in vertices_used][_sage_const_0 ] orbit = [] for sigma in isom_list: if sigma[i-_sage_const_1 ] not in orbit: orbit.append(sigma[i-_sage_const_1 ]) vertices_used.append(sigma[i-_sage_const_1 ]) orbit.sort() orbit_list.append(orbit) return orbit_list for i in range(_sage_const_1 ,nr): for k in range(_sage_const_1 ,nc): if M[i,k][_sage_const_0 ] == _sage_const_2 and M[i,k][_sage_const_1 ] == M[i,k][_sage_const_2 ]: count *= _sage_const_2 L = [] for k in range(_sage_const_1 ,nc): if M[i,k] != _sage_const_0 : if sum(_sage_const_1 for j in range(_sage_const_1 ,nr) if M[j,k] != _sage_const_0 ) == _sage_const_1 : L.append([M[_sage_const_0 ,k],M[i,k]]) count *= aut(L) for j in range(_sage_const_1 ,i): L = [] for k in range(_sage_const_1 ,nc): if M[i,k] != _sage_const_0 and M[j,k] != _sage_const_0 : L.append([M[i,k],M[j,k]]) count *= aut(L) return count
[docs]def graph_list_isomorphisms(G1,G2,only_one=False): if G1.invariant != G2.invariant: return [] else: return list_isomorphisms(G1.M,G2.M,G1.vertex_groupings,G2.vertex_groupings,only_one)
[docs]def list_isomorphisms(M1,M2,group1,group2,only_one=False): # Warning: does not count loops! # If this is too slow, we can probably improve by caching a list of automorphisms and applying those to the first isom found. nr,nc = M1.nrows(),M2.ncols() PermList = [Permutations(range(len(group))) for group in group1] isom_list = [] for sigma_data in itertools.product(*PermList): sigma = [_sage_const_0 for i in range(nr-_sage_const_1 )] for i in range(len(group1)): for j in range(len(group1[i])): sigma[group1[i][j]-_sage_const_1 ] = group2[i][sigma_data[i][j]] good = True for i in range(_sage_const_1 ,nr): ii = sigma[i-_sage_const_1 ] for j in range(_sage_const_1 ,i): jj = sigma[j-_sage_const_1 ] L1 = [] for k in range(_sage_const_1 ,nc): if M1[i,k] != _sage_const_0 and M1[j,k] != _sage_const_0 : L1.append([M1[i,k],M1[j,k]]) L1.sort() L2 = [] for k in range(_sage_const_1 ,nc): if M2[ii,k] != _sage_const_0 and M2[jj,k] != _sage_const_0 : L2.append([M2[ii,k],M2[jj,k]]) L2.sort() if L1 != L2: good = False break if good == False: break if good: cols1 = [[M1[i,j] for i in range(nr)] for j in range(_sage_const_1 ,nc)] cols2 = [[M2[_sage_const_0 ,j]] + [M2[sigma[i-_sage_const_1 ],j] for i in range(_sage_const_1 ,nr)] for j in range(_sage_const_1 ,nc)] edge_group1 = [] edge_group2 = [] used1 = [] for j in range(_sage_const_1 ,nc): if j not in used1: edge_group1.append([]) edge_group2.append([]) for k in range(_sage_const_1 ,nc): if cols1[k-_sage_const_1 ] == cols1[j-_sage_const_1 ]: edge_group1[-_sage_const_1 ].append(k) used1.append(k) if cols2[k-_sage_const_1 ] == cols1[j-_sage_const_1 ]: edge_group2[-_sage_const_1 ].append(k) edge_PermList = [Permutations(range(len(edge_group))) for edge_group in edge_group1] for edge_sigma_data in itertools.product(*edge_PermList): edge_sigma = [_sage_const_0 for i in range(nc-_sage_const_1 )] for i in range(len(edge_group1)): for j in range(len(edge_group1[i])): edge_sigma[edge_group1[i][j]-_sage_const_1 ] = edge_group2[i][edge_sigma_data[i][j]] isom_list.append([sigma,edge_sigma]) if only_one: return isom_list return isom_list
[docs]def aut(L): if len(L) == _sage_const_0 : return _sage_const_1 L.sort() total = _sage_const_1 n = _sage_const_1 last = L[_sage_const_0 ] for l in L[_sage_const_1 :]: if l == last: n += _sage_const_1 else: n = _sage_const_1 total *= n last = l return total
[docs]def degenerate(G_list,moduli_type=MODULI_ST): mod_size = moduli_type + _sage_const_1 if moduli_type == MODULI_SMALL: mod_size = MODULI_SM + _sage_const_1 G_list_new = [[] for i in range(mod_size)] for which_type in range(mod_size): for G in G_list[which_type]: for i in range(_sage_const_1 ,G.num_vertices()+_sage_const_1 ): row = list(G.M[i]) m = row[_sage_const_0 ] + sum(row) if m < _sage_const_4 : continue row1 = [_sage_const_0 for j in range(len(row))] while [_sage_const_2 *x for x in row1] <= row: if row1[_sage_const_0 ] == _sage_const_1 and moduli_type <= MODULI_RT: break if row1[_sage_const_0 ] + sum(row1) >= _sage_const_2 and row1[_sage_const_0 ] + sum(row1) <= m-_sage_const_2 : row2 = [row[j] - row1[j] for j in range(len(row))] G_copy = Graph(G.M) G_copy.split_vertex(i,row1,row2) new_type = which_type if new_type == MODULI_SM: new_type = MODULI_RT if new_type == MODULI_RT and row1[_sage_const_0 ] > _sage_const_0 : new_type = MODULI_CT G_list_new[new_type].append(G_copy) row1[-_sage_const_1 ] += _sage_const_1 for j in range(_sage_const_1 ,len(row)): if row1[-j] <= row[-j]: break row1[-j] = _sage_const_0 row1[-j-_sage_const_1 ] += _sage_const_1 for i in range(mod_size): G_list_new[i] = remove_isomorphic(G_list_new[i]) return G_list_new
[docs]def dim_form(g,n,moduli_type=MODULI_ST): if moduli_type == MODULI_ST: return _sage_const_3 *g-_sage_const_3 +n if moduli_type == MODULI_CT: return _sage_const_2 *g-_sage_const_3 +n if moduli_type == MODULI_RT: if g > _sage_const_0 : return g-_sage_const_2 +n else: return n-_sage_const_3 if moduli_type == MODULI_SM: if n == _sage_const_0 : return g-_sage_const_2 elif g >= _sage_const_1 : return g-_sage_const_1 else: return _sage_const_0 if moduli_type == MODULI_SMALL: return _sage_const_1000 return _sage_const_3 *g-_sage_const_3 +n
[docs]def decorate(G_list,r,moduli_type=MODULI_ST): mod_size = moduli_type + _sage_const_1 if moduli_type == MODULI_SMALL: mod_size = MODULI_SM + _sage_const_1 G_list_new = [[] for i in range(mod_size)] for which_type in range(mod_size): for G in G_list[which_type]: G_deco = [[] for i in range(mod_size)] G.compute_degree_vec() nr,nc = G.M.nrows(),G.M.ncols() two_list = [] one_list = [] for i in range(_sage_const_1 ,nr): for j in range(_sage_const_1 ,nc): if G.M[i,j] == _sage_const_2 : two_list.append([i,j]) elif G.M[i,j] == _sage_const_1 : one_list.append([i,j]) a = nr-_sage_const_1 b = len(two_list) c = len(one_list) dims = [[dim_form(G.M[i+_sage_const_1 ,_sage_const_0 ][_sage_const_0 ], G.degree(i+_sage_const_1 ), mod_type) for i in range(a)] for mod_type in range(mod_size)] for vec in IntegerVectors(r,a+b+c): new_type = which_type if moduli_type > MODULI_SMALL: test_dims = vec[:a] for i in range(b): test_dims[two_list[i][_sage_const_0 ]-_sage_const_1 ] += vec[a+i] for i in range(c): test_dims[one_list[i][_sage_const_0 ]-_sage_const_1 ] += vec[a+b+i] for mod_type in range(which_type,mod_size): for i in range(a): if test_dims[i] > dims[mod_type][i]: new_type = mod_type + _sage_const_1 break if new_type > moduli_type: continue S_list = [] for i in range(a): S_list.append(Partitions(vec[i])) for i in range(a,a+b): S_list.append([[vec[i]-j,j] for j in range(vec[i]/_sage_const_2 + _sage_const_1 )]) S = itertools.product(*S_list) for vec2 in S: G_copy = Graph(G.M) for i in range(a): for j in vec2[i]: G_copy.M[i+_sage_const_1 ,_sage_const_0 ] += X**j for i in range(a,a+b): G_copy.M[two_list[i-a][_sage_const_0 ],two_list[i-a][_sage_const_1 ]] += vec2[i][_sage_const_0 ]*X + vec2[i][_sage_const_1 ]*X**_sage_const_2 for i in range(c): G_copy.M[one_list[i][_sage_const_0 ],one_list[i][_sage_const_1 ]] += vec[i+a+b]*X G_deco[new_type].append(G_copy) for mod_type in range(mod_size): G_list_new[mod_type] += remove_isomorphic(G_deco[mod_type]) return G_list_new
[docs]def remove_isomorphic(G_list): G_list_new = [] inv_dict = {} count = _sage_const_0 for G1 in G_list: G1.compute_invariant() if not inv_dict.has_key(G1.invariant): inv_dict[G1.invariant] = [] good = True for i in inv_dict[G1.invariant]: if graph_isomorphic(G1,G_list_new[i]): good = False break if good: G_list_new.append(G1) inv_dict[G1.invariant].append(count) count += _sage_const_1 return G_list_new
[docs]def num_strata(g,r,markings=(),moduli_type=MODULI_ST): return len(capply(all_strata,g,r,markings,moduli_type))
[docs]def num_pure_strata(g,r,markings=(),moduli_type=MODULI_ST): return len(capply(all_pure_strata,g,r,markings,moduli_type))
[docs]def single_stratum(num,g,r,markings=(),moduli_type=MODULI_ST): return capply(all_strata,g,r,markings,moduli_type)[num]
[docs]def single_pure_stratum(num,g,r,markings=(),moduli_type=MODULI_ST): return capply(all_pure_strata,g,r,markings,moduli_type)[num]
[docs]def autom_count(num,g,r,markings=(),moduli_type=MODULI_ST): return graph_count_automorphisms(single_stratum(num,g,r,markings,moduli_type))
[docs]def pure_strata_autom_count(num,g,r,markings=(),moduli_type=MODULI_ST): return graph_count_automorphisms(single_pure_stratum(num,g,r,markings,moduli_type))
[docs]def unpurify_map(g,r,markings=(),moduli_type=MODULI_ST): unpurify = {} pure_strata = [capply(all_pure_strata,g,r0,markings,moduli_type) for r0 in range(r+_sage_const_1 )] impure_strata = capply(all_strata,g,r,markings,moduli_type) for i in range(len(impure_strata)): G = Graph(impure_strata[i].M) G.purify() r0 = G.num_edges() - len(markings) found = False for j in range(len(pure_strata[r0])): if G.M == pure_strata[r0][j].M: G_key = (r0, j) found = True break if not found: print "ERROR! Purification failed." if not unpurify.has_key(G_key): unpurify[G_key] = [] unpurify[G_key].append(i) return unpurify
[docs]def all_strata(g,r,markings=(),moduli_type=MODULI_ST): mod_size = moduli_type + _sage_const_1 if moduli_type == MODULI_SMALL: mod_size = MODULI_SM + _sage_const_1 big_list = [[] for i in range(mod_size)] for loops in range(g+_sage_const_1 ): if loops == _sage_const_1 and moduli_type <= MODULI_CT: break if loops > r: break for edges in range(r-loops+_sage_const_1 ): if edges == _sage_const_1 and moduli_type <= MODULI_SM: break G = Graph() G.add_vertex(g-loops) for k in range(loops): G.add_edge(_sage_const_1 ,_sage_const_1 ) for k in markings: G.add_edge(_sage_const_1 ,_sage_const_0 ,k) GGG = [[] for i in range(mod_size)] if loops == _sage_const_0 : if edges == _sage_const_0 : GGG[MODULI_SM] = [G] else: GGG[MODULI_RT] = [G] else: GGG[MODULI_ST] = [G] for k in range(edges): GGG = degenerate(GGG,moduli_type) GGG = decorate(GGG,r-loops-edges,moduli_type) for i in range(mod_size): big_list[i] += GGG[i] combined_list = [] for i in range(mod_size): combined_list += big_list[i] for G in combined_list: G.compute_degree_vec() G.set_target_parity() return combined_list
[docs]def all_pure_strata(g,r,markings=(),moduli_type=MODULI_ST): big_list = [[] for i in range(moduli_type+_sage_const_1 )] for loops in range(g+_sage_const_1 ): if loops == _sage_const_1 and moduli_type <= MODULI_CT: break if loops > r: break for edges in range(r-loops,r-loops+_sage_const_1 ): if edges >= _sage_const_1 and moduli_type <= MODULI_SM: break G = Graph() G.add_vertex(g-loops) for k in range(loops): G.add_edge(_sage_const_1 ,_sage_const_1 ) for k in markings: G.add_edge(_sage_const_1 ,_sage_const_0 ,k) G.compute_invariant() GGG = [[] for i in range(moduli_type+_sage_const_1 )] if loops == _sage_const_0 : if edges == _sage_const_0 : GGG[MODULI_SM] = [G] else: GGG[MODULI_RT] = [G] else: GGG[MODULI_ST] = [G] for k in range(edges): GGG = degenerate(GGG,moduli_type) for i in range(moduli_type+_sage_const_1 ): big_list[i] += GGG[i] combined_list = [] for i in range(moduli_type+_sage_const_1 ): combined_list += big_list[i] return combined_list
############################
[docs]def C_coeff(m,term): n = term - floor(m/_sage_const_3 ) if n < _sage_const_0 : return _sage_const_0 if (m % _sage_const_3 ) == _sage_const_0 : return A_list[n] else: return B_list[n]
[docs]def dual_C_coeff(i,j,parity): total = _sage_const_0 k = parity % _sage_const_2 while (floor(k/_sage_const_3 ) <= i): if (k % _sage_const_3 ) == _sage_const_2 : k += _sage_const_2 continue total += (-_sage_const_1 )**(floor(k/_sage_const_3 ))*C_coeff(k,i)*C_coeff(-_sage_const_2 -k,j) k += _sage_const_2 return total
[docs]def poly_to_partition(F): mmm = F.degree() target_partition = [] for i in range(_sage_const_1 ,mmm+_sage_const_1 ): for j in range(F[i]): target_partition.append(i) return tuple(target_partition)
[docs]def kappa_coeff(sigma,kappa_0,target_partition): total = _sage_const_0 num_ones = sum(_sage_const_1 for i in sigma if i == _sage_const_1 ) for i in range(_sage_const_0 ,num_ones+_sage_const_1 ): for injection in Permutations(range(len(target_partition)),len(sigma)-i): term = binomial(num_ones,i)*binomial(kappa_0 + len(target_partition) + i-_sage_const_1 , i)*factorial(i) for j in range(len(sigma)-i): term *= C_coeff(sigma[j+i],target_partition[injection[j]]) for j in range(len(target_partition)): if j in injection: continue term *= C_coeff(_sage_const_0 ,target_partition[j]) total += term total = (-_sage_const_1 )**(len(target_partition)+len(sigma))*total/aut(list(target_partition)) return total
[docs]def FZ_kappa_factor(num,sigma,g,r,markings=(),moduli_type=MODULI_ST): G = single_stratum(num,g,r,markings,moduli_type) L = [] nv = G.num_vertices() for i in range(_sage_const_1 ,nv+_sage_const_1 ): L.append((_sage_const_2 *G.M[i,_sage_const_0 ][_sage_const_0 ]+G.degree(i)-_sage_const_2 ,G.M[i,_sage_const_0 ]-G.M[i,_sage_const_0 ][_sage_const_0 ])) LL = [] tau = [] for i in range(nv): min = -_sage_const_1 for j in range(nv): if (i == _sage_const_0 or L[j] > LL[-_sage_const_1 ] or L[j] == LL[-_sage_const_1 ] and j > tau[-_sage_const_1 ]) and (min == -_sage_const_1 or L[j] < L[min]): min = j tau.append(min) LL.append(L[min]) factor_dict = capply(FZ_kappa_factor2,tuple(LL),sigma) factor_vec = [_sage_const_0 for i in range(_sage_const_1 << nv)] for parity_key in factor_dict.keys(): parity = _sage_const_0 for i in range(nv): if parity_key[i] == _sage_const_1 : parity += _sage_const_1 << tau[i] factor_vec[parity] = factor_dict[parity_key] return factor_vec
[docs]def FZ_marking_factor(num,marking_vec,g,r,markings=(),moduli_type=MODULI_ST): G = single_stratum(num,g,r,markings,moduli_type) nv = G.num_vertices() ne = G.num_edges() num_parities = _sage_const_2 **nv PPP_list = [] for marks in marking_vec: PPP_list.append(Permutations(marks[_sage_const_1 ])) PPP = itertools.product(*PPP_list) marking_factors = [_sage_const_0 for i in range(num_parities)] incident_vertices = [] for mark_type in marking_vec: incident_vertices.append([]) for k in range(_sage_const_1 ,ne+_sage_const_1 ): if G.M[_sage_const_0 ,k] == mark_type[_sage_const_0 ]: for i in range(_sage_const_1 ,nv+_sage_const_1 ): if G.M[i,k] != _sage_const_0 : incident_vertices[-_sage_const_1 ].append((i-_sage_const_1 ,G.M[i,k][_sage_const_1 ])) break for perms in PPP: parity = _sage_const_0 marking_factor = _sage_const_1 for marks_index in range(len(marking_vec)): for count in range(len(incident_vertices[marks_index])): marking_factor *= C_coeff(perms[marks_index][count],incident_vertices[marks_index][count][_sage_const_1 ]) parity ^= (perms[marks_index][count] % _sage_const_2 ) << incident_vertices[marks_index][count][_sage_const_0 ] marking_factors[parity] += marking_factor return marking_factors
[docs]def FZ_kappa_factor2(L,sigma): nv = len(L) mmm = max((_sage_const_0 ,)+sigma) sigma_grouped = [_sage_const_0 for i in range(mmm)] for i in sigma: sigma_grouped[i-_sage_const_1 ] += _sage_const_1 S_list = [] for i in sigma_grouped: S_list.append(IntegerVectors(i,nv)) S = itertools.product(*S_list) kappa_factors = {} for parity in itertools.product(*[(_sage_const_0 ,_sage_const_1 ) for i in range(nv)]): kappa_factors[tuple(parity)] = _sage_const_0 for assignment in S: assigned_sigma = [[] for j in range(nv)] for i in range(mmm): for j in range(nv): for k in range(assignment[i][j]): assigned_sigma[j].append(i+_sage_const_1 ) sigma_auts = _sage_const_1 parity = [_sage_const_0 for i in range(nv)] kappa_factor = _sage_const_1 for j in range(nv): sigma_auts *= aut(assigned_sigma[j]) parity[j] += sum(assigned_sigma[j]) parity[j] %= _sage_const_2 kappa_factor *= capply(kappa_coeff,tuple(assigned_sigma[j]),L[j][_sage_const_0 ],poly_to_partition(L[j][_sage_const_1 ])) kappa_factors[tuple(parity)] += kappa_factor/sigma_auts return kappa_factors
[docs]def FZ_hedge_factor(num,g,r,markings=(),moduli_type=MODULI_ST): G = single_stratum(num,g,r,markings,moduli_type) nv = G.num_vertices() num_parities = _sage_const_2 **nv ne = G.num_edges() edge_list = [] for k in range(_sage_const_1 ,ne+_sage_const_1 ): if G.M[_sage_const_0 ,k] == _sage_const_0 : edge_list.append([k]) for i in range(_sage_const_1 ,nv+_sage_const_1 ): if G.M[i,k] != _sage_const_0 : edge_list[-_sage_const_1 ].append(i) if G.M[i,k][_sage_const_0 ] == _sage_const_2 : edge_list[-_sage_const_1 ].append(i) hedge_factors = [_sage_const_0 for i in range(num_parities)] for edge_parities in itertools.product(*[[_sage_const_0 ,_sage_const_1 ] for i in edge_list]): parity = _sage_const_0 for i in range(len(edge_list)): if edge_parities[i] == _sage_const_1 : parity ^= _sage_const_1 << (edge_list[i][_sage_const_1 ]-_sage_const_1 ) parity ^= _sage_const_1 << (edge_list[i][_sage_const_2 ]-_sage_const_1 ) hedge_factor = _sage_const_1 for i in range(len(edge_list)): if edge_list[i][_sage_const_1 ] == edge_list[i][_sage_const_2 ]: hedge_factor *= capply(dual_C_coeff,G.M[edge_list[i][_sage_const_1 ],edge_list[i][_sage_const_0 ]][_sage_const_1 ],G.M[edge_list[i][_sage_const_1 ],edge_list[i][_sage_const_0 ]][_sage_const_2 ],edge_parities[i] % _sage_const_2 ) else: hedge_factor *= capply(dual_C_coeff,G.M[edge_list[i][_sage_const_1 ],edge_list[i][_sage_const_0 ]][_sage_const_1 ],G.M[edge_list[i][_sage_const_2 ],edge_list[i][_sage_const_0 ]][_sage_const_1 ],edge_parities[i] % _sage_const_2 ) hedge_factors[parity] += hedge_factor return hedge_factors
[docs]def FZ_coeff(num,FZ_param,g,r,markings=(),moduli_type=MODULI_ST): sigma = FZ_param[_sage_const_0 ] marking_vec = FZ_param[_sage_const_1 ] G = single_stratum(num,g,r,markings,moduli_type) nv = G.num_vertices() graph_auts = capply(autom_count,num,g,r,markings,moduli_type) h1_factor = _sage_const_2 **G.h1() num_parities = _sage_const_2 **nv marking_factors = capply(FZ_marking_factor,num,marking_vec,g,r,markings,moduli_type) kappa_factors = capply(FZ_kappa_factor,num,sigma,g,r,markings,moduli_type) hedge_factors = capply(FZ_hedge_factor,num,g,r,markings,moduli_type) total = _sage_const_0 for i in range(num_parities): if marking_factors[i] == _sage_const_0 : continue for j in range(num_parities): total += marking_factors[i]*kappa_factors[j]*hedge_factors[i ^ j ^ G.target_parity] total /= h1_factor*graph_auts return total
[docs]def interior_FZ(g,r,markings=(),moduli_type=MODULI_ST): ngen = num_strata(g,r,markings,moduli_type) #print "%s generators" % ngen relations = [] FZpl = FZ_param_list(_sage_const_3 *r-g-_sage_const_1 ,markings) #print "%s codim 0 relations to compute" % len(FZpl) ccccc = _sage_const_0 for FZ_param in FZpl: # if ccccc % 5 == 0: # print "%s done" % ccccc ccccc += _sage_const_1 relation = [capply(FZ_coeff,i,FZ_param,g,r,markings,moduli_type) for i in range(ngen)] relations.append(relation) return relations
[docs]def possibly_new_FZ(g,r,n=_sage_const_0 ,moduli_type=MODULI_ST): m = _sage_const_3 *r-g-_sage_const_1 -n if m < _sage_const_0 : return [] dprint("Start FZ (%s,%s,%s,%s): %s",g,r,n,moduli_type,floor(get_memory_usage())) markings = tuple([_sage_const_1 for i in range(n)]) ngen = num_strata(g,r,markings,moduli_type) relations = [] for i in range(m+_sage_const_1 ): if m-i % _sage_const_2 == _sage_const_1 : continue for sigma in Partitions(i): if len([j for j in sigma if j%_sage_const_3 != _sage_const_1 ]) > _sage_const_0 : continue if n > _sage_const_0 : FZ_param = (tuple(sigma), ((_sage_const_1 , markings),)) else: FZ_param = (tuple(sigma), ()) relation = [] for j in range(ngen): coeff = capply(FZ_coeff,j,FZ_param,g,r,markings,moduli_type) if coeff != _sage_const_0 : relation.append([j,coeff]) relations.append(relation) dprint("End FZ (%s,%s,%s,%s): %s",g,r,n,moduli_type,floor(get_memory_usage())) return relations
[docs]def boundary_FZ(g,r,markings=(),moduli_type=MODULI_ST): if moduli_type <= MODULI_SM: return [] generators = capply(all_strata,g,r,markings,moduli_type) ngen = len(generators) #print "%s generators" % ngen relations = [] old_count = _sage_const_0 for r0 in range(_sage_const_1 ,r): strata = capply(all_strata,g,r0,markings,moduli_type) for G in strata: vertex_orbits = graph_count_automorphisms(G,True) for i in [orbit[_sage_const_0 ] for orbit in vertex_orbits]: good = True for j in range(G.M.ncols()): if R(G.M[i,j][_sage_const_0 ]) != G.M[i,j]: good = False break if good: g2 = G.M[i,_sage_const_0 ][_sage_const_0 ] if _sage_const_3 *(r-r0) < g2 + _sage_const_1 : continue d = G.degree(i) if dim_form(g2,d,moduli_type) < r-r0: continue strata2 = capply(all_strata,g2,r-r0,tuple(range(_sage_const_1 ,d+_sage_const_1 )),moduli_type) which_gen_list = [-_sage_const_1 for num in range(len(strata2))] for num in range(len(strata2)): G_copy = Graph(G.M) G_copy.replace_vertex_with_graph(i,strata2[num]) which_gen_list[num] = num_of_stratum(G_copy,g,r,markings,moduli_type) rFZpl = reduced_FZ_param_list(G,i,g2,d,_sage_const_3 *(r-r0)-g2-_sage_const_1 ) #print "Computing %s relations to insert at vertex %s into" % (len(rFZpl), i) #print G.M ccccc = _sage_const_0 for FZ_param in rFZpl: relation = [_sage_const_0 for k in range(ngen)] for num in range(len(strata2)): if which_gen_list[num] != -_sage_const_1 : relation[which_gen_list[num]] += capply(FZ_coeff,num,FZ_param,g2,r-r0,tuple(range(_sage_const_1 ,d+_sage_const_1 )),moduli_type) relations.append(relation) ccccc += _sage_const_1 # if ccccc % 5 == 0: # print "%s done" % ccccc return relations
[docs]def list_all_FZ(g,r,markings=(),moduli_type=MODULI_ST): relations = copy(capply(interior_FZ,g,r,markings,moduli_type)) if moduli_type > MODULI_SM: relations += capply(boundary_FZ,g,r,markings,moduli_type) if len(relations) == _sage_const_0 : ngen = num_strata(g,r,markings,moduli_type) relations.append([_sage_const_0 for i in range(ngen)]) return relations
[docs]def reduced_FZ_param_list(G,v,g,d,n): params = FZ_param_list(n,tuple(range(_sage_const_1 ,d+_sage_const_1 ))) graph_params = [] M = matrix(R,_sage_const_2 ,d+_sage_const_1 ) M[_sage_const_0 ,_sage_const_0 ] = -_sage_const_1 for i in range(_sage_const_1 ,d+_sage_const_1 ): M[_sage_const_0 ,i] = i for p in params: G_copy = Graph(G.M) M[_sage_const_1 ,_sage_const_0 ] = -g-_sage_const_1 for j in p[_sage_const_0 ]: M[_sage_const_1 ,_sage_const_0 ] += X**j for i in range(_sage_const_1 ,d+_sage_const_1 ): M[_sage_const_1 ,i] = _sage_const_1 + p[_sage_const_1 ][i-_sage_const_1 ][_sage_const_1 ][_sage_const_0 ]*X G_p = Graph(M) G_copy.replace_vertex_with_graph(v,G_p) graph_params.append([p,G_copy]) params_reduced = [] graphs_seen = [] for x in graph_params: x[_sage_const_1 ].compute_invariant() good = True for GG in graphs_seen: if graph_isomorphic(x[_sage_const_1 ],GG): good = False break if good: graphs_seen.append(x[_sage_const_1 ]) params_reduced.append(x[_sage_const_0 ]) return params_reduced
[docs]def FZ_param_list(n,markings=()): if n < _sage_const_0 : return [] final_list = [] mmm = max((_sage_const_0 ,)+markings) markings_grouped = [_sage_const_0 for i in range(mmm)] for i in markings: markings_grouped[i-_sage_const_1 ] += _sage_const_1 markings_best = [] for i in range(mmm): if markings_grouped[i] > _sage_const_0 : markings_best.append([i+_sage_const_1 ,markings_grouped[i]]) for j in range(n/_sage_const_2 + _sage_const_1 ): for n_vec in IntegerVectors(n-_sage_const_2 *j,_sage_const_1 +len(markings_best)): S_list = [[list(sigma) for sigma in Partitions(n_vec[_sage_const_0 ]).list() if sum(_sage_const_1 for l in sigma if (l % _sage_const_3 ) == _sage_const_2 ) == _sage_const_0 ]] for i in range(len(markings_best)): S_list.append(Partitions(n_vec[i+_sage_const_1 ]+markings_best[i][_sage_const_1 ],length=markings_best[i][_sage_const_1 ]).list()) S_list[-_sage_const_1 ] = [[k - _sage_const_1 for k in sigma] for sigma in S_list[-_sage_const_1 ] if sum(_sage_const_1 for l in sigma if (l % _sage_const_3 ) == _sage_const_0 ) == _sage_const_0 ] for S in itertools.product(*S_list): final_list.append((tuple(S[_sage_const_0 ]),tuple([(markings_best[k][_sage_const_0 ],tuple(S[k+_sage_const_1 ])) for k in range(len(markings_best))]))) return final_list
[docs]def FZ_matrix(g,r,markings=(),moduli_type=MODULI_ST): return matrix(list_all_FZ(g,r,markings,moduli_type))
#####################################
[docs]def contraction_table(g,r,markings=(),moduli_type=MODULI_ST): contraction_dict = {} pure_strata = [capply(all_pure_strata,g,r0,markings,moduli_type) for r0 in range(r+_sage_const_1 )] for r0 in range(r+_sage_const_1 ): for ii in range(len(pure_strata[r0])): G = pure_strata[r0][ii] S = [j for j in range(_sage_const_1 ,G.M.ncols()) if G.M[_sage_const_0 ,j] == _sage_const_0 ] contractions = {} for edge_subset in itertools.product(*[[_sage_const_0 ,_sage_const_1 ] for i in range(r0)]): key = tuple(i for i in range(r0) if edge_subset[i] == _sage_const_0 ) A = [S[i] for i in key] A.reverse() vlist = [[i] for i in range(_sage_const_1 ,G.M.nrows())] elist = [i for i in range(_sage_const_1 ,G.M.ncols())] Gcopy = Graph(G.M) for i in A: vlist,elist = Gcopy.contract(i,vlist,elist) Gcopy.compute_invariant() rnew = r0 - len(A) contraction_result = [] for i in range(len(pure_strata[rnew])): L = graph_list_isomorphisms(pure_strata[rnew][i], Gcopy, True) if len(L) > _sage_const_0 : contraction_result.append((rnew,i)) contraction_result.append(L[_sage_const_0 ]) break contraction_result.append((vlist,elist)) contractions[key] = contraction_result for edge_assignment in itertools.product(*[[_sage_const_0 ,_sage_const_1 ,_sage_const_2 ] for i in range(r0)]): if sum(_sage_const_1 for i in edge_assignment if i == _sage_const_1 ) > r-r0: continue key1 = tuple(i for i in range(r0) if edge_assignment[i] == _sage_const_0 ) B = [S[i] for i in range(r0) if edge_assignment[i] == _sage_const_1 ] key2 = tuple(i for i in range(r0) if edge_assignment[i] == _sage_const_2 ) if key1 > key2: continue contract1 = contractions[key1] contract2 = contractions[key2] dict_key = [contract1[_sage_const_0 ],contract2[_sage_const_0 ]] dict_entry = [contract1[_sage_const_1 ],contract2[_sage_const_1 ]] if dict_key[_sage_const_0 ] > dict_key[_sage_const_1 ]: dict_key.reverse() dict_entry.reverse() dict_entry = [(r0,ii),B,contract2[_sage_const_2 ],contract1[_sage_const_2 ]] + dict_entry else: dict_entry = [(r0,ii),B,contract1[_sage_const_2 ],contract2[_sage_const_2 ]] + dict_entry dict_key = tuple(dict_key) if not contraction_dict.has_key(dict_key): contraction_dict[dict_key] = [] contraction_dict[dict_key].append(dict_entry) if dict_key[_sage_const_0 ] == dict_key[_sage_const_1 ]: contraction_dict[dict_key].append(dict_entry[:_sage_const_2 ] + [dict_entry[_sage_const_3 ],dict_entry[_sage_const_2 ],dict_entry[_sage_const_5 ],dict_entry[_sage_const_4 ]]) return contraction_dict
[docs]def multiply(r1,i1,r2,i2,g,rmax,markings=(),moduli_type=MODULI_ST): unpurify = capply(unpurify_map,g,r1+r2,markings,moduli_type) gens = capply(all_strata,g,r1+r2,markings,moduli_type) ngens = num_strata(g,r1+r2,markings,moduli_type) answer = [_sage_const_0 for i in range(ngens)] pure_strata = [capply(all_pure_strata,g,r,markings,moduli_type) for r in range(rmax+_sage_const_1 )] contraction_dict = capply(contraction_table,g,rmax,markings,moduli_type) G1 = single_stratum(i1,g,r1,markings,moduli_type) G2 = single_stratum(i2,g,r2,markings,moduli_type) G1copy = Graph(G1.M) G2copy = Graph(G2.M) G1copy.purify() G2copy.purify() pure_r1 = G1copy.num_edges() - len(markings) pure_r2 = G2copy.num_edges() - len(markings) found = False for i in range(len(pure_strata[pure_r1])): if G1copy.M == pure_strata[pure_r1][i].M: G1_key = (pure_r1, i) found = True break if not found: print "ERROR! Purification failed." found = False for i in range(len(pure_strata[pure_r2])): if G2copy.M == pure_strata[pure_r2][i].M: G2_key = (pure_r2, i) found = True break if not found: print "ERROR! Purification failed." if G1_key > G2_key: return capply(multiply,r2,i2,r1,i1,g,rmax,markings,moduli_type) if not contraction_dict.has_key((G1_key,G2_key)): return answer for L in contraction_dict[(G1_key,G2_key)]: H = pure_strata[L[_sage_const_0 ][_sage_const_0 ]][L[_sage_const_0 ][_sage_const_1 ]] Hloops = [] if moduli_type > MODULI_CT: for i in range(_sage_const_1 ,H.M.nrows()): for j in range(_sage_const_1 ,H.M.ncols()): if H.M[i,j][_sage_const_0 ] == _sage_const_2 : Hloops.append((i,j)) auts = capply(pure_strata_autom_count,L[_sage_const_0 ][_sage_const_1 ],g,L[_sage_const_0 ][_sage_const_0 ],markings,moduli_type) B = L[_sage_const_1 ] if len(B) == pure_r1 and len(B) == pure_r2: auts *= _sage_const_2 aut_cosets1 = capply(automorphism_cosets,i1,g,r1,markings,moduli_type) aut_cosets2 = capply(automorphism_cosets,i2,g,r2,markings,moduli_type) auts /= aut_cosets1[_sage_const_0 ]*aut_cosets2[_sage_const_0 ] for isom1 in aut_cosets1[_sage_const_1 ]: for isom2 in aut_cosets2[_sage_const_1 ]: Hcopy = Graph(H.M) vmap1 = [_sage_const_0 for i in range(G1.M.nrows())] for i in range(_sage_const_1 ,G1.M.nrows()): vmap1[i] = L[_sage_const_2 ][_sage_const_0 ][L[_sage_const_4 ][_sage_const_0 ][isom1[_sage_const_0 ][i-_sage_const_1 ]-_sage_const_1 ]-_sage_const_1 ] emap1 = [_sage_const_0 for i in range(G1.M.ncols())] for i in range(_sage_const_1 ,G1.M.ncols()): emap1[i] = L[_sage_const_2 ][_sage_const_1 ][L[_sage_const_4 ][_sage_const_1 ][isom1[_sage_const_1 ][i-_sage_const_1 ]-_sage_const_1 ]-_sage_const_1 ] vmap2 = [_sage_const_0 for i in range(G2.M.nrows())] for i in range(_sage_const_1 ,G2.M.nrows()): vmap2[i] = L[_sage_const_3 ][_sage_const_0 ][L[_sage_const_5 ][_sage_const_0 ][isom2[_sage_const_0 ][i-_sage_const_1 ]-_sage_const_1 ]-_sage_const_1 ] emap2 = [_sage_const_0 for i in range(G2.M.ncols())] for i in range(_sage_const_1 ,G2.M.ncols()): emap2[i] = L[_sage_const_3 ][_sage_const_1 ][L[_sage_const_5 ][_sage_const_1 ][isom2[_sage_const_1 ][i-_sage_const_1 ]-_sage_const_1 ]-_sage_const_1 ] psilooplist = [] psiindexlist = [] loop_factor = _sage_const_1 for i in range(_sage_const_1 ,G1.M.nrows()): for j in range(_sage_const_1 ,G1.M.ncols()): if G1.M[i,j][_sage_const_0 ] != _sage_const_0 : if G1.M[i,j][_sage_const_0 ] == _sage_const_1 : if G1.M[i,j][_sage_const_1 ] != _sage_const_0 : jj = emap1[j] for ii in vmap1[i]: if H.M[ii,jj] != _sage_const_0 : Hcopy.M[ii,jj] += G1.M[i,j][_sage_const_1 ]*X break elif G1.M[i,j][_sage_const_1 ] == _sage_const_0 : loop_factor *= _sage_const_2 else: jj = emap1[j] psilooplist.append([[G1.M[i,j][_sage_const_1 ],G1.M[i,j][_sage_const_2 ]],[G1.M[i,j][_sage_const_2 ],G1.M[i,j][_sage_const_1 ]]]) psiindexlist.append([jj]) for ii in vmap1[i]: for k in range(H.M[ii,jj][_sage_const_0 ]): psiindexlist[-_sage_const_1 ].append(ii) for i in range(_sage_const_1 ,G2.M.nrows()): for j in range(_sage_const_1 ,G2.M.ncols()): if G2.M[i,j][_sage_const_0 ] != _sage_const_0 : if G2.M[i,j][_sage_const_0 ] == _sage_const_1 : if G2.M[i,j][_sage_const_1 ] != _sage_const_0 : if G2.M[i,j][_sage_const_0 ] == _sage_const_1 : jj = emap2[j] for ii in vmap2[i]: if H.M[ii,jj] != _sage_const_0 : Hcopy.M[ii,jj] += G2.M[i,j][_sage_const_1 ]*X break elif G2.M[i,j][_sage_const_1 ] == _sage_const_0 : loop_factor *= _sage_const_2 else: jj = emap2[j] psilooplist.append([[G2.M[i,j][_sage_const_1 ],G2.M[i,j][_sage_const_2 ]],[G2.M[i,j][_sage_const_2 ],G2.M[i,j][_sage_const_1 ]]]) psiindexlist.append([jj]) for ii in vmap2[i]: for k in range(H.M[ii,jj][_sage_const_0 ]): psiindexlist[-_sage_const_1 ].append(ii) Klocationlist = [] Kindexlist = [] for i in range(_sage_const_1 ,G1.M.nrows()): for r in range(_sage_const_1 ,rmax+_sage_const_1 ): for k in range(G1.M[i,_sage_const_0 ][r]): Klocationlist.append(vmap1[i]) Kindexlist.append(r) for i in range(_sage_const_1 ,G2.M.nrows()): for r in range(_sage_const_1 ,rmax+_sage_const_1 ): for k in range(G2.M[i,_sage_const_0 ][r]): Klocationlist.append(vmap2[i]) Kindexlist.append(r) psilist = [] for j in B: S = [i for i in range(_sage_const_1 ,H.M.nrows()) if H.M[i,j][_sage_const_0 ] != _sage_const_0 ] if len(S) == _sage_const_2 : psilist.append([[S[_sage_const_0 ],j],[S[_sage_const_1 ],j]]) else: psilooplist.append([[_sage_const_0 ,_sage_const_1 ],[_sage_const_1 ,_sage_const_0 ]]) psiindexlist.append([j,S[_sage_const_0 ],S[_sage_const_0 ]]) for psiloopvals in itertools.product(*psilooplist): for Klocs in itertools.product(*Klocationlist): for psilocs in itertools.product(*psilist): Hcopycopy = Graph(Hcopy.M) for i in range(len(psiindexlist)): Hcopycopy.M[psiindexlist[i][_sage_const_1 ],psiindexlist[i][_sage_const_0 ]] += psiloopvals[i][_sage_const_0 ]*X if psiindexlist[i][_sage_const_1 ] == psiindexlist[i][_sage_const_2 ]: Hcopycopy.M[psiindexlist[i][_sage_const_1 ],psiindexlist[i][_sage_const_0 ]] += psiloopvals[i][_sage_const_1 ]*X**_sage_const_2 else: Hcopycopy.M[psiindexlist[i][_sage_const_2 ],psiindexlist[i][_sage_const_0 ]] += psiloopvals[i][_sage_const_1 ]*X for i in range(len(Kindexlist)): Hcopycopy.M[Klocs[i],_sage_const_0 ] += X**Kindexlist[i] for i in psilocs: Hcopycopy.M[i[_sage_const_0 ],i[_sage_const_1 ]] += X for k in Hloops: if Hcopycopy.M[k][_sage_const_2 ] > Hcopycopy.M[k][_sage_const_1 ]: Hcopycopy.M[k] += (Hcopycopy.M[k][_sage_const_2 ] - Hcopycopy.M[k][_sage_const_1 ])*(X - X**_sage_const_2 ) Hcopycopy.compute_invariant() for which_gen in unpurify[L[_sage_const_0 ]]: if graph_isomorphic(Hcopycopy,gens[which_gen]): answer[which_gen] += (-_sage_const_1 )**len(B)*loop_factor/auts break return answer
[docs]def check_associativity(g,r1,r2,r3,markings=(),moduli_type=MODULI_ST): ngens1 = num_strata(g,r1,markings,moduli_type) ngens2 = num_strata(g,r2,markings,moduli_type) ngens3 = num_strata(g,r3,markings,moduli_type) print "%s*%s*%s = %s associators to compute:" % (ngens1, ngens2, ngens3, ngens1*ngens2*ngens3) count = _sage_const_0 for i1 in range(ngens1): for i2 in range(ngens2): for i3 in range(ngens3): a = capply(multiply,r1,i1,r2,i2,g,r1+r2+r3,markings,moduli_type) answer1 = vector([_sage_const_0 for i in range(num_strata(g,r1+r2+r3,markings,moduli_type))]) for j in range(num_strata(g,r1+r2,markings,moduli_type)): if a[j] == _sage_const_0 : continue answer1 += a[j]*vector(capply(multiply,r1+r2,j,r3,i3,g,r1+r2+r3,markings,moduli_type)) a = capply(multiply,r1,i1,r3,i3,g,r1+r2+r3,markings,moduli_type) answer2 = vector([_sage_const_0 for i in range(num_strata(g,r1+r2+r3,markings,moduli_type))]) for j in range(num_strata(g,r1+r3,markings,moduli_type)): if a[j] == _sage_const_0 : continue answer2 += a[j]*vector(capply(multiply,r1+r3,j,r2,i2,g,r1+r2+r3,markings,moduli_type)) if answer1 != answer2: print "Error: %s %s %s" % (i1,i2,i3) count += _sage_const_1 if count % _sage_const_100 == _sage_const_0 : print "%s done" % count
[docs]def gorenstein_precompute(g,r1,markings=(),moduli_type=MODULI_ST): r3 = dim_form(g,len(markings),moduli_type) r2 = r3-r1 D = capply(all_strata,g,r1,markings,moduli_type) D = capply(all_strata,g,r2,markings,moduli_type) D = capply(contraction_table,g,r3,markings,moduli_type) D = capply(unpurify_map,g,r3,markings,moduli_type)
[docs]def pairing_matrix(g,r1,markings=(),moduli_type=MODULI_ST): r3 = dim_form(g,len(markings),moduli_type) r2 = r3-r1 ngens1 = num_strata(g,r1,markings,moduli_type) ngens2 = num_strata(g,r2,markings,moduli_type) ngens3 = num_strata(g,r3,markings,moduli_type) socle_evaluations = [socle_evaluation(i,g,markings,moduli_type) for i in range(ngens3)] pairings = [[_sage_const_0 for i2 in range(ngens2)] for i1 in range(ngens1)] if r1 == r2: # print "%s*%s/2 = %s pairings to compute:" % (ngens1, ngens1+1, ngens1*(ngens1+1)/2) sym = True else: # print "%s*%s = %s pairings to compute:" % (ngens1, ngens2, ngens1*ngens2) sym = False count = _sage_const_0 for i1 in range(ngens1): for i2 in range(ngens2): if sym and i1 > i2: pairings[i1][i2] = pairings[i2][i1] continue L = capply(multiply,r1,i1,r2,i2,g,r3,markings,moduli_type) pairings[i1][i2] = sum([L[k]*socle_evaluations[k] for k in range(ngens3)]) count += _sage_const_1 # if count % 100 == 0: # print "%s done" % count return pairings
[docs]def pairing_submatrix(S1,S2,g,r1,markings=(),moduli_type=MODULI_ST): r3 = dim_form(g,len(markings),moduli_type) r2 = r3-r1 ngens1 = num_strata(g,r1,markings,moduli_type) ngens2 = num_strata(g,r2,markings,moduli_type) ngens3 = num_strata(g,r3,markings,moduli_type) #print "Computing socle evaluation" socle_evaluations = [socle_evaluation(i,g,markings,moduli_type) for i in range(ngens3)] pairings = [[_sage_const_0 for i2 in S2] for i1 in S1] if r1 == r2 and S1 == S2: # print "%s*%s/2 = %s pairings to compute:" % (len(S1), len(S1)+1, len(S1)*(len(S1)+1)/2) sym = True else: # print "%s*%s = %s pairings to compute:" % (len(S1), len(S2), len(S1)*len(S2)) sym = False count = _sage_const_0 for i1 in range(len(S1)): for i2 in range(len(S2)): if sym and i1 > i2: pairings[i1][i2] = pairings[i2][i1] continue L = capply(multiply,r1,S1[i1],r2,S2[i2],g,r3,markings,moduli_type) pairings[i1][i2] = sum([L[k]*socle_evaluations[k] for k in range(ngens3)]) count += _sage_const_1 # if count % 100 == 0: # print "%s done" % count return pairings
#####################################
[docs]def betti(g,r,marked_points=(),moduli_type=MODULI_ST): """ This function returns the predicted rank of the codimension r grading of the tautological ring of the moduli space of stable genus g curves with marked points labeled by the multiset marked_points. g and r should be nonnegative integers and marked_points should be a tuple of positive integers. The parameter moduli_type determines which moduli space to use: - MODULI_ST: all stable curves (this is the default) - MODULI_CT: curves of compact type - MODULI_RT: curves with rational tails - MODULI_SM: smooth curves EXAMPLES:: sage: from admcycles.DR import betti Check rank R^3(bar{M}_2) = 1:: sage: betti(2,3) 1 Check rank R^2(bar{M}_{2,3}) = 44:: sage: betti(2,2,(1,2,3)) 44 Check rank R^2(bar{M}_{2,3})^{S_3} = 20:: sage: betti(2,2,(1,1,1)) 20 Check rank R^2(bar{M}_{2,3})^{S_2} = 32 (S_2 interchanging markings 1 and 2):: sage: betti(2,2,(1,1,2)) 32 Check rank R^2(M^c_4) = rank R^3(M^c_4) = 6:: sage: from admcycles.DR import MODULI_CT, MODULI_RT, MODULI_SM sage: betti(4,2,(),MODULI_CT) 6 sage: betti(4,3,(),MODULI_CT) 6 Check rank R^8(M^rt_{17,2})^(S_2) < R^9(M^rt_{17,2})^(S_2):: sage: betti(17,8,(1,1),MODULI_RT) # long time 122 sage: betti(17,9,(1,1),MODULI_RT) # long time 123 Check rank R^9(M_{20,1}) < rank R^10(M_{20,1}):: sage: betti(20,9,(1,),MODULI_SM) # long time 75 sage: betti(20,10,(1,),MODULI_SM) # long time 76 """ L = list_all_FZ(g,r,marked_points,moduli_type) L.reverse() return (len(L[_sage_const_0 ]) - compute_rank(L))
[docs]def gorenstein(g,r,marked_points=(),moduli_type=MODULI_ST): """ This function returns the rank of the codimension r grading of the Gorenstein quotient of the tautological ring of the moduli space of genus g curves with marked points labeled by the multiset marked_points. g and r should be nonnegative integers and marked_points should be a tuple of positive integers. The parameter moduli_type determines which moduli space to use: - MODULI_ST: all stable curves (this is the default) - MODULI_CT: curves of compact type - MODULI_RT: curves with rational tails EXAMPLES:: sage: from admcycles.DR import gorenstein Check rank Gor^3(bar{M}_{3}) = 10:: sage: gorenstein(3,3) 10 Check rank Gor^2(bar{M}_{2,2}) = 14:: sage: gorenstein(2,2,(1,2)) 14 Check rank Gor^2(bar{M}_{2,2})^{S_2} = 11:: sage: gorenstein(2,2,(1,1)) 11 Check rank Gor^2(M^c_{4}) = 6:: sage: from admcycles.DR import MODULI_CT, MODULI_RT sage: gorenstein(4,2,(),MODULI_CT) 6 Check rank Gor^4(M^rt_{8,2}) = 22:: sage: gorenstein(8,4,(1,2),MODULI_RT) 22 """ #print "Computing lookup tables" gorenstein_precompute(g,r,marked_points,moduli_type) r3 = dim_form(g,len(marked_points),moduli_type) r2 = r3-r S1 = capply(good_generator_list,g,r,marked_points,moduli_type) S2 = capply(good_generator_list,g,r2,marked_points,moduli_type) M = pairing_submatrix(S1,S2,g,r,marked_points,moduli_type) #print "Computing rank" return matrix(M).rank()
#####################################
[docs]def compute_rank(L): count = _sage_const_0 for i in range(len(L)): S = [j for j in range(len(L[_sage_const_0 ])) if L[i][j] != _sage_const_0 ] #print i,len(S) if len(S) == _sage_const_0 : continue count += _sage_const_1 j = S[_sage_const_0 ] T = [ii for ii in range(i+_sage_const_1 ,len(L)) if L[ii][j] != _sage_const_0 ] for k in S[_sage_const_1 :]: rat = L[i][k]/L[i][j] for ii in T: L[ii][k] -= rat*L[ii][j] for ii in range(i+_sage_const_1 ,len(L)): L[ii][j] = _sage_const_0 return count
[docs]def choose_orders_sparse(D,nrows,ncols): row_nums = [_sage_const_0 for i in range(nrows)] col_nums = [_sage_const_0 for j in range(ncols)] for key in D.keys(): row_nums[key[_sage_const_0 ]] += _sage_const_1 col_nums[key[_sage_const_1 ]] += _sage_const_1 row_order = list(range(nrows)) col_order = list(range(ncols)) row_order.sort(key=lambda x: row_nums[x]) col_order.sort(key=lambda x: col_nums[x]) return row_order,col_order
[docs]def choose_orders(L): rows = len(L) if rows == _sage_const_0 : return [],[] cols = len(L[_sage_const_0 ]) row_nums = [_sage_const_0 for i in range(rows)] col_nums = [_sage_const_0 for j in range(cols)] for i in range(rows): for j in range(cols): if L[i][j] != _sage_const_0 : row_nums[i] += _sage_const_1 col_nums[j] += _sage_const_1 row_order = list(range(rows)) col_order = list(range(cols)) row_order.sort(key=lambda x: row_nums[x]) col_order.sort(key=lambda x: col_nums[x]) return row_order,col_order
[docs]def compute_rank2(L,row_order,col_order): count = _sage_const_0 for irow in range(len(row_order)): i = row_order[irow] S = [j for j in col_order if L[i][j] != _sage_const_0 ] #print i,len(S) if len(S) == _sage_const_0 : continue count += _sage_const_1 j = S[_sage_const_0 ] T = [ii for ii in row_order[irow+_sage_const_1 :] if L[ii][j] != _sage_const_0 ] for k in S[_sage_const_1 :]: rat = L[i][k]/L[i][j] for ii in T: L[ii][k] -= rat*L[ii][j] for ii in T: L[ii][j] = _sage_const_0 return count
####################################
[docs]def socle_evaluation(num,g,markings=(),moduli_type=MODULI_ST): answer = _sage_const_1 G = single_stratum(num,g,dim_form(g,len(markings),moduli_type),markings,moduli_type) for i in range(_sage_const_1 ,G.M.nrows()): g0 = G.M[i,_sage_const_0 ][_sage_const_0 ] psilist = [] for j in range(_sage_const_1 ,G.M.ncols()): if G.M[i,j][_sage_const_0 ] > _sage_const_0 : psilist.append(G.M[i,j][_sage_const_1 ]) if G.M[i,j][_sage_const_0 ] == _sage_const_2 : psilist.append(G.M[i,j][_sage_const_2 ]) n0 = len(psilist) dim0 = dim_form(g0,n0,moduli_type) kappalist = [] for j in range(_sage_const_1 ,dim0+_sage_const_1 ): for k in range(G.M[i,_sage_const_0 ][j]): kappalist.append(j) if sum(psilist)+sum(kappalist) != dim0: print "ERROR: wrong dim" return answer *= socle_formula(g0,psilist,kappalist,moduli_type) return answer
[docs]def socle_formula(g,psilist,kappalist,moduli_type=MODULI_ST): if moduli_type == MODULI_CT or g == _sage_const_0 : return CTconst(g)*CTsum(psilist,kappalist) if moduli_type <= MODULI_SM or moduli_type == MODULI_RT: return RTsum(g,psilist,kappalist) if moduli_type == MODULI_ST: return STsum(psilist,kappalist)
[docs]def setparts_recur(symlist,progress): if len(symlist) == _sage_const_0 : return [progress] l = [] for i in Combinations(symlist[_sage_const_1 :]).list(): j = [symlist[_sage_const_0 ]]+i if len(progress) > _sage_const_0 and j < progress[-_sage_const_1 ]: continue cur = _sage_const_0 new_symlist = [] for k in range(len(symlist)): if cur < len(j) and symlist[k] == j[cur]: cur += _sage_const_1 else: new_symlist.append(symlist[k]) l += setparts_recur(new_symlist,progress+[j]) return l
[docs]def setparts_with_auts(symlist): l = setparts_recur(symlist,[]) a = aut(symlist) ll = [] for i in l: b = aut(i) for j in i: b *= aut(j) ll.append([i,a/b]) return ll
[docs]def setparts(symlist): return setparts_recur(symlist,[])
[docs]def multi(sigma): term = factorial(sum(sigma)) for i in sigma: term /= factorial(i) return term
[docs]def multi2(g, sigma): sigma.sort() if sigma[_sage_const_0 ] == _sage_const_0 : total = _sage_const_0 for i in range(len(sigma)-_sage_const_1 ): sigmacopy = sigma[_sage_const_1 :] if sigmacopy[i] > _sage_const_0 : sigmacopy[i] -= _sage_const_1 total += multi2(g,sigmacopy) return total term = factorial(_sage_const_2 *g-_sage_const_3 +len(sigma)) term *= (_sage_const_2 *g-_sage_const_1 ).multifactorial(_sage_const_2 ) term /= factorial(_sage_const_2 *g-_sage_const_1 ) for i in sigma: term /= (_sage_const_2 *i-_sage_const_1 ).multifactorial(_sage_const_2 ) return term
[docs]def STsum(psilist,kappalist): kappalist.sort() total = _sage_const_0 for i in setparts_with_auts(kappalist): total += (-_sage_const_1 )**(len(i[_sage_const_0 ])) * i[_sage_const_1 ] * STrecur([_sage_const_1 +sum(j) for j in i[_sage_const_0 ]] + psilist) return total*(-_sage_const_1 )**(len(kappalist))
[docs]def RTsum(g,psilist,kappalist): kappalist.sort() total = _sage_const_0 for i in setparts_with_auts(kappalist): total += (-_sage_const_1 )**(len(i[_sage_const_0 ])) * i[_sage_const_1 ] * multi2(g, [_sage_const_1 +sum(j) for j in i[_sage_const_0 ]] + psilist) return total*(-_sage_const_1 )**(len(kappalist))
[docs]def CTsum(psilist,kappalist): kappalist.sort() total = _sage_const_0 for i in setparts_with_auts(kappalist): total += (-_sage_const_1 )**(len(i[_sage_const_0 ])) * i[_sage_const_1 ] * multi([_sage_const_1 +sum(j) for j in i[_sage_const_0 ]] + psilist) return total*(-_sage_const_1 )**(len(kappalist))
[docs]def CTconst(g): return abs((_sage_const_2 **(_sage_const_2 *g-_sage_const_1 )-_sage_const_1 )*bernoulli(_sage_const_2 *g))/(_sage_const_2 **(_sage_const_2 *g-_sage_const_1 )*factorial(_sage_const_2 *g))
[docs]def STrecur(psi): psi.sort() return capply(STrecur_calc,tuple(psi))
[docs]def STrecur_calc(psi): n = len(psi) if n == _sage_const_0 : return _sage_const_1 s = sum(psi) if (s - n) % _sage_const_3 != _sage_const_0 : return _sage_const_0 if psi[_sage_const_0 ] == _sage_const_0 : if s == _sage_const_0 and n == _sage_const_3 : return _sage_const_1 total = _sage_const_0 for i in range(n-_sage_const_1 ): psicopy = list(psi[_sage_const_1 :]) if psicopy[i] > _sage_const_0 : psicopy[i] -= _sage_const_1 total += STrecur(psicopy) return total g = (s - n)/_sage_const_3 + _sage_const_1 d = psi[-_sage_const_1 ] total = _sage_const_0 psicopy = [_sage_const_0 ,_sage_const_0 ,_sage_const_0 ,_sage_const_0 ] + list(psi) psicopy[-_sage_const_1 ] += _sage_const_1 total += (_sage_const_2 *d+_sage_const_3 )/_sage_const_12 * STrecur(psicopy) psicopy = [_sage_const_0 ,_sage_const_0 ,_sage_const_0 ] + list(psi) total -= (_sage_const_2 *g+n-_sage_const_1 )/_sage_const_6 * STrecur(psicopy) for I in Subsets(range(n-_sage_const_1 )): psi3 = [_sage_const_0 ,_sage_const_0 ] + [psi[i] for i in I] x = STrecur(psi3) if x == _sage_const_0 : continue psi1 = [_sage_const_0 ,_sage_const_0 ] + [psi[i] for i in range(n-_sage_const_1 ) if i not in I] + [d+_sage_const_1 ] psi2 = list(psi1[_sage_const_1 :]) psi2[-_sage_const_1 ] = d total += ((_sage_const_2 *d+_sage_const_3 )*STrecur(psi1) - (_sage_const_2 *g+n-_sage_const_1 )*STrecur(psi2))*x total /= (_sage_const_2 *g+n-_sage_const_1 )*(_sage_const_2 *g+n-_sage_const_2 ) return total
[docs]def good_generator_list(g,r,markings=(),moduli_type=MODULI_ST): gens = capply(all_strata,g,r,markings,moduli_type) good_gens = [] ngens = len(gens) for num in range(ngens): G = gens[num] good = True for i in range(_sage_const_1 ,G.M.nrows()): g = G.M[i,_sage_const_0 ][_sage_const_0 ] codim = _sage_const_0 for d in range(_sage_const_1 ,r+_sage_const_1 ): if G.M[i,_sage_const_0 ][d] != _sage_const_0 : if _sage_const_3 *d > g: good = False break codim += d*G.M[i,_sage_const_0 ][d] if not good: break for j in range(_sage_const_1 ,G.M.ncols()): codim += G.M[i,j][_sage_const_1 ] codim += G.M[i,j][_sage_const_2 ] if codim > _sage_const_0 and codim >= g: good = False break if good: good_gens.append(num) return good_gens
[docs]def automorphism_cosets(num,g,r,markings=(),moduli_type=MODULI_ST): G = single_stratum(num,g,r,markings,moduli_type) pureG = Graph(G.M) pureG.purify() pureG.compute_invariant() pure_auts = graph_list_isomorphisms(pureG, pureG) num_pure = len(pure_auts) impure_auts = graph_list_isomorphisms(G, G) num_impure = len(impure_auts) chosen_auts = [] used_auts = [] v = G.num_vertices() e = G.num_edges() for i in range(num_pure): if i not in used_auts: chosen_auts.append(pure_auts[i]) for g in impure_auts: sigma = [[pure_auts[i][_sage_const_0 ][g[_sage_const_0 ][k]-_sage_const_1 ] for k in range(v)],[pure_auts[i][_sage_const_1 ][g[_sage_const_1 ][k]-_sage_const_1 ] for k in range(e)]] for ii in range(num_pure): if pure_auts[ii] == sigma: used_auts.append(ii) break return [num_impure,chosen_auts]
#############################################
[docs]def FZ_rels(g,r,markings=(),moduli_type=MODULI_ST): return span(FZ_matrix(g,r,markings,moduli_type)).basis()
[docs]def goren_rels(g,r,markings=(),moduli_type=MODULI_ST): gorenstein_precompute(g,r,markings,moduli_type) r3 = dim_form(g,len(markings),moduli_type) r2 = r3-r S1 = range(num_strata(g,r,markings,moduli_type)) S2 = capply(good_generator_list,g,r2,markings,moduli_type) M = pairing_submatrix(S1,S2,g,r,markings,moduli_type) return Matrix(M).kernel().basis()
[docs]def kappa_conversion(sigma): answer = [] for spart in setparts_with_auts(list(sigma)): coeff = spart[_sage_const_1 ] poly = R(_sage_const_0 ) for part in spart[_sage_const_0 ]: coeff *= factorial(len(part) - _sage_const_1 ) poly += X**sum(part) answer.append([poly,coeff]) return answer
[docs]def kappa_conversion_inverse(sigma): answer = [] for spart in setparts_with_auts(list(sigma)): coeff = spart[_sage_const_1 ]*(-_sage_const_1 )**(len(sigma)-len(spart[_sage_const_0 ])) poly = R(_sage_const_0 ) for part in spart[_sage_const_0 ]: poly += X**sum(part) answer.append([poly,coeff]) return answer
[docs]def convert_to_monomial_basis(num,g,r,markings=(),moduli_type=MODULI_ST): answer = [] G = single_stratum(num,g,r,markings,moduli_type) genus_vec = [] kappa_vec = [] for i in range(_sage_const_1 ,G.M.nrows()): genus_vec.append(G.M[i,_sage_const_0 ][_sage_const_0 ]) kappa_vec.append([]) for j in range(_sage_const_1 ,r+_sage_const_1 ): for k in range(G.M[i,_sage_const_0 ][j]): kappa_vec[-_sage_const_1 ].append(j) kappa_vec[-_sage_const_1 ] = capply(kappa_conversion,tuple(kappa_vec[-_sage_const_1 ])) for choice in itertools.product(*kappa_vec): coeff = _sage_const_1 GG = Graph(G.M) for i in range(_sage_const_1 ,G.M.nrows()): GG.M[i,_sage_const_0 ] = genus_vec[i-_sage_const_1 ] + choice[i-_sage_const_1 ][_sage_const_0 ] coeff *= choice[i-_sage_const_1 ][_sage_const_1 ] answer.append((num_of_stratum(GG,g,r,markings,moduli_type), coeff)) return answer
[docs]def convert_to_pushforward_basis(num,g,r,markings=(),moduli_type=MODULI_ST): answer = [] G = single_stratum(num,g,r,markings,moduli_type) genus_vec = [] kappa_vec = [] for i in range(_sage_const_1 ,G.M.nrows()): genus_vec.append(G.M[i,_sage_const_0 ][_sage_const_0 ]) kappa_vec.append([]) for j in range(_sage_const_1 ,r+_sage_const_1 ): for k in range(G.M[i,_sage_const_0 ][j]): kappa_vec[-_sage_const_1 ].append(j) kappa_vec[-_sage_const_1 ] = capply(kappa_conversion_inverse,tuple(kappa_vec[-_sage_const_1 ])) for choice in itertools.product(*kappa_vec): coeff = _sage_const_1 GG = Graph(G.M) for i in range(_sage_const_1 ,G.M.nrows()): GG.M[i,_sage_const_0 ] = genus_vec[i-_sage_const_1 ] + choice[i-_sage_const_1 ][_sage_const_0 ] coeff *= choice[i-_sage_const_1 ][_sage_const_1 ] answer.append((num_of_stratum(GG,g,r,markings,moduli_type), coeff)) return answer
[docs]def convert_vector_to_monomial_basis(vec,g,r,markings=(),moduli_type=MODULI_ST): l = len(vec) vec2 = [_sage_const_0 for i in range(l)] for i in range(l): if vec[i] != _sage_const_0 : for x in capply(convert_to_monomial_basis,i,g,r,markings,moduli_type): vec2[x[_sage_const_0 ]] += x[_sage_const_1 ]*vec[i] return vec2
[docs]def convert_vector_to_pushforward_basis(vec,g,r,markings=(),moduli_type=MODULI_ST): l = len(vec) vec2 = [_sage_const_0 for i in range(l)] for i in range(l): if vec[i] != _sage_const_0 : for x in capply(convert_to_pushforward_basis,i,g,r,markings,moduli_type): vec2[x[_sage_const_0 ]] += x[_sage_const_1 ]*vec[i] return vec2
[docs]def kappa_multiple(vec,which_kappa,g,r,n=_sage_const_0 ,moduli_type=MODULI_ST): vec2 = [] for x in vec: for y in capply(single_kappa_multiple,x[_sage_const_0 ],which_kappa,g,r,n,moduli_type): vec2.append([y[_sage_const_0 ], x[_sage_const_1 ]*y[_sage_const_1 ]]) vec2 = simplify_sparse(vec2) return vec2
[docs]def psi_multiple(vec,which_psi,g,r,n=_sage_const_0 ,moduli_type=MODULI_ST): vec2 = [] for x in vec: for y in capply(single_psi_multiple,x[_sage_const_0 ],which_psi,g,r,n,moduli_type): vec2.append([y[_sage_const_0 ], x[_sage_const_1 ]*y[_sage_const_1 ]]) vec2 = simplify_sparse(vec2) return vec2
[docs]def insertion_pullback(vec,g,r,n=_sage_const_0 ,new_mark=_sage_const_1 ,moduli_type=MODULI_ST): vec2 = [] for x in vec: for y in capply(single_insertion_pullback,x[_sage_const_0 ],g,r,n,new_mark,moduli_type): vec2.append([y[_sage_const_0 ], x[_sage_const_1 ]*y[_sage_const_1 ]]) vec2 = simplify_sparse(vec2) return vec2
[docs]def insertion_pullback2(vec,g,r,n=_sage_const_0 ,new_mark=_sage_const_1 ,moduli_type=MODULI_ST): vec2 = [] for x in vec: for y in capply(single_insertion_pullback2,x[_sage_const_0 ],g,r,n,new_mark,moduli_type): vec2.append([y[_sage_const_0 ], x[_sage_const_1 ]*y[_sage_const_1 ]]) vec2 = simplify_sparse(vec2) return vec2
[docs]def strata_invariant_lookup(g,r,markings=(),moduli_type=MODULI_ST): inv_dict = {} L = capply(all_strata,g,r,markings,moduli_type) for i in range(len(L)): if not inv_dict.has_key(L[i].invariant): inv_dict[L[i].invariant] = [] inv_dict[L[i].invariant].append(i) return inv_dict
[docs]def num_of_stratum(G,g,r,markings=(),moduli_type=MODULI_ST): G.compute_invariant() L = capply(all_strata,g,r,markings,moduli_type) LL = capply(strata_invariant_lookup,g,r,markings,moduli_type) # for debugging purposes try: x = LL[G.invariant] except KeyError: print("num_of_stratum(G={}, g={}, r={}, markings={}, moduli_type={}".format( G,g,r,markings,moduli_type)) print("G.invariant = {}".format(G.invariant)) print("LL.keys() = ") for l in LL: print(l) raise if len(x) == _sage_const_1 : return x[_sage_const_0 ] for i in x: if graph_isomorphic(G,L[i]): return i print "ERROR" print (g,r,markings,moduli_type) print G.M
[docs]def single_psi_multiple(num,which_psi,g,r,n=_sage_const_0 ,moduli_type=MODULI_ST): markings = tuple(range(_sage_const_1 ,n+_sage_const_1 )) G = single_stratum(num,g,r,markings,moduli_type) answer = [] for j in range(_sage_const_1 ,G.M.ncols()): if G.M[_sage_const_0 ,j] == which_psi: good_j = j break for i in range(_sage_const_1 ,G.M.nrows()): if G.M[i,good_j] != _sage_const_0 : deg = _sage_const_0 dim_used = _sage_const_0 for j in range(_sage_const_1 ,r+_sage_const_1 ): dim_used += j*G.M[i,_sage_const_0 ][j] for j in range(_sage_const_1 ,G.M.ncols()): dim_used += G.M[i,j][_sage_const_1 ] + G.M[i,j][_sage_const_2 ] deg += G.M[i,j][_sage_const_0 ] if dim_used < dim_form(G.M[i,_sage_const_0 ][_sage_const_0 ],deg,moduli_type): GG = Graph(G.M) GG.M[i,good_j] += X answer.append((num_of_stratum(GG,g,r+_sage_const_1 ,markings,moduli_type),_sage_const_1 )) break return answer
[docs]def single_kappa_multiple(num,which_kappa,g,r,n=_sage_const_0 ,moduli_type=MODULI_ST): markings = tuple(range(_sage_const_1 ,n+_sage_const_1 )) G = single_stratum(num,g,r,markings,moduli_type) answer = [] for i in range(_sage_const_1 ,G.M.nrows()): deg = _sage_const_0 dim_used = _sage_const_0 for j in range(_sage_const_1 ,r+_sage_const_1 ): dim_used += j*G.M[i,_sage_const_0 ][j] for j in range(_sage_const_1 ,G.M.ncols()): dim_used += G.M[i,j][_sage_const_1 ] + G.M[i,j][_sage_const_2 ] deg += G.M[i,j][_sage_const_0 ] if dim_used + which_kappa <= dim_form(G.M[i,_sage_const_0 ][_sage_const_0 ],deg,moduli_type): GG = Graph(G.M) GG.M[i,_sage_const_0 ] += X**which_kappa answer.append((num_of_stratum(GG,g,r+which_kappa,markings,moduli_type),_sage_const_1 )) for j in range(_sage_const_1 ,r+_sage_const_1 ): if G.M[i,_sage_const_0 ][j] > _sage_const_0 : GG = Graph(G.M) GG.M[i,_sage_const_0 ] += X**(j+which_kappa) GG.M[i,_sage_const_0 ] -= X**j answer.append((num_of_stratum(GG,g,r+which_kappa,markings,moduli_type),-G.M[i,_sage_const_0 ][j])) return answer
# Only doing this for markings=range(1,n+1) right now. # Also, this function uses the monomial basis.
[docs]def single_kappa_psi_multiple(num,kappa_partition,psi_exps,g,r,n=_sage_const_0 ,moduli_type=MODULI_ST): markings = tuple(range(_sage_const_1 ,n+_sage_const_1 )) G = single_stratum(num,g,r,markings,moduli_type) GG = Graph(G.M) for j in range(_sage_const_1 ,G.M.ncols()): if GG.M[_sage_const_0 ,j] != _sage_const_0 : for i in range(_sage_const_1 ,G.M.nrows()): if GG.M[i,j] != _sage_const_0 : GG.M[i,j] += psi_exps[GG.M[_sage_const_0 ,j][_sage_const_0 ]-_sage_const_1 ]*X break rnew = r+sum(kappa_partition)+sum(psi_exps) answer = [] kappa_options = [range(_sage_const_1 ,G.M.nrows()) for i in range(len(kappa_partition))] for kappa_distrib in itertools.product(*kappa_options): GGG = Graph(GG.M) for i in range(len(kappa_partition)): GGG.M[kappa_distrib[i],_sage_const_0 ] += X**(kappa_partition[i]) is_bad = False for i in range(_sage_const_1 ,GGG.M.nrows()): deg = _sage_const_0 dim_used = _sage_const_0 for j in range(_sage_const_1 ,rnew+_sage_const_1 ): dim_used += j*GGG.M[i,_sage_const_0 ][j] for j in range(_sage_const_1 ,GGG.M.ncols()): dim_used += GGG.M[i,j][_sage_const_1 ] + GGG.M[i,j][_sage_const_2 ] deg += GGG.M[i,j][_sage_const_0 ] if dim_used > dim_form(GGG.M[i,_sage_const_0 ][_sage_const_0 ],deg,moduli_type): is_bad = True break if is_bad: continue answer.append((num_of_stratum(GGG,g,rnew,markings,moduli_type), _sage_const_1 )) return answer
[docs]def single_insertion_pullback(num,g,r,n=_sage_const_0 ,new_mark=_sage_const_1 ,moduli_type=MODULI_ST): markings = tuple(range(_sage_const_1 ,n+_sage_const_1 )) new_markings = tuple(range(_sage_const_1 ,n+_sage_const_2 )) G = single_stratum(num,g,r,markings,moduli_type) answer = [] for i in range(_sage_const_1 ,G.M.nrows()): GG = Graph(G.M) for j in range(_sage_const_1 ,G.M.ncols()): if GG.M[_sage_const_0 ,j][_sage_const_0 ] >= new_mark: GG.M[_sage_const_0 ,j] += _sage_const_1 GG.add_edge(i,_sage_const_0 ,new_mark) answer.append((num_of_stratum(GG,g,r,new_markings,moduli_type), _sage_const_1 )) for j in range(_sage_const_1 ,r+_sage_const_1 ): for k in range(GG.M[i,_sage_const_0 ][j]): GGG = Graph(GG.M) GGG.M[i,_sage_const_0 ] -= X**j GGG.M[i,-_sage_const_1 ] += j*X answer.append((num_of_stratum(GGG,g,r,new_markings,moduli_type), -_sage_const_1 )) if moduli_type <= MODULI_SM: continue for j in range(_sage_const_1 ,G.M.ncols()): if G.M[i,j][_sage_const_0 ] == _sage_const_1 : if G.M[i,j][_sage_const_1 ] >= _sage_const_1 : x = G.M[i,j][_sage_const_1 ] GGG = Graph(GG.M) row1 = [GG.M[i,k] for k in range(GG.M.ncols())] row2 = [_sage_const_0 for k in range(GG.M.ncols())] row1[j] = _sage_const_0 row1[-_sage_const_1 ] = _sage_const_0 row2[j] = _sage_const_1 row2[-_sage_const_1 ] = _sage_const_1 GGG.split_vertex(i,row1,row2) GGG.M[-_sage_const_2 ,-_sage_const_1 ] += (x-_sage_const_1 )*X answer.append((num_of_stratum(GGG,g,r,new_markings,moduli_type), -_sage_const_1 )) if G.M[i,j][_sage_const_0 ] == _sage_const_2 : if G.M[i,j][_sage_const_1 ] >= _sage_const_1 or G.M[i,j][_sage_const_2 ] >= _sage_const_1 : x = G.M[i,j][_sage_const_1 ] y = G.M[i,j][_sage_const_2 ] row1 = [GG.M[i,k] for k in range(GG.M.ncols())] row2 = [_sage_const_0 for k in range(GG.M.ncols())] row1[j] = _sage_const_0 row1[-_sage_const_1 ] = _sage_const_0 row2[j] = _sage_const_1 row2[-_sage_const_1 ] = _sage_const_1 if y >= _sage_const_1 : row1[j] = _sage_const_1 + x*X GGG = Graph(GG.M) GGG.split_vertex(i,row1,row2) GGG.M[-_sage_const_2 ,-_sage_const_1 ] += (y-_sage_const_1 )*X answer.append((num_of_stratum(GGG,g,r,new_markings,moduli_type), -_sage_const_1 )) if x >= _sage_const_1 : row1[j] = _sage_const_1 + y*X GGG = Graph(GG.M) GGG.split_vertex(i,row1,row2) GGG.M[-_sage_const_2 ,-_sage_const_1 ] += (x-_sage_const_1 )*X answer.append((num_of_stratum(GGG,g,r,new_markings,moduli_type), -_sage_const_1 )) return answer
[docs]def single_insertion_pullback2(num,g,r,n=_sage_const_0 ,new_mark=_sage_const_1 ,moduli_type=MODULI_ST): markings = tuple(range(_sage_const_1 ,n+_sage_const_1 )) new_markings = tuple(range(_sage_const_1 ,n+_sage_const_2 )) G = single_stratum(num,g,r,markings,MODULI_SMALL) answer = [] for i in range(_sage_const_1 ,G.M.nrows()): GG = Graph(G.M) for j in range(_sage_const_1 ,G.M.ncols()): if GG.M[_sage_const_0 ,j][_sage_const_0 ] >= new_mark: GG.M[_sage_const_0 ,j] += _sage_const_1 GG.add_edge(i,_sage_const_0 ,new_mark) answer.append((num_of_stratum(GG,g,r,new_markings,moduli_type), _sage_const_1 )) for j in range(_sage_const_1 ,r+_sage_const_1 ): for k in range(GG.M[i,_sage_const_0 ][j]): GGG = Graph(GG.M) GGG.M[i,_sage_const_0 ] -= X**j GGG.M[i,-_sage_const_1 ] += j*X answer.append((num_of_stratum(GGG,g,r,new_markings,moduli_type), -_sage_const_1 )) if moduli_type <= MODULI_SM: continue for j in range(_sage_const_1 ,G.M.ncols()): if G.M[i,j][_sage_const_0 ] == _sage_const_1 : if G.M[i,j][_sage_const_1 ] >= _sage_const_1 : x = G.M[i,j][_sage_const_1 ] GGG = Graph(GG.M) row1 = [GG.M[i,k] for k in range(GG.M.ncols())] row2 = [_sage_const_0 for k in range(GG.M.ncols())] row1[j] = _sage_const_0 row1[-_sage_const_1 ] = _sage_const_0 row2[j] = _sage_const_1 row2[-_sage_const_1 ] = _sage_const_1 GGG.split_vertex(i,row1,row2) GGG.M[-_sage_const_2 ,-_sage_const_1 ] += (x-_sage_const_1 )*X answer.append((num_of_stratum(GGG,g,r,new_markings,moduli_type), -_sage_const_1 )) return answer
[docs]def num_new_rels(g,r,n=_sage_const_0 ,moduli_type=MODULI_ST): return len(capply(choose_basic_rels,g,r,n,moduli_type))
[docs]def choose_basic_rels(g,r,n=_sage_const_0 ,moduli_type=MODULI_ST): if _sage_const_3 *r < g+n+_sage_const_1 : return [] sym_ngen = num_strata(g,r,tuple([_sage_const_1 for i in range(n)]),moduli_type) if moduli_type == MODULI_SMALL and r > dim_form(g,n,MODULI_SM): sym_possible_rels = [[[i,_sage_const_1 ]] for i in range(sym_ngen)] else: sym_possible_rels = possibly_new_FZ(g,r,n,moduli_type) if len(sym_possible_rels) == _sage_const_0 : return [] dprint("Start basic_rels (%s,%s,%s,%s): %s",g,r,n,moduli_type,floor(get_memory_usage())) previous_rels = derived_rels(g,r,n,moduli_type) nrels = len(previous_rels) dprint("%s gens, %s oldrels",sym_ngen,nrels) D = {} for i in range(nrels): for x in previous_rels[i]: D[i,x[_sage_const_0 ]] = x[_sage_const_1 ] if nrels > _sage_const_0 : row_order,col_order = choose_orders_sparse(D,nrels,sym_ngen) previous_rank = compute_rank_sparse(D,row_order,col_order) dprint("rank %s",previous_rank) else: previous_rank = _sage_const_0 row_order = [] col_order = list(range(sym_ngen)) answer = [] for j in range(len(sym_possible_rels)): for x in sym_possible_rels[j]: D[nrels,x[_sage_const_0 ]] = x[_sage_const_1 ] row_order.append(nrels) nrels += _sage_const_1 if compute_rank_sparse(D,row_order,col_order) > previous_rank: answer.append(unsymmetrize_vec(sym_possible_rels[j],g,r,tuple(range(_sage_const_1 ,n+_sage_const_1 )),moduli_type)) previous_rank += _sage_const_1 dprint("rank %s",previous_rank) #if len(answer) > 0: #print "%s,%s,%s: %s" % (g,r,n,len(answer)) dprint("End basic_rels (%s,%s,%s,%s): %s",g,r,n,moduli_type,floor(get_memory_usage())) dprint("%s,%s,%s,%s: rank %s",g,r,n,moduli_type,sym_ngen-previous_rank) if moduli_type > -_sage_const_1 : dsave("sparse-%s,%s,%s,%s|%s,%s,%s",g,r,n,moduli_type,len(answer),sym_ngen-previous_rank,floor(get_memory_usage())) #if moduli_type >= 0 and sym_ngen-previous_rank != betti(g,r,tuple([1 for i in range(n)]),moduli_type): # dprint("ERROR: %s,%s,%s,%s",g,r,n,moduli_type) # return return answer
[docs]def recursive_betti(g,r,markings=(),moduli_type=MODULI_ST): dprint("Start recursive_betti (%s,%s,%s,%s): %s",g,r,markings,moduli_type,floor(get_memory_usage())) n = len(markings) if r > dim_form(g,n,moduli_type): return _sage_const_0 ngen = num_strata(g,r,markings,moduli_type) dprint("%s gens",ngen) relations = [] partial_sym_map = capply(partial_symmetrize_map,g,r,markings,moduli_type) for rel in capply(choose_basic_rels,g,r,n,moduli_type): rel2 = [] for x in rel: rel2.append([partial_sym_map[x[_sage_const_0 ]], x[_sage_const_1 ]]) rel2 = simplify_sparse(rel2) relations.append(rel2) for rel in capply(interior_derived_rels,g,r,n,moduli_type): rel2 = [] for x in rel: rel2.append([partial_sym_map[x[_sage_const_0 ]], x[_sage_const_1 ]]) rel2 = simplify_sparse(rel2) relations.append(rel2) dprint("%s gens, %s rels so far",ngen,len(relations)) dprint("Middle recursive_betti (%s,%s,%s,%s): %s",g,r,markings,moduli_type,floor(get_memory_usage())) if moduli_type > MODULI_SM: for r0 in range(_sage_const_1 ,r): strata = capply(all_strata,g,r0,markings,moduli_type) for G in strata: vertex_orbits = graph_count_automorphisms(G,True) for i in [orbit[_sage_const_0 ] for orbit in vertex_orbits]: good = True for j in range(G.M.ncols()): if R(G.M[i,j][_sage_const_0 ]) != G.M[i,j]: good = False break if good: g2 = G.M[i,_sage_const_0 ][_sage_const_0 ] if _sage_const_3 *(r-r0) < g2 + _sage_const_1 : continue d = G.degree(i) if dim_form(g2,d,moduli_type) < r-r0: continue strata2 = capply(all_strata,g2,r-r0,tuple(range(_sage_const_1 ,d+_sage_const_1 )),moduli_type) which_gen_list = [-_sage_const_1 for num in range(len(strata2))] for num in range(len(strata2)): G_copy = Graph(G.M) G_copy.replace_vertex_with_graph(i,strata2[num]) which_gen_list[num] = num_of_stratum(G_copy,g,r,markings,moduli_type) rel_list = copy(capply(choose_basic_rels,g2,r-r0,d,moduli_type)) rel_list += capply(interior_derived_rels,g2,r-r0,d,moduli_type) for rel0 in rel_list: relation = [] for x in rel0: num = x[_sage_const_0 ] if which_gen_list[num] != -_sage_const_1 : relation.append([which_gen_list[num], x[_sage_const_1 ]]) relation = simplify_sparse(relation) relations.append(relation) dprint("%s gens, %s rels",ngen,len(relations)) dsave("sparse-%s-gens-%s-rels",ngen,len(relations)) dprint("Middle recursive_betti (%s,%s,%s,%s): %s",g,r,markings,moduli_type,floor(get_memory_usage())) relations = remove_duplicates2(relations) dprint("%s gens, %s distinct rels",ngen,len(relations)) dsave("sparse-%s-distinct-rels",len(relations)) rank = _sage_const_0 D = {} nrels = len(relations) for i in range(nrels): for x in relations[i]: D[i,x[_sage_const_0 ]] = x[_sage_const_1 ] if nrels > _sage_const_0 : row_order,col_order = choose_orders_sparse(D,nrels,ngen) rank = compute_rank_sparse(D,row_order,col_order) dsave("sparse-answer-%s",ngen-rank) return ngen - rank
[docs]def subsequences(seq,l): answer = [] m = len(seq) for vec in IntegerVectors(l,m,max_part=_sage_const_1 ): answer.append([seq[i] for i in range(m) if vec[i] == _sage_const_1 ]) return answer
[docs]def pullback_derived_rels(g,r,n=_sage_const_0 ,moduli_type=MODULI_ST): if r == _sage_const_0 : return [] dprint("Start pullback_derived (%s,%s,%s,%s): %s",g,r,n,moduli_type,floor(get_memory_usage())) answer = [] for n0 in range(n): if dim_form(g,n0,moduli_type) >= r: basic_rels = capply(choose_basic_rels,g,r,n0,moduli_type) for rel in basic_rels: for vec in subsequences(range(_sage_const_1 ,n+_sage_const_1 ),n-n0): rel2 = copy(rel) for i in range(n-n0): rel2 = insertion_pullback(rel2,g,r,n0+i,vec[i],moduli_type) answer.append(rel2) else: basic_rels = capply(choose_basic_rels,g,r,n0,MODULI_SMALL) k = r - dim_form(g,n0,moduli_type) for rel in basic_rels: for vec in subsequences(range(_sage_const_1 ,n0+k-_sage_const_1 +_sage_const_1 ),k-_sage_const_1 ): for vec2 in subsequences(range(_sage_const_1 ,n+_sage_const_1 ),n-n0-k+_sage_const_1 ): rel2 = copy(rel) for i in range(k-_sage_const_1 ): rel2 = insertion_pullback(rel2,g,r,n0+i,vec[i],MODULI_SMALL) rel2 = insertion_pullback2(rel2,g,r,n0+k-_sage_const_1 ,vec2[_sage_const_0 ],moduli_type) for i in range(n-n0-k): rel2 = insertion_pullback(rel2,g,r,n0+k+i,vec2[i+_sage_const_1 ],moduli_type) answer.append(rel2) dprint("End pullback_derived (%s,%s,%s,%s): %s",g,r,n,moduli_type,floor(get_memory_usage())) return answer
[docs]def interior_derived_rels(g,r,n=_sage_const_0 ,moduli_type=MODULI_ST): dprint("Start interior_derived (%s,%s,%s,%s): %s",g,r,n,moduli_type,floor(get_memory_usage())) markings = tuple(range(_sage_const_1 ,n+_sage_const_1 )) answer = copy(capply(pullback_derived_rels,g,r,n,moduli_type)) for r0 in range(r): pullback_rels = copy(capply(choose_basic_rels,g,r0,n,moduli_type)) pullback_rels += capply(pullback_derived_rels,g,r0,n,moduli_type) for rel in pullback_rels: for i in range(r-r0+_sage_const_1 ): for sigma in Partitions(i): for tau in IntegerVectors(r-r0-i,n): rel2 = copy(rel) rcur = r0 for m in range(n): for mm in range(tau[m]): rel2 = psi_multiple(rel2,m+_sage_const_1 ,g,rcur,n,moduli_type) rcur += _sage_const_1 for m in sigma: rel2 = kappa_multiple(rel2,m,g,rcur,n,moduli_type) rcur += m answer.append(rel2) dprint("End interior_derived (%s,%s,%s,%s): %s",g,r,n,moduli_type,floor(get_memory_usage())) return answer
[docs]def symmetrize_map(g,r,markings=(),moduli_type=MODULI_ST): markings2 = tuple([_sage_const_1 for i in markings]) gens = capply(all_strata,g,r,markings,moduli_type) map = [] for G in gens: GG = Graph(G.M) for i in range(_sage_const_1 ,GG.M.ncols()): if GG.M[_sage_const_0 ,i][_sage_const_0 ] > _sage_const_0 : GG.M[_sage_const_0 ,i] = R(_sage_const_1 ) map.append(num_of_stratum(GG,g,r,markings2,moduli_type)) return map
[docs]def partial_symmetrize_map(g,r,markings=(),moduli_type=MODULI_ST): markings1 = tuple(range(_sage_const_1 ,len(markings)+_sage_const_1 )) gens = capply(all_strata,g,r,markings1,moduli_type) map = [] for G in gens: GG = Graph(G.M) for i in range(_sage_const_1 ,GG.M.ncols()): if GG.M[_sage_const_0 ,i][_sage_const_0 ] > _sage_const_0 : GG.M[_sage_const_0 ,i] = R(markings[GG.M[_sage_const_0 ,i][_sage_const_0 ]-_sage_const_1 ]) map.append(num_of_stratum(GG,g,r,markings,moduli_type)) return map
[docs]def unsymmetrize_map(g,r,markings=(),moduli_type=MODULI_ST): markings2 = tuple([_sage_const_1 for i in markings]) sym_map = capply(symmetrize_map,g,r,markings,moduli_type) map = [[] for i in range(num_strata(g,r,markings2,moduli_type))] for i in range(len(sym_map)): map[sym_map[i]].append(i) return map
[docs]def unsymmetrize_vec(vec,g,r,markings=(),moduli_type=MODULI_ST): unsym_map = capply(unsymmetrize_map,g,r,markings,moduli_type) vec2 = [] for x in vec: aut = len(unsym_map[x[_sage_const_0 ]]) for j in unsym_map[x[_sage_const_0 ]]: vec2.append([j, QQ((x[_sage_const_1 ], aut))]) vec2 = simplify_sparse(vec2) return vec2
[docs]def derived_rels(g,r,n=_sage_const_0 ,moduli_type=MODULI_ST): dprint("Start derived (%s,%s,%s,%s): %s",g,r,n,moduli_type,floor(get_memory_usage())) markings = tuple([_sage_const_1 for i in range(n)]) generators = capply(all_strata,g,r,markings,moduli_type) ngen = len(generators) sym_map = capply(symmetrize_map,g,r,tuple(range(_sage_const_1 ,n+_sage_const_1 )),moduli_type) answer = [] for rel in capply(interior_derived_rels,g,r,n,moduli_type): rel2 = [] for x in rel: rel2.append([sym_map[x[_sage_const_0 ]], x[_sage_const_1 ]]) rel2 = simplify_sparse(rel2) answer.append(rel2) if moduli_type <= MODULI_SM: return answer for r0 in range(_sage_const_1 ,r): strata = capply(all_strata,g,r0,markings,moduli_type) for G in strata: vertex_orbits = graph_count_automorphisms(G,True) for i in [orbit[_sage_const_0 ] for orbit in vertex_orbits]: good = True for j in range(G.M.ncols()): if R(G.M[i,j][_sage_const_0 ]) != G.M[i,j]: good = False break if good: g2 = G.M[i,_sage_const_0 ][_sage_const_0 ] if _sage_const_3 *(r-r0) < g2 + _sage_const_1 : continue d = G.degree(i) if dim_form(g2,d,moduli_type) < r-r0: continue strata2 = capply(all_strata,g2,r-r0,tuple(range(_sage_const_1 ,d+_sage_const_1 )),moduli_type) which_gen_list = [-_sage_const_1 for num in range(len(strata2))] for num in range(len(strata2)): G_copy = Graph(G.M) G_copy.replace_vertex_with_graph(i,strata2[num]) which_gen_list[num] = num_of_stratum(G_copy,g,r,markings,moduli_type) rel_list = copy(capply(choose_basic_rels,g2,r-r0,d,moduli_type)) rel_list += capply(interior_derived_rels,g2,r-r0,d,moduli_type) for rel0 in rel_list: relation = [] for x in rel0: num = x[_sage_const_0 ] if which_gen_list[num] != -_sage_const_1 : relation.append([which_gen_list[num], x[_sage_const_1 ]]) relation = simplify_sparse(relation) answer.append(relation) answer = remove_duplicates2(answer) dprint("End derived (%s,%s,%s,%s): %s",g,r,n,moduli_type,floor(get_memory_usage())) return answer
[docs]def remove_duplicates(L): LL = [] for elt in L: duplicate = False for elt2 in LL: if elt2 == elt: duplicate = True break if not duplicate: LL.append(elt) return LL
[docs]def remove_duplicates2(L): L.sort() LL = [] for i in range(len(L)): if i == _sage_const_0 or L[i] != L[i-_sage_const_1 ]: LL.append(L[i]) return LL
[docs]def simplify_sparse(vec): vec.sort() vec2 = [] last_index = None for x in vec: if x[_sage_const_0 ] == last_index: if vec2[-_sage_const_1 ][_sage_const_1 ] == -x[_sage_const_1 ]: vec2 = vec2[:-_sage_const_1 ] last_index = None else: vec2[-_sage_const_1 ][_sage_const_1 ] += x[_sage_const_1 ] else: vec2.append(x) last_index = x[_sage_const_0 ] return vec2
[docs]def compute_rank_sparse(D,row_order,col_order): count = _sage_const_0 nrows = len(row_order) ncols = len(col_order) row_order_rank = [-_sage_const_1 for i in range(nrows)] col_order_rank = [-_sage_const_1 for i in range(ncols)] for i in range(nrows): row_order_rank[row_order[i]] = i for i in range(ncols): col_order_rank[col_order[i]] = i row_contents = [set() for i in range(nrows)] col_contents = [set() for i in range(ncols)] for x in D.keys(): row_contents[x[_sage_const_0 ]].add(x[_sage_const_1 ]) col_contents[x[_sage_const_1 ]].add(x[_sage_const_0 ]) for i in row_order: S = [] for j in row_contents[i]: S.append(j) if len(S) == _sage_const_0 : continue count += _sage_const_1 S.sort(key=lambda x: col_order_rank[x]) j = S[_sage_const_0 ] T = [] for ii in col_contents[j]: if row_order_rank[ii] > row_order_rank[i]: T.append(ii) for k in S[_sage_const_1 :]: rat = QQ((D[i,k], D[i,j])) for ii in T: if not D.has_key((ii,k)): D[ii,k] = _sage_const_0 row_contents[ii].add(k) col_contents[k].add(ii) D[ii,k] -= rat*D[ii,j] if D[ii,k] == _sage_const_0 : D.pop((ii,k)) row_contents[ii].remove(k) col_contents[k].remove(ii) for ii in T: D.pop((ii,j)) row_contents[ii].remove(j) col_contents[j].remove(ii) return count
##########################################
[docs]def veto_for_DR(num,g,r,markings=(),moduli_type=MODULI_ST): G = single_stratum(num,g,r,markings,moduli_type) nr = G.M.nrows() nc = G.M.ncols() marked_vertices = [] for i in range(_sage_const_1 ,nr): if G.M[i,_sage_const_0 ] != R(G.M[i,_sage_const_0 ][_sage_const_0 ]): return true for j in range(_sage_const_1 ,nc): if G.M[i,j][_sage_const_0 ] == _sage_const_2 : return true if G.M[_sage_const_0 ,j] != _sage_const_0 and G.M[i,j] != _sage_const_0 : marked_vertices.append(i) for ii in range(_sage_const_1 ,nr): S = set(marked_vertices) S.add(ii) did_something = true while did_something: did_something = false for i in tuple(S): if i == ii: continue for j in range(_sage_const_1 ,nc): if G.M[i,j] != _sage_const_0 : for i2 in range(_sage_const_1 ,nr): if G.M[i2,j] != _sage_const_0 and i2 not in S: S.add(i2) did_something = true if len(S) < nr-_sage_const_1 : return true return false
[docs]def find_nonsep_pairs(num,g,r,markings=(),moduli_type=MODULI_ST): G = single_stratum(num,g,r,markings,moduli_type) nr = G.M.nrows() nc = G.M.ncols() answer = [] for i1 in range(_sage_const_1 ,nr): for i2 in range(_sage_const_1 ,i1): found_edge = false for j in range(_sage_const_1 ,nc): if G.M[i1,j] != _sage_const_0 and G.M[i2,j] != _sage_const_0 : found_edge = true break if not found_edge: continue S = set([i1]) did_something = true while did_something: did_something = false for i3 in tuple(S): for j in range(_sage_const_1 ,nc): if G.M[i3,j] != _sage_const_0 : for i4 in range(_sage_const_1 ,nr): if G.M[i4,j] != _sage_const_0 and i4 not in S and (i3 != i1 or i4 != i2): S.add(i4) did_something = true if i2 in S: answer.append([i2,i1]) return answer
##### assumes r <= 4 for now
[docs]def DR_coeff_is_known(num,g,r,markings=(),moduli_type=MODULI_ST): return (r <= _sage_const_4 ) # old stuff G = single_stratum(num,g,r,markings,moduli_type) nr = G.M.nrows() nc = G.M.ncols() nonsep_pairs = find_nonsep_pairs(num,g,r,markings,moduli_type) if len(nonsep_pairs) > _sage_const_3 : return false if len(nonsep_pairs) == _sage_const_3 : i1 = nonsep_pairs[_sage_const_0 ][_sage_const_0 ] i2 = nonsep_pairs[_sage_const_0 ][_sage_const_1 ] i3 = nonsep_pairs[_sage_const_2 ][_sage_const_1 ] jlist = [] for j in range(_sage_const_1 ,nc): if len([_sage_const_1 for i in [i1,i2,i3] if G.M[i,j] != _sage_const_0 ]) == _sage_const_2 : jlist.append(j) if len(jlist) > _sage_const_3 : return false for j in jlist: for i in [i1,i2,i3]: if G.M[i,j][_sage_const_1 ] != _sage_const_0 : return false return true # this stuff is old, keeping it around for now just in case for j in range(_sage_const_1 ,nc): ilist = [] for i in range(_sage_const_1 ,nr): if G.M[i,j] != _sage_const_0 : ilist.append(i) if len(ilist) == _sage_const_1 : continue ii1 = ilist[_sage_const_0 ] ii2 = ilist[_sage_const_1 ] jlist = [] for jj in range(_sage_const_1 ,nc): if G.M[ii1,jj] != _sage_const_0 and G.M[ii2,jj] != _sage_const_0 : jlist.append(jj) if len(jlist) == _sage_const_1 or jlist[_sage_const_0 ] != j: continue count = len(jlist) for ii in [ii1,ii2]: for jj in jlist: count += G.M[ii,jj][_sage_const_1 ] if count > _sage_const_3 : return false return true
##### next function doesn't work on arbitrary graphs yet, see above function
[docs]def DR_coeff(num,g,r,n=_sage_const_0 ,dvector=(),moduli_type=MODULI_ST): markings = tuple(range(_sage_const_1 ,n+_sage_const_1 )) G = single_stratum(num,g,r,markings,moduli_type) nr = G.M.nrows() nc = G.M.ncols() answer = QQ((1, capply(autom_count,num,g,r,markings,moduli_type))) def balance(ilist): bal = [_sage_const_0 for i1 in ilist] for j1 in range(_sage_const_1 ,nc): if G.M[_sage_const_0 ,j1] != _sage_const_0 : for i3 in range(_sage_const_1 ,nr): if G.M[i3,j1] != _sage_const_0 : S = set([i3]) did_something = true while did_something: did_something = false for i4 in tuple(S): for j2 in range(_sage_const_1 ,nc): if G.M[i4,j2] != _sage_const_0 : for i5 in range(_sage_const_1 ,nr): if G.M[i5,j2] != _sage_const_0 and i5 not in S and (i4 not in ilist or i5 not in ilist): S.add(i5) did_something = true for kk in range(len(ilist)): if ilist[kk] in S: bal[kk] += dvector[G.M[_sage_const_0 ,j1][_sage_const_0 ] - _sage_const_1 ] return bal def poly4(a,b,c): return (_sage_const_1 /_sage_const_420 )*(-_sage_const_15 *(a**_sage_const_8 +b**_sage_const_8 +c**_sage_const_8 )+_sage_const_40 *(a**_sage_const_7 *b+a**_sage_const_7 *c+b**_sage_const_7 *a+b**_sage_const_7 *c+c**_sage_const_7 *a+c**_sage_const_7 *b)-_sage_const_28 *(a**_sage_const_6 *b**_sage_const_2 +a**_sage_const_6 *c**_sage_const_2 +b**_sage_const_6 *a**_sage_const_2 +b**_sage_const_6 *c**_sage_const_2 +c**_sage_const_6 *a**_sage_const_2 +c**_sage_const_6 *b**_sage_const_2 )-_sage_const_112 *(a**_sage_const_6 *b*c+b**_sage_const_6 *a*c+c**_sage_const_6 *a*b)+_sage_const_84 *(a**_sage_const_5 *b**_sage_const_2 *c+a**_sage_const_5 *c**_sage_const_2 *b+b**_sage_const_5 *a**_sage_const_2 *c+b**_sage_const_5 *c**_sage_const_2 *a+c**_sage_const_5 *a**_sage_const_2 *b+c**_sage_const_5 *b**_sage_const_2 *a)-_sage_const_70 *(a**_sage_const_4 *b**_sage_const_2 *c**_sage_const_2 +b**_sage_const_4 *a**_sage_const_2 *c**_sage_const_2 +c**_sage_const_4 *a**_sage_const_2 *b**_sage_const_2 )) nonsep_pairs = find_nonsep_pairs(num,g,r,markings,moduli_type) if len(nonsep_pairs) == _sage_const_4 : cycle = [nonsep_pairs[_sage_const_0 ][_sage_const_0 ],nonsep_pairs[_sage_const_0 ][_sage_const_1 ],_sage_const_0 ,_sage_const_0 ] for i in range(_sage_const_1 ,_sage_const_4 ): for j in range(_sage_const_2 ): if nonsep_pairs[i][j] == cycle[_sage_const_0 ]: cycle[_sage_const_3 ] = nonsep_pairs[i][_sage_const_1 -j] if nonsep_pairs[i][j] == cycle[_sage_const_1 ]: cycle[_sage_const_2 ] = nonsep_pairs[i][_sage_const_1 -j] dbal = balance(cycle) a = dbal[_sage_const_0 ] b = dbal[_sage_const_0 ] + dbal[_sage_const_1 ] c = dbal[_sage_const_0 ] + dbal[_sage_const_1 ] + dbal[_sage_const_2 ] answer *= poly4(a,b,c) elif len(nonsep_pairs) == _sage_const_3 : i1 = nonsep_pairs[_sage_const_0 ][_sage_const_0 ] i2 = nonsep_pairs[_sage_const_0 ][_sage_const_1 ] i3 = nonsep_pairs[_sage_const_2 ][_sage_const_1 ] cycle = [i1,i2,i3] cycle4 = cycle + [cycle[_sage_const_0 ]] dbal = balance(cycle) dbal4 = dbal + [dbal[_sage_const_0 ]] done = false for k in range(_sage_const_3 ): jlist = [j for j in range(_sage_const_1 ,nc) if (G.M[cycle4[k],j] != _sage_const_0 and G.M[cycle4[k+_sage_const_1 ],j] != _sage_const_0 )] if len(jlist) == _sage_const_2 : answer *= -_sage_const_1 /_sage_const_6 *poly4(_sage_const_0 ,dbal4[k],-dbal4[k+_sage_const_1 ]) done = true elif len(jlist) == _sage_const_1 : if G.M[cycle4[k],jlist[_sage_const_0 ]][_sage_const_1 ] + G.M[cycle4[k+_sage_const_1 ],jlist[_sage_const_0 ]][_sage_const_1 ] > _sage_const_0 : answer *= -_sage_const_1 /_sage_const_2 *poly4(_sage_const_0 ,dbal4[k],-dbal4[k+_sage_const_1 ]) done = true if not done: answer *= (dbal[_sage_const_0 ]**_sage_const_6 + dbal[_sage_const_1 ]**_sage_const_6 + dbal[_sage_const_2 ]**_sage_const_6 - _sage_const_10 *dbal[_sage_const_0 ]**_sage_const_2 *dbal[_sage_const_1 ]**_sage_const_2 *dbal[_sage_const_2 ]**_sage_const_2 )/_sage_const_30 for j in range(_sage_const_1 ,nc): ilist = [] for i in range(_sage_const_1 ,nr): if G.M[i,j] != _sage_const_0 : ilist.append(i) if len(ilist) == _sage_const_1 : x = G.M[ilist[_sage_const_0 ],j][_sage_const_1 ] answer /= factorial(x) answer *= dvector[G.M[_sage_const_0 ,j][_sage_const_0 ] - _sage_const_1 ]**(_sage_const_2 *x) continue if ilist in nonsep_pairs: continue ii1 = ilist[_sage_const_0 ] ii2 = ilist[_sage_const_1 ] jlist = [] for jj in range(_sage_const_1 ,nc): if G.M[ii1,jj] != _sage_const_0 and G.M[ii2,jj] != _sage_const_0 : jlist.append(jj) if len(jlist) == _sage_const_1 : x1 = G.M[ii1,j][_sage_const_1 ] x2 = G.M[ii2,j][_sage_const_1 ] answer *= -_sage_const_1 answer /= factorial(x1) answer /= factorial(x2) answer /= x1+x2+_sage_const_1 answer *= balance([ii1,ii2])[_sage_const_0 ]**(_sage_const_2 *(x1+x2+_sage_const_1 )) continue elif len(jlist) == _sage_const_2 : if jlist[_sage_const_0 ] != j: continue xvec = [G.M[ii1,jlist[_sage_const_0 ]][_sage_const_1 ], G.M[ii1,jlist[_sage_const_1 ]][_sage_const_1 ], G.M[ii2,jlist[_sage_const_0 ]][_sage_const_1 ], G.M[ii2,jlist[_sage_const_1 ]][_sage_const_1 ]] x = sum(xvec) if x == _sage_const_0 : answer *= -_sage_const_1 /_sage_const_6 answer *= balance([ii1,ii2])[_sage_const_0 ]**_sage_const_4 elif x == _sage_const_1 : answer *= -_sage_const_1 /_sage_const_30 answer *= balance([ii1,ii2])[_sage_const_0 ]**_sage_const_6 elif x == _sage_const_2 : if xvec in [[_sage_const_2 ,_sage_const_0 ,_sage_const_0 ,_sage_const_0 ],[_sage_const_0 ,_sage_const_2 ,_sage_const_0 ,_sage_const_0 ],[_sage_const_0 ,_sage_const_0 ,_sage_const_2 ,_sage_const_0 ],[_sage_const_0 ,_sage_const_0 ,_sage_const_0 ,_sage_const_2 ]]: answer *= -_sage_const_1 /_sage_const_168 elif xvec in [[_sage_const_1 ,_sage_const_1 ,_sage_const_0 ,_sage_const_0 ],[_sage_const_0 ,_sage_const_0 ,_sage_const_1 ,_sage_const_1 ],[_sage_const_1 ,_sage_const_0 ,_sage_const_0 ,_sage_const_1 ],[_sage_const_0 ,_sage_const_1 ,_sage_const_1 ,_sage_const_0 ]]: answer *= -_sage_const_1 /_sage_const_280 elif xvec in [[_sage_const_1 ,_sage_const_0 ,_sage_const_1 ,_sage_const_0 ],[_sage_const_0 ,_sage_const_1 ,_sage_const_0 ,_sage_const_1 ]]: answer *= -_sage_const_1 /_sage_const_84 answer *= balance([ii1,ii2])[_sage_const_0 ]**_sage_const_8 continue elif len(jlist) == _sage_const_3 : if jlist[_sage_const_0 ] != j: continue xvec = [G.M[ii1,jlist[_sage_const_0 ]][_sage_const_1 ], G.M[ii1,jlist[_sage_const_1 ]][_sage_const_1 ], G.M[ii1,jlist[_sage_const_2 ]][_sage_const_1 ], G.M[ii2,jlist[_sage_const_0 ]][_sage_const_1 ], G.M[ii2,jlist[_sage_const_1 ]][_sage_const_1 ], G.M[ii2,jlist[_sage_const_2 ]][_sage_const_1 ]] x = sum(xvec) if x == _sage_const_0 : answer *= -_sage_const_1 /_sage_const_90 answer *= balance([ii1,ii2])[_sage_const_0 ]**_sage_const_6 elif x == _sage_const_1 : answer *= -_sage_const_1 /_sage_const_840 answer *= balance([ii1,ii2])[_sage_const_0 ]**_sage_const_8 continue elif len(jlist) == _sage_const_4 : if jlist[_sage_const_0 ] != j: continue answer *= -_sage_const_1 /_sage_const_2520 answer *= balance([ii1,ii2])[_sage_const_0 ]**_sage_const_8 return answer
[docs]def DR_compute_old(g,r,n=_sage_const_0 ,dvector=(),moduli_type=MODULI_ST): answer = [] markings = tuple(range(_sage_const_1 ,n+_sage_const_1 )) for i in range(num_strata(g,r,markings,moduli_type)): if veto_for_DR(i,g,r,markings,moduli_type) or not DR_coeff_is_known(i,g,r,markings,moduli_type): answer.append(_sage_const_0 ) else: answer.append(DR_coeff(i,g,r,n,dvector,moduli_type)) return vector(answer)
[docs]def DR_sym_compute_old(g,r,n,dvector=(),moduli_type=MODULI_ST): unsym_version = DR_compute_old(g,r,n,dvector,moduli_type) sym_map = capply(symmetrize_map,g,r,tuple(range(_sage_const_1 ,n+_sage_const_1 )),moduli_type) ngen = num_strata(g,r,tuple([_sage_const_1 for i in range(n)]),moduli_type) sym_version = [_sage_const_0 for i in range(ngen)] for i in range(len(sym_map)): sym_version[sym_map[i]] += unsym_version[i] return vector(sym_version)
[docs]def DR_uncomputed(g,r,markings,moduli_type=MODULI_ST): #markings = tuple(range(1,n+1)) for i in range(num_strata(g,r,markings,moduli_type)): if not veto_for_DR(i,g,r,markings,moduli_type) and not DR_coeff_is_known(i,g,r,markings,moduli_type): print i print single_stratum(i,g,r,markings,moduli_type).M print "---------------------"
[docs]def constant_cycle(L,g,r,markings=(),moduli_type=MODULI_ST): answer = [_sage_const_0 for i in range(num_strata(g,r,markings,moduli_type))] for i in L: answer[i] += _sage_const_1 /capply(autom_count,i,g,r,markings,moduli_type) return vector(answer)
[docs]def reduce_with_rels(B,vec): vec2 = copy(vec) for row in B: for i in range(len(row)): if row[i] != _sage_const_0 : if vec2[i] != _sage_const_0 : vec2 -= QQ((vec2[i],row[i]))*row break return vec2
[docs]def DR_coeff_setup(num,g,r,n=_sage_const_0 ,dvector=(),kval=_sage_const_0 ,moduli_type=MODULI_ST): markings = tuple(range(_sage_const_1 ,n+_sage_const_1 )) G = single_stratum(num,g,r,markings,moduli_type) nr = G.M.nrows() nc = G.M.ncols() edge_list = [] exp_list = [] scalar_factor = _sage_const_1 /capply(autom_count,num,g,r,markings,moduli_type) for i in range(_sage_const_1 ,nr): for j in range(_sage_const_1 ,G.M[i,_sage_const_0 ].degree()+_sage_const_1 ): scalar_factor /= factorial(G.M[i,_sage_const_0 ][j]) scalar_factor /= factorial(j)**G.M[i,_sage_const_0 ][j] scalar_factor *= (-_sage_const_1 )**G.M[i,_sage_const_0 ][j] scalar_factor *= (kval**_sage_const_2 )**(j*G.M[i,_sage_const_0 ][j]) given_weights = [-kval*(_sage_const_2 *G.M[i+_sage_const_1 ,_sage_const_0 ][_sage_const_0 ] - _sage_const_2 + sum([G.M[i+_sage_const_1 ][j][_sage_const_0 ] for j in range(_sage_const_1 ,nc)])) for i in range(nr-_sage_const_1 )] for j in range(_sage_const_1 ,nc): ilist = [i for i in range(_sage_const_1 ,nr) if G.M[i,j] != _sage_const_0 ] if G.M[_sage_const_0 ,j] == _sage_const_0 : if len(ilist) == _sage_const_1 : i1 = ilist[_sage_const_0 ] i2 = ilist[_sage_const_0 ] exp1 = G.M[i1,j][_sage_const_1 ] exp2 = G.M[i1,j][_sage_const_2 ] else: i1 = ilist[_sage_const_0 ] i2 = ilist[_sage_const_1 ] exp1 = G.M[i1,j][_sage_const_1 ] exp2 = G.M[i2,j][_sage_const_1 ] edge_list.append([i1-_sage_const_1 ,i2-_sage_const_1 ]) exp_list.append(exp1 + exp2 + _sage_const_1 ) scalar_factor /= -factorial(exp1)*factorial(exp2)*(exp1+exp2+_sage_const_1 ) else: exp1 = G.M[ilist[_sage_const_0 ],j][_sage_const_1 ] scalar_factor *= dvector[G.M[_sage_const_0 ,j][_sage_const_0 ]-_sage_const_1 ]**(_sage_const_2 *exp1)/factorial(exp1) given_weights[ilist[_sage_const_0 ] - _sage_const_1 ] += dvector[G.M[_sage_const_0 ,j][_sage_const_0 ] - _sage_const_1 ] return edge_list,exp_list,given_weights,scalar_factor
[docs]def DR_coeff_new(num,g,r,n=_sage_const_0 ,dvector=(),kval=_sage_const_0 ,moduli_type=MODULI_ST): markings = tuple(range(_sage_const_1 ,n+_sage_const_1 )) G = single_stratum(num,g,r,markings,moduli_type) nr = G.M.nrows() nc = G.M.ncols() edge_list,exp_list,given_weights,scalar_factor = DR_coeff_setup(num,g,r,n,dvector,kval,moduli_type) m0 = ceil(sum([abs(i) for i in dvector])/_sage_const_2 ) + g*abs(kval) h0 = nc - nr - n + _sage_const_1 deg = _sage_const_2 *sum(exp_list) mrange = list(range(m0 + _sage_const_1 , m0 + deg + _sage_const_2 )) mvalues = [] for m in mrange: total = _sage_const_0 for weight_data in itertools.product(*[list(range(m)) for i in range(len(edge_list))]): vertex_weights = copy(given_weights) for i in range(len(edge_list)): vertex_weights[edge_list[i][_sage_const_0 ]] += weight_data[i] vertex_weights[edge_list[i][_sage_const_1 ]] -= weight_data[i] if len([i for i in vertex_weights if i % m != _sage_const_0 ]) > _sage_const_0 : continue term = _sage_const_1 for i in range(len(edge_list)): term *= weight_data[i]**(_sage_const_2 *exp_list[i]) total += term mvalues.append(QQ((total,m**h0))) mpoly = interpolate(mrange, mvalues) return mpoly.subs(X = _sage_const_0 )*scalar_factor
[docs]def DR_coeff_setup_m(m,num,g,r,n=_sage_const_0 ,dvector=(),kval=_sage_const_0 ,moduli_type=MODULI_ST): markings = tuple(range(_sage_const_1 ,n+_sage_const_1 )) G = single_stratum(num,g,r,markings,moduli_type) nr = G.M.nrows() nc = G.M.ncols() edge_list = [] exp_list = [] scalar_factor = _sage_const_1 /capply(autom_count,num,g,r,markings,moduli_type) for i in range(_sage_const_1 ,nr): for j in range(_sage_const_1 ,G.M[i,_sage_const_0 ].degree()+_sage_const_1 ): scalar_factor /= factorial(G.M[i,_sage_const_0 ][j]) scalar_factor /= factorial(j)**G.M[i,_sage_const_0 ][j] scalar_factor *= (-_sage_const_1 )**G.M[i,_sage_const_0 ][j] scalar_factor *= (kval**_sage_const_2 - kval*m + m**_sage_const_2 /_sage_const_6 )**(j*G.M[i,_sage_const_0 ][j]) given_weights = [-kval*(_sage_const_2 *G.M[i+_sage_const_1 ,_sage_const_0 ][_sage_const_0 ] - _sage_const_2 + sum([G.M[i+_sage_const_1 ][j][_sage_const_0 ] for j in range(_sage_const_1 ,nc)])) for i in range(nr-_sage_const_1 )] for j in range(_sage_const_1 ,nc): ilist = [i for i in range(_sage_const_1 ,nr) if G.M[i,j] != _sage_const_0 ] if G.M[_sage_const_0 ,j] == _sage_const_0 : if len(ilist) == _sage_const_1 : i1 = ilist[_sage_const_0 ] i2 = ilist[_sage_const_0 ] exp1 = G.M[i1,j][_sage_const_1 ] exp2 = G.M[i1,j][_sage_const_2 ] else: i1 = ilist[_sage_const_0 ] i2 = ilist[_sage_const_1 ] exp1 = G.M[i1,j][_sage_const_1 ] exp2 = G.M[i2,j][_sage_const_1 ] edge_list.append([i1-_sage_const_1 ,i2-_sage_const_1 ]) exp_list.append(exp1 + exp2 + _sage_const_1 ) scalar_factor /= -factorial(exp1)*factorial(exp2)*(exp1+exp2+_sage_const_1 ) else: exp1 = G.M[ilist[_sage_const_0 ],j][_sage_const_1 ] dval = dvector[G.M[_sage_const_0 ,j][_sage_const_0 ]-_sage_const_1 ] if dval < _sage_const_0 : dval = dval + m scalar_factor *= (dval**_sage_const_2 - dval*m + m**_sage_const_2 /_sage_const_6 )**(exp1)/factorial(exp1) given_weights[ilist[_sage_const_0 ] - _sage_const_1 ] += dvector[G.M[_sage_const_0 ,j][_sage_const_0 ] - _sage_const_1 ] return edge_list,exp_list,given_weights,scalar_factor
[docs]def DR_coeff_m(m,num,g,r,n=_sage_const_0 ,dvector=(),kval=_sage_const_0 ,moduli_type=MODULI_ST): markings = tuple(range(_sage_const_1 ,n+_sage_const_1 )) G = single_stratum(num,g,r,markings,moduli_type) nr = G.M.nrows() nc = G.M.ncols() edge_list,exp_list,given_weights,scalar_factor = DR_coeff_setup_m(m,num,g,r,n,dvector,kval,moduli_type) h0 = nc - nr - n + _sage_const_1 total = _sage_const_0 for weight_data in itertools.product(*[list(range(m)) for i in range(len(edge_list))]): vertex_weights = copy(given_weights) for i in range(len(edge_list)): vertex_weights[edge_list[i][_sage_const_0 ]] += weight_data[i] vertex_weights[edge_list[i][_sage_const_1 ]] -= weight_data[i] if len([i for i in vertex_weights if i % m != _sage_const_0 ]) > _sage_const_0 : continue term = _sage_const_1 for i in range(len(edge_list)): dval = weight_data[i] term *= (dval**_sage_const_2 - dval*m + m**_sage_const_2 /_sage_const_6 )**(exp_list[i]) total += term total /= m**h0 total *= scalar_factor return total
[docs]def interpolate(A,B): X=SR.var('X') l = len(A) poly = _sage_const_0 M = [] for i in range(l): M.append(vector([a**i for a in A])) M.append(vector(B)) ker = Matrix(M).kernel().basis()[_sage_const_0 ] for i in range(l): poly += QQ((-ker[i],ker[-_sage_const_1 ]))*X**i return poly
[docs]def DR_compute(g,r,n=_sage_const_0 ,dvector=(),kval=_sage_const_0 ,moduli_type=MODULI_ST): answer = [] markings = tuple(range(_sage_const_1 ,n+_sage_const_1 )) for i in range(num_strata(g,r,markings,moduli_type)): answer.append(DR_coeff_new(i,g,r,n,dvector,kval,moduli_type)) return vector(answer)/_sage_const_2 **r
[docs]def DR_compute_m(m,g,r,n=_sage_const_0 ,dvector=(),kval=_sage_const_0 ,moduli_type=MODULI_ST): answer = [] markings = tuple(range(_sage_const_1 ,n+_sage_const_1 )) for i in range(num_strata(g,r,markings,moduli_type)): answer.append(DR_coeff_m(m,i,g,r,n,dvector,kval,moduli_type)) return vector(answer)
[docs]def DR_sparse(g,r,n=_sage_const_0 ,dvector=(),kval=_sage_const_0 ,moduli_type=MODULI_ST): vec = DR_compute(g,r,n,dvector,kval,moduli_type) sparse_vec = [] for i in range(len(vec)): if vec[i] != _sage_const_0 : sparse_vec.append([i,vec[i]]) return sparse_vec
[docs]def DR_psi_check(g,n,dvector,which_psi): vec = DR_sparse(g,g,n,dvector,_sage_const_0 ) r = g while (r < _sage_const_3 *g-_sage_const_3 +n): vec = psi_multiple(vec,which_psi,g,r,n) r += _sage_const_1 total = _sage_const_0 for pair in vec: total += pair[_sage_const_1 ]*socle_evaluation(pair[_sage_const_0 ],g,tuple(range(_sage_const_1 ,n+_sage_const_1 ))) return total/_sage_const_2 **g
[docs]def DR_reduced(g,dvector=()): n = len(dvector) r = g kval = ZZ(sum(dvector)/(_sage_const_2 *g-_sage_const_2 +n)) vec = DR_compute(g,r,n,dvector,kval) rels = FZ_rels(g,r,tuple(range(_sage_const_1 ,n+_sage_const_1 ))) vec2 = reduce_with_rels(rels,vec) return vec2
[docs]def list_strata(g,r,n=_sage_const_0 ,moduli_type=MODULI_ST): L = all_strata(g,r,tuple(range(_sage_const_1 ,n+_sage_const_1 )),moduli_type) for i in range(len(L)): print "generator %s" % i print L[i].M print "-------------------------"