# Getting started
from admcycles import *

# Entering tautological classes - combinations of divisors
t1=3*sepbdiv(1,(1,2),3,4)-psiclass(4,3,4)^2  # classes live on \Mbar_{3,4}
reset_g_n(2,1)
t2=-1/3*irrbdiv()*lambdaclass(1) # classes live on \Mbar_{2,1}

# Entering tautological classes - listing generators in given degree
list_tautgens(2,0,2) #generators of R^2(\Mbar_{2,0})
L=tautgens(2,0,2);
t3=2*L[3]+L[4]
t3

# Entering a stable graph manually
G = StableGraph([1,1],[[2],[3]],[(2,3)]); G # graph has two genus 1 components, connected by an edge

# Basic operations with tautclasses
s1=psiclass(3,1,3)^2 # psi_3 on \Mbar_{1,3}
s1.forgetful_pushforward([3]) # pi_* (psi_3) under map pi: \Mbar_{1,3} -> \Mbar_{1,2} forgetting marking 3

s2=psiclass(2,1,2)
s2.forgetful_pullback([3]) # pi^* (psi_2) under map pi: \Mbar_{1,3} -> \Mbar_{1,2}

s3=psiclass(2,1,3)*psiclass(3,1,3)^2
s3.evaluate() # compute the degree of the zero cycle s3=psi_2*(psi_3)^2 on \Mbar_{1,3}

s4=psiclass(2,1,2)^2+psiclass(1,1,2)*psiclass(2,1,2)
s4.evaluate()

# Using simplify() to reduce number of terms in tautclass
psisum = psiclass(1,2,1) + 3 * psiclass(1,2,1); psisum
psisimple = psisum.simplify(); psisimple

# A basis of the tautological ring and tautological relations
generating_indices(2,0,2) # gives list of indices of tautgens(2,0,2) forming a basis of R^2(\Mbar_{2,0})
t3.toTautbasis(2,0,2) # expresses class t3 in terms of this basis

# Verifying the divisor relation kappa - psi + delta_0 on \bar M_{1,4}
reset_g_n(1,4)
bgraphs=[bd for bd in list_strata(1,4,1) if bd.numvert()>1]
del0=sum([bd.to_tautclass() for bd in bgraphs]) # sum of boundary classes with separating node
psisum=sum([psiclass(i) for i in range(1,4+1)]) # sum of psi-classes
rel=kappaclass(1)-psisum+del0
rel.is_zero()

# Comparing classes on open subsets of \Mbar_{g,n} using parameter
# moduli = 'st', 'tl', 'ct', 'rt' or 'sm'
kappaclass(1,3,0).toTautbasis(moduli='sm')
lambdaclass(1,3,0).toTautbasis(moduli='sm')
diff = lambdaclass(1,3,0) - (1/12)*kappaclass(1,3,0)
diff.is_zero(moduli='sm')

# Pulling back tautological classes to a boundary divisor
bdry=StableGraph([2,2],[[1],[2]],[(1,2)])
generator=tautgens(4,0,2)[3]; generator
pullback=bdry.boundary_pullback(generator); 
pullback.totensorTautbasis(2)
pullback.totensorTautbasis(2,vecout=true)

# Pushing forward classes from the boundary
B=StableGraph([2,1],[[4,1,2],[3,5]],[(4,5)])
Bclass=B.boundary_pushforward() # class of undecorated boundary divisor
si1=B.boundary_pushforward([fundclass(2,3),-psiclass(2,1,2)]); si1
# si1 is obtained by pushing forward the fundamental class on the genus 2 vertex times -psi_h on the second vertex (where h is the half-edge)
si2=B.boundary_pushforward([-psiclass(1,2,3),fundclass(1,2)]); si2
(Bclass*Bclass-si1-si2).is_zero()

# Double ramification cycles using DR_cycle(g,A)
A=vector((2,4,-6)); B=vector((-3,-1,4))
diff = DR_cycle(1,A)*DR_cycle(1,B)-DR_cycle(1,A)*DR_cycle(1,A+B)
diff.is_zero(moduli='tl') # vanishes by paper [Holmes-Pixton-Schmitt]
diff.is_zero(moduli='st') # does not vanish on locus of all stable curves

# Computing DR-cycles as classes with polynomial coefficients in the input
R.<a1,a2,a3,b1,b2,b3> = PolynomialRing(QQ,6)
A = vector((a1,a2,a3)); B = vector((b1,b2,b3))
diff = DR_cycle(1,A)*DR_cycle(1,B)-DR_cycle(1,A)*DR_cycle(1,A+B)
diff.is_zero(moduli='tl')

# Checking intersection numbers of DR-cycles with lambdaclass from [Buryak-Rossi]
intersect = DR_cycle(1,A)*DR_cycle(1,B)*lambdaclass(1,1,3)
f = intersect.evaluate(); factor(f)
g = f.subs({a3:-a1-a2,b3:-b1-b2}); factor(g)

# Strata of k-differentials using Strataclass(g,k,mu) with mu vector of zero/pole multiplicities
L=Strataclass(2,1,(3,-1)); L.is_zero()
L=Strataclass(2,1,(2,)); (L-Hyperell(2,1)).is_zero()

# Computing Chern classes of R pi_* O(D) for universal curve \C_{g,n} -> \Mbar_{g,n} using generalized_lambda
g=3;n=1; l=1;d=[0];a=[]
s=lambdaclass(2,g,n)
t=generalized_lambda(2,l,d,a,g,n)
(s-t).is_zero()

# Computing the cycle of the hyperelliptic locus in genus 3
H=Hyperell(3,0,0) 
# the cycle of hyperell. curves of genus 3 with 0 marked fixed points of the involution and 0 marked pairs of conjugate points
H.toTautbasis()
# we compare with the known expression H=9*lambda-delta_0-3*delta_1
reset_g_n(3,0)
H2=9*lambdaclass(1)-(1/2)*irrbdiv()-3*sepbdiv(1,())
H2.toTautbasis()

# Specifying Hurwitz cycles
G=PermutationGroup([(1,2)]) # G=Z/2Z
H=HurData(G,[G[1],G[1]]) # this means there are two ramification points with stabilizer generated by G[1], the generator of G

vbeta=Hidentify(2,H,vecout=true) # identify the locus of bielliptic curves (C,p,q) with marked ramification points p,q in \Mbar_{2,2} in terms of the generating set tautgens(2,2,3) of R^3(\Mbar_{2,2})
vector(vbeta)

H2=HurData(G,[G[1],G[1],G[0]])

G=PermutationGroup([(1,2)])
H=HurData(G,[G[1],G[1]])
Biell=Hidentify(2,H,markings=[]) # identify the locus of bielliptic curves C, remembering none of the marked ramification points, inside \Mbar_{2,0}
Biell.toTautbasis(2,0,1) # represent it in the basis of R^1(\Mbar_{2,0})

reset_g_n(2,0)
Biell2=3/4*irrbdiv()+ 3*sepbdiv(1,()) # this is the formula for the bielliptic locus known from the literature
Biell2.toTautbasis(2,0,1)

# Computing the Hurwitz-Hodge integral \int_[\bar B_{2,2,0}] lambda_2
(Biell*lambdaclass(2,2,0)).evaluate() # = 1/48

# Computing Hurwitz-Hodge integral of cyclic triple covers of genus 0 curve against lambda_1
# See [Owens-Somerstep]
G = PermutationGroup([(1,2,3)]); list(G)
H = HurData(G,[G[1],G[1],G[2],G[2]]) #n=2, m=2
t = Hidentify(2,H,markings=[])
(t*lambdaclass(1,2,0)).evaluate()