DR.py

# This file was *autogenerated* from the file DR.sage
from __future__ import absolute_import, print_function
import itertools
from six.moves import range

from copy import copy

from sage.misc.cachefunc import cached_function
from sage.rings.all import PolynomialRing, ZZ
from sage.rings.rational_field import QQ
from sage.arith.misc import factorial, binomial, bernoulli, multinomial
from sage.matrix.constructor import matrix, Matrix
from sage.combinat.all import Combinations, IntegerVectors, Partitions, Permutations, Subsets
from sage.functions.other import floor, ceil
from sage.symbolic.ring import SR
from sage.modules.free_module_element import vector
from sage.modules.free_module import span
from sage.misc.getusage import get_memory_usage

# Sage code by Aaron Pixton ([email protected])
#
# 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)


R = PolynomialRing(ZZ, 1, order='lex', names=('X',))
X = R.gen()
A_list = [ZZ(factorial(6*n))/(factorial(3*n)*factorial(2*n))
          for n in range(100)]
B_list = [ZZ(factorial(6*n+1))/((6*n-1)*factorial(3*n)*factorial(2*n))
          for n in range(100)]
MODULI_SMALL = -1
MODULI_SM = 0
MODULI_RT = 1
MODULI_CT = 2
MODULI_ST = 3
ENABLE_DPRINT = True
ENABLE_DSAVE = False


def dprint(string, *args):
    if ENABLE_DPRINT:
        print(string % args)


def dsave(string, *args):
    if ENABLE_DSAVE:
        save(0, string % args)


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)+1,
                            1, [-1]+genus_list)
        else:
            self.M = matrix(R, 1, 1, -1)

    def num_vertices(self):
        return self.M.nrows() - 1

    def num_edges(self):
        return self.M.ncols() - 1

    def h1(self):
        return self.M.ncols() - self.M.nrows() + 1

    def add_vertex(self, g):
        self.M = self.M.stack(matrix(1, self.M.ncols()))
        self.M[-1, 0] = g

    def add_edge(self, i1, i2, marking=0):
        self.M = self.M.augment(matrix(self.M.nrows(), 1))
        self.M[0, -1] = marking
        if i1 > 0:
            self.M[i1, -1] += 1
        if i2 > 0:
            self.M[i2, -1] += 1

    def del_vertex(self, i):
        self.M = self.M[:i].stack(self.M[(i + 1):])

    def del_edge(self, i):
        if i == self.num_edges():
            self.M = self.M[0:, :i]
        else:
            self.M = self.M[0:, :i].augment(self.M[0:, (i + 1):])

    def compute_degree_vec(self):
        self.degree_vec = [0 for i in range(1, self.M.nrows())]
        for i in range(1, self.M.nrows()):
            for j in range(1, self.M.ncols()):
                self.degree_vec[i - 1] += self.M[i, j][0]

    def degree(self, i):
        return self.degree_vec[i - 1]

    def split_vertex(self, i, row1, row2):
        self.M = self.M.stack(matrix(2, self.M.ncols(), row1+row2))
        self.add_edge(self.M.nrows()-2,
                      self.M.nrows()-1)
        self.del_vertex(i)

    def set_target_parity(self):
        self.target_parity = 0
        for i in range(1, self.M.nrows()):
            local_parity = 1 + self.M[i, 0][0]
            for j in range(1, self.M.ncols()):
                local_parity += self.M[i, j][1] + self.M[i, j][2]
            for j in range(1, self.M[i, 0].degree()+1):
                local_parity += j*self.M[i, 0][j]
            local_parity %= 2
            self.target_parity += (local_parity << (i-1))

    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(1, self.M.ncols()):
            for j in range(self.M[i, k][0]):
                hedge_list.append(k)
        self.del_vertex(i)
        for j in range(G.num_edges() - len(hedge_list)):
            self.add_edge(0, 0)
        for j in range(G.num_vertices()):
            self.add_vertex(G.M[j+1, 0])
        col = ne+1
        for k in range(1, G.M.ncols()):
            if G.M[0, k] > 0:
                mark = ZZ(G.M[0, k])
                for j in range(G.num_vertices()):
                    if self.M[nv+j, hedge_list[mark-1]] == 0:
                        self.M[nv+j, hedge_list[mark-1]] = G.M[j+1, k]
                    elif G.M[j+1, k] != 0:
                        a = self.M[nv+j, hedge_list[mark-1]][1]
                        b = G.M[j+1, k][1]
                        self.M[nv+j, hedge_list[mark-1]] = 2 + \
                            max(a, b)*X + min(a, b)*X**2
            else:
                for j in range(G.num_vertices()):
                    self.M[nv+j, col] = G.M[j+1, k]
                col += 1

    def compute_invariant(self):
        nr, nc = self.M.nrows(), self.M.ncols()
        self.invariant = [[self.M[i, 0], [], [], [
            [] for j in range(1, nr)]] for i in range(1, nr)]
        for k in range(1, nc):
            L = [i for i in range(1, nr)
                 if self.M[i, k] != 0]
            if len(L) == 1:
                if self.M[0, k] != 0:
                    self.invariant[L[0] -
                                   1][2].append((self.M[0, k], self.M[L[0], k]))
                else:
                    self.invariant[L[0]-1][1].append(self.M[L[0], k])
            else:
                self.invariant[L[0]-1][3][L[1]-1].append((self.M[L[0], k],
                                                          self.M[L[1], k]))
                self.invariant[L[1]-1][3][L[0]-1].append((self.M[L[1], k],
                                                          self.M[L[0], k]))
        for i in range(1, nr):
            self.invariant[i - 1][3] = [term for term in self.invariant[i - 1][3]
                                        if len(term)]
            for k in range(len(self.invariant[i - 1][3])):
                self.invariant[i - 1][3][k].sort()
                self.invariant[i -
                               1][3][k] = tuple(self.invariant[i - 1][3][k])
            self.invariant[i - 1][3].sort()
            self.invariant[i - 1][3] = tuple(self.invariant[i - 1][3])
            self.invariant[i - 1][2].sort()
            self.invariant[i - 1][2] = tuple(self.invariant[i - 1][2])
            self.invariant[i - 1][1].sort()
            self.invariant[i - 1][1] = tuple(self.invariant[i - 1][1])
            self.invariant[i - 1] = tuple(self.invariant[i - 1])
        vertex_invariants = [[i, self.invariant[i-1]]
                             for i in range(1, nr)]
        self.invariant.sort()
        self.invariant = tuple(self.invariant)
        vertex_invariants.sort(key=lambda x: x[1])
        self.vertex_groupings = []
        for i in range(nr-1):
            if i == 0 or vertex_invariants[i][1] != vertex_invariants[i-1][1]:
                self.vertex_groupings.append([])
            self.vertex_groupings[-1].append(
                vertex_invariants[i][0])

    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][0])

    def contract(self, i, vlist, elist):
        # assumes graph is undecorated
        if self.M[0, i] != 0:
            print("ERROR: cannot contract a marking")
            return
        S = [row for row in range(
            1, self.M.nrows()) if self.M[row, i] != 0]
        if len(S) == 1:
            self.M[S[0], 0] += 1
            self.del_edge(i)
            elist = elist[:(i-1)] + elist[i:]
        else:
            self.del_edge(i)
            elist = elist[:(i-1)] + elist[i:]
            self.add_vertex(0)
            self.M[-1] += self.M[S[0]]
            self.M[-1] += self.M[S[1]]
            self.del_vertex(S[1])
            self.del_vertex(S[0])
            vlist = vlist[:(S[0]-1)] + vlist[S[0]:(S[1]-1)] + \
                vlist[S[1]:] + [vlist[S[0]-1] + vlist[S[1]-1]]
        return vlist, elist


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)


def isomorphic(M1, M2, group1, group2):
    nr, nc = M1.nrows(), M1.ncols()
    PermList = [Permutations(list(range(len(group)))) for group in group1]
    for sigma_data in itertools.product(*PermList):
        sigma = [0 for i in range(nr - 1)]
        for i in range(len(group1)):
            for j in range(len(group1[i])):
                sigma[group1[i][j] - 1] = group2[i][sigma_data[i][j]]
        good = True
        for i in range(1, nr):
            ii = sigma[i-1]
            for j in range(1, i):
                jj = sigma[j-1]
                L1 = []
                for k in range(1, nc):
                    if M1[i, k] != 0 and M1[j, k] != 0:
                        L1.append([M1[i, k], M1[j, k]])
                L1.sort()
                L2 = []
                for k in range(1, nc):
                    if M2[ii, k] != 0 and M2[jj, k] != 0:
                        L2.append([M2[ii, k], M2[jj, k]])
                L2.sort()
                if L1 != L2:
                    good = False
                    break
            if good is False:
                break
        if good:
            return True
    return False


def graph_count_automorphisms(G, vertex_orbits=False):
    return count_automorphisms(G.M, G.vertex_groupings, vertex_orbits)


def count_automorphisms(M, grouping, vertex_orbits=False):
    nr, nc = M.nrows(), M.ncols()
    count = ZZ.zero()
    PermList = [Permutations(list(range(len(group)))) for group in grouping]
    if vertex_orbits:
        isom_list = []
    for sigma_data in itertools.product(*PermList):
        sigma = [0 for i in range(nr-1)]
        for i in range(len(grouping)):
            for j in range(len(grouping[i])):
                sigma[grouping[i][j]-1] = grouping[i][sigma_data[i][j]]
        good = True
        for i in range(1, nr):
            ii = sigma[i-1]
            for j in range(1, i):
                jj = sigma[j-1]
                L1 = []
                for k in range(1, nc):
                    if M[i, k] != 0 and M[j, k] != 0:
                        L1.append([M[i, k], M[j, k]])
                L1.sort()
                L2 = []
                for k in range(1, nc):
                    if M[ii, k] != 0 and M[jj, k] != 0:
                        L2.append([M[ii, k], M[jj, k]])
                L2.sort()
                if L1 != L2:
                    good = False
                    break
            if good is False:
                break
        if good:
            count += 1
            if vertex_orbits:
                isom_list.append(sigma)

    if vertex_orbits:
        orbit_list = []
        vertices_used = []
        while len(vertices_used) < nr - 1:
            i = [ii for ii in range(1, nr) if ii not in vertices_used][0]
            orbit = []
            for sigma in isom_list:
                if sigma[i - 1] not in orbit:
                    orbit.append(sigma[i - 1])
                    vertices_used.append(sigma[i - 1])
            orbit.sort()
            orbit_list.append(orbit)
        return orbit_list

    for i in range(1, nr):
        for k in range(1, nc):
            if M[i, k][0] == 2 and M[i, k][1] == M[i, k][2]:
                count *= 2
        L = []
        for k in range(1, nc):
            if M[i, k] != 0:
                if sum(1 for j in range(1, nr) if M[j, k] != 0) == 1:
                    L.append([M[0, k], M[i, k]])
        count *= aut(L)

        for j in range(1, i):
            L = []
            for k in range(1, nc):
                if M[i, k] != 0 and M[j, k] != 0:
                    L.append([M[i, k], M[j, k]])
            count *= aut(L)
    return count


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)


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(list(range(len(group)))) for group in group1]
    isom_list = []
    for sigma_data in itertools.product(*PermList):
        sigma = [0 for i in range(nr - 1)]
        for i in range(len(group1)):
            for j in range(len(group1[i])):
                sigma[group1[i][j] - 1] = group2[i][sigma_data[i][j]]
        good = True
        for i in range(1, nr):
            ii = sigma[i-1]
            for j in range(1, i):
                jj = sigma[j-1]
                L1 = []
                for k in range(1, nc):
                    if M1[i, k] != 0 and M1[j, k] != 0:
                        L1.append([M1[i, k], M1[j, k]])
                L1.sort()
                L2 = []
                for k in range(1, nc):
                    if M2[ii, k] != 0 and M2[jj, k] != 0:
                        L2.append([M2[ii, k], M2[jj, k]])
                L2.sort()
                if L1 != L2:
                    good = False
                    break
            if not good:
                break
        if good:
            cols1 = [[M1[i, j] for i in range(nr)] for j in range(1, nc)]
            cols2 = [[M2[0, j]] + [M2[sigma[i - 1], j] for i in range(1, nr)]
                     for j in range(1, nc)]
            edge_group1 = []
            edge_group2 = []
            used1 = []
            for j in range(1, nc):
                if j not in used1:
                    edge_group1.append([])
                    edge_group2.append([])
                    for k in range(1, nc):
                        if cols1[k-1] == cols1[j-1]:
                            edge_group1[-1].append(k)
                            used1.append(k)
                        if cols2[k-1] == cols1[j-1]:
                            edge_group2[-1].append(k)
            edge_PermList = [Permutations(
                list(range(len(edge_group)))) for edge_group in edge_group1]
            for edge_sigma_data in itertools.product(*edge_PermList):
                edge_sigma = [0 for i in range(nc - 1)]
                for i in range(len(edge_group1)):
                    for j in range(len(edge_group1[i])):
                        edge_sigma[edge_group1[i][j] -
                                   1] = edge_group2[i][edge_sigma_data[i][j]]
                isom_list.append([sigma, edge_sigma])
                if only_one:
                    return isom_list
    return isom_list


def aut(L):
    """
    Return the cardinality of the automorphism group of the list ``L``.

    EXAMPLES::

       sage: from admcycles.DR import aut
       sage: aut([])
       1
       sage: aut([4,1,3,2])
       1
       sage: aut([4,5,6,5,4,4,6])
       24
    """
    if not L:
        return ZZ.one()
    L.sort()
    total = ZZ.one()
    n = 1
    last = L[0]
    for l in L[1:]:
        if l == last:
            n += 1
            total *= n
        else:
            n = 1
            last = l
    return total


def degenerate(G_list, moduli_type=MODULI_ST):
    mod_size = moduli_type + 1
    if moduli_type == MODULI_SMALL:
        mod_size = MODULI_SM + 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(1, G.num_vertices()+1):
                row = list(G.M[i])
                m = row[0] + sum(row)
                if m < 4:
                    continue
                row1 = [0 for j in range(len(row))]
                while [2 * x for x in row1] <= row:
                    if row1[0] == 1 and moduli_type <= MODULI_RT:
                        break
                    if row1[0] + sum(row1) >= 2 and row1[0] + sum(row1) <= m - 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[0] > 0:
                            new_type = MODULI_CT
                        G_list_new[new_type].append(G_copy)
                    row1[-1] += 1
                    for j in range(1, len(row)):
                        if row1[-j] <= row[-j]:
                            break
                        row1[-j] = 0
                        row1[-j-1] += 1
    for i in range(mod_size):
        G_list_new[i] = remove_isomorphic(G_list_new[i])
    return G_list_new


def dim_form(g, n, moduli_type=MODULI_ST):
    g = ZZ(g)
    n = ZZ(n)
    if moduli_type == MODULI_ST:
        return 3 * g - 3 + n
    if moduli_type == MODULI_CT:
        return 2 * g - 3 + n
    if moduli_type == MODULI_RT:
        if g > 0:
            return g - 2 + n
        else:
            return n - 3
    if moduli_type == MODULI_SM:
        if n == 0:
            return g - 2
        elif g >= 1:
            return g - 1
        else:
            return ZZ.zero()
    if moduli_type == MODULI_SMALL:
        return ZZ(1000)
    return 3 * g - 3 + n


def decorate(G_list, r, moduli_type=MODULI_ST):
    mod_size = moduli_type + 1
    if moduli_type == MODULI_SMALL:
        mod_size = MODULI_SM + 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(1, nr):
                for j in range(1, nc):
                    if G.M[i, j] == 2:
                        two_list.append([i, j])
                    elif G.M[i, j] == 1:
                        one_list.append([i, j])
            a = nr-1
            b = len(two_list)
            c = len(one_list)
            dims = [[dim_form(G.M[i + 1, 0][0], G.degree(i+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][0] - 1] += vec[a + i]
                    for i in range(c):
                        test_dims[one_list[i][0] - 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 + 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(floor(vec[i]/ZZ(2) + 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 + 1, 0] += X**j
                    for i in range(a, a+b):
                        G_copy.M[two_list[i-a][0], two_list[i-a]
                                 [1]] += vec2[i][0]*X + vec2[i][1]*X**2
                    for i in range(c):
                        G_copy.M[one_list[i][0],
                                 one_list[i][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


def remove_isomorphic(G_list):
    G_list_new = []
    inv_dict = {}
    count = 0
    for G1 in G_list:
        G1.compute_invariant()
        if G1.invariant not in inv_dict:
            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 += 1
    return G_list_new


def num_strata(g, r, markings=(), moduli_type=MODULI_ST):
    return len(all_strata(g, r, markings, moduli_type))


def num_pure_strata(g, r, markings=(), moduli_type=MODULI_ST):
    return len(all_pure_strata(g, r, markings, moduli_type))


def single_stratum(num, g, r, markings=(), moduli_type=MODULI_ST):
    return all_strata(g, r, markings, moduli_type)[num]


def single_pure_stratum(num, g, r, markings=(), moduli_type=MODULI_ST):
    return all_pure_strata(g, r, markings, moduli_type)[num]


@cached_function
def autom_count(num, g, r, markings=(), moduli_type=MODULI_ST):
    return graph_count_automorphisms(single_stratum(num, g, r, markings, moduli_type))


@cached_function
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))


@cached_function
def unpurify_map(g, r, markings=(), moduli_type=MODULI_ST):
    unpurify = {}
    pure_strata = [all_pure_strata(g, r0, markings, moduli_type)
                   for r0 in range(r+1)]
    impure_strata = all_strata(g, r, markings, moduli_type)
    for i, strati in enumerate(impure_strata):
        G = Graph(strati.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 G_key not in unpurify:
            unpurify[G_key] = []
        unpurify[G_key].append(i)
    return unpurify


@cached_function
def all_strata(g, r, markings=(), moduli_type=MODULI_ST):
    mod_size = moduli_type + 1
    if moduli_type == MODULI_SMALL:
        mod_size = MODULI_SM + 1
    big_list = [[] for i in range(mod_size)]
    for loops in range(g+1):
        if loops == 1 and moduli_type <= MODULI_CT:
            break
        if loops > r:
            break
        for edges in range(r-loops+1):
            if edges == 1 and moduli_type <= MODULI_SM:
                break
            G = Graph()
            G.add_vertex(g-loops)
            for k in range(loops):
                G.add_edge(1, 1)
            for k in markings:
                G.add_edge(1, 0, k)
            GGG = [[] for i in range(mod_size)]
            if loops == 0:
                if edges == 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


@cached_function
def all_pure_strata(g, r, markings=(), moduli_type=MODULI_ST):
    big_list = [[] for i in range(moduli_type+1)]
    for loops in range(g+1):
        if loops == 1 and moduli_type <= MODULI_CT:
            break
        if loops > r:
            break
        for edges in range(r-loops, r-loops+1):
            if edges >= 1 and moduli_type <= MODULI_SM:
                break
            G = Graph()
            G.add_vertex(g-loops)
            for k in range(loops):
                G.add_edge(1, 1)
            for k in markings:
                G.add_edge(1, 0, k)
            G.compute_invariant()
            GGG = [[] for i in range(moduli_type+1)]
            if loops == 0:
                if edges == 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+1):
                big_list[i] += GGG[i]
    combined_list = []
    for i in range(moduli_type+1):
        combined_list += big_list[i]
    return combined_list

############################


def C_coeff(m, term):
    n = term - m // 3
    if n < 0:
        return 0
    if m % 3 == 0:
        return A_list[n]
    else:
        return B_list[n]


@cached_function
def dual_C_coeff(i, j, parity):
    total = ZZ.zero()
    k = parity % 2
    while k // 3 <= i:
        if k % 3 != 2:
            total += (-1)**(k // 3) * C_coeff(k, i) * C_coeff(-2 - k, j)
        k += 2
    return total


def poly_to_partition(F):
    mmm = F.degree()
    target_partition = []
    for i in range(1, mmm+1):
        for j in range(F[i]):
            target_partition.append(i)
    return tuple(target_partition)


@cached_function
def kappa_coeff(sigma, kappa_0, target_partition):
    total = ZZ.zero()
    num_ones = sum(1 for i in sigma if i == 1)
    for i in range(0, num_ones+1):
        for injection in Permutations(list(range(len(target_partition))), len(sigma)-i):
            term = binomial(num_ones, i)*binomial(kappa_0 +
                                                  len(target_partition) + i-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(0, target_partition[j])
            total += term
    total = (-1)**(len(target_partition)+len(sigma)) * \
        total/aut(list(target_partition))
    return total


@cached_function
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(1, nv+1):
        L.append((2 * G.M[i, 0][0] + G.degree(i) -
                  2, G.M[i, 0] - G.M[i, 0][0]))
    LL = []
    tau = []
    for i in range(nv):
        mini = -1
        for j in range(nv):
            if (i == 0 or L[j] > LL[-1] or L[j] == LL[-1] and j > tau[-1]) and (mini == -1 or L[j] < L[mini]):
                mini = j
        tau.append(mini)
        LL.append(L[mini])
    factor_dict = FZ_kappa_factor2(tuple(LL), sigma)
    factor_vec = [0 for i in range(1 << nv)]
    for parity_key in factor_dict:
        parity = 0
        for i in range(nv):
            if parity_key[i] == 1:
                parity += 1 << tau[i]
        factor_vec[parity] = factor_dict[parity_key]
    return factor_vec


@cached_function
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 = ZZ(2) ** nv
    PPP_list = []
    for marks in marking_vec:
        PPP_list.append(Permutations(marks[1]))
    PPP = itertools.product(*PPP_list)
    marking_factors = [0 for i in range(num_parities)]
    incident_vertices = []
    for mark_type in marking_vec:
        incident_vertices.append([])
        for k in range(1, ne+1):
            if G.M[0, k] == mark_type[0]:
                for i in range(1, nv+1):
                    if G.M[i, k] != 0:
                        incident_vertices[-1].append((i - 1, G.M[i, k][1]))
                        break
    for perms in PPP:
        parity = 0
        marking_factor = 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][1])
                parity ^= (perms[marks_index][count] %
                           ZZ(2)) << incident_vertices[marks_index][count][0]
        marking_factors[parity] += marking_factor
    return marking_factors


@cached_function
def FZ_kappa_factor2(L, sigma):
    nv = len(L)
    mmm = max((0,)+sigma)
    sigma_grouped = [0 for i in range(mmm)]
    for i in sigma:
        sigma_grouped[i-1] += 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(*[(0, 1) for i in range(nv)]):
        kappa_factors[tuple(parity)] = 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+1)
        sigma_auts = 1
        parity = [0 for i in range(nv)]
        kappa_factor = ZZ.one()
        for j in range(nv):
            sigma_auts *= aut(assigned_sigma[j])
            parity[j] += sum(assigned_sigma[j])
            parity[j] %= ZZ(2)
            kappa_factor *= kappa_coeff(
                tuple(assigned_sigma[j]), L[j][0], poly_to_partition(L[j][1]))
        kappa_factors[tuple(parity)] += kappa_factor/sigma_auts
    return kappa_factors


@cached_function
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 = ZZ(2) ** nv
    ne = G.num_edges()
    edge_list = []
    for k in range(1, ne+1):
        if G.M[0, k] == 0:
            edge_list.append([k])
            for i in range(1, nv+1):
                if G.M[i, k] != 0:
                    edge_list[-1].append(i)
                if G.M[i, k][0] == 2:
                    edge_list[-1].append(i)
    hedge_factors = [0 for i in range(num_parities)]
    for edge_parities in itertools.product(*[[0, 1] for i in edge_list]):
        parity = 0
        for i in range(len(edge_list)):
            if edge_parities[i] == 1:
                parity ^= 1 << (edge_list[i][1]-1)
                parity ^= 1 << (edge_list[i][2]-1)
        hedge_factor = 1
        for i in range(len(edge_list)):
            if edge_list[i][1] == edge_list[i][2]:
                hedge_factor *= dual_C_coeff(G.M[edge_list[i][1], edge_list[i][0]][1],
                                             G.M[edge_list[i][1], edge_list[i][0]][2], edge_parities[i] % ZZ(2))
            else:
                hedge_factor *= dual_C_coeff(G.M[edge_list[i][1], edge_list[i][0]][1],
                                             G.M[edge_list[i][2], edge_list[i][0]][1], edge_parities[i] % ZZ(2))
        hedge_factors[parity] += hedge_factor
    return hedge_factors


@cached_function
def FZ_coeff(num, FZ_param, g, r, markings=(), moduli_type=MODULI_ST):
    sigma = FZ_param[0]
    marking_vec = FZ_param[1]
    G = single_stratum(num, g, r, markings, moduli_type)
    nv = G.num_vertices()
    graph_auts = autom_count(num, g, r, markings, moduli_type)
    h1_factor = ZZ(2) ** G.h1()
    num_parities = ZZ(2) ** nv

    marking_factors = FZ_marking_factor(
        num, marking_vec, g, r, markings, moduli_type)
    kappa_factors = FZ_kappa_factor(num, sigma, g, r, markings, moduli_type)
    hedge_factors = FZ_hedge_factor(num, g, r, markings, moduli_type)

    total = ZZ.zero()
    for i in range(num_parities):
        if marking_factors[i] == 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


@cached_function
def interior_FZ(g, r, markings=(), moduli_type=MODULI_ST):
    ngen = num_strata(g, r, markings, moduli_type)
    relations = []
    FZpl = FZ_param_list(3 * r - g - 1, markings)
    for FZ_param in FZpl:
        relation = [FZ_coeff(i, FZ_param, g, r, markings, moduli_type)
                    for i in range(ngen)]
        relations.append(relation)
    return relations


def possibly_new_FZ(g, r, n=0, moduli_type=MODULI_ST):
    m = 3 * r - g - 1 - n
    if m < 0:
        return []
    dprint("Start FZ (%s,%s,%s,%s): %s", g, r, n, moduli_type,
           floor(get_memory_usage()))
    markings = (1,) * n
    ngen = num_strata(g, r, markings, moduli_type)
    relations = []
    for i in range(m + 1):
        if m - i % 2:
            continue
        for sigma in Partitions(i):
            # TODO: the line above should iterate directly over the set
            # of partitions with all parts of size = 1 (mod 3)
            if any(j % 3 != 1 for j in sigma):
                continue
            if n > 0:
                FZ_param = (tuple(sigma), ((1, markings),))
            else:
                FZ_param = (tuple(sigma), ())
            relation = []
            for j in range(ngen):
                coeff = FZ_coeff(j, FZ_param, g, r, markings, moduli_type)
                if coeff != 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


@cached_function
def boundary_FZ(g, r, markings=(), moduli_type=MODULI_ST):
    if moduli_type <= MODULI_SM:
        return []
    generators = all_strata(g, r, markings, moduli_type)
    ngen = len(generators)
    relations = []
    for r0 in range(1, r):
        strata = all_strata(g, r0, markings, moduli_type)
        for G in strata:
            vertex_orbits = graph_count_automorphisms(G, True)
            for i in [orbit[0] for orbit in vertex_orbits]:
                good = True
                for j in range(G.M.ncols()):
                    if R(G.M[i, j][0]) != G.M[i, j]:
                        good = False
                        break
                if good:
                    g2 = G.M[i, 0][0]
                    if 3 * (r - r0) < g2 + 1:
                        continue
                    d = G.degree(i)
                    if dim_form(g2, d, moduli_type) < r-r0:
                        continue
                    strata2 = all_strata(
                        g2, r-r0, tuple(range(1, d+1)), moduli_type)
                    which_gen_list = [-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, 3 * (r - r0) - g2 - 1)
                    ccccc = 0
                    for FZ_param in rFZpl:
                        relation = [0] * ngen
                        for num in range(len(strata2)):
                            if which_gen_list[num] != -1:
                                relation[which_gen_list[num]] += FZ_coeff(num, FZ_param, g2, r-r0, tuple(
                                    range(1, d+1)), moduli_type)
                        relations.append(relation)
                        ccccc += 1
    return relations


def list_all_FZ(g, r, markings=(), moduli_type=MODULI_ST):
    relations = copy(interior_FZ(g, r, markings, moduli_type))
    if moduli_type > MODULI_SM:
        relations += boundary_FZ(g, r, markings, moduli_type)
    if not relations:
        ngen = num_strata(g, r, markings, moduli_type)
        relations.append([0 for i in range(ngen)])
    return relations


def reduced_FZ_param_list(G, v, g, d, n):
    params = FZ_param_list(n, tuple(range(1, d+1)))
    graph_params = []
    M = matrix(R, 2, d + 1)
    M[0, 0] = -1
    for i in range(1, d+1):
        M[0, i] = i
    for p in params:
        G_copy = Graph(G.M)
        M[1, 0] = -g-1
        for j in p[0]:
            M[1, 0] += X**j
        for i in range(1, d + 1):
            M[1, i] = 1 + p[1][i-1][1][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[1].compute_invariant()
        good = True
        for GG in graphs_seen:
            if graph_isomorphic(x[1], GG):
                good = False
                break
        if good:
            graphs_seen.append(x[1])
            params_reduced.append(x[0])
    return params_reduced


def FZ_param_list(n, markings=()):
    if n < 0:
        return []
    final_list = []
    mmm = max((0,)+markings)
    markings_grouped = [0 for i in range(mmm)]
    for i in markings:
        markings_grouped[i-1] += 1
    markings_best = []
    for i in range(mmm):
        if markings_grouped[i] > 0:
            markings_best.append([i+1, markings_grouped[i]])
    for j in range(floor(n/ZZ(2) + 1)):
        for n_vec in IntegerVectors(n-2 * j, 1 + len(markings_best)):
            S_list = [[list(sigma) for sigma in Partitions(n_vec[0])
                       if not any(l % 3 == 2 for l in sigma)]]
            # TODO: the line above should iterate directly over the set
            # of partitions with no parts of size = 2 (mod 3)
            for i in range(len(markings_best)):
                S_list.append(Partitions(
                    n_vec[i+1]+markings_best[i][1], length=markings_best[i][1]).list())
                S_list[-1] = [[k - 1 for k in sigma] for sigma in S_list[-1]
                              if sum(1 for l in sigma if (l % 3) == 0) == 0]
            for S in itertools.product(*S_list):
                final_list.append((tuple(S[0]), tuple([(markings_best[k][0], tuple(
                    S[k+1])) for k in range(len(markings_best))])))
    return final_list


def FZ_matrix(g, r, markings=(), moduli_type=MODULI_ST):
    return matrix(list_all_FZ(g, r, markings, moduli_type))

#####################################


@cached_function
def contraction_table(g, r, markings=(), moduli_type=MODULI_ST):
    contraction_dict = {}
    pure_strata = [all_pure_strata(g, r0, markings, moduli_type)
                   for r0 in range(r+1)]
    for r0 in range(r+1):
        for ii in range(len(pure_strata[r0])):
            G = pure_strata[r0][ii]
            S = [j for j in range(1, G.M.ncols()) if G.M[0, j] == 0]
            contractions = {}
            for edge_subset in itertools.product(*[[0, 1] for i in range(r0)]):
                key = tuple(i for i in range(
                    r0) if edge_subset[i] == 0)
                A = [S[i] for i in key]
                A.reverse()
                vlist = [[i] for i in range(1, G.M.nrows())]
                elist = list(range(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 L:
                        contraction_result.append((rnew, i))
                        contraction_result.append(L[0])
                        break
                contraction_result.append((vlist, elist))
                contractions[key] = contraction_result

            for edge_assignment in itertools.product(*[[0, 1, 2] for i in range(r0)]):
                if sum(1 for i in edge_assignment if i == 1) > r-r0:
                    continue
                key1 = tuple(i for i in range(
                    r0) if edge_assignment[i] == 0)
                B = [S[i]
                     for i in range(r0) if edge_assignment[i] == 1]
                key2 = tuple(i for i in range(
                    r0) if edge_assignment[i] == 2)
                if key1 > key2:
                    continue
                contract1 = contractions[key1]
                contract2 = contractions[key2]
                dict_key = [contract1[0], contract2[0]]
                dict_entry = [contract1[1], contract2[1]]
                if dict_key[0] > dict_key[1]:
                    dict_key.reverse()
                    dict_entry.reverse()
                    dict_entry = [(r0, ii), B, contract2[2],
                                  contract1[2]] + dict_entry
                else:
                    dict_entry = [(r0, ii), B, contract1[2],
                                  contract2[2]] + dict_entry
                dict_key = tuple(dict_key)
                if dict_key not in contraction_dict:
                    contraction_dict[dict_key] = []
                contraction_dict[dict_key].append(dict_entry)
                if dict_key[0] == dict_key[1]:
                    contraction_dict[dict_key].append(
                        dict_entry[:2] + [dict_entry[3], dict_entry[2], dict_entry[5], dict_entry[4]])
    return contraction_dict


@cached_function
def multiply(r1, i1, r2, i2, g, rmax, markings=(), moduli_type=MODULI_ST):
    unpurify = unpurify_map(g, r1+r2, markings, moduli_type)
    gens = all_strata(g, r1+r2, markings, moduli_type)
    ngens = num_strata(g, r1+r2, markings, moduli_type)
    answer = [0 for i in range(ngens)]
    pure_strata = [all_pure_strata(g, r, markings, moduli_type)
                   for r in range(rmax+1)]
    contraction_dict = 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 multiply(r2, i2, r1, i1, g, rmax, markings, moduli_type)

    if (G1_key, G2_key) not in contraction_dict:
        return answer
    for L in contraction_dict[(G1_key, G2_key)]:
        H = pure_strata[L[0][0]][L[0][1]]
        Hloops = []
        if moduli_type > MODULI_CT:
            for i in range(1, H.M.nrows()):
                for j in range(1, H.M.ncols()):
                    if H.M[i, j][0] == 2:
                        Hloops.append((i, j))
        auts = pure_strata_autom_count(
            L[0][1], g, L[0][0], markings, moduli_type)
        B = L[1]
        if len(B) == pure_r1 and len(B) == pure_r2:
            auts *= 2
        aut_cosets1 = automorphism_cosets(i1, g, r1, markings, moduli_type)
        aut_cosets2 = automorphism_cosets(i2, g, r2, markings, moduli_type)
        auts /= aut_cosets1[0]*aut_cosets2[0]
        for isom1 in aut_cosets1[1]:
            for isom2 in aut_cosets2[1]:
                Hcopy = Graph(H.M)
                vmap1 = [0 for i in range(G1.M.nrows())]
                for i in range(1, G1.M.nrows()):
                    vmap1[i] = L[2][0][L[4][0][isom1[0]
                                               [i-1]-1]-1]
                emap1 = [0 for i in range(G1.M.ncols())]
                for i in range(1, G1.M.ncols()):
                    emap1[i] = L[2][1][L[4][1][isom1[1]
                                               [i-1]-1]-1]
                vmap2 = [0 for i in range(G2.M.nrows())]
                for i in range(1, G2.M.nrows()):
                    vmap2[i] = L[3][0][L[5][0][isom2[0]
                                               [i-1]-1]-1]
                emap2 = [0 for i in range(G2.M.ncols())]
                for i in range(1, G2.M.ncols()):
                    emap2[i] = L[3][1][L[5][1][isom2[1]
                                               [i-1]-1]-1]

                psilooplist = []
                psiindexlist = []
                loop_factor = ZZ.one()
                for i in range(1, G1.M.nrows()):
                    for j in range(1, G1.M.ncols()):
                        if G1.M[i, j][0] != 0:
                            if G1.M[i, j][0] == 1:
                                if G1.M[i, j][1] != 0:
                                    jj = emap1[j]
                                    for ii in vmap1[i]:
                                        if H.M[ii, jj] != 0:
                                            Hcopy.M[ii, jj] += G1.M[i, j][1]*X
                                            break
                            elif G1.M[i, j][1] == 0:
                                loop_factor *= 2
                            else:
                                jj = emap1[j]
                                psilooplist.append([[G1.M[i, j][1], G1.M[i, j][2]], [
                                                   G1.M[i, j][2], G1.M[i, j][1]]])
                                psiindexlist.append([jj])
                                for ii in vmap1[i]:
                                    for k in range(H.M[ii, jj][0]):
                                        psiindexlist[-1].append(ii)
                for i in range(1, G2.M.nrows()):
                    for j in range(1, G2.M.ncols()):
                        if G2.M[i, j][0] != 0:
                            if G2.M[i, j][0] == 1:
                                if G2.M[i, j][1] != 0:
                                    if G2.M[i, j][0] == 1:
                                        jj = emap2[j]
                                        for ii in vmap2[i]:
                                            if H.M[ii, jj] != 0:
                                                Hcopy.M[ii,
                                                        jj] += G2.M[i, j][1]*X
                                                break
                            elif G2.M[i, j][1] == 0:
                                loop_factor *= 2
                            else:
                                jj = emap2[j]
                                psilooplist.append([[G2.M[i, j][1], G2.M[i, j][2]], [
                                                   G2.M[i, j][2], G2.M[i, j][1]]])
                                psiindexlist.append([jj])
                                for ii in vmap2[i]:
                                    for k in range(H.M[ii, jj][0]):
                                        psiindexlist[-1].append(ii)

                Klocationlist = []
                Kindexlist = []
                for i in range(1, G1.M.nrows()):
                    for r in range(1, rmax+1):
                        for k in range(G1.M[i, 0][r]):
                            Klocationlist.append(vmap1[i])
                            Kindexlist.append(r)
                for i in range(1, G2.M.nrows()):
                    for r in range(1, rmax+1):
                        for k in range(G2.M[i, 0][r]):
                            Klocationlist.append(vmap2[i])
                            Kindexlist.append(r)

                psilist = []
                for j in B:
                    S = [i for i in range(1, H.M.nrows()) if H.M[i, j][0] != 0]
                    if len(S) == 2:
                        psilist.append([[S[0], j], [S[1], j]])
                    else:
                        psilooplist.append([[0, 1], [1, 0]])
                        psiindexlist.append([j, S[0], S[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][1],
                                            psiindexlist[i][0]] += psiloopvals[i][0]*X
                                if psiindexlist[i][1] == psiindexlist[i][2]:
                                    Hcopycopy.M[psiindexlist[i][1],
                                                psiindexlist[i][0]] += psiloopvals[i][1]*X**2
                                else:
                                    Hcopycopy.M[psiindexlist[i][2],
                                                psiindexlist[i][0]] += psiloopvals[i][1]*X
                            for i in range(len(Kindexlist)):
                                Hcopycopy.M[Klocs[i],
                                            0] += X**Kindexlist[i]
                            for i in psilocs:
                                Hcopycopy.M[i[0], i[1]] += X
                            for k in Hloops:
                                if Hcopycopy.M[k][2] > Hcopycopy.M[k][1]:
                                    Hcopycopy.M[k] += (Hcopycopy.M[k]
                                                       [2] - Hcopycopy.M[k][1])*(X - X**2)
                            Hcopycopy.compute_invariant()
                            for which_gen in unpurify[L[0]]:
                                if graph_isomorphic(Hcopycopy, gens[which_gen]):
                                    answer[which_gen] += (-1)**len(B) * \
                                        loop_factor/auts
                                    break
    return answer


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 = 0
    for i1 in range(ngens1):
        for i2 in range(ngens2):
            for i3 in range(ngens3):
                a = multiply(r1, i1, r2, i2, g, r1+r2 +
                             r3, markings, moduli_type)
                answer1 = vector(
                    [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] == 0:
                        continue
                    answer1 += a[j]*vector(multiply(r1+r2, j, r3,
                                                    i3, g, r1+r2+r3, markings, moduli_type))
                a = multiply(r1, i1, r3, i3, g, r1+r2 +
                             r3, markings, moduli_type)
                answer2 = vector(
                    [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] == 0:
                        continue
                    answer2 += a[j]*vector(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 += 1
                if count % 100 == 0:
                    print("%s done" % count)


def gorenstein_precompute(g, r1, markings=(), moduli_type=MODULI_ST):
    r3 = dim_form(g, len(markings), moduli_type)
    r2 = r3 - r1
    all_strata(g, r1, markings, moduli_type)
    all_strata(g, r2, markings, moduli_type)
    contraction_table(g, r3, markings, moduli_type)
    unpurify_map(g, r3, markings, moduli_type)


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 = [[0 for i2 in range(ngens2)] for i1 in range(ngens1)]
    sym = bool(r1 == r2)
    for i1 in range(ngens1):
        for i2 in range(ngens2):
            if sym and i1 > i2:
                pairings[i1][i2] = pairings[i2][i1]
                continue
            L = multiply(r1, i1, r2, i2, g, r3, markings, moduli_type)
            pairings[i1][i2] = sum([L[k]*socle_evaluations[k]
                                    for k in range(ngens3)])
    return pairings


@cached_function
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)
    socle_evaluations = [socle_evaluation(
        i, g, markings, moduli_type) for i in range(ngens3)]
    pairings = [[0 for i2 in S2] for i1 in S1]
    sym = bool(r1 == r2 and S1 == S2)
    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 = 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)])
    return pairings

#####################################


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

    We can use this to check that rank R^8(M^rt_{17,2})^(S_2)=122 < R^9(M^rt_{17,2})^(S_2) = 123.
    Indeed, betti(17,8,(1,1),MODULI_RT) = 122 and betti(17,9,(1,1),MODULI_RT)=123
    Similarly, we have rank R^9(M_{20,1}) = 75 < rank R^10(M_{20,1}) = 76
    Indeed, betti(20,9,(1,),MODULI_SM) = 75 and betti(20,10,(1,),MODULI_SM) = 76.
    """
    L = list_all_FZ(g, r, marked_points, moduli_type)
    L.reverse()
    return (len(L[0]) - compute_rank(L))


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
    """
    gorenstein_precompute(g, r, marked_points, moduli_type)
    r3 = dim_form(g, len(marked_points), moduli_type)
    r2 = r3-r
    S1 = good_generator_list(g, r, marked_points, moduli_type)
    S2 = good_generator_list(g, r2, marked_points, moduli_type)
    M = pairing_submatrix(tuple(S1), tuple(S2), g, r,
                          marked_points, moduli_type)
    return matrix(M).rank()

#####################################


def compute_rank(L):
    count = 0
    for i in range(len(L)):
        S = [j for j in range(len(L[0])) if L[i][j] != 0]
        if not S:
            continue
        count += 1
        j = S[0]
        T = [ii for ii in range(i+1, len(L))
             if L[ii][j] != 0]
        for k in S[1:]:
            rat = L[i][k]/L[i][j]
            for ii in T:
                L[ii][k] -= rat*L[ii][j]
        for ii in range(i+1, len(L)):
            L[ii][j] = 0
    return count


def choose_orders_sparse(D, nrows, ncols):
    row_nums = [0 for i in range(nrows)]
    col_nums = [0 for j in range(ncols)]
    for key in D:
        row_nums[key[0]] += 1
        col_nums[key[1]] += 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


def choose_orders(L):
    rows = len(L)
    if rows == 0:
        return [], []
    cols = len(L[0])
    row_nums = [0 for i in range(rows)]
    col_nums = [0 for j in range(cols)]
    for i in range(rows):
        for j in range(cols):
            if L[i][j] != 0:
                row_nums[i] += 1
                col_nums[j] += 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


def compute_rank2(L, row_order, col_order):
    count = 0
    for irow in range(len(row_order)):
        i = row_order[irow]
        S = [j for j in col_order if L[i][j] != 0]
        if not S:
            continue
        count += 1
        j = S[0]
        T = [ii for ii in row_order[irow+1:]
             if L[ii][j] != 0]
        for k in S[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] = 0
    return count

####################################


def socle_evaluation(num, g, markings=(), moduli_type=MODULI_ST):
    answer = 1
    G = single_stratum(num, g, dim_form(
        g, len(markings), moduli_type), markings, moduli_type)
    for i in range(1, G.M.nrows()):
        g0 = G.M[i, 0][0]
        psilist = []
        for j in range(1, G.M.ncols()):
            if G.M[i, j][0] > 0:
                psilist.append(G.M[i, j][1])
                if G.M[i, j][0] == 2:
                    psilist.append(G.M[i, j][2])
        n0 = len(psilist)
        dim0 = dim_form(g0, n0, moduli_type)
        kappalist = []
        for j in range(1, dim0+1):
            for k in range(G.M[i, 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


def socle_formula(g, psilist, kappalist, moduli_type=MODULI_ST):
    if moduli_type == MODULI_CT or g == 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)


def setparts_recur(symlist, progress):
    if not symlist:
        return [progress]
    l = []
    for i in Combinations(symlist[1:]):
        j = [symlist[0]]+i
        if progress and j < progress[-1]:
            continue
        cur = 0
        new_symlist = []
        for k in range(len(symlist)):
            if cur < len(j) and symlist[k] == j[cur]:
                cur += 1
            else:
                new_symlist.append(symlist[k])
        l += setparts_recur(new_symlist, progress+[j])
    return l


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


def multi2(g, sigma):
    sigma.sort()
    if sigma[0] == 0:
        total = ZZ.zero()
        for i in range(len(sigma) - 1):
            sigmacopy = sigma[1:]
            if sigmacopy[i] > 0:
                sigmacopy[i] -= 1
                total += multi2(g, sigmacopy)
        return total
    term = factorial(2 * g - 3 + len(sigma))
    term *= ZZ(2 * g - 1).multifactorial(2)
    term /= factorial(2 * g - 1)
    for i in sigma:
        term /= ZZ(2 * i - 1).multifactorial(2)
    return term


def STsum(psilist, kappalist):
    kappalist.sort()
    total = 0
    for i in setparts_with_auts(kappalist):
        total += (-1)**(len(i[0])) * i[1] * \
            STrecur([1 + sum(j) for j in i[0]] + psilist)
    return total*(-1)**(len(kappalist))


def RTsum(g, psilist, kappalist):
    kappalist.sort()
    total = ZZ.zero()
    for i0, i1 in setparts_with_auts(kappalist):
        total += (-1) ** len(i0) * i1 * \
            multi2(g, [1 + sum(j) for j in i0] + psilist)
    return total * (-1)**(len(kappalist))


def CTsum(psilist, kappalist):
    kappalist.sort()
    total = ZZ.zero()
    for i0, i1 in setparts_with_auts(kappalist):
        total += (-1)**len(i0) * i1 * \
            multinomial([1 + sum(j) for j in i0] + psilist)
    return total * (-1)**len(kappalist)


def CTconst(g):
    return abs((ZZ(2) ** (ZZ(2) * g-1)-1)*bernoulli(ZZ(2) * g))/(ZZ(2) ** (ZZ(2) * g-1)*factorial(ZZ(2) * g))


def STrecur(psi):
    return STrecur_calc(tuple(sorted(psi)))


@cached_function
def STrecur_calc(psi):
    """
    INPUT: a sorted tuple of (nonegative ?) integers

    OUTPUT: a rational
    """
    if not psi:
        return ZZ.one()
    s = sum(psi)
    n = len(psi)
    if (s - n) % 3:
        return ZZ.zero()
    if psi[0] == 0:
        if s == 0 and n == 3:
            return ZZ.one()
        total = ZZ.zero()
        psi_end = list(psi[1:])
        for i in range(n - 1):
            psicopy = list(psi_end)
            if psicopy[i] > 0:
                psicopy[i] -= 1
                total += STrecur(psicopy)
        return total

    g = ZZ(s - n) // 3 + 1  # in ZZ
    d = ZZ(psi[-1])
    total = ZZ.zero()

    psicopy = [0, 0, 0, 0] + list(psi)
    psicopy[-1] += 1
    total += (2 * d + 3) / 12 * STrecur(psicopy)

    psicopy = [0, 0, 0] + list(psi)
    total -= (2 * g + n - 1) / 6 * STrecur(psicopy)

    for I in Subsets(list(range(n - 1))):
        psi3 = [0, 0] + [psi[i] for i in I]
        x = STrecur(psi3)
        if x == 0:
            continue
        psi1 = [0, 0] + [psi[i] for i in range(n - 1) if i not in I] + [d + 1]
        psi2 = list(psi1[1:])
        psi2[-1] = d
        total += ((2 * d + 3)*STrecur(psi1) -
                  (2 * g + n - 1)*STrecur(psi2)) * x
    total /= (2 * g + n - 1) * (2 * g + n - 2)
    return total


@cached_function
def good_generator_list(g, r, markings=(), moduli_type=MODULI_ST):
    gens = 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(1, G.M.nrows()):
            g = G.M[i, 0][0]
            codim = 0
            for d in range(1, r + 1):
                if G.M[i, 0][d] != 0:
                    if 3 * d > g:
                        good = False
                        break
                    codim += d*G.M[i, 0][d]
            if not good:
                break
            for j in range(1, G.M.ncols()):
                codim += G.M[i, j][1]
                codim += G.M[i, j][2]
            if codim > 0 and codim >= g:
                good = False
                break
        if good:
            good_gens.append(num)
    return good_gens


@cached_function
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][0][g[0][k]-1]
                          for k in range(v)], [pure_auts[i][1][g[1][k]-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]

#############################################


def FZ_rels(g, r, markings=(), moduli_type=MODULI_ST):
    return span(FZ_matrix(g, r, markings, moduli_type)).basis()


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 = list(range(num_strata(g, r, markings, moduli_type)))
    S2 = good_generator_list(g, r2, markings, moduli_type)
    M = pairing_submatrix(S1, S2, g, r, markings, moduli_type)
    return Matrix(M).kernel().basis()


@cached_function
def kappa_conversion(sigma):
    answer = []
    for spart in setparts_with_auts(list(sigma)):
        coeff = spart[1]
        poly = R(0)
        for part in spart[0]:
            coeff *= factorial(len(part) - 1)
            poly += X**sum(part)
        answer.append([poly, coeff])
    return answer


@cached_function
def kappa_conversion_inverse(sigma):
    answer = []
    for spart in setparts_with_auts(list(sigma)):
        coeff = spart[1]*(-1)**(len(sigma)-len(spart[0]))
        poly = R(0)
        for part in spart[0]:
            poly += X**sum(part)
        answer.append([poly, coeff])
    return answer


@cached_function
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(1, G.M.nrows()):
        genus_vec.append(G.M[i, 0][0])
        kappa_vec.append([])
        for j in range(1, r+1):
            for k in range(G.M[i, 0][j]):
                kappa_vec[-1].append(j)
        kappa_vec[-1] = kappa_conversion(
            tuple(kappa_vec[-1]))
    for choice in itertools.product(*kappa_vec):
        coeff = 1
        GG = Graph(G.M)
        for i in range(1, G.M.nrows()):
            GG.M[i, 0] = genus_vec[i - 1] + choice[i - 1][0]
            coeff *= choice[i - 1][1]
        answer.append((num_of_stratum(GG, g, r, markings, moduli_type), coeff))
    return answer


@cached_function
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(1, G.M.nrows()):
        genus_vec.append(G.M[i, 0][0])
        kappa_vec.append([])
        for j in range(1, r+1):
            for k in range(G.M[i, 0][j]):
                kappa_vec[-1].append(j)
        kappa_vec[-1] = kappa_conversion_inverse(
            tuple(kappa_vec[-1]))
    for choice in itertools.product(*kappa_vec):
        coeff = 1
        GG = Graph(G.M)
        for i in range(1, G.M.nrows()):
            GG.M[i, 0] = genus_vec[i-1] + \
                choice[i-1][0]
            coeff *= choice[i-1][1]
        answer.append((num_of_stratum(GG, g, r, markings, moduli_type), coeff))
    return answer


def convert_vector_to_monomial_basis(vec, g, r, markings=(), moduli_type=MODULI_ST):
    l = len(vec)
    vec2 = [0 for i in range(l)]
    for i in range(l):
        if vec[i] != 0:
            for x in convert_to_monomial_basis(i, g, r, markings, moduli_type):
                vec2[x[0]] += x[1]*vec[i]
    return vec2


def convert_vector_to_pushforward_basis(vec, g, r, markings=(), moduli_type=MODULI_ST):
    l = len(vec)
    vec2 = [0 for i in range(l)]
    for i in range(l):
        if vec[i] != 0:
            for x in convert_to_pushforward_basis(i, g, r, markings, moduli_type):
                vec2[x[0]] += x[1]*vec[i]
    return vec2


def kappa_multiple(vec, which_kappa, g, r, n=0, moduli_type=MODULI_ST):
    vec2 = []
    for x in vec:
        for y in single_kappa_multiple(x[0], which_kappa, g, r, n, moduli_type):
            vec2.append([y[0], x[1]*y[1]])
    vec2 = simplify_sparse(vec2)
    return vec2


def psi_multiple(vec, which_psi, g, r, n=0, moduli_type=MODULI_ST):
    vec2 = []
    for x in vec:
        for y in single_psi_multiple(x[0], which_psi, g, r, n, moduli_type):
            vec2.append([y[0], x[1]*y[1]])
    vec2 = simplify_sparse(vec2)
    return vec2


def insertion_pullback(vec, g, r, n=0, new_mark=1, moduli_type=MODULI_ST):
    vec2 = []
    for x in vec:
        for y in single_insertion_pullback(x[0], g, r, n, new_mark, moduli_type):
            vec2.append([y[0], x[1]*y[1]])
    vec2 = simplify_sparse(vec2)
    return vec2


def insertion_pullback2(vec, g, r, n=0, new_mark=1, moduli_type=MODULI_ST):
    vec2 = []
    for x in vec:
        for y in single_insertion_pullback2(x[0], g, r, n, new_mark, moduli_type):
            vec2.append([y[0], x[1]*y[1]])
    vec2 = simplify_sparse(vec2)
    return vec2


@cached_function
def strata_invariant_lookup(g, r, markings=(), moduli_type=MODULI_ST):
    inv_dict = {}
    L = all_strata(g, r, markings, moduli_type)
    for i in range(len(L)):
        if L[i].invariant not in inv_dict:
            inv_dict[L[i].invariant] = []
        inv_dict[L[i].invariant].append(i)
    return inv_dict


@cached_function
def num_of_stratum(G, g, r, markings=(), moduli_type=MODULI_ST):
    G.compute_invariant()
    L = all_strata(g, r, markings, moduli_type)
    LL = 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) == 1:
        return x[0]
    for i in x:
        if graph_isomorphic(G, L[i]):
            return i
    print("ERROR")
    print((g, r, markings, moduli_type))
    print(G.M)


@cached_function
def single_psi_multiple(num, which_psi, g, r, n=0, moduli_type=MODULI_ST):
    markings = tuple(range(1, n+1))
    G = single_stratum(num, g, r, markings, moduli_type)
    answer = []
    for j in range(1, G.M.ncols()):
        if G.M[0, j] == which_psi:
            good_j = j
            break
    for i in range(1, G.M.nrows()):
        if G.M[i, good_j] != 0:
            deg = 0
            dim_used = 0
            for j in range(1, r+1):
                dim_used += j*G.M[i, 0][j]
            for j in range(1, G.M.ncols()):
                dim_used += G.M[i, j][1] + G.M[i, j][2]
                deg += G.M[i, j][0]
            if dim_used < dim_form(G.M[i, 0][0], deg, moduli_type):
                GG = Graph(G.M)
                GG.M[i, good_j] += X
                answer.append(
                    (num_of_stratum(GG, g, r+1, markings, moduli_type), 1))
            break
    return answer


@cached_function
def single_kappa_multiple(num, which_kappa, g, r, n=0, moduli_type=MODULI_ST):
    markings = tuple(range(1, n+1))
    G = single_stratum(num, g, r, markings, moduli_type)
    answer = []
    for i in range(1, G.M.nrows()):
        deg = 0
        dim_used = 0
        for j in range(1, r+1):
            dim_used += j*G.M[i, 0][j]
        for j in range(1, G.M.ncols()):
            dim_used += G.M[i, j][1] + G.M[i, j][2]
            deg += G.M[i, j][0]
        if dim_used + which_kappa <= dim_form(G.M[i, 0][0], deg, moduli_type):
            GG = Graph(G.M)
            GG.M[i, 0] += X**which_kappa
            answer.append((num_of_stratum(GG, g, r+which_kappa,
                                          markings, moduli_type), 1))
            for j in range(1, r+1):
                if G.M[i, 0][j] > 0:
                    GG = Graph(G.M)
                    GG.M[i, 0] += X**(j+which_kappa)
                    GG.M[i, 0] -= X**j
                    answer.append((num_of_stratum(
                        GG, g, r+which_kappa, markings, moduli_type), -G.M[i, 0][j]))
    return answer

# Only doing this for markings=list(range(1,n+1)) right now.
# Also, this function uses the monomial basis.


def single_kappa_psi_multiple(num, kappa_partition, psi_exps, g, r, n=0, moduli_type=MODULI_ST):
    markings = tuple(range(1, n+1))
    G = single_stratum(num, g, r, markings, moduli_type)
    GG = Graph(G.M)
    for j in range(1, G.M.ncols()):
        if GG.M[0, j] != 0:
            for i in range(1, G.M.nrows()):
                if GG.M[i, j] != 0:
                    GG.M[i, j] += psi_exps[GG.M[0, j][0]-1]*X
                    break
    rnew = r+sum(kappa_partition)+sum(psi_exps)
    answer = []
    kappa_options = [list(range(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], 0] += X**(kappa_partition[i])
        is_bad = False
        for i in range(1, GGG.M.nrows()):
            deg = 0
            dim_used = 0
            for j in range(1, rnew+1):
                dim_used += j*GGG.M[i, 0][j]
            for j in range(1, GGG.M.ncols()):
                dim_used += GGG.M[i, j][1] + GGG.M[i, j][2]
                deg += GGG.M[i, j][0]
            if dim_used > dim_form(GGG.M[i, 0][0], deg, moduli_type):
                is_bad = True
                break
        if is_bad:
            continue
        answer.append(
            (num_of_stratum(GGG, g, rnew, markings, moduli_type), 1))
    return answer


@cached_function
def single_insertion_pullback(num, g, r, n=0, new_mark=1, moduli_type=MODULI_ST):
    markings = tuple(range(1, n+1))
    new_markings = tuple(range(1, n+2))
    G = single_stratum(num, g, r, markings, moduli_type)
    answer = []
    for i in range(1, G.M.nrows()):
        GG = Graph(G.M)
        for j in range(1, G.M.ncols()):
            if GG.M[0, j][0] >= new_mark:
                GG.M[0, j] += 1
        GG.add_edge(i, 0, new_mark)
        answer.append(
            (num_of_stratum(GG, g, r, new_markings, moduli_type), 1))
        for j in range(1, r+1):
            for k in range(GG.M[i, 0][j]):
                GGG = Graph(GG.M)
                GGG.M[i, 0] -= X**j
                GGG.M[i, -1] += j*X
                answer.append(
                    (num_of_stratum(GGG, g, r, new_markings, moduli_type), -1))
        if moduli_type <= MODULI_SM:
            continue
        for j in range(1, G.M.ncols()):
            if G.M[i, j][0] == 1:
                if G.M[i, j][1] >= 1:
                    x = G.M[i, j][1]
                    GGG = Graph(GG.M)
                    row1 = [GG.M[i, k] for k in range(GG.M.ncols())]
                    row2 = [0 for k in range(GG.M.ncols())]
                    row1[j] = 0
                    row1[-1] = 0
                    row2[j] = 1
                    row2[-1] = 1
                    GGG.split_vertex(i, row1, row2)
                    GGG.M[-2, -
                          1] += (x-1)*X
                    answer.append(
                        (num_of_stratum(GGG, g, r, new_markings, moduli_type), -1))
            if G.M[i, j][0] == 2:
                if G.M[i, j][1] >= 1 or G.M[i, j][2] >= 1:
                    x = G.M[i, j][1]
                    y = G.M[i, j][2]
                    row1 = [GG.M[i, k] for k in range(GG.M.ncols())]
                    row2 = [0 for k in range(GG.M.ncols())]
                    row1[j] = 0
                    row1[-1] = 0
                    row2[j] = 1
                    row2[-1] = 1
                    if y >= 1:
                        row1[j] = 1 + x*X
                        GGG = Graph(GG.M)
                        GGG.split_vertex(i, row1, row2)
                        GGG.M[-2, -
                              1] += (y-1)*X
                        answer.append(
                            (num_of_stratum(GGG, g, r, new_markings, moduli_type), -1))
                    if x >= 1:
                        row1[j] = 1 + y*X
                        GGG = Graph(GG.M)
                        GGG.split_vertex(i, row1, row2)
                        GGG.M[-2, -
                              1] += (x-1)*X
                        answer.append(
                            (num_of_stratum(GGG, g, r, new_markings, moduli_type), -1))
    return answer


@cached_function
def single_insertion_pullback2(num, g, r, n=0, new_mark=1, moduli_type=MODULI_ST):
    markings = tuple(range(1, n+1))
    new_markings = tuple(range(1, n+2))
    G = single_stratum(num, g, r, markings, MODULI_SMALL)
    answer = []
    for i in range(1, G.M.nrows()):
        GG = Graph(G.M)
        for j in range(1, G.M.ncols()):
            if GG.M[0, j][0] >= new_mark:
                GG.M[0, j] += 1
        GG.add_edge(i, 0, new_mark)
        answer.append(
            (num_of_stratum(GG, g, r, new_markings, moduli_type), 1))
        for j in range(1, r+1):
            for k in range(GG.M[i, 0][j]):
                GGG = Graph(GG.M)
                GGG.M[i, 0] -= X**j
                GGG.M[i, -1] += j*X
                answer.append(
                    (num_of_stratum(GGG, g, r, new_markings, moduli_type), -1))
        if moduli_type <= MODULI_SM:
            continue
        for j in range(1, G.M.ncols()):
            if G.M[i, j][0] == 1:
                if G.M[i, j][1] >= 1:
                    x = G.M[i, j][1]
                    GGG = Graph(GG.M)
                    row1 = [GG.M[i, k] for k in range(GG.M.ncols())]
                    row2 = [0 for k in range(GG.M.ncols())]
                    row1[j] = 0
                    row1[-1] = 0
                    row2[j] = 1
                    row2[-1] = 1
                    GGG.split_vertex(i, row1, row2)
                    GGG.M[-2, -
                          1] += (x-1)*X
                    answer.append(
                        (num_of_stratum(GGG, g, r, new_markings, moduli_type), -1))
    return answer


def num_new_rels(g, r, n=0, moduli_type=MODULI_ST):
    return len(choose_basic_rels(g, r, n, moduli_type))


@cached_function
def choose_basic_rels(g, r, n=0, moduli_type=MODULI_ST):
    if 3 * r < g + n + 1:
        return []
    sym_ngen = num_strata(g, r, tuple([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, 1]] for i in range(sym_ngen)]
    else:
        sym_possible_rels = possibly_new_FZ(g, r, n, moduli_type)
    if not sym_possible_rels:
        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[0]] = x[1]
    if nrels > 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 = 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[0]] = x[1]
        row_order.append(nrels)
        nrels += 1
        if compute_rank_sparse(D, row_order, col_order) > previous_rank:
            answer.append(unsymmetrize_vec(sym_possible_rels[j], g, r, tuple(
                range(1, n+1)), moduli_type))
            previous_rank += 1
        dprint("rank %s", previous_rank)
    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 > -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


def recursive_betti(g, r, markings=(), moduli_type=MODULI_ST):
    """
    EXAMPLES::

        sage: from admcycles.DR import recursive_betti
        sage: recursive_betti(2, 3)  # not tested
        ?
    """
    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 0
    ngen = num_strata(g, r, markings, moduli_type)
    dprint("%s gens", ngen)
    relations = []
    partial_sym_map = partial_symmetrize_map(g, r, markings, moduli_type)
    for rel in choose_basic_rels(g, r, n, moduli_type):
        rel2 = []
        for x in rel:
            rel2.append([partial_sym_map[x[0]], x[1]])
        rel2 = simplify_sparse(rel2)
        relations.append(rel2)
    for rel in interior_derived_rels(g, r, n, moduli_type):
        rel2 = []
        for x in rel:
            rel2.append([partial_sym_map[x[0]], x[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(1, r):
            strata = all_strata(g, r0, markings, moduli_type)
            for G in strata:
                vertex_orbits = graph_count_automorphisms(G, True)
                for i in [orbit[0] for orbit in vertex_orbits]:
                    good = True
                    for j in range(G.M.ncols()):
                        if R(G.M[i, j][0]) != G.M[i, j]:
                            good = False
                            break
                    if good:
                        g2 = G.M[i, 0][0]
                        if 3 * (r - r0) < g2 + 1:
                            continue
                        d = G.degree(i)
                        if dim_form(g2, d, moduli_type) < r-r0:
                            continue
                        strata2 = all_strata(
                            g2, r-r0, tuple(range(1, d+1)), moduli_type)
                        which_gen_list = [
                            -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(choose_basic_rels(
                            g2, r-r0, d, moduli_type))
                        rel_list += interior_derived_rels(g2,
                                                          r-r0, d, moduli_type)
                        for rel0 in rel_list:
                            relation = []
                            for x in rel0:
                                num = x[0]
                                if which_gen_list[num] != -1:
                                    relation.append(
                                        [which_gen_list[num], x[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_duplicates(relations)
    nrels = len(relations)
    dprint("%s gens, %s distinct rels", ngen, nrels)
    dsave("sparse-%s-distinct-rels", nrels)
    rank = 0
    D = {}
    for i, reli in enumerate(relations):
        for x in reli:
            D[i, x[0]] = x[1]
    if nrels:
        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


def subsequences(seq, l):
    """
    Return all subsequences of length ``l`` a list ``seq``.

    EXAMPLES::

        sage: from admcycles.DR import subsequences
        sage: subsequences([2,3,5,7],2)
        [[2, 3], [2, 5], [2, 7], [3, 5], [3, 7], [5, 7]]
        sage: subsequences([4,3,2,1],2)
        [[4, 3], [4, 2], [4, 1], [3, 2], [3, 1], [2, 1]]
    """
    S = tuple(range(len(seq)))
    return [[seq[i] for i in subset] for subset in Subsets(S, l)]


@cached_function
def pullback_derived_rels(g, r, n=0, moduli_type=MODULI_ST):
    if r == 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 = choose_basic_rels(g, r, n0, moduli_type)
            for rel in basic_rels:
                for vec in subsequences(list(range(1, n+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 = 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(list(range(1, n0+k-1 + 1)), k-1):
                    for vec2 in subsequences(list(range(1, n+1)), n-n0-k+1):
                        rel2 = copy(rel)
                        for i in range(k-1):
                            rel2 = insertion_pullback(
                                rel2, g, r, n0+i, vec[i], MODULI_SMALL)
                        rel2 = insertion_pullback2(
                            rel2, g, r, n0+k-1, vec2[0], moduli_type)
                        for i in range(n-n0-k):
                            rel2 = insertion_pullback(
                                rel2, g, r, n0+k+i, vec2[i+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


@cached_function
def interior_derived_rels(g, r, n=0, moduli_type=MODULI_ST):
    dprint("Start interior_derived (%s,%s,%s,%s): %s", g,
           r, n, moduli_type, floor(get_memory_usage()))
    answer = copy(pullback_derived_rels(g, r, n, moduli_type))
    for r0 in range(r):
        pullback_rels = copy(choose_basic_rels(g, r0, n, moduli_type))
        pullback_rels += pullback_derived_rels(g, r0, n, moduli_type)
        for rel in pullback_rels:
            for i in range(r-r0+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+1, g, rcur, n, moduli_type)
                                rcur += 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


@cached_function
def symmetrize_map(g, r, markings=(), moduli_type=MODULI_ST):
    markings2 = tuple([1 for i in markings])
    gens = all_strata(g, r, markings, moduli_type)
    map = []
    for G in gens:
        GG = Graph(G.M)
        for i in range(1, GG.M.ncols()):
            if GG.M[0, i][0] > 0:
                GG.M[0, i] = R(1)
        map.append(num_of_stratum(GG, g, r, markings2, moduli_type))
    return map


@cached_function
def partial_symmetrize_map(g, r, markings=(), moduli_type=MODULI_ST):
    markings1 = tuple(range(1, len(markings)+1))
    gens = all_strata(g, r, markings1, moduli_type)
    map = []
    for G in gens:
        GG = Graph(G.M)
        for i in range(1, GG.M.ncols()):
            if GG.M[0, i][0] > 0:
                GG.M[0, i] = R(markings[GG.M[0, i][0]-1])
        map.append(num_of_stratum(GG, g, r, markings, moduli_type))
    return map


@cached_function
def unsymmetrize_map(g, r, markings=(), moduli_type=MODULI_ST):
    markings2 = tuple([1 for i in markings])
    sym_map = 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


def unsymmetrize_vec(vec, g, r, markings=(), moduli_type=MODULI_ST):
    unsym_map = unsymmetrize_map(g, r, markings, moduli_type)
    vec2 = []
    for x in vec:
        aut = len(unsym_map[x[0]])
        for j in unsym_map[x[0]]:
            vec2.append([j, QQ((x[1], aut))])
    vec2 = simplify_sparse(vec2)
    return vec2


def derived_rels(g, r, n=0, moduli_type=MODULI_ST):
    dprint("Start derived (%s,%s,%s,%s): %s", g, r,
           n, moduli_type, floor(get_memory_usage()))
    markings = tuple([1 for _ in range(n)])
    # generators = all_strata(g, r, markings, moduli_type)
    sym_map = symmetrize_map(g, r, tuple(range(1, n+1)), moduli_type)
    answer = []
    for rel in interior_derived_rels(g, r, n, moduli_type):
        rel2 = []
        for x in rel:
            rel2.append([sym_map[x[0]], x[1]])
        rel2 = simplify_sparse(rel2)
        answer.append(rel2)
    if moduli_type <= MODULI_SM:
        return answer
    for r0 in range(1, r):
        strata = all_strata(g, r0, markings, moduli_type)
        for G in strata:
            vertex_orbits = graph_count_automorphisms(G, True)
            for i in [orbit[0] for orbit in vertex_orbits]:
                good = True
                for j in range(G.M.ncols()):
                    if R(G.M[i, j][0]) != G.M[i, j]:
                        good = False
                        break
                if good:
                    g2 = G.M[i, 0][0]
                    if 3 * (r - r0) < g2 + 1:
                        continue
                    d = G.degree(i)
                    if dim_form(g2, d, moduli_type) < r-r0:
                        continue
                    strata2 = all_strata(
                        g2, r-r0, tuple(range(1, d+1)), moduli_type)
                    which_gen_list = [
                        -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(choose_basic_rels(
                        g2, r-r0, d, moduli_type))
                    rel_list += interior_derived_rels(g2, r-r0, d, moduli_type)
                    for rel0 in rel_list:
                        relation = []
                        for x0, x1 in rel0:
                            if which_gen_list[x0] != -1:
                                relation.append([which_gen_list[x0], x1])
                        relation = simplify_sparse(relation)
                        answer.append(relation)
    answer = remove_duplicates(answer)
    dprint("End derived (%s,%s,%s,%s): %s", g, r, n,
           moduli_type, floor(get_memory_usage()))
    return answer


def remove_duplicates(L):
    """
    Remove duplicate elements in a list ``L``.

    One cannot use ``set(L)`` because the elements of ``L`` are not hashable.

    INPUT:

    - ``L`` -- a list

    OUTPUT:

    a list

    EXAMPLES::

       sage: from admcycles.DR import remove_duplicates
       sage: remove_duplicates([4,7,6,4,3,3,4,2,2,1])
       [1, 2, 3, 4, 6, 7]
    """
    if not L:
        return L
    L.sort()
    LL = [L[0]]
    for i, Li in enumerate(L[1:]):
        if Li != L[i]:
            LL.append(Li)
    return LL


def simplify_sparse(vec):
    vec.sort()
    vec2 = []
    last_index = None
    for x in vec:
        if x[0] == last_index:
            if vec2[-1][1] == -x[1]:
                vec2 = vec2[:-1]
                last_index = None
            else:
                vec2[-1][1] += x[1]
        else:
            vec2.append(x)
            last_index = x[0]
    return vec2


def compute_rank_sparse(D, row_order, col_order):
    count = 0
    nrows = len(row_order)
    ncols = len(col_order)
    row_order_rank = [-1 for i in range(nrows)]
    col_order_rank = [-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:
        row_contents[x[0]].add(x[1])
        col_contents[x[1]].add(x[0])

    for i in row_order:
        S = []
        for j in row_contents[i]:
            S.append(j)
        if not S:
            continue
        count += 1
        S.sort(key=lambda x: col_order_rank[x])
        j = S[0]
        T = []
        for ii in col_contents[j]:
            if row_order_rank[ii] > row_order_rank[i]:
                T.append(ii)
        for k in S[1:]:
            rat = QQ((D[i, k], D[i, j]))
            for ii in T:
                if (ii, k) not in D:
                    D[ii, k] = 0
                    row_contents[ii].add(k)
                    col_contents[k].add(ii)
                D[ii, k] -= rat*D[ii, j]
                if D[ii, k] == 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

##########################################


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(1, nr):
        if G.M[i, 0] != R(G.M[i, 0][0]):
            return True
        for j in range(1, nc):
            if G.M[i, j][0] == 2:
                return True
            if G.M[0, j] != 0 and G.M[i, j] != 0:
                marked_vertices.append(i)
    for ii in range(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(1, nc):
                    if G.M[i, j] != 0:
                        for i2 in range(1, nr):
                            if G.M[i2, j] != 0 and i2 not in S:
                                S.add(i2)
                                did_something = True
        if len(S) < nr-1:
            return True
    return False


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(1, nr):
        for i2 in range(1, i1):
            found_edge = False
            for j in range(1, nc):
                if G.M[i1, j] != 0 and G.M[i2, j] != 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(1, nc):
                        if G.M[i3, j] != 0:
                            for i4 in range(1, nr):
                                if G.M[i4, j] != 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


def DR_coeff_is_known(num, g, r, markings=(), moduli_type=MODULI_ST):
    return r <= 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) > 3:
        return False
    if len(nonsep_pairs) == 3:
        i1 = nonsep_pairs[0][0]
        i2 = nonsep_pairs[0][1]
        i3 = nonsep_pairs[2][1]
        jlist = []
        for j in range(1, nc):
            if len([1 for i in [i1, i2, i3] if G.M[i, j] != 0]) == 2:
                jlist.append(j)
        if len(jlist) > 3:
            return False
        for j in jlist:
            for i in [i1, i2, i3]:
                if G.M[i, j][1] != 0:
                    return False

    return True

# this stuff is old, keeping it around for now just in case
    for j in range(1, nc):
        ilist = []
        for i in range(1, nr):
            if G.M[i, j] != 0:
                ilist.append(i)
        if len(ilist) == 1:
            continue
        ii1 = ilist[0]
        ii2 = ilist[1]
        jlist = []
        for jj in range(1, nc):
            if G.M[ii1, jj] != 0 and G.M[ii2, jj] != 0:
                jlist.append(jj)
        if len(jlist) == 1 or jlist[0] != j:
            continue
        count = len(jlist)
        for ii in [ii1, ii2]:
            for jj in jlist:
                count += G.M[ii, jj][1]
        if count > 3:
            return False
    return True

# next function doesn't work on arbitrary graphs yet, see above function


def DR_coeff(num, g, r, n=0, dvector=(), moduli_type=MODULI_ST):
    markings = tuple(range(1, n+1))
    G = single_stratum(num, g, r, markings, moduli_type)
    nr = G.M.nrows()
    nc = G.M.ncols()
    answer = QQ((1, autom_count(num, g, r, markings, moduli_type)))

    def balance(ilist):
        bal = [0 for i1 in ilist]
        for j1 in range(1, nc):
            if G.M[0, j1] != 0:
                for i3 in range(1, nr):
                    if G.M[i3, j1] != 0:
                        S = set([i3])
                        did_something = True
                        while did_something:
                            did_something = False
                            for i4 in tuple(S):
                                for j2 in range(1, nc):
                                    if G.M[i4, j2] != 0:
                                        for i5 in range(1, nr):
                                            if G.M[i5, j2] != 0 and i5 not in S and (i4 not in ilist or i5 not in ilist):
                                                S.add(i5)
                                                did_something = True
                        for kk, ikk in enumerate(ilist):
                            if ikk in S:
                                bal[kk] += dvector[G.M[0, j1][0] - 1]
        return bal

    def poly4(a, b, c):
        return (ZZ.one() / 420)*(-15 * (a**8 + b**8 + c**8)+40 * (a**7 * b+a**7 * c+b**7 * a+b**7 * c+c**7 * a+c**7 * b)-28 * (a**6 * b**2 + a**6 * c**2 + b**6 * a**2 + b**6 * c**2 + c**6 * a**2 + c**6 * b**2)-112 * (a**6 * b*c+b**6 * a*c+c**6 * a*b)+84 * (a**5 * b**2 * c+a**5 * c**2 * b+b**5 * a**2 * c+b**5 * c**2 * a+c**5 * a**2 * b+c**5 * b**2 * a)-70 * (a**4 * b**2 * c**2 + b**4 * a**2 * c**2 + c**4 * a**2 * b**2))

    nonsep_pairs = find_nonsep_pairs(num, g, r, markings, moduli_type)
    if len(nonsep_pairs) == 4:
        cycle = [nonsep_pairs[0][0], nonsep_pairs[0]
                 [1], 0, 0]
        for i in range(1, 4):
            for j in range(2):
                if nonsep_pairs[i][j] == cycle[0]:
                    cycle[3] = nonsep_pairs[i][1 - j]
                if nonsep_pairs[i][j] == cycle[1]:
                    cycle[2] = nonsep_pairs[i][1 - j]
        dbal = balance(cycle)
        a = dbal[0]
        b = dbal[0] + dbal[1]
        c = dbal[0] + dbal[1] + dbal[2]
        answer *= poly4(a, b, c)
    elif len(nonsep_pairs) == 3:
        i1 = nonsep_pairs[0][0]
        i2 = nonsep_pairs[0][1]
        i3 = nonsep_pairs[2][1]
        cycle = [i1, i2, i3]
        cycle4 = cycle + [cycle[0]]
        dbal = balance(cycle)
        dbal4 = dbal + [dbal[0]]
        done = False
        for k in range(3):
            jlist = [j for j in range(1, nc) if (
                G.M[cycle4[k], j] != 0 and G.M[cycle4[k+1], j] != 0)]
            if len(jlist) == 2:
                answer *= -ZZ.one() / 6 * \
                    poly4(0, dbal4[k], -dbal4[k+1])
                done = True
            elif len(jlist) == 1:
                if G.M[cycle4[k], jlist[0]][1] + G.M[cycle4[k+1], jlist[0]][1] > 0:
                    answer *= -ZZ.one() / 2 * \
                        poly4(0, dbal4[k], -dbal4[k+1])
                    done = True
        if not done:
            answer *= (dbal[0]**6 + dbal[1]**6 + dbal[2]**6 -
                       10 * dbal[0]**2 * dbal[1]**2 * dbal[2]**2)/ZZ(30)

    for j in range(1, nc):
        ilist = []
        for i in range(1, nr):
            if G.M[i, j] != 0:
                ilist.append(i)
        if len(ilist) == 1:
            x = G.M[ilist[0], j][1]
            answer /= factorial(x)
            answer *= dvector[G.M[0, j][0] -
                              1]**(ZZ(2) * x)
            continue
        if ilist in nonsep_pairs:
            continue
        ii1 = ilist[0]
        ii2 = ilist[1]
        jlist = []
        for jj in range(1, nc):
            if G.M[ii1, jj] != 0 and G.M[ii2, jj] != 0:
                jlist.append(jj)
        if len(jlist) == 1:
            x1 = G.M[ii1, j][1]
            x2 = G.M[ii2, j][1]
            answer *= -1
            answer /= factorial(x1)
            answer /= factorial(x2)
            answer /= x1+x2+1
            answer *= balance([ii1, ii2])[0]**(ZZ(2) *
                                               (x1+x2+1))
            continue
        elif len(jlist) == 2:
            if jlist[0] != j:
                continue
            xvec = [G.M[ii1, jlist[0]][1], G.M[ii1, jlist[1]][1],
                    G.M[ii2, jlist[0]][1], G.M[ii2, jlist[1]][1]]
            x = sum(xvec)
            if x == 0:
                answer *= -ZZ.one() / 6
                answer *= balance([ii1, ii2])[0]**4
            elif x == 1:
                answer *= -ZZ.one() / 30
                answer *= balance([ii1, ii2])[0]**6
            elif x == 2:
                if xvec in [[2, 0, 0, 0], [0, 2, 0, 0], [0, 0, 2, 0], [0, 0, 0, 2]]:
                    answer *= -ZZ.one() / 168
                elif xvec in [[1, 1, 0, 0], [0, 0, 1, 1], [1, 0, 0, 1], [0, 1, 1, 0]]:
                    answer *= -ZZ.one() / 280
                elif xvec in [[1, 0, 1, 0], [0, 1, 0, 1]]:
                    answer *= -ZZ.one() / 84
                answer *= balance([ii1, ii2])[0]**8
            continue
        elif len(jlist) == 3:
            if jlist[0] != j:
                continue
            xvec = [G.M[ii1, jlist[0]][1], G.M[ii1, jlist[1]][1], G.M[ii1, jlist[2]]
                    [1], G.M[ii2, jlist[0]][1], G.M[ii2, jlist[1]][1], G.M[ii2, jlist[2]][1]]
            x = sum(xvec)
            if x == 0:
                answer *= -ZZ.one() / 90
                answer *= balance([ii1, ii2])[0]**6
            elif x == 1:
                answer *= -ZZ.one() / 840
                answer *= balance([ii1, ii2])[0]**8
            continue
        elif len(jlist) == 4:
            if jlist[0] != j:
                continue
            answer *= -ZZ.one() / 2520
            answer *= balance([ii1, ii2])[0]**8
    return answer


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("---------------------")


def reduce_with_rels(B, vec):
    vec2 = copy(vec)
    for row in B:
        for i in range(len(row)):
            if row[i] != 0:
                if vec2[i] != 0:
                    vec2 -= QQ((vec2[i], row[i]))*row
                break
    return vec2


def DR_coeff_setup(num, g, r, n=0, dvector=(), kval=0, moduli_type=MODULI_ST):
    markings = tuple(range(1, n+1))
    G = single_stratum(num, g, r, markings, moduli_type)
    nr = G.M.nrows()
    nc = G.M.ncols()
    edge_list = []
    exp_list = []
    scalar_factor = ZZ.one() / \
        autom_count(num, g, r, markings, moduli_type)
    for i in range(1, nr):
        for j in range(1, G.M[i, 0].degree()+1):
            scalar_factor /= factorial(G.M[i, 0][j])
            scalar_factor /= factorial(j)**G.M[i, 0][j]
            scalar_factor *= (-1)**G.M[i, 0][j]
            scalar_factor *= (kval**2)**(j * G.M[i, 0][j])
    given_weights = [-kval*(2*G.M[i+1, 0][0] - 2 + sum([G.M[i+1][j][0]
                                                        for j in range(1, nc)])) for i in range(nr-1)]
    for j in range(1, nc):
        ilist = [i for i in range(1, nr)
                 if G.M[i, j] != 0]
        if G.M[0, j] == 0:
            if len(ilist) == 1:
                i1 = ilist[0]
                i2 = ilist[0]
                exp1 = G.M[i1, j][1]
                exp2 = G.M[i1, j][2]
            else:
                i1 = ilist[0]
                i2 = ilist[1]
                exp1 = G.M[i1, j][1]
                exp2 = G.M[i2, j][1]
            edge_list.append([i1-1, i2-1])
            exp_list.append(exp1 + exp2 + 1)
            scalar_factor /= - \
                factorial(exp1)*factorial(exp2)*(exp1+exp2+1)
        else:
            exp1 = G.M[ilist[0], j][1]
            scalar_factor *= dvector[G.M[0, j][0] -
                                     1]**(ZZ(2) * exp1)/factorial(exp1)
            given_weights[ilist[0] -
                          1] += dvector[G.M[0, j][0] - 1]
    return edge_list, exp_list, given_weights, scalar_factor


def DR_coeff_new(num, g, r, n=0, dvector=(), kval=0, moduli_type=MODULI_ST):
    markings = tuple(range(1, n+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])/ZZ(2)) + g*abs(kval)
    h0 = nc - nr - n + 1
    deg = 2 * sum(exp_list)
    mrange = list(range(m0 + 1, m0 + deg + 2))
    mvalues = []
    for m in mrange:
        total = 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][0]] += weight_data[i]
                vertex_weights[edge_list[i][1]] -= weight_data[i]
            if any(i % m for i in vertex_weights):
                continue
            term = 1
            for i in range(len(edge_list)):
                term *= weight_data[i]**(ZZ(2) * exp_list[i])
            total += term
        mvalues.append(QQ((total, m**h0)))
    mpoly = interpolate(mrange, mvalues)
    return mpoly.subs(X=0)*scalar_factor


def DR_coeff_setup_m(m, num, g, r, n=0, dvector=(), kval=0, moduli_type=MODULI_ST):
    markings = tuple(range(1, n+1))
    G = single_stratum(num, g, r, markings, moduli_type)
    nr = G.M.nrows()
    nc = G.M.ncols()
    edge_list = []
    exp_list = []
    scalar_factor = ZZ.one() / \
        autom_count(num, g, r, markings, moduli_type)
    for i in range(1, nr):
        for j in range(1, G.M[i, 0].degree()+1):
            scalar_factor /= factorial(G.M[i, 0][j])
            scalar_factor /= factorial(j)**G.M[i, 0][j]
            scalar_factor *= (-1)**G.M[i, 0][j]
            scalar_factor *= (kval**2 - kval*m + m**2 /
                              ZZ(6))**(j*G.M[i, 0][j])
    given_weights = [-kval*(ZZ(2) * G.M[i+1, 0][0] - ZZ(2) + sum(
        [G.M[i+1][j][0] for j in range(1, nc)])) for i in range(nr-1)]
    for j in range(1, nc):
        ilist = [i for i in range(1, nr)
                 if G.M[i, j] != 0]
        if G.M[0, j] == 0:
            if len(ilist) == 1:
                i1 = ilist[0]
                i2 = ilist[0]
                exp1 = G.M[i1, j][1]
                exp2 = G.M[i1, j][2]
            else:
                i1 = ilist[0]
                i2 = ilist[1]
                exp1 = G.M[i1, j][1]
                exp2 = G.M[i2, j][1]
            edge_list.append([i1-1, i2-1])
            exp_list.append(exp1 + exp2 + 1)
            scalar_factor /= - \
                factorial(exp1)*factorial(exp2)*(exp1+exp2+1)
        else:
            exp1 = G.M[ilist[0], j][1]
            dval = dvector[G.M[0, j][0] - 1]
            if dval < 0:
                dval = dval + m
            scalar_factor *= (dval**2 - dval*m + m**2 /
                              ZZ(6))**(exp1)/factorial(exp1)
            given_weights[ilist[0] -
                          1] += dvector[G.M[0, j][0] - 1]
    return edge_list, exp_list, given_weights, scalar_factor


def DR_coeff_m(m, num, g, r, n=0, dvector=(), kval=0, moduli_type=MODULI_ST):
    markings = tuple(range(1, n+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 + 1
    total = ZZ.zero()
    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][0]] += weight_data[i]
            vertex_weights[edge_list[i][1]] -= weight_data[i]
        if any(i % m for i in vertex_weights):
            continue
        term = 1
        for i in range(len(edge_list)):
            dval = weight_data[i]
            term *= (dval**2 - dval*m + m**2 / ZZ(6))**exp_list[i]
        total += term
    total /= m**h0
    total *= scalar_factor
    return total


def interpolate(A, B):
    X = SR.var('X')
    l = len(A)
    poly = 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()[0]
    for i in range(l):
        poly += (QQ(-ker[i])/QQ(ker[-1]))*X**i
    return poly


def DR_compute(g, r, n=0, dvector=(), kval=0, moduli_type=MODULI_ST):
    answer = []
    markings = tuple(range(1, n+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)/ZZ(2) ** r


def DR_compute_m(m, g, r, n=0, dvector=(), kval=0, moduli_type=MODULI_ST):
    answer = []
    markings = tuple(range(1, n+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)


def DR_sparse(g, r, n=0, dvector=(), kval=0, moduli_type=MODULI_ST):
    """
    Same as :func:`DR_compute` except that it returns the answer as a
    sparse vector.

    EXAMPLES::

        sage: from admcycles.DR import DR_sparse
        sage: DR_sparse(1,1,2,(2,-2))
        [[1, 2], [2, 2], [4, -1/24]]
    """
    vec = DR_compute(g, r, n, dvector, kval, moduli_type)
    return [[i, veci] for i, veci in enumerate(vec) if veci != 0]


def DR_psi_check(g, n, dvector, which_psi):
    vec = DR_sparse(g, g, n, dvector, 0)
    r = g
    while r < 3 * g - 3 + n:
        vec = psi_multiple(vec, which_psi, g, r, n)
        r += 1
    total = ZZ.zero()
    for a, b in vec:
        total += b * socle_evaluation(a, g, tuple(range(1, n + 1)))
    return total / 2**g


def DR_reduced(g, dvector=()):
    """
    Same as :func:`DR_compute` except that it only requires two arguments
    (g and dvector) and simplifies the answer using the 3-spin
    tautological relations.

    EXAMPLES::

        sage: from admcycles.DR import DR_reduced
        sage: DR_reduced(1,(2,-2))
        (0, 0, 0, 4, 1/8)
    """
    g = ZZ(g)
    n = len(dvector)
    r = g
    kval = ZZ(sum(dvector) / (2 * g - 2 + n))
    vec = DR_compute(g, r, n, dvector, kval)
    vec = vector(QQ, (QQ(a) for a in vec))
    rels = FZ_rels(g, r, tuple(range(1, n + 1)))
    return reduce_with_rels(rels, vec)


def list_strata(g, r, n=0, moduli_type=MODULI_ST):
    """
    Displays the list of all strata.

    EXAMPLES::

        sage: from admcycles.DR import list_strata
        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]
        ------------------------------
    """
    L = all_strata(g, r, tuple(range(1, n + 1)), moduli_type)
    for i, Li in enumerate(L):
        print("generator %s" % i)
        print(Li.M)
        print("-" * 30)