WIP: FT State Preparation Through Code Symmetries

pyproject.toml

@@ -261,6 +261,7 @@ ignore = [
     "PLR2004", # Magic value used in comparison
     "PLC0415", # Import should be at top of file
     "S101",    # Use of assert detected
+    "T201",    # avoid removing print statements
 ]
 isort.required-imports = ["from __future__ import annotations"]
 

src/mqt/qecc/analog_information_decoding/simulators/analog_tannergraph_decoding.py

@@ -247,7 +247,7 @@ def run(self, samples: int) -> dict[str, Any]:
                     self.code_params,
                     self.eb_precision,
                 ):
-                    print("Result has converged.")  # noqa: T201
+                    print("Result has converged.")
                     break
 
         x_ler, x_ler_eb, x_wer, x_wer_eb = calculate_error_rates(x_success_cnt, runs, self.code_params)

src/mqt/qecc/analog_information_decoding/simulators/quasi_single_shot_v2.py

@@ -283,6 +283,6 @@ def run(self, samples: int = 1) -> dict[str, Any]:
             if run % self.save_interval == 1:
                 self._save_results(success_cnt, run)
                 if _check_convergence(success_cnt, run, self.code_params, self.eb_precision):
-                    print("Converged")  # noqa: T201
+                    print("Converged")
                     break
         return self._save_results(success_cnt, run)

src/mqt/qecc/analog_information_decoding/simulators/simulation.py

@@ -735,7 +735,7 @@ def run(self, samples: int) -> dict[str, Any]:
                     self.code_params,
                     self.eb_precision,
                 ):
-                    print("Result has converged.")  # noqa: T201
+                    print("Result has converged.")
                     break
 
         x_ler, x_ler_eb, x_wer, x_wer_eb = calculate_error_rates(x_success_cnt, runs, self.code_params)

src/mqt/qecc/circuit_synthesis/permutation_st.py

@@ -0,0 +1,69 @@
+"""This module implements the features of state preparation with qubit permutation."""
+
+from __future__ import annotations
+
+import logging
+from typing import TYPE_CHECKING
+
+import numpy as np
+from qiskit import QuantumCircuit
+
+from mqt.qecc.codes.css_code import CSSCode
+
+if TYPE_CHECKING:  # pragma: no cover
+    import numpy.typing as npt
+
+logger = logging.getLogger(__name__)
+
+
+def get_steane_prep_circ() -> QuantumCircuit:
+    """Function that builds up a state preparation circuit for the Steane code."""
+    circ = QuantumCircuit(7)
+    circ.h(0)
+    circ.h(1)
+    circ.h(3)
+    circ.cx(0, 6)
+    circ.cx(1, 2)
+    circ.cx(3, 5)
+    circ.cx(0, 4)
+    circ.cx(1, 5)
+    circ.cx(0, 2)
+    circ.cx(3, 4)
+    circ.cx(5, 6)
+    return circ
+
+
+def get_steane_prep_circ_twice() -> QuantumCircuit:
+    """Function that builds up a state preparation circuit for the Steane code."""
+    circ = QuantumCircuit(14)
+    circ.h(0)
+    circ.h(1)
+    circ.h(3)
+    circ.cx(0, 6)
+    circ.cx(1, 2)
+    circ.cx(3, 5)
+    circ.cx(0, 4)
+    circ.cx(1, 5)
+    circ.cx(0, 2)
+    circ.cx(3, 4)
+    circ.cx(5, 6)
+    circ.h(7)
+    circ.h(8)
+    circ.h(10)
+    circ.cx(7, 13)
+    circ.cx(8, 9)
+    circ.cx(10, 12)
+    circ.cx(7, 11)
+    circ.cx(8, 12)
+    circ.cx(7, 9)
+    circ.cx(10, 11)
+    circ.cx(12, 13)
+    return circ
+
+
+def get_steane_css_object() -> CSSCode:
+    """Return the Steane Code as CSS object."""
+    distance: int = 3
+    hx: npt.NDArray[np.int8] = np.array([[0, 0, 0, 1, 1, 1, 1], [0, 1, 1, 0, 0, 1, 1], [1, 0, 1, 0, 1, 0, 1]])
+    hz: npt.NDArray[np.int8] = np.array([[0, 0, 0, 1, 1, 1, 1], [0, 1, 1, 0, 0, 1, 1], [1, 0, 1, 0, 1, 0, 1]])
+    return CSSCode(distance, hx, hz, distance, distance)

src/mqt/qecc/circuit_synthesis/simulation.py

@@ -32,6 +32,7 @@ def __init__(
         p: float,
         zero_state: bool = True,
         parallel_gates: bool = True,
+        check_logical_0: bool = False,
     ) -> None:
         """Initialize the simulator.
 
@@ -41,6 +42,7 @@ def __init__(
             p: The error rate.
             zero_state: Whether thezero state is prepared or nor.
             parallel_gates: Whether to allow for parallel execution of gates.
+            check_logical_0: Whether to check flag measurements or the logical state
         """
         if code.Hx is None or code.Hz is None:
             msg = "The code must have both X and Z checks."
@@ -62,6 +64,7 @@ def __init__(
         self.n_measurements = 0
         self.stim_circ = stim.Circuit()
         self.decoder = LutDecoder(code)
+        self.check_logical_0 = check_logical_0
         self.set_p(p)
 
     def set_p(self, p: float) -> None:
@@ -271,7 +274,14 @@ def _simulate_batch(self, shots: int = 1024) -> tuple[int, int]:
 
         # Filter events where the verification circuit flagged
         verification_measurements = self.x_verification_measurements + self.z_verification_measurements
-        index_array = np.where(np.all(detection_events[:, verification_measurements] == 0, axis=1))[0]
+        if self.check_logical_0:
+            # Compute dot products for all rows
+            dot_products = np.dot(detection_events[:, verification_measurements], self.code.Hx.T) % 2
+
+            # Find rows where all dot products are zero
+            index_array = np.where(np.all(dot_products == 0, axis=1))[0]
+        else:
+            index_array = np.where(np.all(detection_events[:, verification_measurements] == 0, axis=1))[0]
         filtered_events = detection_events[index_array].astype(np.int8)
 
         if len(filtered_events) == 0:  # All events were discarded

src/mqt/qecc/circuit_synthesis/state_prep.py

@@ -15,6 +15,7 @@
 from ..codes import InvalidCSSCodeError
 from .synthesis_utils import (
     build_css_circuit_from_cnot_list,
+    get_permutation_group,
     heuristic_gaussian_elimination,
     iterative_search_with_timeout,
     measure_flagged,
@@ -34,6 +35,7 @@
     import numpy.typing as npt
     from qiskit import DAGNode
     from qiskit.quantum_info import PauliList
+    from sympy.combinatorics import Permutation
 
     from ..codes import CSSCode
 
@@ -239,6 +241,71 @@
     return StatePrepCircuit(circ, code, zero_state)
 
 
+def fs_disjunct_spc(
+    code: CSSCode, automorph_generators: list[list[int]], optimize_depth: bool = True, zero_state: bool = True
+) -> None:
+    """Returns dictionary of permutations for the given QEC.
+
+    This function takes a QEC and the corresponding strong automorphism group. It then checks every permutation and
+    returns a dictionary with all permutations as keys and the overlapping propagated errors in the fault set.
+
+    Args:
+        code: The CSS code to prepare the state for.
+        optimize_depth: If True, optimize the depth of the circuit. This may lead to a higher number of CNOTs.
+        zero_state: If True, prepare the +1 eigenstate of the Z basis. If False, prepare the +1 eigenstate of the X basis.
+        automorph_generators: foobar
+    """
+    logger.info("Starting heuristic state preparation.")
+    if code.Hx is None or code.Hz is None:
+        msg = "The code must have both X and Z stabilizers defined."
+        raise InvalidCSSCodeError(msg)
+
+    checks = code.Hx if zero_state else code.Hz
+    assert checks is not None
+    checks, cnots = heuristic_gaussian_elimination(checks, parallel_elimination=optimize_depth)
+    # NOTE: everything above is code duplication from heuristic_prep_circuit()
+
+    permutation_group = get_permutation_group(group_generators=automorph_generators)
+    # BUG: I have those prints since the stateprep circuit below gives another fault set than the one
+    # returned all the way at the bottom without this if section and I do not understand why
+    circ = _build_state_prep_circuit_from_back(checks.copy(), cnots.copy(), zero_state)
+    stateprepcirc = StatePrepCircuit(circ, code, zero_state)
+
+    # NOTE: get fault sets based on code distance
+    stateprepcirc.compute_fault_set()
+    fault_set_1 = stateprepcirc.x_fault_sets[1]
+    if stateprepcirc.code.distance == 5:
+        fault_set_1 = fault_set_1[np.where(fault_set_1.sum(axis=1) > 2)]
+    if stateprepcirc.code.distance == 7:
+        stateprepcirc.compute_fault_set(2)
+        fault_set_2 = stateprepcirc.x_fault_sets[2]
+        fault_set_1 = fault_set_1[np.where(fault_set_1.sum(axis=1) > 3)]
+        fault_set_2 = fault_set_2[np.where(fault_set_2.sum(axis=1) > 3)]
+
+    # HACK: Overlap search has 3 nested for loops, check for improvement
+    overlapping_faults: list[npt.NDArray[np.int8]] = []
+    perm_overlap: dict[Permutation, list[npt.NDArray[np.int8]]] = {}
+    original_faults = fault_set_1
+    print(f"fault set 1: \n{fault_set_1}")
+    if stateprepcirc.code.distance == 7:
+        original_faults = np.vstack(fault_set_1, fault_set_2)
+        print(f"fault set 2: \n{fault_set_2}")
+    print(f"comparison fault set: \n{original_faults}")
+    for perm in permutation_group:
+        perm_fs = fault_set_1[:, perm.array_form]
+        for pf in perm_fs:
+            for ff in original_faults:
+                if code.stabilizer_eq_x_error(pf, ff):
+                    # FIX: without this print statement append gets changed to extend on save
+                    print(f"permutation: {perm}")
+                    print(f"stabilizer equivalent error: \n{pf}(permuted) and \n{ff}(not permuted)")
+                    overlapping_faults.append(pf)
+        perm_overlap[perm] = overlapping_faults
+        overlapping_faults = []
+
+    return perm_overlap
+
+
 def depth_optimal_prep_circuit(
     code: CSSCode,
     zero_state: bool = True,

src/mqt/qecc/circuit_synthesis/synthesis_utils.py

@@ -12,6 +12,7 @@
 from ldpc import mod2
 from qiskit import AncillaRegister, ClassicalRegister, QuantumCircuit
 from stim import Circuit
+from sympy.combinatorics import Permutation, PermutationGroup
 
 if TYPE_CHECKING:  # pragma: no cover
     from collections.abc import Callable
@@ -96,7 +97,7 @@ def heuristic_gaussian_elimination(
         matrix: The matrix to perform Gaussian elimination on.
         parallel_elimination: Whether to prioritize elimination steps that act on disjoint columns.
 
-    returns:
+    Returns:
         The reduced matrix and a list of the elimination steps taken. The elimination steps are represented as tuples of the form (i, j) where i is the column being eliminated with and j is the column being eliminated.
     """
     matrix = matrix.copy()
@@ -152,6 +153,20 @@ def is_reduced() -> bool:
     return matrix, eliminations
 
 
+def get_permutation_group(group_generators: list[list[int]]) -> list[Permutation]:
+    """Based on the generators of the permutation group find the whole permutation group.
+
+    Args:
+        group_generators: A list of generators of the permutation group. Each generator is given as a list of integers describing the permutation. E.g. for a S7 generator: [0, 3, 2, 1, 6, 5, 4]
+
+    Returns:
+        Returns a list of Permutation object coming form sympy.combinatorics.
+    """
+    group_generators = [Permutation(generator) for generator in group_generators]
+    g = PermutationGroup(group_generators)
+    return list(g.generate())
+
+
 def gaussian_elimination_min_column_ops(
     matrix: npt.NDArray[np.int8],
     termination_criteria: Callable[[Any], z3.BoolRef],
@@ -166,7 +181,7 @@ def gaussian_elimination_min_column_ops(
         termination_criteria: A function that takes a boolean matrix as input and returns a Z3 boolean expression that is true if the matrix is considered reduced.
         max_eliminations: The maximum number of eliminations to perform.
 
-    returns:
+    Returns:
         The reduced matrix and a list of the elimination steps taken. The elimination steps are represented as tuples of the form (i, j) where i is the column being eliminated with and j is the column being eliminated.
     """
     n = matrix.shape[1]
@@ -235,7 +250,7 @@ def gaussian_elimination_min_parallel_eliminations(
         termination_criteria: A function that takes a boolean matrix as input and returns a Z3 boolean expression that is true if the matrix is considered reduced.
         max_parallel_steps: The maximum number of parallel elimination steps to perform.
 
-    returns:
+    Returns:
         The reduced matrix and a list of the elimination steps taken. The elimination steps are represented as tuples of the form (i, j) where i is the column being eliminated with and j is the column being eliminated.
     """
     columns = np.array([

src/mqt/qecc/ecc_qiskit_wrapper.py

@@ -81,7 +81,7 @@ def print_simulation_results(result: Result, n_shots: int, threshold_probability
         # Print all results > threshold_probability
         if summarized_counts[result_id] / n_shots > threshold_probability or printed_results == 0:
             result_string = str(result_id)
-            print("State |" + result_string + "> probability " + str(summarized_counts[result_id] / n_shots))  # noqa: T201
+            print("State |" + result_string + "> probability " + str(summarized_counts[result_id] / n_shots))
             printed_results += 1
             if printed_results == 1000:
                 break
@@ -154,7 +154,7 @@ def main() -> None:
     ecc_frequency = args.fq
     ecc_export_filename = args.e
     if forced_simulator is not None and "stabilizer" in forced_simulator and "A" in error_channels:
-        print(  # noqa: T201
+        print(
             'Warning: Non-unitary errors (such as for example amplitude damping ("A")) are not suitable for simulation '
             "with a stabilizer based simulator and may cause an error during the simulation."
         )
@@ -168,7 +168,7 @@ def main() -> None:
     circ = load(open_qasm_file)
 
     if not any(gate.operation.name == "measure" for gate in circ.data):
-        print("Warning: No measurement gates found in the circuit. Adding measurement gates to all qubits.")  # noqa: T201
+        print("Warning: No measurement gates found in the circuit. Adding measurement gates to all qubits.")
         circ.measure_all()
 
     # Initializing the quantum circuit
@@ -178,13 +178,13 @@ def main() -> None:
         circ = loads(result["circ"])
 
     if ecc_export_filename is not None:
-        print("Exporting circuit to: " + str(ecc_export_filename))  # noqa: T201
+        print("Exporting circuit to: " + str(ecc_export_filename))
         with pathlib.Path(ecc_export_filename).open("w", encoding=locale.getpreferredencoding(False)) as f:
             dump(circ, f)
         return
 
     size = circ.num_qubits
-    print(  # noqa: T201
+    print(
         "_____Trying to simulate with "
         + str(error_channels)
         + " (prob="