1. Utilities¶

def density_to_graph(input_matrix, threshold=0.1, validate=True):
    """
    Converts a density matrix state representation to a graph representation via an adjacency matrix.
    It assumes qubit systems.

    :param input_matrix: a density matrix to be converted to a graph
    :type input_matrix: numpy.ndarray
    :param threshold: a minimum threshold value of negativity for the state to assume entangled
    :type threshold: double or float
    :param validate: validate the input density is a graph state if True
    :type validate: bool
    :raise AssertionError: if input_matrix is not a valid graph state when validate is True
    :return: an adjacency matrix representation for a pure state or a list of adjacency matrices for a mixed state
    :rtype: numpy.ndarray or list[(float, numpy.ndarray)]
    """
    if dmf.is_pure(input_matrix):
        graph_adj = _density_to_graph_pure(input_matrix, threshold)
    else:
        # eigen-decomposition
        eigenvalues, eigenvectors = linalg.eigh(input_matrix)
        graph_adj = []
        for i in range(len(eigenvectors)):
            rho_i = eigenvectors[i] @ np.conjugate(eigenvectors[i].T)
            graph_adj.append((eigenvalues, _density_to_graph_pure(rho_i, threshold)))

    if validate:
        assert np.allclose(
            input_matrix, graph_to_density(graph_adj)
        ), "Input matrix is not a graph state."
    return graph_adj

def density_to_stabilizer(input_matrix):
    """
    Convert a density matrix to stabilizer via graph representation

    This works only for graph states, not for general stabilizer states that are not graph states

    :param input_matrix: a density matrix
    :type input_matrix: numpy.ndarray
    :return: a stabilizer representation
    :rtype: list[(float, StabilizerTableau)]
    """
    # first check if the state is pure
    if dmf.is_pure(input_matrix):
        graph = _density_to_graph_pure(input_matrix)
        stabilizer = _graph_to_stabilizer_pure(graph)
        return [(1.0, stabilizer)]
    else:
        # eigen-decomposition
        eigenvalues, eigenvectors = linalg.eigh(input_matrix)
        stabilizer = []
        for i in range(len(eigenvectors)):
            rho_i = eigenvectors[i] @ np.conjugate(eigenvectors[i].T)
            graph = density_to_graph(rho_i)
            tmp_stabilizer = graph_to_stabilizer(graph)
            stabilizer.append((eigenvalues[i], tmp_stabilizer[0][1]))

        return stabilizer

def graph_to_density(input_graph):
    """
    Builds a density matrix representation from a graph (adjacency matrix or networkx.Graph)

    :param input_graph: the graph from which we will build a density matrix
    :type input_graph: networkx.Graph or numpy.ndarray
    :raise TypeError: if input_graph is not of the type of networkx.Graph or a graphiq.backends.graph.state.Graph
        or an adjacency matrix given by a numpy array
    :return: a density matrix as a numpy array
    :rtype: numpy.ndarray
    """

    if isinstance(input_graph, list):
        rho = 0
        for p_i, graph_i in input_graph:
            state_i = _graph_to_density_pure(graph_i)
            rho += p_i * state_i
        return rho
    elif isinstance(input_graph, (nx.Graph, np.ndarray)):
        return _graph_to_density_pure(input_graph)

    else:
        raise TypeError(
            "Input graph must be NetworkX graph or adjacency matrix "
            "using numpy.ndarray or mixed state representation for graph. "
        )

def graph_to_stabilizer(input_graph):
    """
    Convert a graph to stabilizer represented by StabilizerTableau

    :param input_graph: the input graph to be converted to the stabilizer
    :type input_graph: networkX.Graph or numpy.ndarray
    :raise TypeError: if input_graph is not of the type of networkx.Graph or
            an adjacency matrix given by a numpy array or a list
    :return: two binary matrices representing the stabilizer generators
    :rtype: list[(float, StabilizerTableau)]
    """
    if isinstance(input_graph, list):
        results = []
        for p_i, graph_i in input_graph:
            tableau_i = _graph_to_stabilizer_pure(graph_i)
            results.append((p_i, tableau_i))
        return results
    else:
        tableau = _graph_to_stabilizer_pure(input_graph)
        return [(1.0, tableau)]

def mixed_graph_equivalency(graph1, graph2):
    """
    Identify if two mixed state graph representations are the same.

    Assuming no repetitive graph in each graph list. Two graphs are equal (not isomorphism).

    :param graph1: the first graph
    :type graph1: list((float, networkX.Graph)) or networkX.Graph
    :param graph2: the second graph
    :type graph2: list((float, networkX.Graph)) or networkX.Graph
    :return: True if graph1 is the same as graph2; False otherwise
    :rtype: bool
    """
    if isinstance(graph1, list) and len(graph1) == 1:
        graph1 = graph1[0][1]
    if isinstance(graph2, list) and len(graph2) == 1:
        graph2 = graph2[0][1]
    if isinstance(graph1, list) and isinstance(graph2, list):
        if len(graph1) == len(graph2):
            graph2_copy = copy.deepcopy(graph2)
            for p_i, g_i in graph1:
                for q_i, t_i in graph2_copy:
                    if np.equal(p_i, q_i) and nx.utils.graphs_equal(g_i, t_i):
                        graph2_copy.remove((q_i, t_i))
                        break
            return len(graph2_copy) == 0
        else:
            return False
    elif isinstance(graph1, nx.Graph) and isinstance(graph2, nx.Graph):
        return nx.utils.graphs_equal(graph1, graph2)
    else:
        return False

def mixed_stabilizer_equivalency(stab1, stab2):
    """
     Identify if two mixed state stabilizer representations are the same.

    :param stab1: the first stabilizer state
    :type stab1: list[(float, StabilizerTableau)] or StabilizerTableau
    :param stab2: the second stabilizer state
    :type stab2: list[(float, StabilizerTableau)] or StabilizerTableau
    :return: if stab1 is the same as stab2
    :rtype: bool
    """
    if isinstance(stab1, list) and len(stab1) == 1:
        stab1 = stab1[0][1]
    if isinstance(stab2, list) and len(stab2) == 1:
        stab2 = stab2[0][1]
    if isinstance(stab1, list) and isinstance(stab2, list):
        if len(stab1) == len(stab2):
            stab2_copy = copy.deepcopy(stab2)
            for p_i, s_i in stab1:
                for q_i, t_i in stab2_copy:
                    if np.equal(p_i, q_i) and s_i == t_i:
                        stab2_copy.remove((q_i, t_i))
                        break
            return len(stab2_copy) == 0
        else:
            return False
    elif isinstance(stab1, StabilizerTableau) and isinstance(stab2, StabilizerTableau):
        return stab1 == stab2
    else:
        return False

def stabilizer_to_density(input_stabilizer):
    """
    Convert a stabilizer to a density matrix

    :param input_stabilizer: the stabilizer representation
    :type input_stabilizer: list[(float, StabilizerTableau)] or StabilizerTableau
    :return: a density matrix
    :rtype: numpy.ndarray
    """
    if isinstance(input_stabilizer, list):
        if isinstance(input_stabilizer[0], str):
            x_matrix, z_matrix = sfu.string_to_symplectic(input_stabilizer)
            stab = StabilizerTableau([x_matrix, z_matrix])
            return _stabilizer_to_density_pure(stab)
        elif isinstance(input_stabilizer[0], tuple):
            rho = 0

            for p_i, tableau_i in input_stabilizer:
                rho += p_i * _stabilizer_to_density_pure(tableau_i)
        else:
            raise ValueError(
                "Invalid stabilizer input: use either StabilizerTableau or a mixed state representation."
            )
    elif isinstance(input_stabilizer, StabilizerTableau):
        return _stabilizer_to_density_pure(input_stabilizer)
    else:
        raise ValueError(
            "Invalid stabilizer input: use either StabilizerTableau or a mixed state representation."
        )

def stabilizer_to_graph(input_stabilizer, validate=True):
    """
    Convert a stabilizer to graph.

    :param input_stabilizer: the stabilizer representation in terms of (probability, StabilizerTableau)
    :type input_stabilizer: list[(float, StabilizerTableau)] or StabilizerTableau
    :param validate: validate if input_stabilizer is a graph state when enabled
    :type validate: bool
    :raises AssertionError: if input stabilizer is not a graph state when validate is True
    :return: a graph representation
    :rtype: [(float, numpy.ndarray)]
    """
    graph_list = []
    if isinstance(input_stabilizer, list):
        for each_stabilizer in input_stabilizer:
            tableau = each_stabilizer[1]
            graph = _graph_finder(tableau.x_matrix, tableau.z_matrix)
            graph_list.append((each_stabilizer[0], graph))
    elif isinstance(input_stabilizer, StabilizerTableau):
        tableau = input_stabilizer
        graph = _graph_finder(tableau.x_matrix, tableau.z_matrix)
        graph_list.append((1.0, graph))
    else:
        raise ValueError(
            "Please input the stabilizer representation in terms of probability, X matrix, Z matrix"
        )

    if validate:
        assert mixed_stabilizer_equivalency(
            input_stabilizer, graph_to_stabilizer(graph_list)
        ), "Input stabilizer is not a graph state."

    return graph_list

def state_to_graph(state):
    """
    A helper function to turn any valid representation into a graph. It also returns the StabilizerTableau corresponding
     to the initial state, and the gate sequence needed to convert the state into a graph-state.

    :param state: the state to be converted to graph
    :type state: StabilizerTableau or CliffordTableau or networkX.Graph or np.ndarray
    :return: tuple (graph, input state's stabilizer tableau, gate_list)
    :rtype: tuple (networkX.Graph, StabilizerTableau, list)
    """
    state = copy.deepcopy(state)
    # returns graph, input state's tableau, gate_list
    if isinstance(state, nx.Graph):
        tab = get_stabilizer_tableau_from_graph(state)
        return state, tab, []
    elif isinstance(state, np.ndarray):
        try:
            graph = nx.from_numpy_array(state)
            tab = get_stabilizer_tableau_from_graph(graph)
            return graph, [], tab
        except:
            raise ValueError(
                "the input numpy array is not a valid adjacency matrix, try fixing it or using other valid input types"
            )

    elif isinstance(state, CliffordTableau):
        z_matrix = state.stabilizer_z
        x_matrix = state.stabilizer_x
        tab = state.to_stabilizer()
    elif isinstance(state, StabilizerTableau):
        z_matrix = state.z_matrix
        x_matrix = state.x_matrix
        tab = state
    else:
        raise ValueError(
            "input data should either be a adjacency matrix, graph, Clifford or Stabilizer tableau"
        )
    graph, (h_pos, p_dag_pos) = _graph_finder(x_matrix, z_matrix, get_ops_data=True)
    gate_list = [("H", pos) for pos in h_pos] + [("P_dag", pos) for pos in p_dag_pos]

    # phase correction; adding Z gates at the end to make the phase of the transformed state equal to an ideal graph
    g_tab = get_stabilizer_tableau_from_graph(graph)
    phase_correction = _phase_correction(tab, g_tab, gate_list)

    gate_list += phase_correction
    return graph, tab, gate_list

def is_lc_equivalent(adj_matrix1, adj_matrix2, mode="deterministic", seed=0):
    r"""
    Determines whether two graph states are local-Clifford equivalent or not, given the adjacency matrices of the two.
    It takes two adjacency matrices as input and returns a numpy.ndarray containing $n$ ($2 \times 2$ arrays
    = clifford operations on each qubit.

    :param adj_matrix1: the adjacency matrix of the first graph. This is equal to the binary matrix for representing
        Pauli Z part of the symplectic binary representation of the stabilizer generators
    :type adj_matrix1: numpy.ndarray
    :param adj_matrix2: the adjacency matrix of the second graph
    :type adj_matrix2: numpy.ndarray
    :param mode: the chosen mode for finding solutions. It can be either 'deterministic' (default) or 'random'.
    :type mode: str
    :param seed: an optional input to set the random seed for the random search approach
    :type seed: int
    :raises AssertionError: if the number of rows in the row reduced matrix is less than the rank of coefficient matrix
        or if the number of linearly dependent columns is not equal to :math:`4n - rank`
        (for :math:`n` being the number of nodes in the graph)
    :return: If a solution is found, returns an array of single-qubit Clifford :math:`2 \\times 2` matrices
        in the symplectic formalism. If not, graphs are not LC equivalent and returns None.
    :rtype: bool, numpy.ndarray or None
    """

    n_nodes = np.shape(adj_matrix1)[0]
    # get the coefficient matrix for the system of linear equations
    coeff_matrix = _coeff_maker(adj_matrix1, adj_matrix2)

    # row reduction applied
    reduced_coeff_matrix, _, last_nonzero_row_index = slinalg.row_reduction(
        coeff_matrix, coeff_matrix * 0
    )

    rank = last_nonzero_row_index + 1
    # check for rank to see how many independent equations are there = rank of the matrix
    if rank >= 4 * n_nodes:
        # Those two graph states are not LC equivalent for sure
        return False, None
    # update the matrix to remove zero rows
    reduced_coeff_matrix = np.array([row for row in reduced_coeff_matrix if row.any()])
    assert (
        np.shape(reduced_coeff_matrix)[0] == rank
    ), "The number of remaining rows is less than the rank!"
    # rank = np.shape(reduced_coeff_matrix)[0]

    # finding linearly dependent columns
    col_list = _col_finder(reduced_coeff_matrix)
    length = len(col_list)
    assert length == 4 * n_nodes - rank, "column list is not correct"

    # if solution basis' length, which is the dimension of the solution basis, is less than 4 then we should check every
    # possible vector in solution basis:
    solution_basis = _solution_basis_finder(reduced_coeff_matrix, col_list)
    if len(solution_basis) < 5:
        basis_dimension = len(solution_basis)
        all_solutions = [*range(2**basis_dimension)]
        all_solutions = [list(format(i, f"0{basis_dimension}b")) for i in all_solutions]
        all_solutions = np.array(all_solutions).astype(int)
        all_solutions = all_solutions.T
        solution_basis = np.transpose(solution_basis, axes=(1, 2, 0)).reshape(
            4 * n_nodes, basis_dimension
        )
        all_solutions = (solution_basis @ all_solutions) % 2
        for solution in all_solutions.T:
            if _is_valid_clifford(solution.reshape(4 * n_nodes, 1)):
                # convert the solution to an array of n * (2 X 2) matrices
                valid_solution = solution.reshape(4 * n_nodes, 1)
                return True, valid_solution.reshape(n_nodes, 2, 2)
        return False, None

    if mode == "random":
        # Use random mode to get the fast convergence for large states
        # Check for random solutions 1000 times
        rand_solution = _random_checker(
            reduced_coeff_matrix, col_list, trial_count=1000, seed=seed
        )
        if isinstance(rand_solution, np.ndarray):
            # random result
            return True, rand_solution.reshape(n_nodes, 2, 2)
        else:
            return False, None

    elif mode == "deterministic":
        basis = _solution_basis_finder(reduced_coeff_matrix, col_list)
        possible_combinations = combinations(basis, 2)
        solution_set = []
        for basis_combination in possible_combinations:
            possible_solution = basis_combination[0] + basis_combination[1]
            solution_set.append(possible_solution % 2)
        for solution in solution_set:
            if _is_valid_clifford(solution):
                # convert the solution to an array of n * (2 X 2) matrices
                return True, solution.reshape(n_nodes, 2, 2)

        # states are NOT LC equivalent
        return False, None
    else:
        raise ValueError(
            'The mode should be either "random" or "deterministic" (default)'
        )

def iso_equal_check(graph1, graph2):
    """
    Checks if the graph1 is local-Clifford equivalent to any graph that is isomorphic to graph2

    :param graph1: initial graph
    :type graph1: networkx.Graph
    :param graph2: target graph
    :type graph2: networkx.Graph
    :return: If equivalent (True,the graph that graph1 is equivalent to) and if not, (False, graph1 itself)
    :rtype: (bool, networkx.Graph)
    """
    iso_graphs_g2 = iso_graph_finder(graph2)
    adj_matrix1 = nx.to_numpy_array(graph1)
    iso_g2_adj_matrices = [nx.to_numpy_array(graph) for graph in iso_graphs_g2]
    for adj_matrix2 in iso_g2_adj_matrices:
        success, solution = is_lc_equivalent(adj_matrix1, adj_matrix2)
        if isinstance(solution, np.ndarray):
            return True, nx.to_networkx_graph(adj_matrix2)

    return False, graph1

def iso_graph_finder(input_graph):
    """
    Generates the list of all graphs that are isomorphic to the input graph G.
    Scales with n! and faces runtime or memory issues for large graphs.

    :param input_graph: input graph
    :type input_graph: networkx.Graph
    :return: the list of graphs that are isomorphic to G
    :rtype: list[networkx.Graph]
    """
    iso_graphs = []
    list_nodes = sorted(input_graph)
    n_nodes = len(list_nodes)
    all_permu = list(permutations(list_nodes, len(list_nodes)))

    for each_permu in all_permu:
        adj_matrix = np.zeros((n_nodes, n_nodes))
        map_dict = dict(zip(list_nodes, each_permu))
        g_copy = nx.relabel_nodes(input_graph, map_dict, copy=True)
        for node_id1 in list_nodes:
            for node_id2 in list(g_copy.neighbors(node_id1)):
                adj_matrix[node_id1, node_id2] = 1

        g_copy = nx.to_networkx_graph(adj_matrix)
        iso_graphs.append(g_copy)

    return iso_graphs

def lc_graph_operations(adj_matrix, solution):
    """
    Finds a list of local complementations needed to be applied on each node of the given graph
    to transform it into the target graph given the Clifford transformation matrix,
    which is the output of the solver function.

    :param adj_matrix: The adjacency matrix of the first graph. This is equal to the binary matrix for
        representing Pauli Z part of the symplectic binary representation of the stabilizer generators
    :type adj_matrix: numpy.ndarray
    :param solution: an array of single-qubit Clifford :math:`2 \\times 2` matrices in the symplectic formalism
    :type solution: numpy.ndarray
    :return: a list of the names of the operations that need to be applied on each qubit in the correct order.
    :rtype: list[str]
    """
    # takes an adjacency matrix and the solution (Clifford operation) and returns the list of local complementations
    # needed for graph transformation.
    r_matrix = _R_matrix(adj_matrix, solution)
    n = np.shape(r_matrix)[0]
    singles_list = []
    doubles_list = []

    while _condition(r_matrix):
        r_matrix, s_list = _singles(r_matrix)
        singles_list.extend(s_list)

    while not np.array_equal(r_matrix, np.eye(n)):
        r_matrix, d_list = _doubles(r_matrix)
        doubles_list.extend(d_list)
    g_list = singles_list
    for i in doubles_list:
        g_list.append(i[0])
        g_list.append(i[1])
        g_list.append(i[0])
    return g_list

def local_clifford_ops(solution):
    """
    Finds a list of operators needed to be applied on each qubit of the first graph to transform in to the second,
    given the Clifford transformation matrix, which is the output of the solver function.

    :param solution: an array of single-qubit Clifford :math:`2 \\times 2` matrices in the symplectic formalism
    :type solution: numpy.ndarray
    :return: a list of the names of the operations that need to be applied on each qubit in the correct order.
    :rtype: list[str]
    """
    # The order of the operations is the same as the qubits' labels in the graphs

    # allowed operations on single qubits in binary symplectic representation
    identity = np.array([[1, 0], [0, 1]])
    hadamard = np.array([[0, 1], [1, 0]])
    p = np.array([[1, 1], [0, 1]])
    ph = np.array([[1, 1], [1, 0]])
    hp_dagger = np.array([[0, 1], [1, 1]])
    php = np.array([[1, 0], [1, 1]])

    ops_list = [identity, hadamard, p, ph, hp_dagger, php]
    ops_list_str = ["I", "H", "P", "P H", "H P_dag", "P H P"]
    ops_dict = zip(list(range(len(ops_list))), ops_list_str)
    ops_dict = dict(ops_dict)
    ops_names = []
    for i in solution:
        for j in range(len(ops_list)):
            if np.array_equal(i, ops_list[j]):
                ops_names.append(ops_dict[j])
    return ops_names

def local_comp_graph(input_graph, node_id):
    """
    Applies a local complementation on the node (specified by node_id) of the input graph input_graph
    and returns the resulting graph.

    :param input_graph: input graph
    :type input_graph: networkx.Graph
    :param node_id: the index of the node to apply tho local complementation on
    :type node_id: int
    :raises AssertionError: if the node index node_id is not in the graph
    :return: a new transformed graph
    :rtype: networkx.Graph
    """
    # TODO: integrate this with the state representation functions
    n_nodes = input_graph.number_of_nodes()
    assert n_nodes > node_id >= 0, "node index is not in the graph"
    adj_matrix = nx.to_numpy_array(input_graph).astype(int)
    identity = np.eye(n_nodes, n_nodes)
    gamma_matrix = np.zeros((n_nodes, n_nodes))

    # gamma has only a single 1 element on position = diag(i)
    gamma_matrix[node_id, node_id] = 1

    new_adj_matrix = (
        adj_matrix
        @ (
            gamma_matrix @ adj_matrix
            + adj_matrix[node_id, node_id] * gamma_matrix
            + identity
        )
        % 2
    ) % 2
    for j in range(n_nodes):
        new_adj_matrix[j, j] = 0
    new_graph = nx.to_networkx_graph(new_adj_matrix)

    return new_graph

def automorph_check(adj1, labels_arr):
    """
    Given an initial adjacency matrix (of a graph), and an array of new permutations of the nodes' labels, this function
    return an array of distinct adjacency matrices resulted from those permutations.
    Note that the order will not remain the same as the order of the input labels.

    :param adj1: initial adjacency matrix
    :type adj1: numpy.ndarray
    :param labels_arr:  an array of new labeling
    :type labels_arr: numpy.ndarray
    :return: An array of adjacency matrices
    :rtype: numpy.ndarray
    """
    # uses set to remove redundancies
    # the set includes the adjacency matrices in flatten form so that they can turn into tuples to be members of a set
    adj_set = {tuple(adj1.astype(int).flatten())}
    adj_list = [adj1]
    n_node = adj1.shape[0]
    for label in labels_arr:
        new_adj = relabel(adj1, label)
        adj_set.add(tuple(new_adj.flatten()))
    # remove the initial adjacency matrix from the set to just keep the new ones, since the original adj1 is already in
    # the final list of adjacencies as the first element
    adj_set.remove(tuple(adj1.astype(int).flatten()))
    for flat_adj in adj_set:
        remade_adj = np.array(flat_adj)
        adj_list.append(remade_adj.reshape(n_node, n_node))
    return np.array(adj_list)

def check_isomorphism(graph, g_list, _only_auto=False):
    """
    check if the provided graph is isomorphic to any graph in the g_list

    :param graph: graph to check
    :type graph: nx.Graph
    :param g_list: graph list to check against
    :type g_list: list
    :param _only_auto: if set to true, instead of isomorphism, only automorphism will be checked
    :type _only_auto: bool
    :return: True of False, if an isomorphism case was detected.
    :rtype: bool
    """
    iso = False
    check = _equal_graphs if _only_auto else nx.is_isomorphic
    for g in g_list:
        if check(graph, g):
            iso = True
            break
    return iso

def emitter_sorted(adj_arr):
    """
    Takes an iterable of adjacency matrices (ideally output of the <iso_finder> function) as input and sort them by number of
    emitters needed to generate those graph states.

    :param adj_arr: An array of adjacency matrices which can be the output of the <iso_finder> function
    :type adj_arr: numpy.ndarray
    :return: a list of tuples (adj matrices, n_emitter) sorted by number of emitters
    :rtype: list
    """
    tuple_list = []
    for adj in adj_arr:
        g = nx.from_numpy_array(adj)
        n_emit = height_max(graph=g)
        tuple_list.append((adj, n_emit))
    sorted_list = sorted(tuple_list, key=lambda x: x[1])
    return sorted_list

def get_relabel_map(g1, g2):
    """
    Finds the map between nodes of g1 and g2 if they are isomorphic.

    :param g1: networkx graph or adjacency matrix
    :type g1: nx.Graph or np.ndarray
    :param g2: networkx graph or adjacency matrix
    :type g2: nx.Graph or np.ndarray
    :return: a dictionary mapping each node of g1 to a node of g2. The map may not be unique.
    :rtype: dict
    """
    # g1 and g2 can be adj matrices or nx graphs.
    if not isinstance(g1, nx.Graph):
        g1 = nx.from_numpy_array(g1)
        assert isinstance(g1, nx.Graph)
    if not isinstance(g2, nx.Graph):
        g2 = nx.from_numpy_array(g2)
        assert isinstance(g2, nx.Graph)
    if np.array_equal(nx.to_numpy_array(g1), nx.to_numpy_array(g2)):
        return {**{-1: "self"}, **dict(zip(g1.nodes(), g2.nodes()))}
    GM = isomorphism.GraphMatcher(g1, g2)
    assert GM.is_isomorphic()
    return GM.mapping

def iso_finder(
    adj_matrix,
    n_iso,
    rel_inc_thresh=0.2,
    allow_exhaustive=True,
    sort_emit=False,
    label_map=False,
    thresh=None,
    seed=None,
):
    """
    The function permutes the labels of the vertices of a graph to get n_iso distinct isomorphic graphs. The original
    graph will also be returned as the first element of the list and counted toward the number of found cases.
    The maximum number of possible distinct cases may be less than n_iso. The graph G with the nodes relabeled using
    consecutive integers.

    :param adj_matrix: initial adjacency matrix or graph
    :type adj_matrix: numpy.ndarray
    :param n_iso: number of the isomorphic graphs to look for (including the original graph)
    :type n_iso: int
    :param rel_inc_thresh: a threshold value between 0 and 1. The closer to 0 the closer we get to an exhaustive search.
    :type rel_inc_thresh: float
    :param allow_exhaustive: whether exhaustive search over all possible permutations is allowed or not. The runtime may
    become too long if this parameter is True.
    :type allow_exhaustive: bool
    :param sort_emit: if True the outcome is sorted by the number of emitters needed to generate each graph
    :type sort_emit: bool
    :param label_map: if True, the function also returns a list of dictionaries, each representing the label mapping
     between the original input graph and the one in the returned list.
    :type label_map: bool
    :param thresh: a threshold value that determines how many random trials is performed in search for new permutations.
    The default is 5 times the number of needed permutations (=n_iso).
    :type thresh: int
    :param seed: the seed for random sampling of labels
    :type seed: int
    :return: An array of adjacency matrices. The length of the array may be less than n_iso since the maximum number of
    isomorphic graphs for the particular input may be reached. Or the function iteration is interrupted due to threshold
    values. If the label_map is True, function returns a tuple, the first item being the array of adjacency matrices and
     the second one a list of dictionaries that are label maps between the original and relabeled graphs.
    :rtype: numpy.ndarray or tuple (numpy.ndarray, list[dict])
    """
    n_node = adj_matrix.shape[0]
    n_max = np.math.factorial(n_node)
    n_label = n_iso
    labels_arr = _label_finder(n_label, n_node, seed=seed, thresh=thresh)
    adj_arr = automorph_check(adj_matrix, labels_arr)
    if len(adj_arr) >= n_iso:
        adj_arr = adj_arr[:n_iso]
        return adj_arr
    else:
        rel_inc = 1  # relative increase ratio of the number of non-redundant cases as the number of labeling increase
        all_checked = False
        while len(adj_arr) < n_iso and rel_inc > rel_inc_thresh and not all_checked:
            success_ratio = len(adj_arr) / n_label if len(adj_arr) != 0 else 0.5
            add_n = int(n_label / success_ratio) + 1 - n_label
            n_label = int(n_label / success_ratio) + 1
            if n_label > n_max:
                n_label = n_max
                all_checked = True
                if (not allow_exhaustive) and success_ratio < 0.5:
                    warnings.warn(
                        f"Only {len(adj_arr)} isomorphic graphs were found due to high symmetry."
                        f"This may not be the maximum value. Allow for exhaustive search to check"
                    )
                    return adj_arr

            # update the seed used in _add_labels to get new values when we repeat it in the loop
            seed = seed + 1 if (seed is not None) else None
            labels_arr = _add_labels(
                labels_arr,
                add_n,
                exhaustive=allow_exhaustive * all_checked,
                seed=seed,
                thresh=thresh,
            )
            n1 = len(adj_arr)
            adj_arr = automorph_check(adj_matrix, labels_arr)
            n2 = len(adj_arr)
            if n2 >= n1:
                # update increase ratio and relative increase ratio
                inc_ratio = n2 / (n1 + 1)  # plus one is to avoid division by zero
                rel_inc = inc_ratio * success_ratio
        if all_checked:
            warnings.warn(f"Maximum of {n2} possible isomorphic graphs exist")
        elif rel_inc < rel_inc_thresh:
            warnings.warn(
                f"Only {n2} isomorphic graph were found. Consider decreasing rel_inc_thresh to possibly "
                f"get more."
            )
        else:
            pass
        if sort_emit:
            adj_arr = np.array([x[0] for x in emitter_sorted(adj_arr[:n_iso])])
        if label_map:
            mapping = []
            for new_adj in adj_arr:
                mapping.append(get_relabel_map(adj_matrix, new_adj))
            return adj_arr[:n_iso], mapping
        return adj_arr[:n_iso]

def lc_orbit_finder(
    graph: nx.Graph,
    comp_depth=None,
    orbit_size_thresh=None,
    with_iso=False,
    rand=False,
    rep_allowed=False,
):
    """
    Given a graph this functions tries all possible local-complementation sequences of length up to comp_depth to
    come up with new distinct graphs in the orbit of the input graph. The comp_depth determines the maximum depth of the
     orbit explored.

    :param graph: original graph
    :type graph: nx.Graph
    :param comp_depth: the maximum length of the sequence of local-complementations applied on the graph; if None,
                        continue till the required number of graphs are found, or no new graphs are found.
    :type comp_depth: int
    :param orbit_size_thresh: sets a limit on the maximum number of orbit graphs to look for
    :type orbit_size_thresh: int
    :param with_iso: if true, isomorph graphs will be kept in the orbit
    :type with_iso: bool
    :param rep_allowed: if true the orbit finder does not check for iso or auto morphism; useful for large graphs.
    :type rand: bool
    :return: list of distinct graphs in the orbit of original graph
    :rtype: list[nx.Graph]
    """
    orbit_list = [graph.copy()]
    node_list = [*graph.nodes]
    new_g = graph
    if rand:
        for i in range(min(10, len(graph))):
            node_x = np.random.randint(0, len(graph))
            new_g = local_comp_graph(new_g, node_x)
        orbit_list = [new_g]
    if orbit_size_thresh == 1:
        return orbit_list
    new_graphs = 1
    i = 0
    if comp_depth is None:
        cond = lambda x: True
    else:
        cond = lambda x: bool(x < comp_depth)

    while cond(i):
        len_before = len(orbit_list)
        # iterate over the new graphs appended to the end of the orbit list
        for graph in orbit_list[-new_graphs:]:
            if rand:
                np.random.shuffle(node_list)
            for node in node_list:
                if graph.degree(node) > 1:
                    g_lc = local_comp_graph(graph, node)
                    if not rep_allowed:
                        if not check_isomorphism(g_lc, orbit_list, _only_auto=with_iso):
                            orbit_list.append(g_lc)
                    else:
                        orbit_list.append(g_lc)
                if (
                    orbit_size_thresh is not None
                    and len(orbit_list) >= orbit_size_thresh
                ):
                    return orbit_list[:orbit_size_thresh]
                if (
                    rand and len(orbit_list) > len_before
                ):  # in random case we only keep 1 new graph
                    break
        # orbit_list = remove_iso(orbit_list)
        len_after = len(orbit_list)
        new_graphs = len_after - len_before
        # print("new graphs", new_graphs)
        if new_graphs == 0:
            break
        i += 1
    return orbit_list