import sympy as sp
import numpy as np
import functools as ft
import pickle
import multiprocessing as mp
mp.set_start_method('fork') # 'spawn' might be better on some platforms

from sage.graphs.graph_generators import graphs

import sys

# helper functions before main routine
def coeff_func_for_part(the_graph,box_size,partition,ext_ext_edges=[],ext_int_edges=[],int_int_edges=[],num_int_vert=0):
    """Constructs a coefficient function for a given partition
    
    Constructs a function that returns the coefficient of tuples in the `box_size`-th box space
    of the spin planar algebra, specified by the partition.
    
    Parameters
    ----------
    box_size : int
    partition : list of sets of integers 0 .. box_size-1. You can omit singleton sets, so that
    [{1,2}] and [{0},{1,2}] would have the same effect.
    ext_ext_edges : list of pairs (a, b)
    """

    def coeff_func(*variables):
        if len(variables) == 1 and variables[0] == "diag":
            return diag_for_4_box_elem(partition,ext_ext_edges,ext_int_edges,int_int_edges,num_int_vert)
        if len(variables) != box_size:
            raise Exception("number of variables does not match the box size.")
        for part in partition:
            vars_in_part = [variables[i] for i in range(box_size) if i in part]
            if len(set(vars_in_part)) > 1:
                return 0
        if any(True for edge in ext_ext_edges if not variables[edge[0]] in the_graph.neighbors(variables[edge[1]])):
            return 0
        if num_int_vert == 0:
            return 1
        int_vert_candidates = []
        for i in range(num_int_vert):
            ext_vert_neighs = [set(the_graph.neighbors(variables[edge[0]])) for edge in ext_int_edges if edge[1] == i]
            cand_set = ft.reduce(lambda x, y: x.intersection(y), ext_vert_neighs, set(range(the_graph.order())))
            int_vert_candidates.append(cand_set)
        def branch_candidates_for(prev_verts):
            i = len(prev_verts)
            prev_int_verts_conn_by_edges_0 = {edge[0] for edge in int_int_edges if edge[1] == i and edge[0] < i}
            prev_int_verts_conn_by_edges_1 = {edge[1] for edge in int_int_edges if edge[0] == i and edge[1] < i}
            prev_int_verts_conn_by_edges = prev_int_verts_conn_by_edges_0.union(prev_int_verts_conn_by_edges_1)
            nvc = [v for v in int_vert_candidates[i] if not any(True for j in prev_int_verts_conn_by_edges if not v in the_graph.neighbors(prev_verts[j]))]
            if len(prev_verts) == num_int_vert - 1:
                return [prev_verts + [v] for v in nvc]
            else:
                return [x for v in nvc for x in branch_candidates_for(prev_verts + [v])]
        valid_int_conf = branch_candidates_for([])
        return len(valid_int_conf)
    return coeff_func
    
def paths_for_4_box_elem(partition,ext_ext_edges=[],ext_int_edges=[],int_int_edges=[],num_int_vert=0):
    output =""
    if num_int_vert == 1:
        output += r'\node (inn_v0) at (0,0) [inn-vert] {};' + '\n'
    elif num_int_vert == 2:
        output += r'\node (inn_v0) at (-0.3,0) [inn-vert] {};' + '\n'
        output += r'\node (inn_v1) at (0.3,0) [inn-vert] {};' + '\n'
    elif num_int_vert == 3:
        output += r'\node (inn_v0) at (-0.4,0.3) [inn-vert] {};' + '\n'
        output += r'\node (inn_v1) at (0.4,0.3) [inn-vert] {};' + '\n'
        output += r'\node (inn_v2) at (0,-0.4) [inn-vert] {};' + '\n'
    elif num_int_vert == 4:
        output += r'\node (inn_v0) at (0,0.5) [inn-vert] {};' + '\n'
        output += r'\node (inn_v1) at (0.5,0) [inn-vert] {};' + '\n'
        output += r'\node (inn_v2) at (0,-0.5) [inn-vert] {};' + '\n'
        output += r'\node (inn_v3) at (-0.5,0) [inn-vert] {};' + '\n'
    if any({0,1}.issubset(s) for s in partition):
        output += r'\draw [thick,double] (corner0) -- (corner1);' + '\n'
    if any({1,2}.issubset(s) for s in partition):
        output += r'\draw [thick,double] (corner1) -- (corner2);' + '\n'
    if any({2,3}.issubset(s) for s in partition):
        output += r'\draw [thick,double] (corner2) -- (corner3);' + '\n'
    if any({3,0}.issubset(s) for s in partition):
        output += r'\draw [thick,double] (corner3) -- (corner0);' + '\n'
    if not any({0,1}.issubset(s) for s in partition) and not any({3,0}.issubset(s) for s in partition) and any({0,2}.issubset(s) for s in partition):
        output += r'\draw [thick,double] (corner0) -- (corner2);' + '\n'    
    if not any({0,1}.issubset(s) for s in partition) and not any({1,2}.issubset(s) for s in partition) and any({1,3}.issubset(s) for s in partition):
        output += r'\draw [thick,double] (corner1) -- (corner3);' + '\n'    
    for ee_edge in ext_ext_edges:
        output += r'\draw [] (corner%d) -- (corner%d);' % ee_edge + '\n'
    for ei_edge in ext_int_edges:
        output += r'\draw [] (corner%d) -- (inn_v%d);' % ei_edge + '\n'
    for ii_edge in int_int_edges:
        output += r'\draw [] (inn_v%d) -- (inn_v%d);' % ii_edge + '\n'
    return output

def diag_for_4_box_elem(partition,ext_ext_edges=[],ext_int_edges=[],int_int_edges=[],num_int_vert=0):
    output = r'\begin{tikzpicture}' + '\n'
    output += r'\draw[bgrect] (-1.3,-1.3) rectangle (1.3,1.3);' + '\n'
    output += r'\node (corner0) at (-1,1) [ext-vert] {};' + '\n'
    output += r'\node (corner1) at (1,1) [ext-vert] {};' + '\n'
    output += r'\node (corner2) at (1,-1) [ext-vert] {};' + '\n'
    output += r'\node (corner3) at (-1,-1) [ext-vert] {};' + '\n'
    output += paths_for_4_box_elem(partition,ext_ext_edges,ext_int_edges,int_int_edges,num_int_vert)
    output += r'\end{tikzpicture}' + '\n'
    return output
 
def rot_ind(x,step,length):
    y = x + step
    if y < length:
        return y
    return y-length

def rot_partition(partition,step,length):
    """Rotate partitions
    
    Parameters
    ==========
    partition : collection of sets of indexes for external leaves, singletons can be omitted
    step : int, for step of rotations
    length : int, period of rotation
    """
    return [set([rot_ind(x,step,length) for x in block]) for block in partition]

def rot_e_e_edges(edges,step,length):
    """Rotated edges between external leaves
    
    Parameters
    ==========
    edges : collection of pairs of indexes for external leaves
    step : int, for step of rotations
    length : int, period of rotation
    """
    return [(rot_ind(edge[0],step,length),rot_ind(edge[1],step,length)) for edge in edges]

def rot_e_i_edges(edges,step,length):
    """Rotated edges between external leaves and internal vertices
    
    Parameters
    ==========
    edges : collection of pairs (index for external leaves, index for interval vertices)
    step : int, for step of rotations
    length : int, period of rotation
    """
    return [(rot_ind(edge[0],step,length),edge[1]) for edge in edges]

def coeff_funcs_parts_e_e_rotated_3(the_graph,partition,ee_edges):
    """Rotated coefficient functions based on edges between external leaves
    
    Parameters
    ==========
    partition : collection of sets of indexes for external leaves, singletons can be omitted
    ee_edges : collection of pairs of indexes for external leaves
    """
    return [coeff_func_for_part(the_graph,3,rot_partition(partition,i,3),rot_e_e_edges(ee_edges,i,3)) for i in range(3)]
def coeff_funcs_parts_e_e_e_i_rotated_3(the_graph,partition,ee_edges,ei_edges,num_int_e):
    """Rotated coefficient functions based on edges between external leaves,
    and those between external leaves and internal vertices
    
    Parameters
    ==========
    partition : collection of sets of indexes for external leaves, singletons can be omitted
    ee_edges : collection of pairs of indexes for external leaves
    ei_edges : collection of pairs (index for external leaves, index for interval vertices)
    num_int_e : int, number of internal vertices
    """
    return [coeff_func_for_part(the_graph,3,rot_partition(partition,i,3),rot_e_e_edges(ee_edges,i,3),rot_e_i_edges(ei_edges,i,3),[],num_int_e) for i in range(3)]

def coeff_funcs_parts_e_e_rotated_4(the_graph,partition,ee_edges):
    return [coeff_func_for_part(the_graph,4,rot_partition(partition,i,4),rot_e_e_edges(ee_edges,i,4)) for i in range(4)]
def coeff_funcs_parts_e_e_e_i_rotated_4(the_graph,partition,ee_edges,ei_edges,num_int_e):
    return [coeff_func_for_part(the_graph,4,rot_partition(partition,i,4),rot_e_e_edges(ee_edges,i,4),rot_e_i_edges(ei_edges,i,4),[],num_int_e) for i in range(4)]
def coeff_funcs_parts_e_e_e_i_i_i_rotated_4(the_graph,partition,ee_edges,ei_edges,ii_edges,num_int_e):
    return [coeff_func_for_part(the_graph,4,rot_partition(partition,i,4),rot_e_e_edges(ee_edges,i,4),rot_e_i_edges(ei_edges,i,4),ii_edges,num_int_e) for i in range(4)]

## main routine
# the_graph = graphs.OrthogonalPolarGraph(6,2,'-')
def compute(the_graph):
    print(f"We are working with {the_graph}.")
    num_vert = the_graph.order()
    # find orbits of vertices
    vert_orbs = the_graph.automorphism_group(return_group=False, orbits=True)
    if len(vert_orbs) > 1:
        print("The graph is not vertex transitive.")
    sgl_orb_rep = [(orb[0],) for orb in vert_orbs]
    # double orbits
    dbl_vert_orb_rep = []
    for sgl_orb in sgl_orb_rep:
        singleton = [[sgl_orb[0]]]
        scd_coord_orbits = the_graph.automorphism_group(partition=singleton+[list(set(the_graph.vertices()).difference(set(sgl_orb)))],orbits=True,return_group=False)
        for orb in scd_coord_orbits:
            dbl_orb = sgl_orb + (orb[0],)
            dbl_vert_orb_rep += [dbl_orb]
    print(f"There are {len(dbl_vert_orb_rep)} orbits in the product of vertices.")
    # triple orbits
    trip_vert_orb_rep = []
    for dbl_orb in dbl_vert_orb_rep:
        singletons = [[x] for x in set(dbl_orb)]
        trip_coord_orbits = the_graph.automorphism_group(partition=singletons+[list(set(the_graph.vertices()).difference(set(dbl_orb)))],orbits=True,return_group=False)
        for orb in trip_coord_orbits:
            trip_orb = dbl_orb + (orb[0],)
            trip_vert_orb_rep += [trip_orb]
    print(f"There are {len(trip_vert_orb_rep)} orbits in the triple product of vertices.")
    # quadruple orbits
    quad_vert_orb_rep = []
    quad_vert_orb_fourth_cand = {}
    for trip_orb in trip_vert_orb_rep:
        singletons = [[x] for x in set(trip_orb)]
        fourth_coord_orbits = the_graph.automorphism_group(partition=singletons+[list(set(the_graph.vertices()).difference(set(trip_orb)))],orbits=True,return_group=False)
        for orb in fourth_coord_orbits:
            quad_orb = trip_orb + (orb[0],)
            quad_vert_orb_rep += [quad_orb]
            quad_vert_orb_fourth_cand[quad_orb] = set(orb)
    print(f"There are {len(quad_vert_orb_rep)} orbits in the quadruple product of vertices.")
    
    # TL diagrams
    nc_partitions_4 = [[]] + [rot_partition([{0,1}],i,4) for i in range(4)] + [[{1,3}],[{0,2}]] + [rot_partition([{0,1,2}],i,4) for i in range(4)] + [[{0,1,2,3}],[{0,1},{2,3}],[{0,3},{1,2}]]
    func_list_4_box_sp = [coeff_func_for_part(the_graph,4,this_nc_part) for this_nc_part in nc_partitions_4]
    # elements with one internal edge
    func_list_4_box_sp += coeff_funcs_parts_e_e_rotated_4(the_graph,[],[(0,1)])
    
    func_list_4_box_sp += [coeff_func_for_part(the_graph,4,[],[(0,2)]),coeff_func_for_part(the_graph,4,[],[(1,3)])]
    
    func_list_4_box_sp += coeff_funcs_parts_e_e_rotated_4(the_graph,[{0,1}],[(0,3)])
    
    func_list_4_box_sp += coeff_funcs_parts_e_e_rotated_4(the_graph,[{0,1}],[(0,2)])
    
    func_list_4_box_sp += coeff_funcs_parts_e_e_rotated_4(the_graph,[{0,1}],[(2,3)])
    
    func_list_4_box_sp += coeff_funcs_parts_e_e_rotated_4(the_graph,[{1,3}],[(0,1)])
    
    func_list_4_box_sp += coeff_funcs_parts_e_e_rotated_4(the_graph,[{1,2,3}],[(0,1)])
    
    func_list_4_box_sp += [coeff_func_for_part(the_graph,4,[{0,1},{2,3}],[(1,2)]),coeff_func_for_part(the_graph,4,[{0,3},{1,2}],[(0,1)])]
    
    # two internal edges
    func_list_4_box_sp += coeff_funcs_parts_e_e_rotated_4(the_graph,[],[(0,1),(1,2)])
    
    func_list_4_box_sp += coeff_funcs_parts_e_e_rotated_4(the_graph,[],[(0,1),(0,2)])
    
    func_list_4_box_sp += coeff_funcs_parts_e_e_rotated_4(the_graph,[],[(0,1),(1,3)])
    
    func_list_4_box_sp += [coeff_func_for_part(the_graph,4,[],[(0,1),(2,3)]),coeff_func_for_part(the_graph,4,[],[(0,3),(1,2)])]
    
    func_list_4_box_sp += coeff_funcs_parts_e_e_rotated_4(the_graph,[{0,1}],[(0,3),(0,2)])
    
    func_list_4_box_sp += coeff_funcs_parts_e_e_rotated_4(the_graph,[{0,1}],[(0,3),(2,3)])
    
    func_list_4_box_sp += coeff_funcs_parts_e_e_rotated_4(the_graph,[{0,1}],[(0,2),(2,3)])
    
    func_list_4_box_sp += [coeff_func_for_part(the_graph,4,[{0,2}],[(1,2),(2,3)]),coeff_func_for_part(the_graph,4,[{1,3}],[(0,1),(1,2)])]
    
    # three internal edges
    func_list_4_box_sp += coeff_funcs_parts_e_e_rotated_4(the_graph,[],[(0,1),(0,3),(1,3)])
    
    func_list_4_box_sp += coeff_funcs_parts_e_e_rotated_4(the_graph,[],[(0,3),(2,3),(1,2)])
    
    func_list_4_box_sp += coeff_funcs_parts_e_e_rotated_4(the_graph,[],[(0,1),(0,2),(0,3)])
    
    func_list_4_box_sp += coeff_funcs_parts_e_e_rotated_4(the_graph,[{0,1}],[(0,3),(0,2),(2,3)])
    
    func_list_4_box_sp += [coeff_func_for_part(the_graph,4,[],[(0,2),(0,3),(1,2)]),coeff_func_for_part(the_graph,4,[],[(0,1),(1,3),(2,3)])]
    
    func_list_4_box_sp += [coeff_func_for_part(the_graph,4,[],[(1,3),(0,3),(1,2)]),coeff_func_for_part(the_graph,4,[],[(0,1),(0,2),(2,3)])]
    
    # four internal edges
    func_list_4_box_sp += coeff_funcs_parts_e_e_rotated_4(the_graph,[],[(0,3),(2,3),(1,2),(1,3)])
    
    func_list_4_box_sp += coeff_funcs_parts_e_e_rotated_4(the_graph,[],[(0,3),(2,3),(1,2),(0,2)])
    
    func_list_4_box_sp += [coeff_func_for_part(the_graph,4,[],[(0,1),(1,2),(2,3),(0,3)]),coeff_func_for_part(the_graph,4,[],[],[(0,0),(1,0),(2,0),(3,0)],[],1)]
    
    # five internal edges
    func_list_4_box_sp += [coeff_func_for_part(the_graph,4,[],[(0,1),(1,2),(2,3),(0,3),(1,3)]),coeff_func_for_part(the_graph,4,[],[(0,1),(1,2),(2,3),(0,3),(0,2)])]
    
    func_list_4_box_sp += coeff_funcs_parts_e_e_e_i_rotated_4(the_graph,[],[(0,1)],[(0,0),(1,0),(2,0),(3,0)],1)
    
    # six internal edges
    func_list_4_box_sp += coeff_funcs_parts_e_e_e_i_rotated_4(the_graph,[],[(0,1),(1,2)],[(0,0),(1,0),(2,0),(3,0)],1)
    
    func_list_4_box_sp += [coeff_func_for_part(the_graph,4,[],[(0,1),(2,3)],[(0,0),(1,0),(2,0),(3,0)],[],1),coeff_func_for_part(the_graph,4,[],[(0,3),(1,2)],[(0,0),(1,0),(2,0),(3,0)],[],1)]
    
    # seven internal edges
    func_list_4_box_sp += coeff_funcs_parts_e_e_e_i_rotated_4(the_graph,[],[(0,1),(1,2),(2,3)],[(0,0),(1,0),(2,0),(3,0)],1)
    
    func_list_4_box_sp += [coeff_func_for_part(the_graph,4,[],[],[(0,0),(0,1),(1,0),(2,0),(2,1),(3,1)],[(0,1)],2),coeff_func_for_part(the_graph,4,[],[],[(0,0),(1,0),(1,1),(2,1),(3,0),(3,1)],[(0,1)],2)]
    
    # eight internal edges
    func_list_4_box_sp += [coeff_func_for_part(the_graph,4,[],[(0,1),(1,2),(2,3),(3,0)],[(0,0),(1,0),(2,0),(3,0)],[],1)]
    
    func_list_4_box_sp += coeff_funcs_parts_e_e_e_i_i_i_rotated_4(the_graph,[],[(0,1)],[(0,0),(0,1),(1,0),(2,0),(2,1),(3,1)],[(0,1)],2)
    
    func_list_4_box_sp += coeff_funcs_parts_e_e_e_i_i_i_rotated_4(the_graph,[],[(0,1)],[(0,0),(1,0),(1,1),(2,1),(3,0),(3,1)],[(0,1)],2)
    
    # nine internal edges
    func_list_4_box_sp += coeff_funcs_parts_e_e_e_i_i_i_rotated_4(the_graph,[],[],[(0,0),(1,1),(2,1),(2,2),(3,0),(3,2)],[(0,1),(1,2),(0,2)],3)
    
    func_list_4_box_sp += coeff_funcs_parts_e_e_e_i_i_i_rotated_4(the_graph,[],[(0,1),(1,2)],[(0,0),(0,1),(1,0),(2,0),(2,1),(3,1)],[(0,1)],2)
    
    func_list_4_box_sp += coeff_funcs_parts_e_e_e_i_i_i_rotated_4(the_graph,[],[(0,1),(1,2)],[(0,0),(1,0),(1,1),(2,1),(3,0),(3,1)],[(0,1)],2)
    
    func_list_4_box_sp += [coeff_func_for_part(the_graph,4,[],[(0,1),(2,3)],[(0,0),(0,1),(1,0),(2,0),(2,1),(3,1)],[(0,1)],2),coeff_func_for_part(the_graph,4,[],[(0,1),(2,3)],[(0,0),(1,0),(1,1),(2,1),(3,0),(3,1)],[(0,1)],2)]
    
    func_list_4_box_sp += [coeff_func_for_part(the_graph,4,[],[(1,2),(0,3)],[(0,0),(0,1),(1,0),(2,0),(2,1),(3,1)],[(0,1)],2),coeff_func_for_part(the_graph,4,[],[(1,2),(0,3)],[(0,0),(1,0),(1,1),(2,1),(3,0),(3,1)],[(0,1)],2)]
    
    # ten internal edges
    func_list_4_box_sp += coeff_funcs_parts_e_e_e_i_i_i_rotated_4(the_graph,[],[(0,1),(1,2),(2,3)],[(0,0),(0,1),(1,0),(2,0),(2,1),(3,1)],[(0,1)],2)
    
    func_list_4_box_sp += coeff_funcs_parts_e_e_e_i_i_i_rotated_4(the_graph,[],[(0,1),(1,2),(2,3)],[(0,0),(1,0),(1,1),(2,1),(3,0),(3,1)],[(0,1)],2)
    
    func_list_4_box_sp += [coeff_func_for_part(the_graph,4,[],[],[(0,0),(0,1),(0,2),(1,0),(2,0),(2,1),(2,2),(3,2)],[(0,1),(1,2)],3),coeff_func_for_part(the_graph,4,[],[],[(0,0),(1,0),(1,1),(1,2),(2,2),(3,0),(3,1),(3,2)],[(0,1),(1,2)],3)]
    
    func_list_4_box_sp += coeff_funcs_parts_e_e_e_i_i_i_rotated_4(the_graph,[],[],[(0,0),(1,0),(1,1),(1,2),(2,2),(3,0),(3,2)],[(0,1),(1,2),(0,2)],3)
    
    func_list_4_box_sp += coeff_funcs_parts_e_e_e_i_i_i_rotated_4(the_graph,[],[],[(0,0),(1,0),(1,1),(2,1),(2,2),(3,0),(3,2)],[(0,1),(1,2),(0,2)],3)
    
    func_list_4_box_sp += coeff_funcs_parts_e_e_e_i_i_i_rotated_4(the_graph,[],[],[(0,0),(0,1),(1,1),(2,1),(2,2),(3,0),(3,2)],[(0,1),(1,2),(0,2)],3)
    
    # eleven internal edges
    func_list_4_box_sp  += [coeff_func_for_part(the_graph,4,[],[(0,1),(1,2),(2,3),(3,0)],[(0,0),(0,1),(1,0),(2,0),(2,1),(3,1)],[(0,1)],2),coeff_func_for_part(the_graph,4,[],[(0,1),(1,2),(2,3),(3,0)],[(0,0),(1,0),(1,1),(2,1),(3,0),(3,1)],[(0,1)],2)]
    
    func_list_4_box_sp += coeff_funcs_parts_e_e_e_i_i_i_rotated_4(the_graph,[],[(0,1)],[(0,0),(0,1),(0,2),(1,0),(2,0),(2,1),(2,2),(3,2)],[(0,1),(1,2)],3)
    
    func_list_4_box_sp += coeff_funcs_parts_e_e_e_i_i_i_rotated_4(the_graph,[],[(0,1)],[(0,0),(1,0),(1,1),(1,2),(2,2),(3,0),(3,1),(3,2)],[(0,1),(1,2)],3)
    
    func_list_4_box_sp += coeff_funcs_parts_e_e_e_i_i_i_rotated_4(the_graph,[],[(0,1)],[(0,0),(1,0),(1,1),(1,2),(2,2),(3,0),(3,2)],[(0,1),(1,2),(0,2)],3)
    
    func_list_4_box_sp += coeff_funcs_parts_e_e_e_i_i_i_rotated_4(the_graph,[],[(0,3)],[(0,0),(1,0),(1,1),(1,2),(2,2),(3,0),(3,2)],[(0,1),(1,2),(0,2)],3)
    
    func_list_4_box_sp += coeff_funcs_parts_e_e_e_i_i_i_rotated_4(the_graph,[],[(1,2)],[(0,0),(1,0),(1,1),(1,2),(2,2),(3,0),(3,2)],[(0,1),(1,2),(0,2)],3)
    
    func_list_4_box_sp += coeff_funcs_parts_e_e_e_i_i_i_rotated_4(the_graph,[],[(2,3)],[(0,0),(1,0),(1,1),(1,2),(2,2),(3,0),(3,2)],[(0,1),(1,2),(0,2)],3)
    
    # twelve internal edges
    func_list_4_box_sp += coeff_funcs_parts_e_e_e_i_i_i_rotated_4(the_graph,[],[(0,1)],[(0,0),(1,0),(1,1),(1,2),(2,2),(3,0),(3,2)],[(0,1),(1,2),(0,2)],3)
    
    func_list_4_box_sp += [coeff_func_for_part(the_graph,4,[],[],[(0,0),(0,3),(1,0),(1,1),(2,1),(2,2),(3,2),(3,3)],[(0,1),(1,2),(2,3),(3,0)],4)]
    
    func_list_4_box_sp += [coeff_func_for_part(the_graph,4,[],[],[(0,0),(0,1),(1,0),(1,2),(2,2),(2,3),(3,1),(3,3)],[(0,1),(0,2),(1,3),(2,3)],4)]
    
    # thirteen internal edges
    func_list_4_box_sp += coeff_funcs_parts_e_e_e_i_i_i_rotated_4(the_graph,[],[(0,1)],[(0,0),(0,3),(1,0),(1,1),(2,1),(2,2),(3,2),(3,3)],[(0,1),(1,2),(2,3),(3,0)],4)
    
    func_list_4_box_sp += [coeff_func_for_part(the_graph,4,[],[],[(0,1),(1,0),(1,2),(2,2),(2,3),(3,1),(3,3)],[(0,1),(0,2),(1,2),(1,3),(2,3)],4)]
    
    func_list_4_box_sp += [coeff_func_for_part(the_graph,4,[],[],[(0,1),(1,0),(1,2),(2,2),(2,3),(3,1),(3,3)],[(0,1),(0,2),(1,2),(1,3)],4)]
    
    func_list_4_box_sp += [coeff_func_for_part(the_graph,4,[],[],[(0,0),(0,1),(1,0),(1,2),(2,2),(2,3),(3,1),(3,3)],[(0,2),(1,2),(1,3),(2,3)],4)]
    
    # more
#    func_list_4_box_sp += [coeff_func_for_part(the_graph,4,[],[],[(0,0),(0,3),(1,0),(1,1),(2,7),(2,8),(3,8),(3,9)],[(0,1),(1,2),(2,3),(3,0),(1,4),(2,4),(2,6),(3,6),(4,7),(4,5),(6,5),(6,9),(5,7),(5,9),(7,8),(8,9)],10)]
#    func_list_4_box_sp += [coeff_func_for_part(the_graph,4,[],[],[(0,0),(0,3),(1,0),(1,1),(2,7),(2,8),(2,4),(3,8),(3,9)],[(0,1),(1,2),(2,3),(3,0),(1,4),(2,4),(2,6),(3,6),(4,7),(4,5),(6,5),(6,9),(5,7),(5,9),(7,8),(8,9)],10)]
#    func_list_4_box_sp += [coeff_func_for_part(the_graph,4,[],[],[(0,0),(0,3),(1,6),(1,7),(2,7),(2,8),(3,2),(3,3)],[(0,1),(1,2),(2,3),(3,0),(1,4),(1,5),(0,4),(2,5),(6,4),(8,5),(9,4),(9,5),(6,7),(7,8),(8,9),(9,6)],10)]
#    func_list_4_box_sp += [coeff_func_for_part(the_graph,4,[],[(0,1)],[(0,0),(0,3),(1,0),(1,1),(2,7),(2,8),(3,8),(3,9)],[(0,1),(1,2),(2,3),(3,0),(1,4),(2,4),(2,6),(3,6),(4,7),(4,5),(6,5),(6,9),(5,7),(5,9),(7,8),(8,9)],10)]

    def myfunc3(x):
        return [func(x[0],x[1],x[2],x[3]) for func in func_list_4_box_sp]
    
    if __name__ == '__main__':
        with mp.Pool() as p:
            Ap_array = p.map(myfunc3, quad_vert_orb_rep)
    else:
        print(f"The script is not run as main, we resort to single-core computation.")
        Ap_array = list(map(myfunc3, quad_vert_orb_rep))
    
    Ap_rank = np.linalg.matrix_rank(np.array(Ap_array))
    print(f"We now have rank {Ap_rank}, spanned by {len(Ap_array[0])} elements, in the 4-box space.")
    if Ap_rank == len(quad_vert_orb_rep):
        print("We found enough elements in the 4-box space to conclude that the quantum automorphism groups is classical.")
    else:
        print("We cannot conclude that the quantum automorphism group is trivial.")