def statistic(D):
    H = D.heights()
    T = [H[i+1] for i in range(len(D)-1) if D[i] != D[i+1]]
    return T[(len(T)-1)//2]

def statistic(D):
    H = D.heights()
    T := [H[i+1] for i in range(len(D)-1) if D[i] != D[i+1]]
    return T[(len(T)-1)//2]

@cached_function
def _level_tuples(n):
	return [tuple(B.tamari_sorting_tuple()[0]) for B in BinaryTrees(n)]

def statistic(tree):
	n = tree.node_number() 
	if n <= 1:                       # C_0 = C_1 = 1; the tree is alone on its level
	  return 1
	target = tuple(tree.tamari_sorting_tuple()[0])
	level = _level_tuples(n)
	return sum(1 for row in level if all(a <= b for a, b in zip(row, target)))

def statistic(tree):
    return @cached_function
def _level_tuples(n):
	return [tuple(B.tamari_sorting_tuple()[0]) for B in BinaryTrees(n)]

def statistic(tree):
	n = tree.node_number() 
	if n <= 1:                       # C_0 = C_1 = 1; the tree is alone on its level
	  return 1
	target = tuplen(tree.tamari_smaller())orting_tuple()[0])
	level = _level_tuples(n)
	return sum(1 for row in level if all(a <= b for a, b in zip(row, target)))

_PATTERN = BinaryTree('[[.,.],[.,[[.,.],.]]]')

def _occ(t, p):
    # A leaf '.' in the pattern is a wildcard; internal nodes must match
    # the left/right branching exactly and contiguously.
    if p.is_empty():
        return True
    if t.is_empty():
        return False
    return _occ(t[0], p[0]) and _occ(t[1], p[1])

def statistic(t):
    # Count the nodes of t at which the pattern sits contiguously.
    if t.is_empty():
        return 0
    return (1 if _occ(t, _PATTERN) else 0) + statistic(t[0]) + statistic(t[1])

# (* in Mathematica *)
# tree = {{{{}, {}}, {{}, {}}}, {{{}, {}}, {{}, {}}}};
# Count[tree, {{{___}, {___}}, {{___}, {{{___}, {___}}, {___}}}}, {0, Infinity}]

#CodeLanguage:_PATTERN = BinaryTree('[[.,.],[.,[[.,.],.]]]')

def _occ(t, p):
    # A leaf '.' in the pattern is a wildcard; internal nodes must match
    # the left/right branching exactly and contiguously.
    if p.is_empty():
        return True
    if t.is_empty():
        return False
    return _occ(t[0], p[0]) and _occ(t[1], p[1])

def statistic(t):
    # Count the nodes of t at which the pattern sits contiguously.
    if t.is_empty():
        return 0
    return (1 if _occ(t, _PATTERN) else 0) + statistic(t[0]) + statistic(t[1])

# (* in Mathematica
 *)
# tree = {{{{}, {}}, {{}, {}}}, {{{}, {}}, {{}, {}}}};
# Count[tree, {{{___}, {___}}, {{___}, {{{___}, {___}}, {___}}}}, {0, Infinity}]

_PATTERN = BinaryTree('[.,[[[[.,.],.],.],.]]')

def _occ(t, p):
    # A leaf '.' in the pattern is a wildcard; internal nodes must match
    # the left/right branching exactly and contiguously.
    if p.is_empty():
        return True
    if t.is_empty():
        return False
    return _occ(t[0], p[0]) and _occ(t[1], p[1])

def statistic(t):
    # Count the nodes of t at which the pattern sits contiguously.
    if t.is_empty():
        return 0
    return (1 if _occ(t, _PATTERN) else 0) + statistic(t[0]) + statistic(t[1])

# (* in Mathematica *)
# tree = {{{{}, {}}, {{}, {}}}, {{{}, {}}, {{}, {}}}};
# Count[tree, {{___}, {{{{{___}, {___}}, {___}}, {___}}, {___}}}, {0, Infinity}]

_PATTERN = BinaryTree('[.,[[[[.,.],.],.],.]]')

def _occ(t, p):
    # A leaf '.' in the pattern is a wildcard; internal nodes must match
    # the left/right branching exactly and contiguously.
    if p.is_empty():
        return True
    if t.is_empty():
        return False
    return _occ(t[0], p[0]) and _occ(t[1], p[1])

def statistic(t):
    # Count the nodes of t at which the pattern sits contiguously.
    if t.is_empty():
        return 0
    return (1 if _occ(t, _PATTERN) else 0) + statistic(t[0]) + statistic(t[1])

# (* in Mathematica *)
# tree = {{{{}, {}}, {{}, {}}}, {{{}, {}}, {{}, {}}}};
# Count[tree, {{___}, {{{{{___}, {___}}, {___}}, {___}}, {___}}}, {0, Infinity}]

_PATTERN = BinaryTree('[.,[[.,.],[[.,.],.]]]')

def _occ(t, p):
    # A leaf '.' in the pattern is a wildcard; internal nodes must match
    # the left/right branching exactly and contiguously.
    if p.is_empty():
        return True
    if t.is_empty():
        return False
    return _occ(t[0], p[0]) and _occ(t[1], p[1])

def statistic(t):
    # Count the nodes of t at which the pattern sits contiguously.
    if t.is_empty():
        return 0
    return (1 if _occ(t, _PATTERN) else 0) + statistic(t[0]) + statistic(t[1])

# (* in Mathematica *)
# tree = {{{{}, {}}, {{}, {}}}, {{{}, {}}, {{}, {}}}};
# Count[tree, {{___}, {{{___}, {___}}, {{{___}, {___}}, {___}}}}, {0, Infinity}]

_PATTERN = BinaryTree('[.,[[.,.],[[.,.],.]]]')

def _occ(t, p):
    # A leaf '.' in the pattern is a wildcard; internal nodes must match
    # the left/right branching exactly and contiguously.
    if p.is_empty():
        return True
    if t.is_empty():
        return False
    return _occ(t[0], p[0]) and _occ(t[1], p[1])

def statistic(t):
    # Count the nodes of t at which the pattern sits contiguously.
    if t.is_empty():
        return 0
    return (1 if _occ(t, _PATTERN) else 0) + statistic(t[0]) + statistic(t[1])

# (* in Mathematica *)
# tree = {{{{}, {}}, {{}, {}}}, {{{}, {}}, {{}, {}}}};
# Count[tree, {{___}, {{{___}, {___}}, {{{___}, {___}}, {___}}}}, {0, Infinity}]

_PATTERN = BinaryTree('[.,[.,[[[.,.],.],.]]]')

def _occ(t, p):
    # A leaf '.' in the pattern is a wildcard; internal nodes must match
    # the left/right branching exactly and contiguously.
    if p.is_empty():
        return True
    if t.is_empty():
        return False
    return _occ(t[0], p[0]) and _occ(t[1], p[1])

def statistic(t):
    # Count the nodes of t at which the pattern sits contiguously.
    if t.is_empty():
        return 0
    return (1 if _occ(t, _PATTERN) else 0) + statistic(t[0]) + statistic(t[1])

# (* in Mathematica *)
# tree = {{{{}, {}}, {{}, {}}}, {{{}, {}}, {{}, {}}}};
# Count[tree, {{___}, {{___}, {{{{___}, {___}}, {___}}, {___}}}}, {0, Infinity}]

_PATTERN = BinaryTree('[.,[.,[[[.,.],.],.]]]')

def _occ(t, p):
    # A leaf '.' in the pattern is a wildcard; internal nodes must match
    # the left/right branching exactly and contiguously.
    if p.is_empty():
        return True
    if t.is_empty():
        return False
    return _occ(t[0], p[0]) and _occ(t[1], p[1])

def statistic(t):
    # Count the nodes of t at which the pattern sits contiguously.
    if t.is_empty():
        return 0
    return (1 if _occ(t, _PATTERN) else 0) + statistic(t[0]) + statistic(t[1])

# (* in Mathematica *)
# tree = {{{{}, {}}, {{}, {}}}, {{{}, {}}, {{}, {}}}};
# Count[tree, {{___}, {{___}, {{{{___}, {___}}, {___}}, {___}}}}, {0, Infinity}]

_PATTERN = BinaryTree('[.,[.,[[.,[.,.]],.]]]')

def _occ(t, p):
    # A leaf '.' in the pattern is a wildcard; internal nodes must match
    # the left/right branching exactly and contiguously.
    if p.is_empty():
        return True
    if t.is_empty():
        return False
    return _occ(t[0], p[0]) and _occ(t[1], p[1])

def statistic(t):
    # Count the nodes of t at which the pattern sits contiguously.
    if t.is_empty():
        return 0
    return (1 if _occ(t, _PATTERN) else 0) + statistic(t[0]) + statistic(t[1])

# --- original FindStat Code (Mathematica), commented out ---
# (* in Mathematica *)
# tree = {{{{}, {}}, {{}, {}}}, {{{}, {}}, {{}, {}}}};
# Count[tree, {{___}, {{___}, {{{___}, {{___}, {___}}}, {___}}}}, {0, Infinity}]

_PATTERN = BinaryTree('[.,[.,[[.,[.,.]],.]]]')

def _occ(t, p):
    # A leaf '.' in the pattern is a wildcard; internal nodes must match
    # the left/right branching exactly and contiguously.
    if p.is_empty():
        return True
    if t.is_empty():
        return False
    return _occ(t[0], p[0]) and _occ(t[1], p[1])

def statistic(t):
    # Count the nodes of t at which the pattern sits contiguously.
    if t.is_empty():
        return 0
    return (1 if _occ(t, _PATTERN) else 0) + statistic(t[0]) + statistic(t[1])

# --- original FindStat Code (Mathematica), commented out ---
# (* in Mathematica *)
# tree = {{{{}, {}}, {{}, {}}}, {{{}, {}}, {{}, {}}}};
# Count[tree, {{___}, {{___}, {{{___}, {{___}, {___}}}, {___}}}}, {0, Infinity}]

_PATTERN = BinaryTree('[.,[.,[.,[[.,.],.]]]]')

def _occ(t, p):
    # A leaf '.' in the pattern is a wildcard; internal nodes must match
    # the left/right branching exactly and contiguously.
    if p.is_empty():
        return True
    if t.is_empty():
        return False
    return _occ(t[0], p[0]) and _occ(t[1], p[1])

def statistic(t):
    # Count the nodes of t at which the pattern sits contiguously.
    if t.is_empty():
        return 0
    return (1 if _occ(t, _PATTERN) else 0) + statistic(t[0]) + statistic(t[1])

# --- original FindStat Code (Mathematica), commented out ---
# (* in Mathematica *)
# tree = {{{{}, {}}, {{}, {}}}, {{{}, {}}, {{}, {}}}};
# Count[tree, {{___}, {{___}, {{___}, {{{___}, {___}}, {___}}}}}, {0, Infinity}]

_PATTERN = BinaryTree('[.,[.,[.,[[.,.],.]]]]')

def _occ(t, p):
    # A leaf '.' in the pattern is a wildcard; internal nodes must match
    # the left/right branching exactly and contiguously.
    if p.is_empty():
        return True
    if t.is_empty():
        return False
    return _occ(t[0], p[0]) and _occ(t[1], p[1])

def statistic(t):
    # Count the nodes of t at which the pattern sits contiguously.
    if t.is_empty():
        return 0
    return (1 if _occ(t, _PATTERN) else 0) + statistic(t[0]) + statistic(t[1])

# --- original FindStat Code (Mathematica), commented out ---
# (* in Mathematica *)
# tree = {{{{}, {}}, {{}, {}}}, {{{}, {}}, {{}, {}}}};
# Count[tree, {{___}, {{___}, {{___}, {{{___}, {___}}, {___}}}}}, {0, Infinity}]

_PATTERN = BinaryTree('[.,[.,[.,[.,[.,.]]]]]')

def _occ(t, p):
    # A leaf '.' in the pattern is a wildcard; internal nodes must match
    # the left/right branching exactly and contiguously.
    if p.is_empty():
        return True
    if t.is_empty():
        return False
    return _occ(t[0], p[0]) and _occ(t[1], p[1])

def statistic(t):
    # Count the nodes of t at which the pattern sits contiguously.
    if t.is_empty():
        return 0
    return (1 if _occ(t, _PATTERN) else 0) + statistic(t[0]) + statistic(t[1])

# (* in Mathematica *)
# tree = {{{{}, {}}, {{}, {}}}, {{{}, {}}, {{}, {}}}};
# Count[tree, {{___}, {{___}, {{___}, {{___}, {{___}, {___}}}}}}, {0, Infinity}]

_PATTERN = BinaryTree('[.,[.,[.,[.,[.,.]]]]]')

def _occ(t, p):
    # A leaf '.' in the pattern is a wildcard; internal nodes must match
    # the left/right branching exactly and contiguously.
    if p.is_empty():
        return True
    if t.is_empty():
        return False
    return _occ(t[0], p[0]) and _occ(t[1], p[1])

def statistic(t):
    # Count the nodes of t at which the pattern sits contiguously.
    if t.is_empty():
        return 0
    return (1 if _occ(t, _PATTERN) else 0) + statistic(t[0]) + statistic(t[1])

# (* in Mathematica *)
# tree = {{{{}, {}}, {{}, {}}}, {{{}, {}}, {{}, {}}}};
# Count[tree, {{___}, {{___}, {{___}, {{___}, {{___}, {___}}}}}}, {0, Infinity}]

_PATTERN = BinaryTree('[.,[[[.,.],.],.]]')

def _occ(t, p):
    # A leaf '.' in the pattern is a wildcard; internal nodes must match
    # the left/right branching exactly and contiguously.
    if p.is_empty():
        return True
    if t.is_empty():
        return False
    return _occ(t[0], p[0]) and _occ(t[1], p[1])

def statistic(t):
    # Count the nodes of t at which the pattern sits contiguously.
    if t.is_empty():
        return 0
    return (1 if _occ(t, _PATTERN) else 0) + statistic(t[0]) + statistic(t[1])

# (* in Mathematica *)
# tree = {{{{}, {}}, {{}, {}}}, {{{}, {}}, {{}, {}}}};
# Count[tree, {{___}, {{{{___}, {___}}, {___}}, {___}}}, {0, Infinity}]

_PATTERN = BinaryTree('[.,[[[.,.],.],.]]')

def _occ(t, p):
    # A leaf '.' in the pattern is a wildcard; internal nodes must match
    # the left/right branching exactly and contiguously.
    if p.is_empty():
        return True
    if t.is_empty():
        return False
    return _occ(t[0], p[0]) and _occ(t[1], p[1])

def statistic(t):
    # Count the nodes of t at which the pattern sits contiguously.
    if t.is_empty():
        return 0
    return (1 if _occ(t, _PATTERN) else 0) + statistic(t[0]) + statistic(t[1])

# (* in Mathematica *)
# tree = {{{{}, {}}, {{}, {}}}, {{{}, {}}, {{}, {}}}};
# Count[tree, {{___}, {{{{___}, {___}}, {___}}, {___}}}, {0, Infinity}]

_PATTERN = BinaryTree('[.,[.,[[.,.],.]]]')

def _occ(t, p):
    # A leaf '.' in the pattern is a wildcard; internal nodes must match
    # the left/right branching exactly and contiguously.
    if p.is_empty():
        return True
    if t.is_empty():
        return False
    return _occ(t[0], p[0]) and _occ(t[1], p[1])

def statistic(t):
    # Count the nodes of t at which the pattern sits contiguously.
    if t.is_empty():
        return 0
    return (1 if _occ(t, _PATTERN) else 0) + statistic(t[0]) + statistic(t[1])

# (* in Mathematica *)
# tree = {{{{}, {}}, {{}, {}}}, {{{}, {}}, {{}, {}}}};
# Count[tree, {{___}, {{___}, {{{___}, {___}}, {___}}}}, {0, Infinity}]

_PATTERN = BinaryTree('[.,[.,[[.,.],.]]]')

def _occ(t, p):
    # A leaf '.' in the pattern is a wildcard; internal nodes must match
    # the left/right branching exactly and contiguously.
    if p.is_empty():
        return True
    if t.is_empty():
        return False
    return _occ(t[0], p[0]) and _occ(t[1], p[1])

def statistic(t):
    # Count the nodes of t at which the pattern sits contiguously.
    if t.is_empty():
        return 0
    return (1 if _occ(t, _PATTERN) else 0) + statistic(t[0]) + statistic(t[1])

# (* in Mathematica *)
# tree = {{{{}, {}}, {{}, {}}}, {{{}, {}}, {{}, {}}}};
# Count[tree, {{___}, {{___}, {{{___}, {___}}, {___}}}}, {0, Infinity}]

_PATTERN = BinaryTree('[.,[.,[.,[.,.]]]]')

def _occ(t, p):
    # A leaf '.' in the pattern is a wildcard; internal nodes must match
    # the left/right branching exactly and contiguously.
    if p.is_empty():
        return True
    if t.is_empty():
        return False
    return _occ(t[0], p[0]) and _occ(t[1], p[1])

def statistic(t):
    # Count the nodes of t at which the pattern sits contiguously.
    if t.is_empty():
        return 0
    return (1 if _occ(t, _PATTERN) else 0) + statistic(t[0]) + statistic(t[1])

# (* in Mathematica *)
# tree = {{{{}, {}}, {{}, {}}}, {{{}, {}}, {{}, {}}}};
# Count[tree, {{___}, {{___}, {{___}, {{___}, {___}}}}}, {0, Infinity}]

_PATTERN = BinaryTree('[.,[.,[.,[.,.]]]]')

def _occ(t, p):
    # A leaf '.' in the pattern is a wildcard; internal nodes must match
    # the left/right branching exactly and contiguously.
    if p.is_empty():
        return True
    if t.is_empty():
        return False
    return _occ(t[0], p[0]) and _occ(t[1], p[1])

def statistic(t):
    # Count the nodes of t at which the pattern sits contiguously.
    if t.is_empty():
        return 0
    return (1 if _occ(t, _PATTERN) else 0) + statistic(t[0]) + statistic(t[1])

# (* in Mathematica *)
# tree = {{{{}, {}}, {{}, {}}}, {{{}, {}}, {{}, {}}}};
# Count[tree, {{___}, {{___}, {{___}, {{___}, {___}}}}}, {0, Infinity}]

_PATTERN = BinaryTree('[.,[.,[.,.]]]')

def _occ(t, p):
    # A leaf '.' in the pattern is a wildcard; internal nodes must match
    # the left/right branching exactly and contiguously.
    if p.is_empty():
        return True
    if t.is_empty():
        return False
    return _occ(t[0], p[0]) and _occ(t[1], p[1])

def statistic(t):
    # Count the nodes of t at which the pattern sits contiguously.
    if t.is_empty():
        return 0
    return (1 if _occ(t, _PATTERN) else 0) + statistic(t[0]) + statistic(t[1])

# (* in Mathematica *)
# tree = {{{{}, {}}, {{}, {}}}, {{{}, {}}, {{}, {}}}};
# Count[tree, {{___}, {{___}, {{___}, {___}}}}, {0, Infinity}]

_PATTERN = BinaryTree('[.,[.,[.,.]]]')

def _occ(t, p):
    # A leaf '.' in the pattern is a wildcard; internal nodes must match
    # the left/right branching exactly and contiguously.
    if p.is_empty():
        return True
    if t.is_empty():
        return False
    return _occ(t[0], p[0]) and _occ(t[1], p[1])

def statistic(t):
    # Count the nodes of t at which the pattern sits contiguously.
    if t.is_empty():
        return 0
    return (1 if _occ(t, _PATTERN) else 0) + statistic(t[0]) + statistic(t[1])

# (* in Mathematica *)
# tree = {{{{}, {}}, {{}, {}}}, {{{}, {}}, {{}, {}}}};
# Count[tree, {{___}, {{___}, {{___}, {___}}}}, {0, Infinity}]

def statistic(C):
    W = CoxeterGroup(C, implementation='coxeter3')
    r = []
    for u in W:
        U = (W(v) for v in W.bruhat_interval(u, W.long_element()))
        next(U)
        for v in U:
            ldiff = v.length()-u.length()-1
            if is_even(ldiff):
                p = W.kazhdan_lusztig_polynomial(u, v)
                r.append(p[ldiff//2])
    return max(r)

def statistic(C):
    W = CoxeterGroup(C, implementation='coxeter3')
    r = []
    for u in W:
	        U = (W(v) for v in W.bruhat_interval(u, W.long_element()))
        next(U)
        for v in U:
            ldiff = v.length()-u.length()-1
            if is_even(ldiff):
                p = W.kazhdan_lusztig_polynomial(u, v)
                r.append(p[ldiff//2])
    return max(r)

def statistic(C):
    """                                                                                                                                               
    sage: statistic(CartanType(["A", 4]))                                                                                                             
    482                                                                                                                                               
    """
    W = CoxeterGroup(C, implementation='coxeter3')
    r = 0
    for u in W:
        U = (W(v) for v in W.bruhat_interval(u, W.long_element()))
        next(U)
        for v in U:
            ldiff = v.length()-u.length()-1
            if is_even(ldiff):
                p = W.kazhdan_lusztig_polynomial(u, v)
                if p[ldiff//2] != 0:
                    r += 1
    return r

def statistic(C):
    """                                                                                                                                               
    sage: statistic(CartanType(["A", 4]))                                                                                                             
    482                                                                                                                                               
    """
    W = CoxeterGroup(C, implementation='coxeter3')
    r = 0
    for u in W:
	        U = (W(v) for v in W.bruhat_interval(u, W.long_element()))
	        next(U)
	        for v in U:
            ldiff = v.length()-u.length()-1
            if is_even(ldiff):
		                p = W.kazhdan_lusztig_polynomial(u, v)
		                if p[ldiff//2] != 0:
                    r += 1
    return r

# code for St000482 taken from https://findstat.org/StatisticsDatabase/St00482/

import random
def Z_game(g,B,ban=[]):
    """
    Input:
        g: a simple graph
        B: a set of initial blue vertices
        ban: a set of banned vertices
    Output:
        return the derived set under the regular CCR-Z.
        Note: vertices in ban cannot make a force, but are still white neighbors if white.
    """    
    V=g.vertices();
    white_neighbors={}; #a dictionary with the structure {v: list of white neighbors}
    white_numbers={}; #a dictionary with the structure {v: number of white neighbors}
    for v in V:
        nbh=g.neighbors(v);
        for b in B:
            try:
                nbh.remove(b);
            except ValueError:
                pass;
        white_neighbors[v]=nbh;
        white_numbers[v]=len(nbh);
    queue=copy(B); #queue stores list of vertices that can possibly make a force
    derived_set=copy(B); #derived_set stores the set of blue vertices
    whole_loop=True;
    while whole_loop: #keep searching if queue!=[]
        try:
            v=queue[0];
            queue.remove(v);
            if v not in ban and white_numbers[v]==1:
                u=white_neighbors[v][0]; #the only white neighbor
                derived_set.append(u); #make the force
                #update white_numbers, white_neighbors, and queue
                if white_numbers[u]==1: 
                    queue.append(u);
                u_nbr=g.neighbors(u);
                for w in u_nbr:
                    white_neighbors[w].remove(u);
                    white_numbers[w]+=-1;
                    if w in derived_set and white_numbers[w]==1:
                        queue.append(w);
        except IndexError:
            whole_loop=False;
    return derived_set;

def St000482(g):
    """
    Input:
        g: a simple graph
    Output:
        return the value of the conventional zero forcing number Z(g).
    """
    V=g.vertices();
    n=g.order();
    lbd=-1;
    ubd=n;
    while ubd-lbd>=2: #apply random algorithm;
        guess=random.choice(range(lbd+1,ubd)); #take an interior point in [lbd,ubd];
        found=False;
        for sub in Combinations(V,guess):
            if n==len(Z_game(g,sub)):
                   found=True;
                   break;
        if found:
            ubd=guess;
        else:
            lbd=guess;
    return ubd;

@cached_function
def statistic(G):
    if G.size() == 0:
        return 0  # other definitions possible

    lower_bound = G.pathwidth()
    upper_bound = St000482(G)  # zero forcing number
    if upper_bound == lower_bound:
        return upper_bound

    for e in G.complement().edges(labels=False):
        H = G.copy(immutable=False)
        H.add_edge(e)
        upper_bound = min(upper_bound, statistic(H.copy(immutable=True)))
        if upper_bound == lower_bound:
            return upper_bound

    return upper_bound

# code for St000482 taken from https://findstat.org/StatisticsDatabase/St00482/

import random
def Z_game(g,B,ban=[]):
    """
    Input:
        g: a simple graph
        B: a set of initial blue vertices
        ban: a set of banned vertices
    Output:
        return the derived set under the regular CCR-Z.
        Note: vertices in ban cannot make a force, but are still white neighbors if white.
    """    
    V=g.vertices();
    white_neighbors={}; #a dictionary with the structure {v: list of white neighbors}
    white_numbers={}; #a dictionary with the structure {v: number of white neighbors}
    for v in V:
        nbh=g.neighbors(v);
        for b in B:
            try:
                nbh.remove(b);
            except ValueError:
                pass;
        white_neighbors[v]=nbh;
        white_numbers[v]=len(nbh);
    queue=copy(B); #queue stores list of vertices that can possibly make a force
    derived_set=copy(B); #derived_set stores the set of blue vertices
    whole_loop=True;
    while whole_loop: #keep searching if queue!=[]
        try:
            v=queue[0];
            queue.remove(v);
            if v not in ban and white_numbers[v]==1:
                u=white_neighbors[v][0]; #the only white neighbor
                derived_set.append(u); #make the force
                #update white_numbers, white_neighbors, and queue
                if white_numbers[u]==1: 
                    queue.append(u);
                u_nbr=g.neighbors(u);
                for w in u_nbr:
                    white_neighbors[w].remove(u);
                    white_numbers[w]+=-1;
                    if w in derived_set and white_numbers[w]==1:
                        queue.append(w);
        except IndexError:
            whole_loop=False;
    return derived_set;

def St000482(g):
    """
    Input:
        g: a simple graph
    Output:
        return the value of the conventional zero forcing number Z(g).
    """
    V=g.vertices();
    n=g.order();
    lbd=-1;
    ubd=n;
    while ubd-lbd>=2: #apply random algorithm;
        guess=random.choice(range(lbd+1,ubd)); #take an interior point in [lbd,ubd];
        found=False;
        for sub in Combinations(V,guess):
            if n==len(Z_game(g,sub)):
                   found=True;
                   break;
        if found:
            ubd=guess;
        else:
            lbd=guess;
    return ubd;

@cached_function
def statistic(G):
    if G.size() == 0:
        return 0  # other definitions possible

    lower_bound = G.pathwidth()
    upper_bound = St000482(G)  # zero forcing number
    if upper_bound == lower_bound:
        return upper_bound

    for e in G.complement().edges(labels=False):
        H = G.copy(immutable=False)
        H.add_edge(e)
        upper_bound = min(upper_bound, statistic(H.copy(immutable=True)))
        if upper_bound == lower_bound:
            return upper_bound

    return upper_bound

from sage.combinat.cyclic_sieving_phenomenon import orbit_decomposition
def toggle_labelling(G, pi, i, j):
    if G.has_edge(pi.index(i), pi.index(j)):
        return pi
    sigma = [j if e == i else i if e == j else e for e in pi]
    return Permutation(sigma)

def toric_promotion_labelling(G, pi):
    n = G.num_verts()
    assert set(G.vertices()) == set(range(n))
    for i in range(1, n):
        pi = toggle_labelling(G, pi, i, i+1)
    return toggle_labelling(G, pi, n, 1)

def toric_promotion_labelling_orbits(G):
    G = G.canonical_label().copy(immutable=True)
    return toric_promotion_labelling_orbits_aux(G)

def toric_promotion_labelling_orbits_aux(G):
    n = G.num_verts()
    return orbit_decomposition(Permutations(n),
                               lambda pi: toric_promotion_labelling(G, pi))

def statistic(G):
    return lcm(len(o) for o in toric_promotion_labelling_orbits(G))

from sage.combinat.cyclic_sieving_phenomenon import orbit_decomposition
def toggle_labelling(G, pi, i, j):
    if G.has_edge(pi.index(i), pi.index(j)):
        return pi
    sigma = [j if e == i else i if e == j else e for e in pi]
    return Permutation(sigma)

def toric_promotion_labelling(G, pi):
    n = G.num_verts()
    assert set(G.vertices()) == set(range(n))
    for i in range(1, n):
        pi = toggle_labelling(G, pi, i, i+1)
    return toggle_labelling(G, pi, n, 1)

def toric_promotion_labelling_orbits(G):
    G = G.canonical_label().copy(immutable=True)
    return toric_promotion_labelling_orbits_aux(G)

@cached_function
def toric_promotion_labelling_orbits_aux(G):
    n = G.num_verts()
    return orbit_decomposition(Permutations(n),
                               lambda pi: toric_promotion_labelling(G, pi))

def statistic(G):
    return lcm(len(o) for o in toric_promotion_labelling_orbits(G))

from sage.combinat.cyclic_sieving_phenomenon import orbit_decomposition
def toggle_labelling(G, pi, i, j):
    if G.has_edge(pi.index(i), pi.index(j)):
        return pi
    sigma = [j if e == i else i if e == j else e for e in pi]
    return Permutation(sigma)

def promotion_labelling(G, pi):
    n = G.num_verts()
    assert set(G.vertices()) == set(range(n))
    for i in range(1, n):
        pi = toggle_labelling(G, pi, i, i+1)
    return pi

def promotion_labelling_orbits(G):
    G = G.canonical_label().copy(immutable=True)
    return promotion_labelling_orbits_aux(G)

def promotion_labelling_orbits_aux(G):
    n = G.num_verts()
    return orbit_decomposition(Permutations(n),
                               lambda pi: promotion_labelling(G, pi))

def statistic(G):
    return len(promotion_labelling_orbits(G))

from sage.combinat.cyclic_sieving_phenomenon import orbit_decomposition
def toggle_labelling(G, pi, i, j):
    if G.has_edge(pi.index(i), pi.index(j)):
        return pi
    sigma = [j if e == i else i if e == j else e for e in pi]
    return Permutation(sigma)

def promotion_labelling(G, pi):
    n = G.num_verts()
    assert set(G.vertices()) == set(range(n))
    for i in range(1, n):
        pi = toggle_labelling(G, pi, i, i+1)
    return pi

def promotion_labelling_orbits(G):
    G = G.canonical_label().copy(immutable=True)
    return promotion_labelling_orbits_aux(G)

@cached_function
def promotion_labelling_orbits_aux(G):
    n = G.num_verts()
    return orbit_decomposition(Permutations(n),
                               lambda pi: promotion_labelling(G, pi))

def statistic(G):
    return len(promotion_labelling_orbits(G))

from sage.combinat.cyclic_sieving_phenomenon import orbit_decomposition
def toggle_labelling(G, pi, i, j):
    if G.has_edge(pi.index(i), pi.index(j)):
        return pi
    sigma = [j if e == i else i if e == j else e for e in pi]
    return Permutation(sigma)

def toric_promotion_labelling(G, pi):
    n = G.num_verts()
    assert set(G.vertices()) == set(range(n))
    for i in range(1, n):
        pi = toggle_labelling(G, pi, i, i+1)
    return toggle_labelling(G, pi, n, 1)

def toric_promotion_labelling_orbits(G):
    G = G.canonical_label().copy(immutable=True)
    return toric_promotion_labelling_orbits_aux(G)

def toric_promotion_labelling_orbits_aux(G):
    n = G.num_verts()
    return orbit_decomposition(Permutations(n),
                               lambda pi: toric_promotion_labelling(G, pi))

def statistic(G):
    return len(toric_promotion_labelling_orbits(G))

from sage.combinat.cyclic_sieving_phenomenon import orbit_decomposition
def toggle_labelling(G, pi, i, j):
    if G.has_edge(pi.index(i), pi.index(j)):
        return pi
    sigma = [j if e == i else i if e == j else e for e in pi]
    return Permutation(sigma)

def toric_promotion_labelling(G, pi):
    n = G.num_verts()
    assert set(G.vertices()) == set(range(n))
    for i in range(1, n):
        pi = toggle_labelling(G, pi, i, i+1)
    return toggle_labelling(G, pi, n, 1)

def toric_promotion_labelling_orbits(G):
    G = G.canonical_label().copy(immutable=True)
    return toric_promotion_labelling_orbits_aux(G)

@cached_function
def toric_promotion_labelling_orbits_aux(G):
    n = G.num_verts()
    return orbit_decomposition(Permutations(n),
                               lambda pi: toric_promotion_labelling(G, pi))

def statistic(G):
    return len(toric_promotion_labelling_orbits(G))

@cached_function
def statistic(G):
    if G.is_planar():
        return 0
    bound = G.size()
    for e in G.edges(labels=False):
        H = G.copy(immutable=False)
        H.delete_edge(e)
        bound = min(bound, 1+statistic(H.canonical_label().copy(immutable=True)))
    return bound

# alternative slower code
# def statistic(G):
#     E = G.edges(labels=False)
#     m = len(E)
#     for sublist in reversed(Subsets(E).list()):
#         if Graph(list(sublist)).is_planar():
#             return m-len(sublist)

@cached_function
def statistic(G):
    if G.is_planar():
        return 0
    bound = G.size()
    for e in G.edges(labels=False):
        H = G.copy(immutable=False)
        H.delete_edge(e)
        bound = min(bound, 1+statistic(H.canonical_label().copy(immutable=True)))
    return bound

# alternative slower code
# def statistic(G):
#     E = G.edges(labels=False)
#     m = len(E)
#     for sublist in reversed(Subsets(E).list()):
#         if Graph(list(sublist)).is_planar():
#             return m-len(sublist)

# Pure combinatorial: count indecomposable reflexive (2-syzygy) modules
# with dominant dimension >= 1.
# For Nakayama algebras, domdim >= 1 means the module is torsionless (a 1st syzygy).
# Every 2-syzygy is automatically a 1st syzygy, so the domdim filter is redundant.
# An indecomposable M(i, l) is a 2-syzygy iff:
#   - M(i, l) is projective (l = K[i]), OR
#   - ∃ M(j, m) with Ω²(M(j, m)) = M(i, l)

def kupisch(D):
    DR = D.reverse()
    H = DR.heights()
    return [1 + H[i] for i, s in enumerate(DR) if s == 0] + [1]

def statistic(D):
    K = kupisch(D)
    n = len(K)

    def omega(i, l):
        if l == K[i]:
            return None
        ni, nl = i + l, K[i] - l
        if ni >= n or nl <= 0:
            return None
        return (ni, nl)

    syzygy2_modules = set()

    # All projective modules are reflexive (2-syzygies for all n)
    for i in range(n):
        syzygy2_modules.add((i, K[i]))

    # Find non-projective indecomposables that appear as Ω² of something
    for j in range(n):
        for l in range(1, K[j] + 1):
            cur = (j, l)
            for _ in range(2):
                if cur is None:
                    break
                cur = omega(cur[0], cur[1])
            if cur is not None:
                syzygy2_modules.add(cur)

    return ZZ(len(syzygy2_modules))

DeclareOperation("IsNthSyzygy",[IsList]);

InstallMethod(IsNthSyzygy, "for a representation of a quiver", [IsList],0,function(LIST)

local M, n, f, N, i, h,W;

M:=LIST[1];

n:=LIST[2];

N:=DualOfModule(NthSyzygy(DualOfModule(M),n));

W:=NthSyzygy(N,n);

if IsDirectSummand(M,W)=true
then
return(1);
else return(0);
fi;
end);


DeclareOperation("NumberOfReflexiveModules",[IsList# Pure combinatorial: count indecomposable reflexive (2-syzygy) modules
# with dominant dimension >= 1.
# For Nakayama algebras, domdim >= 1 means the module is torsionless (a 1st syzygy).
# Every 2-syzygy is automatically a 1st syzygy, so the domdim filter is redundant.
# An indecomposable M(i, l) is a 2-syzygy iff:
#   - M(i, l) is projective (l = K[i]), OR
#   - ∃ M(j, m) with Ω²(M(j, m)) = M(i, l)

def kupisch(D):
    DR = D.reverse()
    H = DR.heights()
    return [1 + H[i] for i, s in enumerate(DR) if s == 0] + [1]

def statistic(D):
    K = kupisch(D)
    n = len(K)

    def omega(i, l):
        if l == K[i]:
            return None
        ni, nl = i + l, K[i] - l
        if ni >= n or nl <= 0:
            return None
        return (ni, nl)

    syzygy2_modules = set()

    # All projective modules are reflexive (2-syzygies for all n)
    for i in range(n):
        syzygy2_modules.add((i, K[i]);)

InstallMethod(NumberOfReflexiveModules, "for a representation of a quiver", [IsList],0,function(LIST)

local M, n, f, N, i, h,W;

L:=LIST[1];

A:=NakayamaAlgebra(L,GF(3));
L:=ARQuiver([A,1000])[2];
LL1:=Filtered(L,x->DominantDimensionOfModule(x,30)>=1);
LL2:=Filtered(LL1,x->IsNthSyzygy([x,2])=1);
return(Size(LL2));


end);    # Find non-projective indecomposables that appear as Ω² of something
    for j in range(n):
        for l in range(1, K[j] + 1):
            cur = (j, l)
            for _ in range(2):
                if cur is None:
                    break
                cur = omega(cur[0], cur[1])
            if cur is not None:
                syzygy2_modules.add(cur)

    return ZZ(len(syzygy2_modules))

def standard_monomial(D, low=None, east=1):                      
    n = len(D)                         
    A = [ i+1 for i,s in enumerate(D) if s == 1-east ]
    if low is None:
        low = [1 .. len(A)]
    low = [ i for i in [1 .. n] if i not in low ]
    vals = []
    j = 0
    blocker = low + [Infinity]
    for i in [1 .. n]:
        if i not in A:
            if i < blocker[0]:
                j += 1
            else:
                blocker.pop(0)
        else:
            if j > 0:
                j -= 1
                vals.append(i)
                blocker.pop(0)
    return tuple(vals)

def statistic(D):
    return len(standard_monomial(D, east=0))

def standard_monomial(D, low=None, east=1):                      
    n = len(D)                         
    A = [ i+1 for i,s in enumerate(D) if s == 1-east ]
    low = [ i+1 for i,s in enumerate(low) if s == 1-east ]
    if low is None:
        low = [1 .. len(A)]
    low = [ i for i in [1 .. n] if i not in low ]
    vals = []
    j = 0
    blocker = low + [Infinity]
    for i in [1 .. n]:
        if i not in A:
            if i < blocker[0]:
                j += 1
            else:
                blocker.pop(0)
        else:
            if j > 0:
                j -= 1
                vals.append(i)
                blocker.pop(0)
    return tuple(vals)

def statistic(D):
    return len(standard_monomial(D, east=0))

def expand(M, ac, bd):
    a, c = ac
    b, d = bd
    m = list(M)
    m.remove(frozenset(ac))
    m.remove(frozenset(bd))
    m1 = m + [frozenset([a, b]), frozenset([c, d])]
    m2 = m + [frozenset([a, d]), frozenset([b, c])]
    return PerfectMatching(m1), PerfectMatching(m2)

def expansion(M):
    input = {M: 1}
    output = {}
    while input:
        m, e = input.popitem()
        try:
            ac, bd = next(m.crossings_iterator())
        except StopIteration:
            output[m] = output.get(m, 0) + e
        else:
            m1, m2 = expand(m, ac, bd)
            input[m1] = input.get(m1, 0) + e
            input[m2] = input.get(m2, 0) + e
    return output

def statistic(M):
    NN = PerfectMatching([[2*i+1,2*(i+1)] for i in range(len(M))])
    return expansion(M).get(NN, 0)

def expand(M, ac, bd):
    a, c = ac
    b, d = bd
    m = list(M)
    m.remove(frozenset(ac))
    m.remove(frozenset(bd))
    m1 = m + [frozenset([a, b]), frozenset([c, d])]
    m2 = m + [frozenset([a, d]), frozenset([b, c])]
    return PerfectMatching(m1), PerfectMatching(m2)

def expansion(M):
    input = {M: 1}
    output = {}
    while input:
        m, e = input.popitem()
        try:
            ac, bd = next(m.crossings_iterator())
        except StopIteration:
            output[m] = output.get(m, 0) + e
        else:
            m1, m2 = expand(m, ac, bd)
            input[m1] = input.get(m1, 0) + e
            input[m2] = input.get(m2, 0) + e
    return output

def statistic(M):
    NN = PerfectMatching([[2*i+1,2*(i+1)] for i in range(len(M))])
    return expansion(M).get(non_nesting_non_crossing(len(M))NN, 0)

@cached_function
def statistic(w):
    def children(m):
        for (a,b),(c,d) in m.crossings():
            m_new = list(m)
            m_new.remove(frozenset((a,b)))
            m_new.remove(frozenset((c,d)))
            A, B = min(a,b), max(a,b)
            C, D = min(c,d), max(c,d)
            if C < A:
                (A,B),(C,D) = (C,D),(A,B)
            yield PerfectMatching(m_new + [frozenset((A,C)), frozenset((B,D))])
            yield PerfectMatching(m_new + [frozenset((A,D)), frozenset((B,C))])

    w = PerfectMatching(sorted([sorted(a) for a in w]))
    l = [statistic(v) for v in children(w)]
    if len(l) == 0:
        return 0
    else:
        return 1+min(l)

@cached_function
def statistic(w):
    def children(m):
        for (a,b),(c,d) in m.crossings():
            m_new = list(m)
            m_new.remove(frozenset((a,b)))
            m_new.remove(frozenset((c,d)))
            A, B = min(a,b), max(a,b)
            C, D = min(c,d), max(c,d)
            if C < A:
                (A,B),(C,D) = (C,D),(A,B)
            yield PerfectMatching(m_new + [frozenset((A,C)), frozenset((B,D))])
            yield PerfectMatching(m_new + [frozenset((A,D)), frozenset((B,C))])

    w = PerfectMatching(sorted([sorted(a) for a in w]))
    l = [statistic(v) for v in children(w)]
    if len(l) == 0:
        return 0
    else:
        return 1+min(l)

@cached_function
def statistic(w):
    def children(m):
        for (a,b),(c,d) in m.nestings():
            m_new = list(m)
            m_new.remove(frozenset((a,b)))
            m_new.remove(frozenset((c,d)))
            A = min(a,b)
            B = max(a,b)
            C = min(c,d)
            D = max(c,d)
            if C < A:
                (A,B),(C,D) = (C,D),(A,B)
            yield PerfectMatching(m_new + [frozenset((A,C)), frozenset((B,D))])

    l = [statistic(v) for v in children(w)]
    i = 0
    while i in l:
        i += 1
    return i

@cached_function
def statistic(w):
    def children(m):
        for (a,b),(c,d) in m.nestings():
            m_new = list(m)
            m_new.remove(frozenset((a,b)))
            m_new.remove(frozenset((c,d)))
            A = min(a,b)
            B = max(a,b)
            C = min(c,d)
            D = max(c,d)
            if C < A:
                (A,B),(C,D) = (C,D),(A,B)
            yield PerfectMatching(m_new + [frozenset((A,C)), frozenset((B,D))])

    l = [statistic(v) for v in children(w)]
    i = 0
    while i in l:
        i += 1
    return i

def statistic(T):
    S = [[e for e in row] for row in T]
    count = 0
    for i, row in enumerate(S):
        for j, e in enumerate(row):
            S[i][j] = e+1
            try:
                _ = SemistandardTableau(S)
            except Exception:
                count += 1
            S[i][j] = e
    return count

def statistic(T):
    S = [[e for e in row] for row in T]
    count = 0
    for i, row in enumerate(S):
        for j, e in enumerate(row):
            S[i][j] = e+1
            try:
                _ = SemistandardTableau(S)
            except StandardErrorException:
                count += 1
            S[i][j] = e
    return count

def statistic(T):
    S = [[e for e in row] for row in T]
    count = 0
    for i, row in enumerate(S):
        for j, e in enumerate(row):
            S[i][j] = e-1
            try:
                _ = SemistandardTableau(S)
            except Exception:
                count += 1
            S[i][j] = e
    return count

def statistic(T):
    S = [[e for e in row] for row in T]
    count = 0
    for i, row in enumerate(S):
        for j, e in enumerate(row):
            S[i][j] = e-1
            try:
                _ = SemistandardTableau(S)
            except StandardErrorException:
                count += 1
            S[i][j] = e
    return count

def statistic(T):
    f = 0
    segments = {}
    for r in range(len(T)):
        for c in range(len(T[r])):
            for j in range(c+1):
                if T[r][j] != r+1:
                    if (r,T[r][j]) not in segments.keys():
                        segments[(r,T[r][j])] = j
    L = segments.items()
    for s in segments.items():
        if s[0][0] != len(T)-1 and s[1] == len(T[s[0][0]+1]) and T[s[0][0]+1][-1] <= s[0][1]:
            f += 1
        if s[0][0] == len(T)-1 and s[1] == 0:
            f += 1
        else:
            for t in L:
                if s[0][0]+1 == t[0][0] and s[1] == t[1]:
                    if s[1] >= 1 and T[s[0][0]+1][s[1]] != T[s[0][0]+1][s[1]-1]:
                        f += 1
                    if s[1] < 1 and T[s[0][0]+1][s[1]] != s[0][0]+2:
                        f += 1
    return f

def statistic(T):
    f = 0
    segments = {}
    for r in range(len(T)):
        for c in range(len(T[r])):
            for j in range(c+1):
                if T[r][j] != r+1:
                    if (r,T[r][j]) not in segments.keys():
                        segments[(r,T[r][j])] = j
    L = segments.items()
    for s in segments.iteritems():
        if s[0][0] != len(T)-1 and s[1] == len(T[s[0][0]+1]) and T[s[0][0]+1][-1] <= s[0][1]:
            f += 1
        if s[0][0] == len(T)-1 and s[1] == 0:
            f += 1
        else:
            for t in L:
                if s[0][0]+1 == t[0][0] and s[1] == t[1]:
                    if s[1] >= 1 and T[s[0][0]+1][s[1]] != T[s[0][0]+1][s[1]-1]:
                        f += 1
                    if s[1] < 1 and T[s[0][0]+1][s[1]] != s[0][0]+2:
                        f += 1
    return f

def statistic(P):
    return len(poset_endomorphisms(P))

def poset_endomorphisms(P):
    P = P.relabel()
    r = P.cardinality()
    S = cartesian_product([range(r)]*r)
    return [pi for pi in S if P.is_poset_morphism(lambda i: pi[i], P)]

# (* Mathematica code using depth first recursion *)

# endoCount::usage := "If rel is a reflexive relation on 1,2,..,Max[rel], then \
# endoCount[rel] is the number of self maps f for which {f[a],f[b]} is in the \
# relation whenever {a,b} is in the relation.”

# endoCount[rel_] := morphismCount[rel, rel]

# morphismCount[rel1_List, rel2_List] :=
# Module[{max1 = Max[rel1], max2 = Max[rel2], down, checkdown, num, ans},
#  Do[down[i] = Select[rel1, Max[#] == i &], {i, max1}];
#  checkdown[f_List] := AllTrue[down[Length[f]], MemberQ[rel2, f[[#1]]] &];
#  num[{}] := Sum[num[{i}], {i, max2}];
#  num[f_List] :=
#   If[checkdown[f],
#    If[Length[f] == max1, 1, Sum[num[Append[f, i]], {i, max2}]], 0];
#  ans = num[{}]; Clear[checkdown, down]; ans]

# (* example *)

# endoCount[{{1, 1}, {1, 2}, {1, 10}, {2, 2}, {3, 2}, {3, 3}, {3,
#   4}, {4, 4}, {5, 4}, {5, 5}, {5, 6}, {6, 6}, {7, 6}, {7, 7}, {7,
#   8}, {8, 8}, {9, 8}, {9, 9}, {9, 10}, {10, 10}}]

# (* gives 10030 *)

# (* Sanity checks are to be added. For instance we must have 
#   Union @@ rel == Range[Max[rel]]
# *)

def statistic(P):
    return len(poset_endomorphisms(P))

def poset_endomorphisms(P):
    P = P.relabel()
    r = P.cardinality()
    S = cartesian_product([range(r)]*r)
    return [pi for pi in S if P.is_poset_morphism(lambda i: pi[i], P)]

# (* Mathematica code using depth first recursion *)

# endoCount::usage := "If rel is a reflexive relation on 1,2,..,Max[rel], then \
# endoCount[rel] is the number of self maps f for which {f[a],f[b]} is in the \
# relation whenever {a,b} is in the relation.”

# endoCount[rel_] := morphismCount[rel, rel]

# morphismCount[rel1_List, rel2_List] :=
# Module[{max1 = Max[rel1], max2 = Max[rel2], down, checkdown, num, ans},
#  Do[down[i] = Select[rel1, Max[#] == i &], {i, max1}];
#  checkdown[f_List] := AllTrue[down[Length[f]], MemberQ[rel2, f[[#1]]] &];
#  num[{}] := Sum[num[{i}], {i, max2}];
#  num[f_List] :=
#   If[checkdown[f],
#    If[Length[f] == max1, 1, Sum[num[Append[f, i]], {i, max2}]], 0];
#  ans = num[{}]; Clear[checkdown, down]; ans]

# (* example *)

# endoCount[{{1, 1}, {1, 2}, {1, 10}, {2, 2}, {3, 2}, {3, 3}, {3,
#   4}, {4, 4}, {5, 4}, {5, 5}, {5, 6}, {6, 6}, {7, 6}, {7, 7}, {7,
#   8}, {8, 8}, {9, 8}, {9, 9}, {9, 10}, {10, 10}}]

# (* gives 10030 *)

# (* Sanity checks are to be added. For instance we must have 
#   Union @@ rel == Range[Max[rel]]
*)# *)

def filling(w):
    w = list(w)
    f = list(w)
    for i in range(len(w)-1):
        if w[i:i+2] == [0, 1]:
            f[i:i+2] = [1, 0]
    if w[len(w)-1] == 0:
        f[len(w)-1] = 1
    return Words([0,1])(f)

def degree(w):
    d = 0
    while w.count(0):
        w = filling(w)
        d += 1
    return d

def statistic(w): 
    return degree(w)

def filling(w):
    w = list(w)
    f = list(w)
    for i in range(len(w)-1):
        if w[i:i+2] == [0, 1]:
            f[i:i+2] = [1, 0]
    if w[len(w)-1] == 0:
        f[len(w)-1] = 1
    return Words([0,1])(f)

def statisticdegree(w):
    d = 0
    while w.count(0):
        w = filling(w)
        d += 1
    return d

def statistic(w): 
    return degree(w)

from collections import defaultdict

def word(pi):
    w = [0]*(len(pi)-1)
    for i in pi.descents(from_zero=True):
        w[i] = 1
    return Words([0,1])(w)

@cached_function
def preimages(n):
    D = defaultdict(int)
    for pi in Permutations(n):
        D[word(pi)] += 1
    return D

def statistic(word):
    return preimages(len(word)+Integer(1))[word]

from collections import defaultdict

def word(pi):
    w = [0]*(len(pi)-1)
    for i in pi.descents(from_zero=True):
        w[i] = 1
    return Words([0,1])(w)

@cached_function
def preimages(n):
    D = defaultdict(int)
    for pi in Permutations(n):
        D[word(pi)] += 1
    return D

def statistic(word):
    return preimages(len(word)+Integer(1))[word]

def statistic(B):
    if B == BinaryTree():
        return 0
    if B == BinaryTree([]):
        return 1
    return sum(statistic(S) for S in B)

# alternative implementation
# BinaryTree -> Bool
# def is_leaf(b):
#     return b == BinaryTree('.')

# # BinaryTree -> Int
# def leaf_count(b):
#     if is_leaf(b):
#         return 0
#     if is_leaf(b[0]) and is_leaf(b[1]):
#         return 1
#     return 0

# # (Num a) => a -> a -> a
# def add(a, b):
#     return a + b

# # (a -> a -> a) -> (BinaryTree -> a) -> BinaryTree -> a
# def fold_tree(f, g, b):
#     if is_leaf(b):
#         return g(b)
#     else:
#         return f(g(b), f(fold_tree(f, g, b[0]), fold_tree(f, g, b[1])))

# for n in [0..7]:
#     for b in BinaryTrees(n):
#         indiv_count = fold_tree(add, leaf_count, b)
#         print b, '=>', indiv_count

def statistic(B):
    if B == BinaryTree():
        return 0
    if B == BinaryTree([]):
        return 1
    return sum(statistic(S) for S in B)

# alternative implementation
# BinaryTree -> Bool
# def is_leaf(b):
#     return b == BinaryTree('.')

# # BinaryTree -> Int
# def leaf_count(b):
#     if is_leaf(b):
#         return 0
#     if is_leaf(b[0]) and is_leaf(b[1]):
#         return 1
#     return 0

# # (Num a) => a -> a -> a
# def add(a, b):
#     return a + b

# # (a -> a -> a) -> (BinaryTree -> a) -> BinaryTree -> a
# def fold_tree(f, g, b):
#     if is_leaf(b):
#         return g(b)
#     else:
#         return f(g(b), f(fold_tree(f, g, b[0]), fold_tree(f, g, b[1])))

# for n in [0..7]:
#     for b in BinaryTrees(n):
#         indiv_count = fold_tree(add, leaf_count, b)
#         print b, '=>', indiv_count

def statistic(G):
    from sage.graphs.convexity_properties import ConvexityProperties
    return sum(1 for V in subsets(G) if set(ConvexityProperties(G).hull(V)) != set(V))

def statistic(G):
    from sage.graphs.convexity_properties import ConvexityProperties
    return sum(1 for V in subsets(G) if set(ConvexityProperties(G).hull(V)) != set(V))

def statistic(G):
    from sage.graphs.convexity_properties import ConvexityProperties
    return sum(1 for V in subsets(G) if set(ConvexityProperties(G).hull(V)) == set(V))

def statistic(G):
    from sage.graphs.convexity_properties import ConvexityProperties
    return sum(1 for V in subsets(G) if set(ConvexityProperties(G).hull(V)) == set(V))

def graph_polynomial(G):
    G.relabel(inplace=False)
    n = G.num_verts()
    R = PolynomialRing(ZZ, "x", n)
    X = R.gens()
    p = R(1)
    for a, b in G.edges(labels=False):
        if a > b:
            a, b = b, a
        p *= X[a] - X[b]
    return p

def statistic(G):
    """
    sage: [statistic(graphs.CompleteBipartiteGraph(n,n)) for n in range(1, 5)]
    [2, 2, 3, 3]
    
    sage: lG = [(graphs.CycleGraph(2*n+1).cartesian_product(graphs.PathGraph(k))) for k in range(1,3) for n in range(1,3)]
    sage: [statistic(G) for G in lG if G.num_verts() <= 6]
    [3, 3, 3]
    """
    lm = graph_polynomial(G).monomials()
    return min(max(m.degrees()) for m in lm)+1

def graph_polynomial(G):
    G.relabel(inplace=False)
    n = G.num_verts()
    R = PolynomialRing(ZZ, "x", n)
    X = R.gens()
    p = R(1)
    for a, b in G.edges(labels=False):
        if a > b:
            a, b = b, a
        p *= X[a] - X[b]
    return p

def statistic(G):
    """
    sage: [statistic(graphs.CompleteBipartiteGraph(n,n)) for n in range(1, 5)]
    [2, 2, 3, 3]
    
    sage: lG = [(graphs.CycleGraph(2*n+1).cartesian_product(graphs.PathGraph(k))) for k in range(1,3) for n in range(1,3)]
    sage: [statistic(G) for G in lG if G.num_verts() <= 6]
    [3, 3, 3]
    """
    lm = graph_polynomial(G).monomials()
    return min(max(m.degrees()) for m in lm)+1

def is_dominating(G,V):        
    return set(G.vertices()) == set(V).union(*[G.neighbors(v) for v in V])

def minimal_dominating_sets(G):
    Vs = list(Subsets(G.vertices()))
    mdoms = []
    while Vs:
        V = Vs.pop(0)
        if is_dominating(G,V):
            mdoms.append(V)
            Vs = list(filter( lambda X: not X.issuperset(V), Vs))
    return mdoms

def statistic(G):
    return len(minimal_dominating_sets(G))

def is_dominating(G,V):        
    return set(G.vertices()) == set(V).union(*[G.neighbors(v) for v in V])

def minimal_dominating_sets(G):
    Vs = list(Subsets(G.vertices()))
    mdoms = []
    while Vs:
        V = Vs.pop(0)
        if is_dominating(G,V):
            mdoms.append(V)
            Vs = list(filter( lambda X: not X.issuperset(V), Vs))
    return mdoms

def statistic(G):
    return len(minimal_dominating_sets(G))

def statistic(G):
    """Return the smallest number such that Alice has a winning strategy.

    EXAMPLES:

    From http://www.dlsu.edu.ph/conferences/dlsu_research_congress/2014/_pdf/proceedings/TPHS-I-011-FT.pdf::

        sage: [statistic(graphs.PathGraph(r)) for r in range(2,12)]
        [2, 2, 3, 3, 3, 3, 3, 3, 3, 3]
    
        sage: [statistic(graphs.CycleGraph(r)) for r in range(2,10)]
        [2, 3, 3, 3, 3, 3, 3, 3]
    
        sage: [statistic(graphs.StarGraph(r)) for r in range(2,10)]
        [2, 2, 2, 2, 2, 2, 2, 2]
    
        sage: [statistic(graphs.CompleteGraph(r)) for r in range(2,10)]
        [2, 3, 4, 5, 6, 7, 8, 9]

    Proposition 4 of the reference seems to miss some corner cases.
    W(3) and W(4) are complete graphs, so Proposition 5 applies.  

    However, W(5) and W(7) seem to be special::

        sage: [statistic(graphs.WheelGraph(r)) for r in range(5,12)] 
        [3, 4, 3, 4, 4, 4, 4]
    
        sage: [[statistic(graphs.CompleteBipartiteGraph(r, s)) for r in range(1,6)] for s in range(1,6)]
        [[2, 2, 2, 2, 2],
        [2, 3, 3, 3, 3],
        [2, 3, 3, 3, 3],
        [2, 3, 3, 3, 3],
        [2, 3, 3, 3, 3]]

    The following contradicts Theorem 1, but is in agreement with
    http://ssltest.cs.umd.edu/~gasarch/TOPICS/graphcolgame/tree3.pdf,
    Theorem 3::

        sage: statistic(graphs.PetersenGraph())
        4

    from wikipedia::

        sage: [statistic(graphs.CompleteGraph(2).cartesian_product(graphs.PathGraph(r))) for r in range(1,6)]
        [2, 3, 3, 4, 4]
    
        sage: [statistic(graphs.CompleteGraph(2).cartesian_product(graphs.CycleGraph(r))) for r in range(1,5)]
        [2, 3, 4, 4]
    
        sage: [statistic(graphs.CompleteGraph(2).cartesian_product(graphs.CompleteGraph(r))) for r in range(1,5)]
        [2, 3, 4, 5]

    """
    G = G.relabel(immutable=True, inplace=False)
    for k in range(G.chromatic_number(), max(G.degree())+2):
        if Alice_wins(G, k):
            return k

def normalize_colours(D):
    """
    From left to right, use small colours first.
    """
    pi = [None]*(max(cc for cc in D if cc is not None)+1)
    i = 0
    for e in D:
        if e is not None and pi[e] is None:
            pi[e] = i
            i += 1
    return tuple([None if e is None else pi[e] for e in D])

@cached_function
def Alice_wins(G, k, C=None):
    """Return whether C is a winning position for Alice, who wants to
    properly color the graph.

    Expect that the vertices are labelled 0,1,...

    """
    def children(D):
        if all(c is None for c in D):
            colours = 1
        else:
            colours = min(k, max(cc for cc in D if cc is not None)+2)
        for v, c in enumerate(D):
            if c is None:
                for d in range(colours):
                    if not any(D[u] == d for u in G.neighbor_iterator(v)):
                        yield D[:v] + tuple([d]) + D[v+1:]

    if C is None:
        C = tuple([None]*G.num_verts())

    uncoloured_vertices = C.count(None)
    if uncoloured_vertices == 0:
        return True

    # let Alice colour a vertex
    for D in children(C):
        # if C was missing precisely one vertex, and we could
        # colour it, Alice wins
        if uncoloured_vertices == 1:
            return True
        # Alice wins if there is a D such that all moves of Bob make Alice win
        has_colouring = False
        for E in children(D):
            has_colouring = True
            if not Alice_wins(G, k, normalize_colours(E)):
                break
        else:
            if has_colouring:
                return True
    return False

def statistic(G):
    """Return the smallest number such that Alice has a winning strategy.

    EXAMPLES:

    From http://www.dlsu.edu.ph/conferences/dlsu_research_congress/2014/_pdf/proceedings/TPHS-I-011-FT.pdf::

        sage: [statistic(graphs.PathGraph(r)) for r in range(2,12)]
        [2, 2, 3, 3, 3, 3, 3, 3, 3, 3]
    
        sage: [statistic(graphs.CycleGraph(r)) for r in range(2,10)]
        [2, 3, 3, 3, 3, 3, 3, 3]
    
        sage: [statistic(graphs.StarGraph(r)) for r in range(2,10)]
        [2, 2, 2, 2, 2, 2, 2, 2]
    
        sage: [statistic(graphs.CompleteGraph(r)) for r in range(2,10)]
        [2, 3, 4, 5, 6, 7, 8, 9]

    Proposition 4 of the reference seems to miss some corner cases.
    W(3) and W(4) are complete graphs, so Proposition 5 applies.  

    However, W(5) and W(7) seem to be special::

        sage: [statistic(graphs.WheelGraph(r)) for r in range(5,12)] 
        [3, 4, 3, 4, 4, 4, 4]
    
        sage: [[statistic(graphs.CompleteBipartiteGraph(r, s)) for r in range(1,6)] for s in range(1,6)]
        [[2, 2, 2, 2, 2],
        [2, 3, 3, 3, 3],
        [2, 3, 3, 3, 3],
        [2, 3, 3, 3, 3],
        [2, 3, 3, 3, 3]]

    The following contradicts Theorem 1, but is in agreement with
    http://ssltest.cs.umd.edu/~gasarch/TOPICS/graphcolgame/tree3.pdf,
    Theorem 3::

        sage: statistic(graphs.PetersenGraph())
        4

    from wikipedia::

        sage: [statistic(graphs.CompleteGraph(2).cartesian_product(graphs.PathGraph(r))) for r in range(1,6)]
        [2, 3, 3, 4, 4]
    
        sage: [statistic(graphs.CompleteGraph(2).cartesian_product(graphs.CycleGraph(r))) for r in range(1,5)]
        [2, 3, 4, 4]
    
        sage: [statistic(graphs.CompleteGraph(2).cartesian_product(graphs.CompleteGraph(r))) for r in range(1,5)]
        [2, 3, 4, 5]

    """
    G = G.relabel(immutable=True, inplace=False)
    for k in range(G.chromatic_number(), max(G.degree())+2):
        if Alice_wins(G, k):
            return k

def normalize_colours(D):
    """
    From left to right, use small colours first.
    """
    pi = [None]*(max(cc for cc in D if cc is not None)+1)
    i = 0
    for e in D:
        if e is not None and pi[e] is None:
            pi[e] = i
            i += 1
    return tuple([None if e is None else pi[e] for e in D])

@cached_function
def Alice_wins(G, k, C=None):
    """Return whether C is a winning position for Alice, who wants to
    properly color the graph.

    Expect that the vertices are labelled 0,1,...

    """
    def children(D):
        if all(c is None for c in D):
            colours = 1
        else:
            colours = min(k, max(cc for cc in D if cc is not None) + +2)
        for v, c in enumerate(D):
            if c is None:
                for d in range(colours):
                    if not any(D[u] == d for u in G.neighbor_iterator(v)):
                        yield D[:v] + tuple([d]) + D[v+1:]

    if C is None:
        C = tuple([None]*G.num_verts())

    uncoloured_vertices = C.count(None)
    if uncoloured_vertices == 0:
        return True

    # let Alice colour a vertex
    for D in children(C):
        # if C was missing precisely one vertex, and we could
        # colour it, Alice wins
        if uncoloured_vertices == 1:
            return True
        # Alice wins if there is a D such that all moves of Bob make Alice win
        has_colouring = False
        for E in children(D):
            has_colouring = True
            if not Alice_wins(G, k, normalize_colours(E)):
                break
        else:
            if has_colouring:
                return True
    return False

def statistic(ct):
    from sage.combinat.crystals.kirillov_reshetikhin import KirillovReshetikhinCrystal
    from sage.combinat.crystals.highest_weight_crystals import HighestWeightCrystal
    if ct.letter in ['A','B','C','D']:
        return crystals.Letters(ct).cardinality()
    elif ct.letter == 'E':
        if ct.rank() == 6:
            B = KirillovReshetikhinCrystal(['E',6,1], 1,1)
            return B.cardinality()
        elif ct.rank() == 7:
            La = ct.root_system().weight_lattice().fundamental_weight(1)     
            T = HighestWeightCrystal(La)
            return T.cardinality()
        elif ct.rank() == 8:
            RC = RiggedConfigurations(['E',8,1], [[1,1]])
            return RC.cardinality()
    elif ct.letter == 'F' and ct.rank() == 4:
        RC = RiggedConfigurations(['F',4,1], [[1,1]])
        return RC.cardinality()

def statistic(ct):
    from sage.combinat.crystals.kirillov_reshetikhin import KirillovReshetikhinCrystal
    from sage.combinat.crystals.highest_weight_crystals import HighestWeightCrystal
    if ct.letter in ['A','B','C','D']:
        return crystals.Letters(ct).cardinality()
    elif ct.letter == 'E':
        if ct.rank() == 6:
            B = KirillovReshetikhinCrystal(['E',6,1], 1,1)
            return B.cardinality()
        elif ct.rank() == 7:
            La = Cct.root_system().weight_lattice().fundamental_weight(1)     
            T = HighestWeightCrystal(La)
            return T.cardinality()
        elif ct.rank() == 8:
            RC = RiggedConfigurations(['E',8,1], [[1,1]])
            return RCT.cardinality()
    elif ct.letter == 'F' and ct.rank() == 4:
        RC = RiggedConfigurations(['F',4,1], [[1,1]])
        return RC.cardinality()

def statistic(p):
    return depth_of_1(set_partition_to_tree(p))

def depth_of_1(T):
    if T.label() == 1:
        return 0
    if T.node_number() == 1:
        return infinity
    else:
        return 1 + min(depth_of_1(C) for C in T)

def set_partition_to_tree(p):
    # the children of the smallest label are the largest remaining
    # element and its partner
    trees = [LabelledRootedTree([LabelledRootedTree([], label=e) for e in b]) for b in p]
    max_label = p.size()+1-len(p) # last labelled node
    while len(trees) > 1:
        max_label += 1
        # find tree with smallest child and all children smaller than max_label
        A = sorted([T for T in trees if max(C.label() for C in T) < max_label],
                   key = lambda T: min(C.label() for C in T))[0]
        trees.remove(A)
        # give it's root node the new label
        A = LabelledRootedTree(A, label=max_label)
        # find tree with child having label max_label
        B = next(T for T in trees
             if any(C.label() == max_label for C in T))
        trees.remove(B)

        # replace the child of B carrying max_label with A
        C = LabelledRootedTree([A] + [T for T in B if T.label() != max_label])

        trees.append(C)

    return trees[0]    

def statistic(p):
    return depth_of_1(set_partition_to_tree(p))

def depth_of_1(T):
    if T.label() == 1:
        return 0
    if T.node_number() == 1:
        return infinity
    else:
        return 1 + min(depth_of_1(C) for C in T)

def set_partition_to_tree(p):
    # the children of the smallest label are the largest remaining
    # element and its partner
    trees = [LabelledRootedTree([LabelledRootedTree([], label=e) for e in b]) for b in p]
    max_label = p.size()+1-len(p) # last labelled node
    while len(trees) > 1:
        max_label += 1
        # find tree with smallest child and all children smaller than max_label
        A = sorted([T for T in trees if max(C.label() for C in T) < max_label],
                   key = lambda T: min(C.label() for C in T))[0]
        trees.remove(A)
        # give it's root node the new label
        A = LabelledRootedTree(A, label=max_label)
        # find tree with child having label max_label
        B = next(T for T in trees
             if any(C.label() == max_label for C in T)).next()
        trees.remove(B)

        # replace the child of B carrying max_label with A
        C = LabelledRootedTree([A] + [T for T in B if T.label() != max_label])

        trees.append(C)

    return trees[0]    

def statistic(m):
    return largest_label_with_larger_sister(matching_to_tree(m))

def largest_label_with_larger_sister(T):
    T0 = T[0].node_number()
    T1 = T[1].node_number()
    if T0 == T1 == 1:
        return min(T[0].label(), T[1].label())
    elif T0 == 1 and T1 > 1:
        return largest_label_with_larger_sister(T[1])
    elif T0 > 1 and T0 == 1:
        return largest_label_with_larger_sister(T[0])
    elif T0 > 1 and T1 > 1:
        return max(largest_label_with_larger_sister(T[0]),
                   largest_label_with_larger_sister(T[1]))

def matching_to_tree(m):
    """
    INPUT:

    - m, a PerfectMatching on {1,...,2n}.

    OUTPUT:

    a decreasingly labelled, unordered full binary tree with n+1 leaves.

    EXAMPLES::

        sage: m = PerfectMatching([(1,4),(2,9),(3,10),(5,7),(6,8),(11,12)])
        sage: ascii_art(matching_to_tree(m))
    """
    # the children of the smallest label are the largest remaining
    # element and its partner
    trees = [LabelledRootedTree([LabelledRootedTree([], label=i),
                                 LabelledRootedTree([], label=j)]) for i, j in m]
    max_label = m.size()//2+1 # last labelled node

    while len(trees) > 1:
        max_label += 1
        # find tree with smallest child and both children smaller than max_label
        A = sorted([T for T in trees if max(T[0].label(), T[1].label()) < max_label],
                   key = lambda T: min(T[0].label(), T[1].label()))[0]
        trees.remove(A)
        # give it's root node the new label
        A = LabelledRootedTree(A, label=max_label)
        # find tree with child having label max_label
        B = next(T for T in trees
             if T[0].label() == max_label or T[1].label() == max_label)
        trees.remove(B)
        # replace B with [B[0], A] or [B[1], A]
        if B[0].label() == max_label:
            C = LabelledRootedTree([A, B[1]])
        else:
            C = LabelledRootedTree([A, B[0]])

        trees.append(C)

    return trees[0]

def statistic(m):
    return largest_label_with_larger_sister(matching_to_tree(m))

def largest_label_with_larger_sister(T):
    T0 = T[0].node_number()
    T1 = T[1].node_number()
    if T0 == T1 == 1:
        return min(T[0].label(), T[1].label())
    elif T0 == 1 and T1 > 1:
        return largest_label_with_larger_sister(T[1])
    elif T0 > 1 and T0 == 1:
        return largest_label_with_larger_sister(T[0])
    elif T0 > 1 and T1 > 1:
        return max(largest_label_with_larger_sister(T[0]),
                   largest_label_with_larger_sister(T[1]))

def matching_to_tree(m):
    """
    INPUT:

    - m, a PerfectMatching on {1,...,2n}.

    OUTPUT:

    a decreasingly labelled, unordered full binary tree with n+1 leaves.

    EXAMPLES::

        sage: m = PerfectMatching([(1,4),(2,9),(3,10),(5,7),(6,8),(11,12)])
        sage: ascii_art(matching_to_tree(m))
    """
    # the children of the smallest label are the largest remaining
    # element and its partner
    trees = [LabelledRootedTree([LabelledRootedTree([], label=i),
                                 LabelledRootedTree([], label=j)]) for i, j in m]
    max_label = m.size()//2+1 # last labelled node

    while len(trees) > 1:
        max_label += 1
        # find tree with smallest child and both children smaller than max_label
        A = sorted([T for T in trees if max(T[0].label(), T[1].label()) < max_label],
                   key = lambda T: min(T[0].label(), T[1].label()))[0]
        trees.remove(A)
        # give it's root node the new label
        A = LabelledRootedTree(A, label=max_label)
        # find tree with child having label max_label
        B = next(T for T in trees
             if T[0].label() == max_label or T[1].label() == max_label).next()
        trees.remove(B)
        # replace B with [B[0], A] or [B[1], A]
        if B[0].label() == max_label:
            C = LabelledRootedTree([A, B[1]])
        else:
            C = LabelledRootedTree([A, B[0]])

        trees.append(C)

    return trees[0]

def statistic(m):
    return max(sister_of_1(matching_to_tree(m)).leaf_labels())

def sister_of_1(T):
    if T[0].label() == 1:
        return T[1]
    if T[1].label() == 1:
        return T[0]
    if 1 in T[0].leaf_labels():
        return sister_of_1(T[0])
    if 1 in T[1].leaf_labels():
        return sister_of_1(T[1])

def matching_to_tree(m):
    """
    INPUT:

    - m, a PerfectMatching on {1,...,2n}.

    OUTPUT:

    a decreasingly labelled, unordered full binary tree with n+1 leaves.

    EXAMPLES::

        sage: m = PerfectMatching([(1,4),(2,9),(3,10),(5,7),(6,8),(11,12)])
        sage: ascii_art(matching_to_tree(m))
             ____None___
            /          /
          _11__      _12_
         /    /     /   /
        3   _10_   6   8_
           /   /      / /
          2   9_     1 4
             / /
            5 7

    """
    # the children of the smallest label are the largest remaining
    # element and its partner
    trees = [LabelledRootedTree([LabelledRootedTree([], label=i),
                                 LabelledRootedTree([], label=j)]) for i, j in m]
    max_label = m.size()//2+1 # last labelled node

    while len(trees) > 1:
        max_label += 1
        # find tree with smallest child and both children smaller than max_label
        A = sorted([T for T in trees if max(T[0].label(), T[1].label()) < max_label],
                   key = lambda T: min(T[0].label(), T[1].label()))[0]
        trees.remove(A)
        # give it's root node the new label
        A = LabelledRootedTree(A, label=max_label)
        # find tree with child having label max_label
        B = next(T for T in trees
             if T[0].label() == max_label or T[1].label() == max_label)
        trees.remove(B)
        # replace B with [B[0], A] or [B[1], A]
        if B[0].label() == max_label:
            C = LabelledRootedTree([A, B[1]])
        else:
            C = LabelledRootedTree([A, B[0]])

        trees.append(C)

    return trees[0]

def statistic(m):
    return max(sister_of_1(matching_to_tree(m)).leaf_labels())

def sister_of_1(T):
    if T[0].label() == 1:
        return T[1]
    if T[1].label() == 1:
        return T[0]
    if 1 in T[0].leaf_labels():
        return sister_of_1(T[0])
    if 1 in T[1].leaf_labels():
        return sister_of_1(T[1])

def matching_to_tree(m):
    """
    INPUT:

    - m, a PerfectMatching on {1,...,2n}.

    OUTPUT:

    a decreasingly labelled, unordered full binary tree with n+1 leaves.

    EXAMPLES::

        sage: m = PerfectMatching([(1,4),(2,9),(3,10),(5,7),(6,8),(11,12)])
        sage: ascii_art(matching_to_tree(m))
             ____None___
            /          /
          _11__      _12_
         /    /     /   /
        3   _10_   6   8_
           /   /      / /
          2   9_     1 4
             / /
            5 7

    """
    # the children of the smallest label are the largest remaining
    # element and its partner
    trees = [LabelledRootedTree([LabelledRootedTree([], label=i),
                                 LabelledRootedTree([], label=j)]) for i, j in m]
    max_label = m.size()//2+1 # last labelled node

    while len(trees) > 1:
        max_label += 1
        # find tree with smallest child and both children smaller than max_label
        A = sorted([T for T in trees if max(T[0].label(), T[1].label()) < max_label],
                   key = lambda T: min(T[0].label(), T[1].label()))[0]
        trees.remove(A)
        # give it's root node the new label
        A = LabelledRootedTree(A, label=max_label)
        # find tree with child having label max_label
        B = next(T for T in trees
             if T[0].label() == max_label or T[1].label() == max_label).next()
        trees.remove(B)
        # replace B with [B[0], A] or [B[1], A]
        if B[0].label() == max_label:
            C = LabelledRootedTree([A, B[1]])
        else:
            C = LabelledRootedTree([A, B[0]])

        trees.append(C)

    return trees[0]

def statistic(m):
    return min(sister_of_1(matching_to_tree(m)).leaf_labels())

def sister_of_1(T):
    if T[0].label() == 1:
        return T[1]
    if T[1].label() == 1:
        return T[0]
    if 1 in T[0].leaf_labels():
        return sister_of_1(T[0])
    if 1 in T[1].leaf_labels():
        return sister_of_1(T[1])

def matching_to_tree(m):
    """
    INPUT:

    - m, a PerfectMatching on {1,...,2n}.

    OUTPUT:

    a decreasingly labelled, unordered full binary tree with n+1 leaves.

    EXAMPLES::

        sage: m = PerfectMatching([(1,4),(2,9),(3,10),(5,7),(6,8),(11,12)])
        sage: ascii_art(matching_to_tree(m))
             ____None___
            /          /
          _11__      _12_
         /    /     /   /
        3   _10_   6   8_
           /   /      / /
          2   9_     1 4
             / /
            5 7

    """
    # the children of the smallest label are the largest remaining
    # element and its partner
    trees = [LabelledRootedTree([LabelledRootedTree([], label=i),
                                 LabelledRootedTree([], label=j)]) for i, j in m]
    max_label = m.size()//2+1 # last labelled node

    while len(trees) > 1:
        max_label += 1
        # find tree with smallest child and both children smaller than max_label
        A = sorted([T for T in trees if max(T[0].label(), T[1].label()) < max_label],
                   key = lambda T: min(T[0].label(), T[1].label()))[0]
        trees.remove(A)
        # give it's root node the new label
        A = LabelledRootedTree(A, label=max_label)
        # find tree with child having label max_label
        B = next(T for T in trees
             if T[0].label() == max_label or T[1].label() == max_label)
        trees.remove(B)
        # replace B with [B[0], A] or [B[1], A]
        if B[0].label() == max_label:
            C = LabelledRootedTree([A, B[1]])
        else:
            C = LabelledRootedTree([A, B[0]])

        trees.append(C)

    return trees[0]

def statistic(m):
    return min(sister_of_1(matching_to_tree(m)).leaf_labels())

def sister_of_1(T):
    if T[0].label() == 1:
        return T[1]
    if T[1].label() == 1:
        return T[0]
    if 1 in T[0].leaf_labels():
        return sister_of_1(T[0])
    if 1 in T[1].leaf_labels():
        return sister_of_1(T[1])

def matching_to_tree(m):
    """
    INPUT:

    - m, a PerfectMatching on {1,...,2n}.

    OUTPUT:

    a decreasingly labelled, unordered full binary tree with n+1 leaves.

    EXAMPLES::

        sage: m = PerfectMatching([(1,4),(2,9),(3,10),(5,7),(6,8),(11,12)])
        sage: ascii_art(matching_to_tree(m))
             ____None___
            /          /
          _11__      _12_
         /    /     /   /
        3   _10_   6   8_
           /   /      / /
          2   9_     1 4
             / /
            5 7

    """
    # the children of the smallest label are the largest remaining
    # element and its partner
    trees = [LabelledRootedTree([LabelledRootedTree([], label=i),
                                 LabelledRootedTree([], label=j)]) for i, j in m]
    max_label = m.size()//2+1 # last labelled node

    while len(trees) > 1:
        max_label += 1
        # find tree with smallest child and both children smaller than max_label
        A = sorted([T for T in trees if max(T[0].label(), T[1].label()) < max_label],
                   key = lambda T: min(T[0].label(), T[1].label()))[0]
        trees.remove(A)
        # give it's root node the new label
        A = LabelledRootedTree(A, label=max_label)
        # find tree with child having label max_label
        B = next(T for T in trees
             if T[0].label() == max_label or T[1].label() == max_label).next()
        trees.remove(B)
        # replace B with [B[0], A] or [B[1], A]
        if B[0].label() == max_label:
            C = LabelledRootedTree([A, B[1]])
        else:
            C = LabelledRootedTree([A, B[0]])

        trees.append(C)

    return trees[0]

def statistic(m):
    return leaves_of_sister_of_1(matching_to_tree(m))

def leaves_of_sister_of_1(T):
    if len(T) == 0:
        if T.label() == 1:
            return 0
        return 1
    if T[0].label() == 1:
        return (T[1].node_number()+1)/2
    if T[1].label() == 1:
        return (T[0].node_number()+1)/2
    if 1 in T[0].leaf_labels():
        return leaves_of_sister_of_1(T[0])
    if 1 in T[1].leaf_labels():
        return leaves_of_sister_of_1(T[1])

def matching_to_tree(m):
    """
    INPUT:

    - m, a PerfectMatching on {1,...,2n}.

    OUTPUT:

    a decreasingly labelled, unordered full binary tree with n+1 leaves.

    EXAMPLES::

        sage: m = PerfectMatching([(1,4),(2,9),(3,10),(5,7),(6,8),(11,12)])
        sage: ascii_art(matching_to_tree(m))
             ____None___
            /          /
          _11__      _12_
         /    /     /   /
        3   _10_   6   8_
           /   /      / /
          2   9_     1 4
             / /
            5 7

    """
    # the children of the smallest label are the largest remaining
    # element and its partner
    trees = [LabelledRootedTree([LabelledRootedTree([], label=i),
                                 LabelledRootedTree([], label=j)]) for i, j in m]
    max_label = m.size()//2+1 # last labelled node

    while len(trees) > 1:
        max_label += 1
        # find tree with smallest child and both children smaller than max_label
        A = sorted([T for T in trees if max(T[0].label(), T[1].label()) < max_label],
                   key = lambda T: min(T[0].label(), T[1].label()))[0]
        trees.remove(A)
        # give it's root node the new label
        A = LabelledRootedTree(A, label=max_label)
        # find tree with child having label max_label
        B = next(T for T in trees
             if T[0].label() == max_label or T[1].label() == max_label)
        trees.remove(B)
        # replace B with [B[0], A] or [B[1], A]
        if B[0].label() == max_label:
            C = LabelledRootedTree([A, B[1]])
        else:
            C = LabelledRootedTree([A, B[0]])

        trees.append(C)

    return trees[0]

def statistic(m):
    return leaves_of_sister_of_1(matching_to_tree(m))

def leaves_of_sister_of_1(T):
    if len(T) == 0:
        if T.label() == 1:
            return 0
        return 1
    if T[0].label() == 1:
        return (T[1].node_number()+1)/2
    if T[1].label() == 1:
        return (T[0].node_number()+1)/2
    if 1 in T[0].leaf_labels():
        return leaves_of_sister_of_1(T[0])
    if 1 in T[1].leaf_labels():
        return leaves_of_sister_of_1(T[1])

def matching_to_tree(m):
    """
    INPUT:

    - m, a PerfectMatching on {1,...,2n}.

    OUTPUT:

    a decreasingly labelled, unordered full binary tree with n+1 leaves.

    EXAMPLES::

        sage: m = PerfectMatching([(1,4),(2,9),(3,10),(5,7),(6,8),(11,12)])
        sage: ascii_art(matching_to_tree(m))
             ____None___
            /          /
          _11__      _12_
         /    /     /   /
        3   _10_   6   8_
           /   /      / /
          2   9_     1 4
             / /
            5 7

    """
    # the children of the smallest label are the largest remaining
    # element and its partner
    trees = [LabelledRootedTree([LabelledRootedTree([], label=i),
                                 LabelledRootedTree([], label=j)]) for i, j in m]
    max_label = m.size()//2+1 # last labelled node

    while len(trees) > 1:
        max_label += 1
        # find tree with smallest child and both children smaller than max_label
        A = sorted([T for T in trees if max(T[0].label(), T[1].label()) < max_label],
                   key = lambda T: min(T[0].label(), T[1].label()))[0]
        trees.remove(A)
        # give it's root node the new label
        A = LabelledRootedTree(A, label=max_label)
        # find tree with child having label max_label
        B = next(T for T in trees
             if T[0].label() == max_label or T[1].label() == max_label).next()
        trees.remove(B)
        # replace B with [B[0], A] or [B[1], A]
        if B[0].label() == max_label:
            C = LabelledRootedTree([A, B[1]])
        else:
            C = LabelledRootedTree([A, B[0]])

        trees.append(C)

    return trees[0]

def statistic(m):
    return trivial_trees_on_path_to_1(matching_to_tree(m))

def trivial_trees_on_path_to_1(T):
    if T.node_number() == 1:
        return 0
    if 1 in T[0].leaf_labels():
        if T[1].node_number() == 1:
            return 1 + trivial_trees_on_path_to_1(T[0])
        else:
            return trivial_trees_on_path_to_1(T[0])
    else:
        if T[0].node_number() == 1:
            return 1 + trivial_trees_on_path_to_1(T[1])
        else:
            return trivial_trees_on_path_to_1(T[1])

def matching_to_tree(m):
    """
    INPUT:

    - m, a PerfectMatching on {1,...,2n}.

    OUTPUT:

    a decreasingly labelled, unordered full binary tree with n+1 leaves.

    EXAMPLES::

        sage: m = PerfectMatching([(1,4),(2,9),(3,10),(5,7),(6,8),(11,12)])
        sage: ascii_art(matching_to_tree(m))
             ____None___
            /          /
          _11__      _12_
         /    /     /   /
        3   _10_   6   8_
           /   /      / /
          2   9_     1 4
             / /
            5 7

    """
    # the children of the smallest label are the largest remaining
    # element and its partner
    trees = [LabelledRootedTree([LabelledRootedTree([], label=i),
                                 LabelledRootedTree([], label=j)]) for i, j in m]
    max_label = m.size()//2+1 # last labelled node

    while len(trees) > 1:
        max_label += 1
        # find tree with smallest child and both children smaller than max_label
        A = sorted([T for T in trees if max(T[0].label(), T[1].label()) < max_label],
                   key = lambda T: min(T[0].label(), T[1].label()))[0]
        trees.remove(A)
        # give it's root node the new label
        A = LabelledRootedTree(A, label=max_label)
        # find tree with child having label max_label
        B = next(T for T in trees
             if T[0].label() == max_label or T[1].label() == max_label)
        trees.remove(B)
        # replace B with [B[0], A] or [B[1], A]
        if B[0].label() == max_label:
            C = LabelledRootedTree([A, B[1]])
        else:
            C = LabelledRootedTree([A, B[0]])

        trees.append(C)

    return trees[0]

def statistic(m):
    return trivial_trees_on_path_to_1(matching_to_tree(m))

def trivial_trees_on_path_to_1(T):
    if T.node_number() == 1:
        return 0
    if 1 in T[0].leaf_labels():
        if T[1].node_number() == 1:
            return 1 + trivial_trees_on_path_to_1(T[0])
        else:
            return trivial_trees_on_path_to_1(T[0])
    else:
        if T[0].node_number() == 1:
            return 1 + trivial_trees_on_path_to_1(T[1])
        else:
            return trivial_trees_on_path_to_1(T[1])

def matching_to_tree(m):
    """
    INPUT:

    - m, a PerfectMatching on {1,...,2n}.

    OUTPUT:

    a decreasingly labelled, unordered full binary tree with n+1 leaves.

    EXAMPLES::

        sage: m = PerfectMatching([(1,4),(2,9),(3,10),(5,7),(6,8),(11,12)])
        sage: ascii_art(matching_to_tree(m))
             ____None___
            /          /
          _11__      _12_
         /    /     /   /
        3   _10_   6   8_
           /   /      / /
          2   9_     1 4
             / /
            5 7

    """
    # the children of the smallest label are the largest remaining
    # element and its partner
    trees = [LabelledRootedTree([LabelledRootedTree([], label=i),
                                 LabelledRootedTree([], label=j)]) for i, j in m]
    max_label = m.size()//2+1 # last labelled node

    while len(trees) > 1:
        max_label += 1
        # find tree with smallest child and both children smaller than max_label
        A = sorted([T for T in trees if max(T[0].label(), T[1].label()) < max_label],
                   key = lambda T: min(T[0].label(), T[1].label()))[0]
        trees.remove(A)
        # give it's root node the new label
        A = LabelledRootedTree(A, label=max_label)
        # find tree with child having label max_label
        B = next(T for T in trees
             if T[0].label() == max_label or T[1].label() == max_label).next()
        trees.remove(B)
        # replace B with [B[0], A] or [B[1], A]
        if B[0].label() == max_label:
            C = LabelledRootedTree([A, B[1]])
        else:
            C = LabelledRootedTree([A, B[0]])

        trees.append(C)

    return trees[0]

def statistic(m):
    T = matching_to_tree(m)
    if 1 in T[0].leaf_labels():
        return min(T[1].leaf_labels())
    return min(T[0].leaf_labels())

def matching_to_tree(m):
    """
    INPUT:

    - m, a PerfectMatching on {1,...,2n}.

    OUTPUT:

    a decreasingly labelled, unordered full binary tree with n+1 leaves.

    EXAMPLES::

        sage: m = PerfectMatching([(1,4),(2,9),(3,10),(5,7),(6,8),(11,12)])
        sage: ascii_art(matching_to_tree(m))
    """
    # the children of the smallest label are the largest remaining
    # element and its partner
    trees = [LabelledRootedTree([LabelledRootedTree([], label=i),
                                 LabelledRootedTree([], label=j)]) for i, j in m]
    max_label = m.size()//2+1 # last labelled node

    while len(trees) > 1:
        max_label += 1
        # find tree with smallest child and both children smaller than max_label
        A = sorted([T for T in trees if max(T[0].label(), T[1].label()) < max_label],
                   key = lambda T: min(T[0].label(), T[1].label()))[0]
        trees.remove(A)
        # give it's root node the new label
        A = LabelledRootedTree(A, label=max_label)
        # find tree with child having label max_label
        B = next(T for T in trees
             if T[0].label() == max_label or T[1].label() == max_label)
        trees.remove(B)
        # replace B with [B[0], A] or [B[1], A]
        if B[0].label() == max_label:
            C = LabelledRootedTree([A, B[1]])
        else:
            C = LabelledRootedTree([A, B[0]])

        trees.append(C)

    return trees[0]

def statistic(m):
    T = matching_to_tree(m)
    if 1 in T[0].leaf_labels():
        return min(T[1].leaf_labels())
    return min(T[0].leaf_labels())

def matching_to_tree(m):
    """
    INPUT:

    - m, a PerfectMatching on {1,...,2n}.

    OUTPUT:

    a decreasingly labelled, unordered full binary tree with n+1 leaves.

    EXAMPLES::

        sage: m = PerfectMatching([(1,4),(2,9),(3,10),(5,7),(6,8),(11,12)])
        sage: ascii_art(matching_to_tree(m))
    """
    # the children of the smallest label are the largest remaining
    # element and its partner
    trees = [LabelledRootedTree([LabelledRootedTree([], label=i),
                                 LabelledRootedTree([], label=j)]) for i, j in m]
    max_label = m.size()//2+1 # last labelled node

    while len(trees) > 1:
        max_label += 1
        # find tree with smallest child and both children smaller than max_label
        A = sorted([T for T in trees if max(T[0].label(), T[1].label()) < max_label],
                   key = lambda T: min(T[0].label(), T[1].label()))[0]
        trees.remove(A)
        # give it's root node the new label
        A = LabelledRootedTree(A, label=max_label)
        # find tree with child having label max_label
        B = next(T for T in trees
             if T[0].label() == max_label or T[1].label() == max_label).next()
        trees.remove(B)
        # replace B with [B[0], A] or [B[1], A]
        if B[0].label() == max_label:
            C = LabelledRootedTree([A, B[1]])
        else:
            C = LabelledRootedTree([A, B[0]])

        trees.append(C)

    return trees[0]

def statistic(m):
    T = matching_to_tree(m)
    if 1 in T[0].leaf_labels():
        return len(T[0].leaf_labels())
    return len(T[1].leaf_labels())

def matching_to_tree(m):
    """
    INPUT:

    - m, a PerfectMatching on {1,...,2n}.

    OUTPUT:

    a decreasingly labelled, unordered full binary tree with n+1 leaves.

    EXAMPLES::

        sage: m = PerfectMatching([(1,4),(2,9),(3,10),(5,7),(6,8),(11,12)])
        sage: ascii_art(matching_to_tree(m))
    """
    # the children of the smallest label are the largest remaining
    # element and its partner
    trees = [LabelledRootedTree([LabelledRootedTree([], label=i),
                                 LabelledRootedTree([], label=j)]) for i, j in m]
    max_label = m.size()//2+1 # last labelled node

    while len(trees) > 1:
        max_label += 1
        # find tree with smallest child and both children smaller than max_label
        A = sorted([T for T in trees if max(T[0].label(), T[1].label()) < max_label],
                   key = lambda T: min(T[0].label(), T[1].label()))[0]
        trees.remove(A)
        # give it's root node the new label
        A = LabelledRootedTree(A, label=max_label)
        # find tree with child having label max_label
        B = next(T for T in trees
             if T[0].label() == max_label or T[1].label() == max_label)
        trees.remove(B)
        # replace B with [B[0], A] or [B[1], A]
        if B[0].label() == max_label:
            C = LabelledRootedTree([A, B[1]])
        else:
            C = LabelledRootedTree([A, B[0]])

        trees.append(C)

    return trees[0]

def statistic(m):
    T = matching_to_tree(m)
    if 1 in T[0].leaf_labels():
        return len(T[0].leaf_labels())
    return len(T[1].leaf_labels())

def matching_to_tree(m):
    """
    INPUT:

    - m, a PerfectMatching on {1,...,2n}.

    OUTPUT:

    a decreasingly labelled, unordered full binary tree with n+1 leaves.

    EXAMPLES::

        sage: m = PerfectMatching([(1,4),(2,9),(3,10),(5,7),(6,8),(11,12)])
        sage: ascii_art(matching_to_tree(m))
    """
    # the children of the smallest label are the largest remaining
    # element and its partner
    trees = [LabelledRootedTree([LabelledRootedTree([], label=i),
                                 LabelledRootedTree([], label=j)]) for i, j in m]
    max_label = m.size()//2+1 # last labelled node

    while len(trees) > 1:
        max_label += 1
        # find tree with smallest child and both children smaller than max_label
        A = sorted([T for T in trees if max(T[0].label(), T[1].label()) < max_label],
                   key = lambda T: min(T[0].label(), T[1].label()))[0]
        trees.remove(A)
        # give it's root node the new label
        A = LabelledRootedTree(A, label=max_label)
        # find tree with child having label max_label
        B = next(T for T in trees
             if T[0].label() == max_label or T[1].label() == max_label).next()
        trees.remove(B)
        # replace B with [B[0], A] or [B[1], A]
        if B[0].label() == max_label:
            C = LabelledRootedTree([A, B[1]])
        else:
            C = LabelledRootedTree([A, B[0]])

        trees.append(C)

    return trees[0]

def statistic(m):
    T = matching_to_tree(m)
    if 1 in T[0].leaf_labels():
        return len(T[1].leaf_labels())
    return len(T[0].leaf_labels())

def matching_to_tree(m):
    """
    INPUT:

    - m, a PerfectMatching on {1,...,2n}.

    OUTPUT:

    a decreasingly labelled, unordered full binary tree with n+1 leaves.

    EXAMPLES::

        sage: m = PerfectMatching([(1,4),(2,9),(3,10),(5,7),(6,8),(11,12)])
        sage: ascii_art(matching_to_tree(m))
    """
    # the children of the smallest label are the largest remaining
    # element and its partner
    trees = [LabelledRootedTree([LabelledRootedTree([], label=i),
                                 LabelledRootedTree([], label=j)]) for i, j in m]
    max_label = m.size()//2+1 # last labelled node

    while len(trees) > 1:
        max_label += 1
        # find tree with smallest child and both children smaller than max_label
        A = sorted([T for T in trees if max(T[0].label(), T[1].label()) < max_label],
                   key = lambda T: min(T[0].label(), T[1].label()))[0]
        trees.remove(A)
        # give it's root node the new label
        A = LabelledRootedTree(A, label=max_label)
        # find tree with child having label max_label
        B = next(T for T in trees
             if T[0].label() == max_label or T[1].label() == max_label)
        trees.remove(B)
        # replace B with [B[0], A] or [B[1], A]
        if B[0].label() == max_label:
            C = LabelledRootedTree([A, B[1]])
        else:
            C = LabelledRootedTree([A, B[0]])

        trees.append(C)

    return trees[0]

def statistic(m):
    T = matching_to_tree(m)
    if 1 in T[0].leaf_labels():
        return len(T[1].leaf_labels())
    return len(T[0].leaf_labels())

def matching_to_tree(m):
    """
    INPUT:

    - m, a PerfectMatching on {1,...,2n}.

    OUTPUT:

    a decreasingly labelled, unordered full binary tree with n+1 leaves.

    EXAMPLES::

        sage: m = PerfectMatching([(1,4),(2,9),(3,10),(5,7),(6,8),(11,12)])
        sage: ascii_art(matching_to_tree(m))
    """
    # the children of the smallest label are the largest remaining
    # element and its partner
    trees = [LabelledRootedTree([LabelledRootedTree([], label=i),
                                 LabelledRootedTree([], label=j)]) for i, j in m]
    max_label = m.size()//2+1 # last labelled node

    while len(trees) > 1:
        max_label += 1
        # find tree with smallest child and both children smaller than max_label
        A = sorted([T for T in trees if max(T[0].label(), T[1].label()) < max_label],
                   key = lambda T: min(T[0].label(), T[1].label()))[0]
        trees.remove(A)
        # give it's root node the new label
        A = LabelledRootedTree(A, label=max_label)
        # find tree with child having label max_label
        B = next(T for T in trees
             if T[0].label() == max_label or T[1].label() == max_label).next()
        trees.remove(B)
        # replace B with [B[0], A] or [B[1], A]
        if B[0].label() == max_label:
            C = LabelledRootedTree([A, B[1]])
        else:
            C = LabelledRootedTree([A, B[0]])

        trees.append(C)

    return trees[0]

def statistic(m):
    def aux(T):
        if len(T) == 0:
            return 0
        return 1 + min(aux(T[0]), aux(T[1]))
    return aux(matching_to_tree(m))

def matching_to_tree(m):
    """
    INPUT:

    - m, a PerfectMatching on {1,...,2n}.

    OUTPUT:

    a decreasingly labelled, unordered full binary tree with n+1 leaves.

    EXAMPLES::

        sage: m = PerfectMatching([(1,4),(2,9),(3,10),(5,7),(6,8),(11,12)])
        sage: ascii_art(matching_to_tree(m))
    """
    # the children of the smallest label are the largest remaining
    # element and its partner
    trees = [LabelledRootedTree([LabelledRootedTree([], label=i),
                                 LabelledRootedTree([], label=j)]) for i, j in m]
    max_label = m.size()//2+1 # last labelled node

    while len(trees) > 1:
        max_label += 1
        # find tree with smallest child and both children smaller than max_label
        A = sorted([T for T in trees if max(T[0].label(), T[1].label()) < max_label],
                   key = lambda T: min(T[0].label(), T[1].label()))[0]
        trees.remove(A)
        # give it's root node the new label
        A = LabelledRootedTree(A, label=max_label)
        # find tree with child having label max_label
        B = next(T for T in trees
             if T[0].label() == max_label or T[1].label() == max_label)
        trees.remove(B)
        # replace B with [B[0], A] or [B[1], A]
        if B[0].label() == max_label:
            C = LabelledRootedTree([A, B[1]])
        else:
            C = LabelledRootedTree([A, B[0]])

        trees.append(C)

    return trees[0]

def statistic(m):
    def aux(T):
        if len(T) == 0:
            return 0
        return 1 + min(aux(T[0]), aux(T[1]))
    return aux(matching_to_tree(m))

def matching_to_tree(m):
    """
    INPUT:

    - m, a PerfectMatching on {1,...,2n}.

    OUTPUT:

    a decreasingly labelled, unordered full binary tree with n+1 leaves.

    EXAMPLES::

        sage: m = PerfectMatching([(1,4),(2,9),(3,10),(5,7),(6,8),(11,12)])
        sage: ascii_art(matching_to_tree(m))
    """
    # the children of the smallest label are the largest remaining
    # element and its partner
    trees = [LabelledRootedTree([LabelledRootedTree([], label=i),
                                 LabelledRootedTree([], label=j)]) for i, j in m]
    max_label = m.size()//2+1 # last labelled node

    while len(trees) > 1:
        max_label += 1
        # find tree with smallest child and both children smaller than max_label
        A = sorted([T for T in trees if max(T[0].label(), T[1].label()) < max_label],
                   key = lambda T: min(T[0].label(), T[1].label()))[0]
        trees.remove(A)
        # give it's root node the new label
        A = LabelledRootedTree(A, label=max_label)
        # find tree with child having label max_label
        B = next(T for T in trees
             if T[0].label() == max_label or T[1].label() == max_label).next()
        trees.remove(B)
        # replace B with [B[0], A] or [B[1], A]
        if B[0].label() == max_label:
            C = LabelledRootedTree([A, B[1]])
        else:
            C = LabelledRootedTree([A, B[0]])

        trees.append(C)

    return trees[0]

gap('LoadPackage("QPA");')

import tempfile as _tf, os as _os
_gap_code = r"""
DeclareOperation("NthRadical",[IsList]);
InstallMethod(NthRadical, "for a representation of a quiver", [IsList], 0, function(LIST)
    local M, N, f, h, i, n;
    M := LIST[1];
    n := LIST[2];
    if n = 0 then
        return(IdentityMapping(M));
    else
        f := RadicalOfModuleInclusion(M);
        N := Source(f);
        for i in [1..n-1] do
            h := RadicalOfModuleInclusion(N);
            N := Source(h);
            f := \*(h,f);
        od;
        return(f);
    fi;
end);

DeclareOperation("ARQuiverNak",[IsList]);
InstallMethod(ARQuiverNak, "for a representation of a quiver", [IsList], 0, function(LIST)
    local A, UU, i, injA, j;
    A := LIST[1];
    injA := IndecInjectiveModules(A);
    UU := [];
    for i in injA do
        for j in [0..Dimension(i)-1] do
            Append(UU,[Source(NthRadical([i,j]))]);
        od;
    od;
    return(UU);
end);

DeclareOperation("numbersimplespdatleast3",[IsList]);

InstallMethod(numbersimplespdatleast3, "for a representation of a quiver", [IsList],0,function(LIST)
    local A, UU, simA;
    A := LIST[1];
    simA := SimpleModules(A);
    UU := Filtered(simA,x->ProjDimensionOfModule(x,30)>=3);
    return(Size(UU));
end);

DeclareOperation("IsNtorsionfree",[IsList]);

InstallMethod(IsNtorsionfree, "for a representation of a quiver", [IsList],0,function(LIST)
    local A, CoRegA, M, i, n, temm23;
    A := LIST[1];
    M := LIST[2];

    n := LIST[3];

    CoRegA := DirectSumOfQPAModules(IndecInjectiveModules(A));
    temm23 := [];
    for i in [0..n-1] do Append(temm23,[Size(ExtOverAlgebra(NthSyzygy(CoRegA,i),DTr(M))[2])]);
    od;
    return(Sum(temm23));
end);

DeclareOperation("topreflexive",[IsList]);

InstallMethod(topreflexive, "for a representation of a quiver", [IsList],0,function(LIST)
    local A, LL, LL2, UU, UU2, x;
    A := LIST[1];
    LL := ARQuiverNak([A]);
    LL2 := Filtered(LL,x->IsProjectiveModule(x)=false and IsNtorsionfree([A,x,2])=0);
    UU := [];
    for x in LL2 do Append(UU,[DimensionVector(TopOfModule(x))]);
    od;
    UU2 := Set(UU);
    return(Size(UU2));
end);
"""
with _tf.NamedTemporaryFile(mode="w", suffix=".g", delete=False, dir="/tmp") as _f:
    _f.write('LoadPackage("QPA");;\n')
    _f.write(_gap_code)
    _tmp = _f.name
gap.eval('Read("' + _tmp + '");')
_os.unlink(_tmp)

def kupisch(D):
    DR = D.reverse()
    H = DR.heights()
    return [1 + H[i] for i, s in enumerate(DR) if s == 0] + [1]

def statistic(D):
    K = kupisch(D)
    A = gap.NakayamaAlgebra(gap.GF(3), K)
    return ZZ(gap.topreflexive([A]))

N_vertices = 13 # the maximal number of vertices we consider
N_Ramsey = 7 # all graphs with Ramsey number at most N_Ramsey
statistic_dict = dict()
def statistic(G):
    return statistic_dict.get(G.canonical_label().copy(immutable=True))

"""
The Ramsey number of a graph.

This is the smallest integer $n$ such that every two-colouring of the
$n$ vertices of the complete graph $K_n$ contains a (not necessarily
induced) monochromatic copy of the given graph. [1]

Thus, the Ramsey number of the complete graph $K_n$ is the ordinary Ramsey number $R(n,n)$.  Very few of these numbers are known. [2,3]

[1] Chvatál,Rödl,Szemerédi,Trotter, The Ramsey number of a graph with bounded maximum degree. doi:10.1016/0095-8956(83)90037-0

[2] Radziszowski, Stanisław P. "Small Ramsey numbers." Electron. J. Combin 1.7 (1994).

[3] [[wikipedia:Ramsey's theorem#Ramsey numbers]]

[4] Hendry, G. R. T. Ramsey numbers for graphs with five vertices

"""
def check_Ramsey(n, already_found=[]):
    r"""
    Colour the complete graph on n vertices with two colours, and
    return the subgraphs which appear in all colourings.

    EXAMPLES::

        sage: L = dict()
        sage: n = 2; L[n] = check_Ramsey(n)
        sage: n = 3; L[n] = check_Ramsey(n, sum(L.values(), []))
        sage: n = 4; L[n] = check_Ramsey(n, sum(L.values(), []))
        sage: n = 5; L[n] = check_Ramsey(n, sum(L.values(), []))
        sage: n = 6; L[n] = check_Ramsey(n, sum(L.values(), []))
        sage: n = 7; L[n] = check_Ramsey(n, sum(L.values(), []))

        sage: n=6; graphics_array([H.plot() for H in L[n] if H.is_connected()])
    """
    def has_copy(H1, H2, H):
        r"""
        Return True if there is a copy of H in H1=s or H2=E\s.
        """
        S1 = H1.subgraph_search(H)
        if not S1 is None:
            return True
        S2 = H2.subgraph_search(H)
        if not S2 is None:
            return True
        return False

    candidates = [H.canonical_label() for k in range(n+1) for H in graphs(k)]
    candidates = [H for H in candidates if H not in already_found]
    G = graphs.CompleteGraph(n)
    E = Set(G.edges(labels=False))
    V = G.vertices()
    for size in range(1+len(E)//2):
        print("check subgraphs of size", size)
        tested = []
        for s in E.subsets(size):
            H1 = Graph(n)
            H1.add_edges(s)
            H1 = H1.canonical_label()
            H2 = Graph(n)
            H2.add_edges(E.difference(s))
            H2 = H2.canonical_label()
            if (H1, H2) not in tested:
                tested += [(H1, H2)]
                candidates = [H for H in candidates if has_copy(H1, H2, H)]

    return candidates

def add_to_dict(G, v):
    H = G.canonical_label().copy(immutable=True)
    if H in statistic_dict:
        assert statistic_dict[H] == v, "The graph %s should have Ramsey number %s, got %s instead"%(repr(H), statistic_dict[H], v)
        print("%s already known to have Ramsey value %s"%(repr(H), v))
    else:
        statistic_dict[H] = v

Ramsey_small = dict()
for n in range(N_Ramsey+1):
    Ramsey_small[n] = check_Ramsey(n, sum(Ramsey_small.values(), []))

for v, lG in Ramsey_small.items():
    for G in lG:
        add_to_dict(G, v)

# table IIIa
Ramsey_almost_complete = dict()
Ramsey_almost_complete[3] = 3
Ramsey_almost_complete[4] = 10
Ramsey_almost_complete[5] = 22

for n, v in Ramsey_almost_complete.items():
    G = graphs.CompleteGraph(n)
    G.delete_edge(G.edges()[0])
    add_to_dict(G, v)

# G8 - G20 are from Table 1 in
# Hendry, G. R. T. Ramsey numbers for graphs with five vertices

G8 = Graph([(1,2),(2,3),(3,4),(4,5),(3,5)]).canonical_label().copy(immutable=True)
G9 = Graph([(1,2),(2,3),(3,4),(4,5),(2,4)]).canonical_label().copy(immutable=True)
G10 = Graph([(1,2),(1,3),(1,4),(1,5),(2,3)]).canonical_label().copy(immutable=True)
G12 = Graph([(1,2),(2,3),(3,4),(4,5),(1,3),(3,5)]).canonical_label().copy(immutable=True)
G13 = Graph([(1,2),(1,3),(1,4),(1,5),(2,3),(3,4)]).canonical_label().copy(immutable=True)
G14 = Graph([(1,2),(1,3),(1,4),(2,5),(2,3),(3,4)]).canonical_label().copy(immutable=True)
G15 = Graph([(1,2),(2,3),(3,4),(4,5),(5,1),(1,3)]).canonical_label().copy(immutable=True)
G16 = Graph([(1,2),(2,3),(3,4),(4,5),(5,1),(1,3),(1,4)]).canonical_label().copy(immutable=True)
G17 = Graph([(1,2),(2,3),(3,4),(4,5),(5,1),(1,3),(2,4)]).canonical_label().copy(immutable=True)
G19 = Graph([(1,2),(1,3),(1,4),(2,3),(2,4),(3,4),(1,5)]).canonical_label().copy(immutable=True)
G20 = Graph([(1,2),(1,3),(1,4),(2,3),(2,4),(3,4),(1,5),(2,5)]).canonical_label().copy(immutable=True)

Ramsey_Hendry = [(G8, 9), (G9, 9), (G10, 9), (G12, 9), (G13, 10), (G14, 10), (G15, 9), (G16, 10), (G17, 10), (G19, 18), (G20, 18)]

for G, v in Ramsey_Hendry:
    add_to_dict(G, v)

# table IVa and IVb
Ramsey_bipartite = dict()
Ramsey_bipartite[3, 3] = 18
Ramsey_bipartite[2, 2] = 6
Ramsey_bipartite[2, 3] = 10
Ramsey_bipartite[2, 4] = 14
Ramsey_bipartite[2, 5] = 18
Ramsey_bipartite[2, 6] = 21
Ramsey_bipartite[2, 7] = 26
Ramsey_bipartite[2, 8] = 30
Ramsey_bipartite[2, 9] = 33
Ramsey_bipartite[2, 10] = 38
Ramsey_bipartite[2, 11] = 42
Ramsey_bipartite[2, 12] = 46
Ramsey_bipartite[2, 13] = 50
Ramsey_bipartite[2, 14] = 54
Ramsey_bipartite[2, 15] = 57
Ramsey_bipartite[2, 16] = 62

for (n,m), v in Ramsey_bipartite.items():
    add_to_dict(graphs.CompleteBipartiteGraph(n,m), v)

# section 3.3.2 a
def Ramsey_star(n):
    m = n-1
    assert m > 0
    if is_even(m):
        return 2*m-1
    else:
        return 2*m

for n in range(2, N_vertices+1):
    add_to_dict(graphs.StarGraph(n-1), Ramsey_star(n))

# section 3.3.2, n
def Ramsey_m_4_star(n):
    assert n % 4 == 0, "The number of vertices, %s, should be divisible by 4"%n
    m = n//4
    assert m >= 2, "There must be at least 2 copies of the star, but there are only %s"%m
    return 5*m-1

for n in range(8, N_vertices+1, 4):
    add_to_dict(sum([graphs.StarGraph(3) for k in range(n//4)], Graph()), Ramsey_m_4_star(n))

# section 4.1, a,b,c
def Ramsey_cycle(n):
    assert n > 2
    if is_odd(n) and n > 3:
        return 2*n-1
    elif is_even(n) and n > 4:
        return n-1+(n//2)
    elif n == 3 or n == 4:
        return 6

for n in range(3, N_vertices+1):
    add_to_dict(graphs.CycleGraph(n), Ramsey_cycle(n))

# section 4.1, e
def Ramsey_m_triangles(n):
    assert n % 3 == 0, "The number of vertices, %s, should be divisible by 3"%n
    m = n // 3
    assert m >= 2, "There must be at least 2 copies of the triangle, but there are only %s"%m
    return 5*m

for n in range(6, N_vertices+1, 3):
    add_to_dict(sum([graphs.CycleGraph(3) for k in range(n//3)], Graph()), Ramsey_m_triangles(n))

# section 4.1, f
def Ramsey_m_squares(n):
    assert n % 4 == 0, "The number of vertices, %s, should be divisible by 4"%n
    m = n // 4
    assert m >= 2, "There must be at least 2 copies of the square, but there are only %s"%m
    return 6*m-1

for n in range(8, N_vertices+1, 4):
    add_to_dict(sum([graphs.CycleGraph(4) for k in range(n//4)], Graph()), Ramsey_m_squares(n))

# section 5.1
def Ramsey_paths(m):
    assert m >= 2
    return m + (m//2) - 1

for n in range(2, N_vertices+1):
    add_to_dict(graphs.PathGraph(n), Ramsey_paths(n))

# section 5.13 b
def Ramsey_union_K2(m):
    assert is_even(m) and m >= 2
    return 3*(m//2) - 1

for n in range(2, N_vertices+1, 2):
    add_to_dict(sum([graphs.CompleteGraph(2) for k in range(n//2)], Graph()), Ramsey_union_K2(n))

Ramsey_wheel = dict()
# table VIII
Ramsey_wheel[3] = 6
Ramsey_wheel[4] = 18
Ramsey_wheel[5] = 15
Ramsey_wheel[6] = 17
Ramsey_wheel[7] = 19

for n, v in Ramsey_wheel.items():
    add_to_dict(graphs.WheelGraph(n), v)

Ramsey_book = dict()
# table IX
Ramsey_book[1] = 6
Ramsey_book[2] = 10
Ramsey_book[3] = 14
Ramsey_book[4] = 18
Ramsey_book[5] = 21
Ramsey_book[6] = 26

for n, v in Ramsey_book.items():
    add_to_dict(Graph(1).join(graphs.StarGraph(n)), v)

# R.J. Faudree and R.H. Schelp, Ramsey Numbers for All Linear Forests, Discrete Mathematics, 16 (1976) 149-155.

def Ramsey_linear_forest(n, j):
    assert is_even(n-j)
    return n + (n-j)//2 - 1

for n in range(2, N_vertices+1):
    for la in Partitions(n):
        if min(la) > 1:
            G = sum([graphs.PathGraph(p) for p in la], Graph())
            add_to_dict(G, Ramsey_linear_forest(n, sum(1 for p in la if is_odd(p))))

# J.W. Grossman, The Ramsey Numbers of the Union of Two Stars, Utilitas Mathematica, 16 (1979) 271-279.

def Ramsey_two_stars(n, m):
    assert n >= m, "Ramsey_two_stars only valid for n >= m, but n = %s and m = %s"%(n,m)
    return max(n+2*m, 2*n+1, n+m+3)

for n in range(2, N_vertices+1):
    for m in range(2,(n-1)//2+1):
        G = graphs.StarGraph(m-1) + graphs.StarGraph(n-m-1)
        add_to_dict(G, Ramsey_two_stars(n-m-1, m-1))

# Yu, P., & Li, Y. (2016). All Ramsey numbers for brooms in graphs. The Electronic Journal of Combinatorics, 23(3), 3-29.
def Ramsey_broom(k, l):
    assert k >= 2
    n = k+l
    if l == 1 or l == 2:
        return Ramsey_star(n)
    if l == 3:
        return Ramsey_star(n-1)
    if l >= 2*k - 1:
        return n + (l+1)//2 - 1
    if 4 <= l <= 2*k - 2:
        return 2*n - 2*((l+1)//2) - 1

def BroomGraph(k, l):
    G = graphs.StarGraph(k)
    assert G.degree(0) == k
    G.add_path([0] + [k+1+i for i in range(l-1)])
    G.layout("tree", save_pos=True)
    return G

for n in range(2, N_vertices+1):
    for k in range(2,n-1):
        G = BroomGraph(k, n-k)
        add_to_dict(G, Ramsey_broom(k, n-k))

# Jerrold W. Grossman, Frank Harary, Maria Klawe, Generalized Ramsey theory for graphs, X: double stars
        
def Ramsey_double_star(k, l):
    assert k >= 2
    n = k+l+2
    if is_odd(k) and l <= 2:
        return max(2*k+1, k+2*l+2)
    if (is_even(k) or l >= 3) and (k^2 <= 2*l or k >= 3*l):
        return max(2*k+2, k+2*l+2)
        
def DoubleStarGraph(k, l):
    assert k >= l, "DoubleStarGraph only defined for k >= l"
    G = graphs.StarGraph(k)
    H = G.disjoint_union(graphs.StarGraph(l), "pairs")
    H.add_edge((0,0),(1,0))
    return H

for n in range(2, N_vertices+1):
    for l in range(2,n//2):
        k = n-l-2
        G = DoubleStarGraph(k, l)
        assert G.num_verts() == n
        if (is_odd(k) and l <= 2) or ((is_even(k) or l >= 3) and (k^2 <= 2*l or k >= 3*l)):
            add_to_dict(G, Ramsey_double_star(k, l))


######################################################################
# finally, add isolated vertices
for G, v in list(statistic_dict.items()):
    for k in range(N_vertices-G.num_verts()):
        add_to_dict(G.disjoint_union(Graph(k)), max(G.num_verts()+k, v))

N_vertices = 13 # the maximal number of vertices we consider
N_Ramsey = 7 # all graphs with Ramsey number at most N_Ramsey
statistic_dict = dict()
def statistic(G):
    return statistic_dict.get(G.canonical_label().copy(immutable=True))

"""
The Ramsey number of a graph.

This is the smallest integer n such that every two-colouring of the
n vertices of the complete graph K_n contains a (not necessarily
induced) monochromatic copy of the given graph. [1]

Thus, the Ramsey number of the complete graph K_n is the ordinary Ramsey number R(n,n).  Very few of these numbers are known. [2,3]

[1] Chvatál,Rödl,Szemerédi,Trotter, The Ramsey number of a graph with bounded maximum degree. doi:10.1016/0095-8956(83)90037-0

[2] Radziszowski, Stanisław P. "Small Ramsey numbers." Electron. J. Combin 1.7 (1994).

[3] [[wikipedia:Ramsey's theorem#Ramsey numbers]]

[4] Hendry, G. R. T. Ramsey numbers for graphs with five vertices

"""
def check_Ramsey(n, already_found=[]):
    r"""
    Colour the complete graph on n vertices with two colours, and
    return the subgraphs which appear in all colourings.

    EXAMPLES::

        sage: L = dict()
        sage: n = 2; L[n] = check_Ramsey(n)
        sage: n = 3; L[n] = check_Ramsey(n, sum(L.values(), []))
        sage: n = 4; L[n] = check_Ramsey(n, sum(L.values(), []))
        sage: n = 5; L[n] = check_Ramsey(n, sum(L.values(), []))
        sage: n = 6; L[n] = check_Ramsey(n, sum(L.values(), []))
        sage: n = 7; L[n] = check_Ramsey(n, sum(L.values(), []))

        sage: n=6; graphics_array([H.plot() for H in L[n] if H.is_connected()])
    """
    def has_copy(H1, H2, H):
        r"""
        Return True if there is a copy of H in H1=s or H2=E\s.
        """
        S1 = H1.subgraph_search(H)
        if not S1 is None:
            return True
        S2 = H2.subgraph_search(H)
        if not S2 is None:
            return True
        return False

    candidates = [H.canonical_label() for k in range(n+1) for H in graphs(k)]
    candidates = [H for H in candidates if H not in already_found]
    G = graphs.CompleteGraph(n)
    E = Set(G.edges(labels=False))
    V = G.vertices()
    for size in range(1+len(E)//2):
        print("check subgraphs of size", size)
        tested = []
        for s in E.subsets(size):
            H1 = Graph(n)
            H1.add_edges(s)
            H1 = H1.canonical_label()
            H2 = Graph(n)
            H2.add_edges(E.difference(s))
            H2 = H2.canonical_label()
            if (H1, H2) not in tested:
                tested += [(H1, H2)]
                candidates = [H for H in candidates if has_copy(H1, H2, H)]

    return candidates

def add_to_dict(G, v):
    H = G.canonical_label().copy(immutable=True)
    if H in statistic_dict:
        assert statistic_dict[H] == v, "The graph %s should have Ramsey number %s, got %s instead"%(repr(H), statistic_dict[H], v)
        print("%s already known to have Ramsey value %s"%(repr(H), v))
    else:
        statistic_dict[H] = v

Ramsey_small = dict()
for n in range(N_Ramsey+1):
    Ramsey_small[n] = check_Ramsey(n, sum(Ramsey_small.values(), []))

for v, lG in Ramsey_small.items():
    for G in lG:
        add_to_dict(G, v)

# table IIIa
Ramsey_almost_complete = dict()
Ramsey_almost_complete[3] = 3
Ramsey_almost_complete[4] = 10
Ramsey_almost_complete[5] = 22

for n, v in Ramsey_almost_complete.items():
    G = graphs.CompleteGraph(n)
    G.delete_edge(G.edges()[0])
    add_to_dict(G, v)

# G8 - G20 are from Table 1 in
# Hendry, G. R. T. Ramsey numbers for graphs with five vertices

G8 = Graph([(1,2),(2,3),(3,4),(4,5),(3,5)]).canonical_label().copy(immutable=True)
G9 = Graph([(1,2),(2,3),(3,4),(4,5),(2,4)]).canonical_label().copy(immutable=True)
G10 = Graph([(1,2),(1,3),(1,4),(1,5),(2,3)]).canonical_label().copy(immutable=True)
G12 = Graph([(1,2),(2,3),(3,4),(4,5),(1,3),(3,5)]).canonical_label().copy(immutable=True)
G13 = Graph([(1,2),(1,3),(1,4),(1,5),(2,3),(3,4)]).canonical_label().copy(immutable=True)
G14 = Graph([(1,2),(1,3),(1,4),(2,5),(2,3),(3,4)]).canonical_label().copy(immutable=True)
G15 = Graph([(1,2),(2,3),(3,4),(4,5),(5,1),(1,3)]).canonical_label().copy(immutable=True)
G16 = Graph([(1,2),(2,3),(3,4),(4,5),(5,1),(1,3),(1,4)]).canonical_label().copy(immutable=True)
G17 = Graph([(1,2),(2,3),(3,4),(4,5),(5,1),(1,3),(2,4)]).canonical_label().copy(immutable=True)
G19 = Graph([(1,2),(1,3),(1,4),(2,3),(2,4),(3,4),(1,5)]).canonical_label().copy(immutable=True)
G20 = Graph([(1,2),(1,3),(1,4),(2,3),(2,4),(3,4),(1,5),(2,5)]).canonical_label().copy(immutable=True)

Ramsey_Hendry = [(G8, 9), (G9, 9), (G10, 9), (G12, 9), (G13, 10), (G14, 10), (G15, 9), (G16, 10), (G17, 10), (G19, 18), (G20, 18)]

for G, v in Ramsey_Hendry:
    add_to_dict(G, v)

# table IVa and IVb
Ramsey_bipartite = dict()
Ramsey_bipartite[3, 3] = 18
Ramsey_bipartite[2, 2] = 6
Ramsey_bipartite[2, 3] = 10
Ramsey_bipartite[2, 4] = 14
Ramsey_bipartite[2, 5] = 18
Ramsey_bipartite[2, 6] = 21
Ramsey_bipartite[2, 7] = 26
Ramsey_bipartite[2, 8] = 30
Ramsey_bipartite[2, 9] = 33
Ramsey_bipartite[2, 10] = 38
Ramsey_bipartite[2, 11] = 42
Ramsey_bipartite[2, 12] = 46
Ramsey_bipartite[2, 13] = 50
Ramsey_bipartite[2, 14] = 54
Ramsey_bipartite[2, 15] = 57
Ramsey_bipartite[2, 16] = 62

for (n,m), v in Ramsey_bipartite.items():
    add_to_dict(graphs.CompleteBipartiteGraph(n,m), v)

# section 3.3.2 a
def Ramsey_star(n):
    m = n-1
    assert m > 0
    if is_even(m):
        return 2*m-1
    else:
        return 2*m

for n in range(2, N_vertices+1):
    add_to_dict(graphs.StarGraph(n-1), Ramsey_star(n))

# section 3.3.2, n
def Ramsey_m_4_star(n):
    assert n % 4 == 0, "The number of vertices, %s, should be divisible by 4"%n
    m = n//4
    assert m >= 2, "There must be at least 2 copies of the star, but there are only %s"%m
    return 5*m-1

for n in range(8, N_vertices+1, 4):
    add_to_dict(sum([graphs.StarGraph(3) for k in range(n//4)], Graph()), Ramsey_m_4_star(n))

# section 4.1, a,b,c
def Ramsey_cycle(n):
    assert n > 2
    if is_odd(n) and n > 3:
        return 2*n-1
    elif is_even(n) and n > 4:
        return n-1+(n//2)
    elif n == 3 or n == 4:
        return 6

for n in range(3, N_vertices+1):
    add_to_dict(graphs.CycleGraph(n), Ramsey_cycle(n))

# section 4.1, e
def Ramsey_m_triangles(n):
    assert n % 3 == 0, "The number of vertices, %s, should be divisible by 3"%n
    m = n // 3
    assert m >= 2, "There must be at least 2 copies of the triangle, but there are only %s"%m
    return 5*m

for n in range(6, N_vertices+1, 3):
    add_to_dict(sum([graphs.CycleGraph(3) for k in range(n//3)], Graph()), Ramsey_m_triangles(n))

# section 4.1, f
def Ramsey_m_squares(n):
    assert n % 4 == 0, "The number of vertices, %s, should be divisible by 4"%n
    m = n // 4
    assert m >= 2, "There must be at least 2 copies of the square, but there are only %s"%m
    return 6*m-1

for n in range(8, N_vertices+1, 4):
    add_to_dict(sum([graphs.CycleGraph(4) for k in range(n//4)], Graph()), Ramsey_m_squares(n))

# section 5.1
def Ramsey_paths(m):
    assert m >= 2
    return m + (m//2) - 1

for n in range(2, N_vertices+1):
    add_to_dict(graphs.PathGraph(n), Ramsey_paths(n))

# section 5.13 b
def Ramsey_union_K2(m):
    assert is_even(m) and m >= 2
    return 3*(m//2) - 1

for n in range(2, N_vertices+1, 2):
    add_to_dict(sum([graphs.CompleteGraph(2) for k in range(n//2)], Graph()), Ramsey_union_K2(n))

Ramsey_wheel = dict()
# table VIII
Ramsey_wheel[3] = 6
Ramsey_wheel[4] = 18
Ramsey_wheel[5] = 15
Ramsey_wheel[6] = 17
Ramsey_wheel[7] = 19

for n, v in Ramsey_wheel.items():
    add_to_dict(graphs.WheelGraph(n), v)

Ramsey_book = dict()
# table IX
Ramsey_book[1] = 6
Ramsey_book[2] = 10
Ramsey_book[3] = 14
Ramsey_book[4] = 18
Ramsey_book[5] = 21
Ramsey_book[6] = 26

for n, v in Ramsey_book.items():
    add_to_dict(Graph(1).join(graphs.StarGraph(n)), v)

# R.J. Faudree and R.H. Schelp, Ramsey Numbers for All Linear Forests, Discrete Mathematics, 16 (1976) 149-155.

def Ramsey_linear_forest(n, j):
    assert is_even(n-j)
    return n + (n-j)//2 - 1

for n in range(2, N_vertices+1):
    for la in Partitions(n):
        if min(la) > 1:
            G = sum([graphs.PathGraph(p) for p in la], Graph())
            add_to_dict(G, Ramsey_linear_forest(n, sum(1 for p in la if is_odd(p))))

# J.W. Grossman, The Ramsey Numbers of the Union of Two Stars, Utilitas Mathematica, 16 (1979) 271-279.

def Ramsey_two_stars(n, m):
    assert n >= m, "Ramsey_two_stars only valid for n >= m, but n = %s and m = %s"%(n,m)
    return max(n+2*m, 2*n+1, n+m+3)

for n in range(2, N_vertices+1):
    for m in range(2,(n-1)//2+1):
        G = graphs.StarGraph(m-1) + graphs.StarGraph(n-m-1)
        add_to_dict(G, Ramsey_two_stars(n-m-1, m-1))

# Yu, P., & Li, Y. (2016). All Ramsey numbers for brooms in graphs. The Electronic Journal of Combinatorics, 23(3), 3-29.
def Ramsey_broom(k, l):
    assert k >= 2
    n = k+l
    if l == 1 or l == 2:
        return Ramsey_star(n)
    if l == 3:
        return Ramsey_star(n-1)
    if l >= 2*k - 1:
        return n + (l+1)//2 - 1
    if 4 <= l <= 2*k - 2:
        return 2*n - 2*((l+1)//2) - 1

def BroomGraph(k, l):
    G = graphs.StarGraph(k)
    assert G.degree(0) == k
    G.add_path([0] + [k+1+i for i in range(l-1)])
    G.layout("tree", save_pos=True)
    return G

for n in range(2, N_vertices+1):
    for k in range(2,n-1):
        G = BroomGraph(k, n-k)
        add_to_dict(G, Ramsey_broom(k, n-k))

# Jerrold W. Grossman, Frank Harary, Maria Klawe, Generalized Ramsey theory for graphs, X: double stars
        
def Ramsey_double_star(k, l):
    assert k >= 2
    n = k+l+2
    if is_odd(k) and l <= 2:
        return max(2*k+1, k+2*l+2)
    if (is_even(k) or l >= 3) and (k^2 <= 2*l or k >= 3*l):
        return max(2*k+2, k+2*l+2)
        
def DoubleStarGraph(k, l):
    assert k >= l, "DoubleStarGraph only defined for k >= l"
    G = graphs.StarGraph(k)
    H = G.disjoint_union(graphs.StarGraph(l), "pairs")
    H.add_edge((0,0),(1,0))
    return H

for n in range(2, N_vertices+1):
    for l in range(2,n//2):
        k = n-l-2
        G = DoubleStarGraph(k, l)
        assert G.num_verts() == n
        if (is_odd(k) and l <= 2) or ((is_even(k) or l >= 3) and (k^2 <= 2*l or k >= 3*l)):
            add_to_dict(G, Ramsey_double_star(k, l))


######################################################################
# finally, add isolated vertices
for G, v in list(statistic_dict.items()):
    for k in range(N_vertices-G.num_verts()):
        add_to_dict(G.disjoint_union(Graph(k)), max(G.num_verts()+k, v))

def statistic(pi):
    i = 1
    pos = 0
    count = 0
    while i < len(pi):
        if pi[pos] == i:
            i += 1
        elif pi[pos] > i:
            count += 1
        pos = (pos+1)%len(pi)
    return count
def statistic(pi):
    n = len(pi)
    pinv = pi.inverse()
    return sum( n-i for i in [1 .. n-1] if pinv(i) > pinv(i+1) )

def statistic(pi):
    i = 1
    pos = 0
    count = 0
    while i < len(pi):
        if pi[pos] == i:
            i += 1
        elif pi[pos] > i:
            count += 1
        pos = (pos+1)%len(pi)
    return count
def statistic(pi):
    n = len(pi)
    pinv = pi.inverse()
    return sum( n-i for i in [1 .. n-1] if pinv(i) > pinv(i+1) )

def statistic(x):
    xw = 1
    lvl = 1
    for step in x:
        if step == 0:
            lvl += 1
        else:
            xw *= lvl
    return xw

def dyck_box_weightstatistic(x):
    xw = 1
    lvl = 1
    for step in x:
        if step == 0:
            lvl += 1
        else:
            xw *= lvl
    return xw

def statistic(p):
    return sum( i+1 for i in range(len(p)-1) if p[i] == p[i+1] == 0 )

def count_prime_leaves(m):
    from sage.graphs.graph_decompositions.modular_decomposition import NodeType
    if m.label() == NodeType.PRIME:
        return len(m.leaf_labels())
    return sum(count_prime_leaves(c) for c in m)

def statistic(G):
    m = G.modular_decomposition(algorithm='habib', style='tree')
    return count_prime_leaves(m)

def count_prime_leaves(m):
    from sage.graphs.graph_decompositions.modular_decomposition import NodeType
def count_prime_leaves(m):
    if m.label() == NodeType.PRIME:
        return len(m.leaf_labels())
    return sum(count_prime_leaves(c) for c in m)

def statistic(G):
    m = G.modular_decomposition(algorithm='habib', style='tree')
    return count_prime_leaves(m)

def statistic(g):
    return sum(h.szeged_index() for h in g.connected_components_subgraphs())

def statistic(g):
    return gsum(h.szeged_index() for h in g.connected_components_subgraphs())

def statistic(C):
    from sage.combinat.q_analogues import q_int
    q = polygen(ZZ, 'q')
    return max(prod(q_int(d, q) for d in WeylGroup(C).degrees()).list())

def statistic(C):
    from sage.combinat.q_analogues import q_int
    q = polygen(ZZ, 'q')
    return max(prod(q_int(d, q) for d in WeylGroup(C).degrees()).list())

def statistic(D):
    D = DyckWord(D)
    n = len(D)    
    tunnels = D.tunnels()
    t = [(i, j) for (i, j) in tunnels if i + j == n]
    res = 0
    for (a1, b1) in t:
        for (a2, b2) in t:
            if a1 < a2 < b2 < b1:
                res += 1
    return res