def statistic(G):
    return sum(cop_number_connected(S) for S in G.connected_components_subgraphs())

def cop_number_connected(S):
    n = S.num_verts()
    if n == 1:
        return 1
    else:
        for s in range(1, n):
            if cop_number_s(S, s):
                return s
        return n

def power_vertex_to_subset_aux(N):
    """
    INPUT:

    - an integer
    - a pair, where the first is an acceptable input and the second is an integer
    """
    if isinstance(N, (int, Integer)):
        return [N]

    return power_vertex_to_subset_aux(N[0]) + [N[1]]

def power_vertex_to_subset(N):
    return set(power_vertex_to_subset_aux(N))

def cop_number_s(G, s):
    """
    G must be connected and have at least 2 vertices.
    """
    assert G.is_connected(), "G must be connected"
    P = G
    for i in range(s-1):
        P = P.strong_product(G)

    V = set(G.vertices())

    psi = {T: V.difference(*[G.neighbors(v)+[v] for v in power_vertex_to_subset(T)]) for T in P}
    changed = True
    while changed:
        changed = False
        for (T1, T2) in P.edges(labels=False):
            newpsiT1 = psi[T1].intersection(set(sum([G.neighbors(v)+[v] for v in psi[T2]], [])))
            if len(newpsiT1) == 0:
                return True
            newpsiT1 = set(newpsiT1)
            if newpsiT1 != psi[T1]:
                changed = True
            newpsiT2 = psi[T2].intersection(set(sum([G.neighbors(v)+[v] for v in psi[T1]], [])))
            if len(newpsiT2) == 0:
                return True
            newpsiT2 = set(newpsiT2)
            if newpsiT2 != psi[T2]:
                changed = True
            psi[T1] = newpsiT1
            psi[T2] = newpsiT2
    return False