1. Solver¶

class TimeReversedSolver(SolverBase):
    """
    This solver finds a quantum circuit that produces the target state in a time-reversed fashion.
    """

    name = "time-reversed"

    one_qubit_ops = list(ops.one_qubit_cliffords())

    def __init__(
        self,
        target,
        metric: MetricBase,
        compiler: CompilerBase,
        circuit: CircuitDAG = None,
        io: IO = None,
        noise_model_mapping=None,
        *args,
        **kwargs,
    ):
        """
        Initialize a TimeReversedSolver object

        :param target: target quantum state
        :type target: QuantumState
        :param metric: a metric to be evaluated for the final state and circuit
        :type metric: MetricBase
        :param compiler: a backend compiler to be used for compilation of the quantum circuit
        :type compiler: CompilerBase
        :param circuit: an initial quantum circuit
        :type circuit: CircuitDAG
        :param io: input/output object for saving logs, intermediate results, circuits, etc.
        :type io: IO
        :param noise_model_mapping: a dictionary that associates each operation type to a noise model
        :type noise_model_mapping: dict
        """

        super().__init__(target, metric, compiler, circuit, io, *args, **kwargs)
        if target.rep_type != "s":
            target.convert_representation("s")
        tableau = target.rep_data.data.to_stabilizer()
        self.n_emitter = self.determine_n_emitters(tableau)
        self.n_photon = tableau.n_qubits
        self.noise_simulation = True

        if noise_model_mapping is None:
            noise_model_mapping = {"e": dict(), "p": dict(), "ee": dict(), "ep": dict()}
            self.noise_simulation = False
        elif type(noise_model_mapping) is not dict:
            raise TypeError(
                f"Datatype {type(noise_model_mapping)} is not a valid noise_model_mapping. "
                f"noise_model_mapping should be a dict or None"
            )
        self.noise_model_mapping = noise_model_mapping

    def solve(self):
        """
        The main function for the solver

        :return: nothing
        :rtype: None
        """

        self.compiler.noise_simulation = self.noise_simulation
        # initialize the circuit
        circuit = CircuitDAG(
            n_emitter=self.n_emitter, n_photon=self.n_photon, n_classical=1
        )

        stabilizer_tableau = self.target.rep_data.data.to_stabilizer()

        # add emitter qubits in |0> state
        for _ in range(self.n_emitter):
            stabilizer_tableau = sfs.insert_qubit(
                stabilizer_tableau, stabilizer_tableau.n_qubits
            )

        # main loop
        for j in range(self.n_photon, 0, -1):
            # transform the stabilizers into echelon gauge
            stabilizer_tableau = sfs.rref(stabilizer_tableau)
            height_list = height.height_func_list(
                stabilizer_tableau.x_matrix, stabilizer_tableau.z_matrix
            )

            height_list = [0] + height_list
            if height_list[j] < height_list[j - 1]:
                # apply time-reversed measurement and update the tableau
                self._time_reversed_measurement(circuit, stabilizer_tableau, j - 1)
                stabilizer_tableau = sfs.rref(stabilizer_tableau)

            # apply photon-absorption and update the stabilizer tableau
            self._add_photon_absorption(circuit, stabilizer_tableau, j - 1)

        stabilizer_tableau = sfs.rref(stabilizer_tableau)

        # make sure that all photonic qubits are in the |0> state
        assert np.array_equal(
            stabilizer_tableau.x_matrix[0 : self.n_photon, 0 : self.n_photon],
            np.zeros((self.n_photon, self.n_photon)),
        )
        assert np.array_equal(
            stabilizer_tableau.z_matrix[0 : self.n_photon, 0 : self.n_photon],
            np.eye(self.n_photon),
        )

        # transform all the emitter qubits to the |0> state
        _, inverse_circuit = sfs.inverse_circuit(stabilizer_tableau.copy())

        self._add_gates_from_str(circuit, stabilizer_tableau, inverse_circuit)

        # correct the phase of generators that act only on emitter qubits
        for i in range(self.n_emitter):
            if stabilizer_tableau.phase[self.n_photon + i] == 1:
                transform.x_gate(stabilizer_tableau, self.n_photon + i)
                self._add_one_qubit_gate(circuit, [ops.SigmaX], self.n_photon + i)

        # final circuit
        circuit.validate()

        compiled_state = self.compiler.compile(circuit)

        compiled_state.partial_trace(
            keep=[*range(self.n_photon)], dims=(self.n_photon + self.n_emitter) * [2]
        )
        # evaluate the metric
        score = self.metric.evaluate(compiled_state, circuit)
        self.result = (score, circuit.copy())

    @staticmethod
    def determine_n_emitters(tableau):
        """
        Determine the minimum number of emitters needed

        :param tableau: a tableau that represents the stabilizer state
        :type tableau: StabilizerTableau
        :return: the minimum number of emitters required to generate the state
        :rtype: int
        """

        tableau = sfs.rref(tableau)
        height_list = height.height_func_list(tableau.x_matrix, tableau.z_matrix)
        return max(height_list)

    def _time_reversed_measurement(self, circuit, tableau, photon_index):
        """
        Apply the time-reversed measurement

        :param circuit: a quantum circuit
        :type circuit: CircuitDAG
        :param tableau: a stabilizer tableau
        :type tableau: StabilizerTableau
        :param photon_index: index of the photon for the conditional X gate
        :type photon_index: int
        :return: nothing
        :rtype: None
        """

        # find the generators that act only on the emitters
        sum_table = (
            tableau.x_matrix[:, 0 : self.n_photon]
            + tableau.z_matrix[:, 0 : self.n_photon]
        )
        possible_generators = np.where(~sum_table.any(axis=1))[0]

        assert len(possible_generators) > 0
        # In Li et al.'s implementation, it always picks the first one
        generator_index = possible_generators[0]

        # find the first nontrivial emitter position in this generator
        emitter_indices = self._find_emitter_indices(tableau, generator_index)
        emitter_index = int(emitter_indices[0])

        self._single_out_emitter(circuit, tableau, generator_index, emitter_index)

        transform.hadamard_gate(tableau, self.n_photon + emitter_index)
        self._add_measurement_cnot_and_reset(circuit, emitter_index, photon_index)

        # transform this generator so that it acts nontrivially on only one emitter and that Pauli is X
        transform.cnot_gate(tableau, self.n_photon + emitter_index, photon_index)

    def _add_photon_absorption(self, circuit, tableau, photon_index):
        """
        Add a photon absorption (time-reversed photon emission) event

        :param circuit: a quantum circuit
        :type circuit: CircuitDAG
        :param tableau: a stabilizer tableau
        :type tableau: StabilizerTableau
        :param photon_index: index of photon to be absorbed
        :type photon_index: int
        :return: nothing
        :rtype: None
        """

        for i in range(tableau.n_qubits - 1, -1, -1):
            if height.leftmost_nontrivial_index(tableau, i) == photon_index:
                generator_index = i
                break

        gate_list = self._change_pauli_type(tableau, generator_index, photon_index, "z")
        self._add_one_qubit_gate(circuit, gate_list, photon_index)

        emitter_indices = self._find_emitter_indices(tableau, generator_index)
        emitter_index = int(emitter_indices[0])
        for i in range(self.n_emitter):
            # transform this generator so that it has only Z on emitter registers
            gate_list = self._change_pauli_type(
                tableau,
                generator_index,
                self.n_photon + i,
                "z",
            )
            self._add_one_qubit_gate(circuit, gate_list, self.n_photon + i)

        # now this generator contains only ZZ on emitter qubits
        # transform this generator so that for emitter qubits, there is only one Z.
        self._transform_generator_emitters(
            circuit, tableau, generator_index, emitter_index
        )

        # if the phase of this generator is -1, flip it by adding X gate to the emitter qubit
        if tableau.phase[generator_index] == 1:
            transform.x_gate(tableau, self.n_photon + emitter_index)
            self._add_one_qubit_gate(
                circuit, [ops.SigmaX], self.n_photon + emitter_index
            )

        # emit (absorb) this photon by the chosen emitter
        self._add_emitter_photon_cnot(circuit, emitter_index, photon_index)
        transform.cnot_gate(tableau, self.n_photon + emitter_index, photon_index)

        # transform all generators (except the chosen generator) with Z on this photonic register
        # by eliminating Z
        pauli_z_list = sfs.one_pauli_type_finder(
            tableau.x_matrix, tableau.z_matrix, [0, photon_index], "z"
        )

        pauli_z_list = np.setdiff1d(pauli_z_list, [generator_index])
        for i in pauli_z_list:
            # use left multiplication instead of right multiplication
            sfs.tab_row_sum(tableau, generator_index, i)

    def _add_gates_from_str(self, circuit, tableau, gate_str_list):
        """
        Add gates to disentangle all emitter qubits. This is used in the last step.

        :param circuit: a quantum circuit
        :type circuit: CircuitDAG
        :param tableau: a stabilizer tableau
        :type tableau: StabilizerTableau
        :param gate_str_list: a list of gates to be applied
        :type gate_str_list: list[(str, int) or (str, int, int)]
        :return: nothing
        :rtype: None
        """

        for gate in gate_str_list:
            if gate[0] == "H":
                self._add_one_qubit_gate(circuit, [ops.Hadamard], gate[1])
                transform.hadamard_gate(tableau, gate[1])
            elif gate[0] == "P":
                # add the inverse of the phase gate
                self._add_one_qubit_gate(circuit, [ops.SigmaZ, ops.Phase], gate[1])
                transform.phase_gate(tableau, gate[1])
            elif gate[0] == "X":
                self._add_one_qubit_gate(circuit, [ops.SigmaX], gate[1])
                transform.x_gate(tableau, gate[1])
            elif gate[0] == "CNOT":
                self._add_one_emitter_cnot(
                    circuit, gate[1] - self.n_photon, gate[2] - self.n_photon
                )
                transform.cnot_gate(tableau, gate[1], gate[2])
            elif gate[0] == "CZ":
                self._add_one_qubit_gate(circuit, [ops.Hadamard], gate[2])
                transform.hadamard_gate(tableau, gate[2])
                self._add_one_emitter_cnot(
                    circuit, gate[1] - self.n_photon, gate[2] - self.n_photon
                )
                transform.cnot_gate(tableau, gate[1], gate[2])
                self._add_one_qubit_gate(circuit, [ops.Hadamard], gate[2])
                transform.hadamard_gate(tableau, gate[2])
            else:
                raise ValueError("Invalid gate in the list.")

    def _change_pauli_type(self, tableau, row, column, result="z"):
        """
        Transform the Pauli of the given qubit (column) in the given generator (row)
        so that it becomes the Pauli that is specified by the input parameter (result).

        :param tableau: a stabilizer tableau of the state
        :type tableau: StabilizerTableau
        :param row: the generator index in the tableau
        :type row: int
        :param column: the qubit index in the tableau
        :type column: int
        :param result: the Pauli type to be transformed to
        :type result: str
        :return: a list of one-qubit gates to be added to the circuit
        :rtype: list[ops.OperationBase]
        """

        gate_list = []
        if tableau.x_matrix[row, column] == 1 and tableau.z_matrix[row, column] == 0:
            # current is X
            if result.lower() == "z":
                transform.hadamard_gate(tableau, column)
                gate_list.append(ops.Hadamard)
            elif result.lower() == "y":
                transform.phase_gate(tableau, column)
                # add the inverse of phase gate to the gate list
                gate_list.append(ops.SigmaZ)
                gate_list.append(ops.Phase)

        elif tableau.x_matrix[row, column] == 1 and tableau.z_matrix[row, column] == 1:
            # current is Y
            if result.lower() == "z":
                transform.phase_dagger_gate(tableau, column)
                transform.hadamard_gate(tableau, column)

                # add the inverse of phase dagger gate and inverse of Hadamard gate to the gate list
                gate_list.append(ops.Phase)
                gate_list.append(ops.Hadamard)
            elif result.lower() == "x":
                transform.phase_dagger_gate(tableau, column)

                # add the inverse of phase dagger gate to the gate list
                gate_list.append(ops.Phase)

        if tableau.x_matrix[row, column] == 0 and tableau.z_matrix[row, column] == 1:
            # current is Z
            if result.lower() == "x":
                transform.hadamard_gate(tableau, column)
                gate_list.append(ops.Hadamard)
            elif result.lower() == "y":
                transform.hadamard_gate(tableau, column)
                transform.phase_gate(tableau, column)

                # add the inverse of Hadamard gate and the inverse of phase gate to the gate list
                gate_list.append(ops.Hadamard)
                gate_list.append(ops.SigmaZ)
                gate_list.append(ops.Phase)

        return gate_list

    def _add_one_qubit_gate(self, circuit, gate_list, index):
        """
        Add a one-qubit gate to the circuit

        :param circuit: a quantum circuit
        :type circuit: CircuitDAG
        :param gate_list: a list of one-qubit gates to be added to the circuit
        :type gate_list: list[ops.OperationBase]
        :param index: the qubit position where this one-qubit gate is applied
        :type index: int
        :return: nothing
        :rtype: None
        """

        if index >= self.n_photon:
            reg_type = "e"
            reg = index - self.n_photon
        else:
            reg_type = "p"
            reg = index

        edge = circuit.dag.out_edges(nbunch=f"{reg_type}{reg}_in", keys=True)
        edge = list(edge)[0]
        next_node = circuit.dag.nodes[edge[1]]
        next_op = next_node["op"]
        if isinstance(next_op, ops.OneQubitGateWrapper):
            # if two OneQubitGateWrapper gates are next to each other, combine them
            gate_list = next_op.operations + gate_list
            gate_list = ops.simplify_local_clifford(gate_list)
            if gate_list == [ops.Identity, ops.Identity]:
                circuit.remove_op(edge[1])
                return
        else:
            # simplify the gate to be one of the 24 local Clifford gates
            gate_list = ops.simplify_local_clifford(gate_list)
            if gate_list == [ops.Identity, ops.Identity]:
                return
        if reg_type == "e":
            noise = self._wrap_noise(gate_list, self.noise_model_mapping["e"])
        else:
            noise = self._wrap_noise(gate_list, self.noise_model_mapping["p"])
        gate = ops.OneQubitGateWrapper(
            gate_list,
            reg_type=reg_type,
            register=reg,
            noise=noise,
        )
        if isinstance(next_op, ops.OneQubitGateWrapper):
            circuit.replace_op(edge[1], gate)
        else:
            circuit.insert_at(gate, [edge])

    def _add_one_emitter_cnot(self, circuit, control_emitter, target_emitter):
        """
        Add a CNOT between two emitter qubits

        :param circuit: a quantum circuit
        :type circuit: CircuitDAG
        :param control_emitter: register index of the control emitter
        :type control_emitter: int
        :param target_emitter: register index of the target emitter
        :type target_emitter: int
        :return: nothing
        :rtype: None
        """

        control_edge = circuit.dag.out_edges(nbunch=f"e{control_emitter}_in", keys=True)
        edge0 = list(control_edge)[0]
        target_edge = circuit.dag.out_edges(nbunch=f"e{target_emitter}_in", keys=True)
        edge1 = list(target_edge)[0]
        gate = ops.CNOT(
            control=circuit.dag.edges[edge0]["reg"],
            control_type="e",
            target=circuit.dag.edges[edge1]["reg"],
            target_type="e",
            noise=self._identify_noise(ops.CNOT, self.noise_model_mapping["ee"]),
        )
        circuit.insert_at(gate, [edge0, edge1])

    def _transform_generator_emitters(
        self, circuit, tableau, generator_index, emitter_index
    ):
        """
        Transform a given generator such that it contains only one Pauli Z on emitter qubits

        :param circuit: a quantum circuit
        :type circuit: CircuitDAG
        :param tableau: a stabilizer tableau
        :type tableau: StabilizerTableau
        :param generator_index: the index of stabilizer generator in the tableau
        :type generator_index: int
        :param emitter_index: the index of chosen emitter register
        :type emitter_index: int
        :return: nothing
        :rtype: None
        """

        if self.n_emitter == 1:
            return

        # only Z
        assert not np.any(tableau.x_matrix[generator_index])
        emitters = np.nonzero(tableau.z_matrix[generator_index, self.n_photon :])[0]
        target_emitter = emitter_index
        rest_emitters = np.setdiff1d(emitters, [target_emitter])
        for control_emitter in rest_emitters:
            # find two emitters with ZZ in the generator specified by generator_index
            # add CNOT to remove one Z
            self._add_one_emitter_cnot(circuit, control_emitter, target_emitter)
            transform.cnot_gate(
                tableau,
                self.n_photon + control_emitter,
                self.n_photon + target_emitter,
            )

    def _transform_generator_emitters_advanced(self, circuit, tableau, generator_index):
        """
        Transform a given generator such that it contains only one Pauli Z on emitter qubits
        This function will be used when we can select which emitter for the shortest circuit depth.

        :param circuit: a quantum circuit
        :type circuit: CircuitDAG
        :param tableau: a stabilizer tableau
        :type tableau: StabilizerTableau
        :param generator_index: the index of stabilizer generator in the tableau
        :type generator_index: int
        :return: the index of emitter chosen for the shortest circuit depth
        :rtype: int
        """

        # wait for the circuit class refactoring
        emitter_depth = circuit.register_depth["e"]
        # allow selection of emitters
        if self.n_emitter == 1:
            return 0

        assert not np.any(tableau.x_matrix[generator_index])
        emitters = np.nonzero(tableau.z_matrix[generator_index, self.n_photon :])[0]
        while len(emitters) > 1:
            # find two emitters with ZZ in the generator specified by generator_index
            # if multiple, find these two with the shortest circuit depths
            min_index = np.argmin(emitter_depth[emitters])
            target_emitter = int(emitters[min_index])
            rest_emitters = np.setdiff1d(emitters, [target_emitter])
            min_index = np.argmin(emitter_depth[rest_emitters])
            control_emitter = int(rest_emitters[min_index])
            self._add_one_emitter_cnot(circuit, control_emitter, target_emitter)
            transform.cnot_gate(
                tableau,
                self.n_photon + control_emitter,
                self.n_photon + target_emitter,
            )

            emitters = np.nonzero(tableau.z_matrix[generator_index, self.n_photon :])[0]
        return int(emitters[0])

    def _add_emitter_photon_cnot(self, circuit, emitter_index, photon_index):
        """
        Add a CNOT gate between two edges

        :param circuit: a quantum circuit
        :type circuit: CircuitDAG
        :return: nothing
        :rtype: None
        """

        emitter_edge = circuit.dag.out_edges(nbunch=f"e{emitter_index}_in", keys=True)
        edge0 = list(emitter_edge)[0]
        photonic_edge = circuit.dag.out_edges(nbunch=f"p{photon_index}_in", keys=True)
        edge1 = list(photonic_edge)[0]

        gate = ops.CNOT(
            control=emitter_index,
            control_type="e",
            target=photon_index,
            target_type="p",
            noise=self._identify_noise(ops.CNOT, self.noise_model_mapping["ep"]),
        )
        gate.add_labels("Fixed")
        circuit.insert_at(gate, [edge0, edge1])

    def _find_emitter_indices(self, tableau, generator_index):
        """
        Find nontrivial emitter positions in this generator

        :param tableau: a stabilizer tableau
        :type tableau: StabilizerTableau
        :param generator_index: the index of the generator
        :type generator_index: int
        :return: nontrivial emitter positions in this generator
        :rtype: numpy.ndarray
        """

        generator_sum = (
            tableau.x_matrix[generator_index, self.n_photon :]
            + tableau.z_matrix[generator_index, self.n_photon :]
        )
        return np.nonzero(generator_sum)[0]

    def _single_out_emitter(self, circuit, tableau, generator_index, emitter_index):
        """
        Turn the chosen generator into a single Pauli Z. This function is used in the time-reversed measurement.

        :param circuit: a quantum circuit
        :type circuit: CircuitDAG
        :param tableau: a stabilizer tableau
        :type tableau: StabilizerTableau
        :param generator_index: the index of stabilizer generator in the tableau
        :type generator_index: int
        :param emitter_index: the index of emitter
        :type emitter_index: int
        :return: nothing
        :rtype: None
        """

        # first transform this generator so that it only contains Zs.
        for i in range(self.n_emitter):
            gate_list = self._change_pauli_type(
                tableau, generator_index, self.n_photon + i, "z"
            )
            if len(gate_list) > 0:
                self._add_one_qubit_gate(circuit, gate_list, self.n_photon + i)

        # disentangle the chosen emitter from the rest
        self._transform_generator_emitters(
            circuit, tableau, generator_index, emitter_index
        )

        if tableau.phase[generator_index] == 1:
            transform.x_gate(tableau, self.n_photon + emitter_index)
            self._add_one_qubit_gate(
                circuit, [ops.SigmaX], self.n_photon + emitter_index
            )

    def _add_measurement_cnot_and_reset(self, circuit, emitter_index, photon_index):
        """
        Add a MeasurementCNOTandReset operation from an emitter qubit to a photonic qubit such that no consecutive
        MeasurementCNOTReset is allowed. This operation cannot be added before the photonic qubit is initialized.

        :param circuit: a quantum circuit
        :type circuit: CircuitDAG
        :return: nothing
        :rtype: None
        """

        emitter_edge = circuit.dag.out_edges(nbunch=f"e{emitter_index}_in", keys=True)
        edge0 = list(emitter_edge)[0]
        photonic_edge = circuit.dag.out_edges(nbunch=f"p{photon_index}_in", keys=True)
        edge1 = list(photonic_edge)[0]

        gate = ops.MeasurementCNOTandReset(
            control=circuit.dag.edges[edge0]["reg"],
            control_type="e",
            target=circuit.dag.edges[edge1]["reg"],
            target_type="p",
            noise=self._identify_noise(
                ops.MeasurementCNOTandReset, self.noise_model_mapping["ep"]
            ),
        )
        gate.add_labels("Fixed")
        circuit.insert_at(gate, [edge0, edge1])

    @property
    def solver_info(self):
        """
        Return the solver setting

        :return: attributes of the solver
        :rtype: dict
        """

        def op_names(op_list):
            op_name = []
            for op_val in op_list:
                if isinstance(op_val, list):
                    op_name.append([op.__name__ for op in op_val])
                else:
                    op_name.append(op_val.__name__)
            return op_name

        return {
            "solver name": self.name,
            "seed": self.last_seed,
            "One-qubit ops": op_names(self.one_qubit_ops),
        }

class AlternateTargetSolver:
    def __init__(
        self,
        target: nx.Graph or QuantumState = None,
        metric: MetricBase = Infidelity,
        compiler: CompilerBase = StabilizerCompiler(),
        noise_compiler: CompilerBase = DensityMatrixCompiler(),
        io: IO = None,
        noise_model_mapping=None,
        solver_setting=None,
        seed=None,
    ):
        """
        Constructor for AlternateGraphSolver

        :param target: target graph state
        :type target: QuantumState or nx.Graph
        :param metric: a metric
        :type metric: MetricBase
        :param compiler: a compiler
        :type compiler: CompilerBase
        :param noise_compiler: compiler for noise simulation
        :type noise_compiler: CompilerBase
        :param io: input/output
        :type io: IO
        :param noise_model_mapping:
        :type noise_model_mapping: dict
        :param solver_setting:
        :type solver_setting:
        :param seed: a random seed
        :type seed: int
        """
        if solver_setting is None:
            solver_setting = AlternateTargetSolverSetting(
                n_iso_graphs=1, n_lc_graphs=1, graph_metric=pre.graph_metric_lists[0]
            )
        if noise_model_mapping is None:
            noise_model_mapping = {"e": dict(), "p": dict(), "ee": dict(), "ep": dict()}
            self.noise_simulation = False
            self.monte_carlo = False
        elif noise_model_mapping == "depolarizing":
            self.noise_simulation = True
            self.depolarizing_rate = solver_setting.depolarizing_rate
            if solver_setting.monte_carlo:
                self.monte_carlo = True
                self.mc_params = solver_setting.monte_carlo_params
                self.mc_params["map"] = self.mc_depol()
            else:
                self.monte_carlo = False
                noise_model_mapping = self.depol_noise_map()
        elif type(noise_model_mapping) is not dict:
            raise TypeError(
                f"Datatype {type(noise_model_mapping)} is not a valid noise_model_mapping. "
                f"noise_model_mapping should be a dict or None"
            )
        else:
            self.noise_simulation = True
            self.monte_carlo = False
            if solver_setting.monte_carlo:
                self.monte_carlo = True
                self.mc_params = solver_setting.monte_carlo_params

        self.noise_model_mapping = noise_model_mapping

        if isinstance(target, nx.Graph):
            self.target_graph = target
            self.target = QuantumState(target, rep_type="g")
        elif isinstance(target, QuantumState):
            if target.mixed:
                raise ValueError(
                    "AlternateTargetSolver does not support mixed states as its target."
                )
            self.target = target
            if target.rep_type == "s":
                target_graph = stabilizer_to_graph(
                    target.rep_data.tableau.to_stabilizer()
                )
                self.target_graph = target_graph[0][1]
            elif target.rep_type == "dm":
                target_graph = density_to_graph(target.rep_data.data)
                self.target_graph = nx.from_numpy_array(target_graph)

            elif target.rep_type == "g":
                self.target_graph = target.rep_data.data
            else:
                raise ValueError("Wrong representation data type for QuantumState")

        self.metric = metric
        self.noise_compiler = noise_compiler
        self.compiler = compiler
        self.solver_setting = solver_setting
        self.seed = seed
        self.result = None
        self.mc_list = []

    def solve(self):
        """
        Finds alternative circuits to generate the target graph or a relabeled version of it by searching through
        alternative isomorphic and LC-equivalent graphs. Notice that the returned circuits generate the target state
        up to relabeling if this feature is enabled in the setting. Otherwise, user's exact target graph is produced.

        :return: a dictionary where keys are the circuits and values are themselves dictionaries containing the LC graph
         used as the intermediate states and relabeling map between the actual target and the graph the circuit
         generates. {circuit: {'g':graph, 'map': relabel_map}}
        :rtype: dict
        """
        setting = self.solver_setting
        n_iso = setting.n_iso_graphs
        n_lc = setting.n_lc_graphs
        adj_matrix = nx.to_numpy_array(self.target_graph)

        iso_adjs = iso_finder(
            adj_matrix,
            n_iso,
            rel_inc_thresh=setting.rel_inc_thresh,
            allow_exhaustive=setting.allow_exhaustive,
            sort_emit=setting.sort_emitter,
            label_map=setting.label_map,
            thresh=setting.iso_thresh,
            seed=self.seed,
        )
        results_list = []
        mc_list = []
        # repeater graph state check
        adj_rgs = nx.to_numpy_array(
            repeater_graph_states(int(len(self.target_graph) / 2))
        )
        target_is_rgs = True if np.array_equal(adj_rgs, adj_matrix) else False
        # list of tuples [(circuit, dict={'g': graph used to find circuit, 'map': relabel map with target})]
        iso_graphs = [nx.from_numpy_array(adj) for adj in iso_adjs]
        for iso_graph in iso_graphs:
            if setting.lc_method == "rgs":
                lc_graphs = rgs_orbit_finder(iso_graph)[:n_lc]
            elif setting.lc_method == "linear":
                lc_graphs = linear_partial_orbit(iso_graph)[:n_lc]
            elif setting.lc_method == "depth_first":
                lc_graphs = depth_first_orbit(iso_graph)[:n_lc]
            elif setting.lc_method == "lc_with_iso":
                lc_graphs = lc_orbit_finder(
                    iso_graph,
                    comp_depth=setting.lc_orbit_depth,
                    orbit_size_thresh=n_lc,
                    with_iso=True,
                )
            elif setting.lc_method == "random":
                lc_graphs = lc_orbit_finder(
                    iso_graph,
                    comp_depth=setting.lc_orbit_depth,
                    orbit_size_thresh=n_lc,
                    rand=True,
                )
            elif setting.lc_method == "random_with_iso":
                lc_graphs = lc_orbit_finder(
                    iso_graph,
                    comp_depth=setting.lc_orbit_depth,
                    orbit_size_thresh=n_lc,
                    with_iso=True,
                    rand=True,
                )
            elif setting.lc_method == "random_with_rep":
                lc_graphs = lc_orbit_finder(
                    iso_graph,
                    comp_depth=setting.lc_orbit_depth,
                    orbit_size_thresh=n_lc,
                    with_iso=True,
                    rand=True,
                    rep_allowed=True,
                )
            elif setting.lc_method is None:
                lc_graphs = lc_orbit_finder(
                    iso_graph, comp_depth=setting.lc_orbit_depth, orbit_size_thresh=n_lc
                )
            else:
                raise ValueError(
                    "LC method is not valid. Please set it to a valid string or None"
                )
            lc_circ_list = []
            lc_score_list = []
            rmap = get_relabel_map(self.target_graph, iso_graph)

            for lc_graph in lc_graphs:
                circuit = graph_to_circ(lc_graph)
                success, conversion_gates = slc.lc_check(
                    lc_graph, iso_graph, validate=True
                )
                try:
                    assert success, "LC graphs are not actually LC equivalent!"
                    slc.state_converter_circuit(lc_graph, iso_graph, validate=True)
                except:
                    raise UserWarning("LC conversion failed")
                if not lc_graph.adj == iso_graph.adj:
                    conversion_ops = slc.str_to_op(conversion_gates)
                else:
                    conversion_ops = []
                for op in conversion_ops:
                    circuit.add(op)
                if self.noise_simulation:
                    if self.monte_carlo:
                        # generate a unique seed for the monte carlo object based on user input seeds
                        if self.seed is not None and self.mc_params["seed"] is not None:
                            rng = np.random.default_rng(
                                self.seed * self.mc_params["seed"]
                            )
                        else:
                            # seedless randomness if at least one of the user's seeds is None
                            rng = np.random.default_rng()
                        n_mc = max(
                            int(self.mc_params["n_sample"] or 0),
                            int(self.mc_params["n_single"] or 0)
                            * int(self.mc_params["n_parallel"] or 0),
                        )
                        mc_seed = rng.integers(
                            low=1, high=int(10e6 * max(n_mc, 2)), size=1
                        )
                        # end of random seed generation
                        mc = mcn.MonteCarloNoise(
                            circuit,
                            self.mc_params["n_sample"],
                            self.mc_params["map"],
                            self.mc_params["compiler"],
                            mc_seed,
                        )
                        mc_list.append(mc)
                        if self.mc_params["n_parallel"] is not None:
                            n_total = (
                                self.mc_params["n_parallel"]
                                * self.mc_params["n_single"]
                            )
                            assert (
                                n_total > 0
                            ), "n_single and n_parallel both must be integers > 1 or None"
                            self.solver_setting.monte_carlo_params["n_sample"] = n_total
                            # multicore parallel processing
                            # ray.init()
                            noise_score = mcn.parallel_monte_carlo(
                                mc,
                                self.mc_params["n_parallel"],
                                self.mc_params["n_single"],
                            )
                            # ray.shutdown()

                        else:
                            noise_score = mc.run()
                    else:
                        circuit = circuit.assign_noise(self.noise_model_mapping)
                        noise_score = self.noise_score(circuit, rmap)
                    lc_score_list.append(noise_score)
                else:
                    lc_score_list.append(0.0001)

                lc_circ_list.append(circuit)

            for i, circ in enumerate(lc_circ_list):
                results_list.append(
                    (circ, {"g": lc_graphs[i], "score": lc_score_list[i], "map": rmap})
                )
            # iso_lc_circuit_dict[get_relabel_map(self.target, iso_graph)] = dict(zip(lc_graphs, lc_circ_list))
        # remove redundant auto-morph graphs
        adj_list = [nx.to_numpy_array(result[1]["g"]) for result in results_list]
        set_list = []
        for i in range(len(adj_list)):
            already_found = False
            for s in set_list:
                if i in s:
                    already_found = True
                    break
            if not already_found:
                s = {i}
                for j in range(i + 1, len(adj_list)):
                    if np.array_equal(adj_list[i], adj_list[j]):
                        s.add(j)
                set_list.append(s)
        # set_list now contains the equivalent group of graphs in the result's dict
        # remove redundant items of the dict
        redundant_indices = [index for s in set_list for index in list(s)[1:]]
        redundant_indices.sort()
        for index in redundant_indices[::-1]:
            del results_list[index]

        # results setter
        self.mc_list = mc_list
        circ_list = [x[0] for x in results_list]
        properties = list(results_list[0][1].keys()) if results_list else None
        self.result = SolverResult(circ_list, properties)
        # [r[1][p] for p in properties for r in results_list]
        for p in properties:
            self.result[p] = [r[1][p] for r in results_list]

        return results_list

    def noise_score(self, circuit, relabel_map):
        """


        :param circuit:
        :type circuit:
        :param relabel_map:
        :type relabel_map:
        :return:
        :rtype:
        """
        graph = self.target_graph
        old_adj = nx.to_numpy_array(graph)
        n = graph.number_of_nodes()
        # relabel the target graph
        new_labels = np.array([relabel_map[i] for i in range(n)])
        new_adj = relabel(old_adj, new_labels)
        graph = nx.from_numpy_array(new_adj)

        c_tableau = get_clifford_tableau_from_graph(graph)
        ideal_state = QuantumState(c_tableau, rep_type="s")
        if isinstance(self.noise_compiler, DensityMatrixCompiler):
            ideal_state = QuantumState(graph_to_density(graph), rep_type="dm")
        metric = Infidelity(ideal_state)

        compiler = self.noise_compiler
        compiler.noise_simulation = True
        compiled_state = compiler.compile(circuit=circuit)
        # trace out emitter qubits
        compiled_state.partial_trace(
            keep=list(range(circuit.n_photons)),
            dims=(circuit.n_photons + circuit.n_emitters) * [2],
        )

        # evaluate the metric
        score = metric.evaluate(compiled_state, circuit)
        score = 0.00001 if (np.isclose(score, 0.0) and score != 0.0) else score
        return score

    def depol_noise_map(self):
        """


        :return:
        :rtype:
        """

        rate = self.depolarizing_rate
        dep_noise_model_mapping = dict()
        dep_noise_model_mapping["e"] = {
            "SigmaX": nm.DepolarizingNoise(rate),
            "SigmaY": nm.DepolarizingNoise(rate),
            "SigmaZ": nm.DepolarizingNoise(rate),
            "Phase": nm.DepolarizingNoise(rate),
            "PhaseDagger": nm.DepolarizingNoise(rate),
            "Hadamard": nm.DepolarizingNoise(rate),
        }
        dep_noise_model_mapping["p"] = {}  # dep_noise_model_mapping["e"]
        # dep_noise_model_mapping["ee"] = {}
        # dep_noise_model_mapping["ep"] = {}
        dep_noise_model_mapping["ee"] = {"CNOT": nm.DepolarizingNoise(rate)}
        dep_noise_model_mapping["ep"] = {"CNOT": nm.DepolarizingNoise(rate)}
        return dep_noise_model_mapping

    def mc_depol(self):
        """
        Returns a Monte-Carlo noise map for depolarizing noise. Currently only emitter gates are noisy.
        :return: mcn.McNoiseMap
        :rtype: mcn.McNoiseMap
        """
        rate = self.depolarizing_rate / 3
        mc_noise = mcn.McNoiseMap()
        for gate_type in [
            "Hadamard",
            "Phase",
            "PhaseDagger",
            "SigmaX",
            "SigmaY",
            "SigmaZ",
        ]:
            mc_noise.add_gate_noise(
                "e",
                gate_type,
                [
                    (nm.PauliError("X"), rate),
                    (nm.PauliError("Y"), rate),
                    (nm.PauliError("Z"), rate),
                ],
            )

        mc_noise.add_gate_noise(
            "ee",
            "CNOT",
            [
                (nm.PauliError("X"), rate),
                (nm.PauliError("Y"), rate),
                (nm.PauliError("Z"), rate),
            ],
        )
        mc_noise.add_gate_noise(
            "ep",
            "CNOT",
            [
                (nm.PauliError("X"), rate),
                (nm.PauliError("Y"), rate),
                (nm.PauliError("Z"), rate),
            ],
        )
        return mc_noise