Repository: QCHackers/tqec (python)

Issue:
Integrate automatic detector computation into `tqec`
**Is your feature request related to a problem? Please describe.**
#249 aims at including in the code base a generic way to automatically compute detectors from a given quantum circuit checking some pre-conditions. The implementation in #249 is voluntarily disconnected from the main `tqec` library, to be able to extract it as an independent library if we think this might be useful to the QEC community. But this disconnection implies that, for the moment, automatic detector computation does not use the additional knowledge provided by the `tqec` library.

In particular, if #249 is merged without significant changes, the workflow that will need to be followed is the following:
1. Create the necessary plaquettes and templates.
2. Assemble the created plaquettes and templates into a structure representing the QEC computation.
3. For each code distance `d` and noise model `e`:
    1. Create the quantum circuit representing the QEC computation.
    2. Compute detectors over the **whole quantum circuit** and add the computed detectors to the circuit.

Even though detector computation is quite efficient most of the time, the computation might get very expensive on large code distance. 

**Describe the solution you'd like**
In practice, it might not be needed to perform detector computation at step 3.2. over the whole circuit. Using knowledge on the code locality (e.g., detectors can only include measurements that are within a ball of radius `X` in space around a central point, and maximum `Y` temporal steps before the current step), we might be able to extract from the representation obtained after step 2. a limited number of sub-circuits (subsets of templates with plaquettes) that:
1. are centered on a specific plaquette (i.e., one index particular location in a template that will be replaced with a plaquette),
2. for each sub-circuit, only include plaquettes that are at a distance `X` from the center plaquette,
3. cover each plaquette in the template.

Basically, finding local (both in time and in space) regions that should contain at least one detector involving a defined central measurement, computing the detector for each of these regions, and re-using the computation when the same local region is found multiple times in the QEC computation.


Pull request:
[WIP] Automatic and spatially-local detector computation
This PR is based on #318.

This is still a work in progress as some changes from #305 are needed to complete this PR.

In the end, this PR should provide all the necessary tools to efficiently compute detectors when making the translation `BlockGraph` -> `cirq.Circuit` (see #298).

Fixes #259 

Patch:
diff --git a/src/tqec/block/block.py b/src/tqec/block/block.py
index 762979b9..c018fa15 100644
--- a/src/tqec/block/block.py
+++ b/src/tqec/block/block.py
@@ -13,7 +13,9 @@ from typing_extensions import override
 from tqec.block.enums import BlockDimension
 from tqec.circuit.circuit import generate_circuit
 from tqec.circuit.operations.measurement import Measurement, RepeatedMeasurement
+from tqec.circuit.operations.operation import Detector
 from tqec.circuit.schedule import ScheduledCircuit, merge_scheduled_circuits
+from tqec.detectors.computation import compute_detectors_in_last_timestep
 from tqec.exceptions import TQECException
 from tqec.plaquette.library.empty import empty_square_plaquette
 from tqec.plaquette.plaquette import Plaquette, Plaquettes
@@ -239,15 +241,13 @@ class StandardComputationBlock(ComputationBlock):
 
     @override
     def instantiate(self) -> cirq.Circuit:
-        circuit = generate_circuit(self.template, self.initial_plaquettes.collection)
+        circuit = generate_circuit(self.template, self.initial_plaquettes)
         if self.repeating_plaquettes is not None:
             repetitions = self.repeating_plaquettes.number_of_rounds(self.template.k)
             plaquettes = self.repeating_plaquettes.plaquettes
-            inner_circuit = generate_circuit(
-                self.template, plaquettes.collection
-            ).freeze()
+            inner_circuit = generate_circuit(self.template, plaquettes).freeze()
             circuit += cirq.CircuitOperation(inner_circuit, repetitions=repetitions)
-        circuit += generate_circuit(self.template, self.final_plaquettes.collection)
+        circuit += generate_circuit(self.template, self.final_plaquettes)
         return circuit
 
     @override
@@ -330,6 +330,34 @@ class StandardComputationBlock(ComputationBlock):
                 measurements.extend(m.measurements())
         return measurements
 
+    def compute_detectors(
+        self, subtemplate_radius_trial_range: range = range(1, 2)
+    ) -> list[list[Detector]]:
+        initial_detectors = compute_detectors_in_last_timestep(
+            self.template, (self.initial_plaquettes,), subtemplate_radius_trial_range
+        )
+        repeating_detectors: list[Detector]
+        final_detectors: list[Detector]
+        if self.repeating_plaquettes is not None:
+            repeating_detectors = compute_detectors_in_last_timestep(
+                self.template,
+                (self.initial_plaquettes, self.repeating_plaquettes.plaquettes),
+                subtemplate_radius_trial_range,
+            )
+            final_detectors = compute_detectors_in_last_timestep(
+                self.template,
+                (self.repeating_plaquettes.plaquettes, self.final_plaquettes),
+                subtemplate_radius_trial_range,
+            )
+        else:
+            repeating_detectors = []
+            final_detectors = compute_detectors_in_last_timestep(
+                self.template,
+                (self.initial_plaquettes, self.final_plaquettes),
+                subtemplate_radius_trial_range,
+            )
+        return [initial_detectors, repeating_detectors, final_detectors]
+
 
 @dataclass
 class Computation:
diff --git a/src/tqec/block/library/open.py b/src/tqec/block/library/open.py
index fbf3cddf..ffd08ea3 100644
--- a/src/tqec/block/library/open.py
+++ b/src/tqec/block/library/open.py
@@ -23,7 +23,6 @@ from tqec.templates.qubit import (
     QubitTemplate,
     QubitVerticalBorders,
 )
-from tqec.templates.scale import LinearFunction
 
 
 def _with_resets_on_data_qubits_on_side(
diff --git a/src/tqec/circuit/circuit.py b/src/tqec/circuit/circuit.py
index 6f25930c..ee852a09 100644
--- a/src/tqec/circuit/circuit.py
+++ b/src/tqec/circuit/circuit.py
@@ -5,17 +5,19 @@ from collections import defaultdict
 from copy import deepcopy
 
 import cirq
+import numpy
+import numpy.typing as npt
 
 from tqec.circuit.schedule import ScheduledCircuit, merge_scheduled_circuits
 from tqec.exceptions import TQECException
-from tqec.plaquette.plaquette import Plaquette
+from tqec.plaquette.plaquette import Plaquette, Plaquettes
 from tqec.position import Displacement
 from tqec.templates.base import Template
 
 
 def generate_circuit(
     template: Template,
-    plaquettes: list[Plaquette] | dict[int, Plaquette] | defaultdict[int, Plaquette],
+    plaquettes: Plaquettes,
 ) -> cirq.Circuit:
     """Generate a quantum circuit from a template and its plaquettes.
 
@@ -49,38 +51,85 @@ def generate_circuit(
     """
     # Check that the user gave enough plaquettes.
     if (
-        not isinstance(plaquettes, defaultdict)
+        not plaquettes.has_default
         and len(plaquettes) != template.expected_plaquettes_number
     ):
         raise TQECException(
             f"{len(plaquettes)} plaquettes have been provided, but "
             f"{template.expected_plaquettes_number} were expected."
         )
-    if isinstance(plaquettes, ty.Mapping) and 0 in plaquettes:
-        raise TQECException(
-            "If using a dictionary, the input plaquettes parameter should not "
-            f"contain the entry 0. Found a value ({plaquettes[0]}) at entry 0."
-        )
-
-    # If plaquettes are given as a list, make that a dict to simplify the following operations
-    if isinstance(plaquettes, list):
-        plaquettes = {i + 1: plaquette for i, plaquette in enumerate(plaquettes)}
 
     # instantiate the template with the appropriate plaquette indices.
     # Index 0 is "no plaquette" by convention and should not be included here.
     _indices = list(range(1, template.expected_plaquettes_number + 1))
     template_plaquettes = template.instantiate(_indices)
     increments = template.get_increments()
+
+    return generate_circuit_from_instantiation(
+        template_plaquettes, plaquettes, increments
+    )
+
+
+def generate_circuit_from_instantiation(
+    plaquette_array: npt.NDArray[numpy.int_],
+    plaquettes: Plaquettes,
+    increments: Displacement,
+) -> cirq.Circuit:
+    """Generate a quantum circuit from an array of plaquette indices and the
+    associated plaquettes.
+
+    This is one of the core methods of the `tqec` package. It generates a quantum circuit
+    from a spatial description of where the plaquettes should be located as well as the
+    actual plaquettes used.
+
+    This function requires that a few pre-conditions on the inputs are met:
+    1. the number of plaquettes provided should match the number of plaquettes required by
+    the provided template.
+    2. all the provided plaquettes should be implemented on cirq.GridQubit instances **only**.
+
+    If any of the above pre-conditions is not met, the inputs are considered invalid, in which
+    case this function **might** raise an error.
+
+    Args:
+        plaquette_array: an array of indices referencing Plaquette instances in the
+            `plaquettes` argument.
+        plaquettes: description of the computation that should happen at
+            different time-slices of the quantum error correction experiment (or
+            at least part of it). If provided as a dictionary, plaquettes should be
+            1-indexed (i.e., ``0 not in plaquettes`` should be ``True``).
+        increments: the displacement between each plaquette origin.
+
+    Returns:
+        a cirq.Circuit instance implementing the (part of) quantum error
+        correction experiment represented by the provided inputs.
+
+    Raises:
+        TQECException: if any index in `plaquette_array` is not correctly associated to a
+            plaquette in `plaquettes` or if `plaquettes` is provided as a dictionary and
+            `0 in plaquettes`.
+    """
+    # Check that the user gave enough plaquettes.
+    indices = numpy.unique(plaquette_array)
+    # Remove the first 0 entry in indices if present
+    if indices[0] == 0:
+        indices = indices[1:]
+    expected_plaquettes_number = len(indices)
+    if not plaquettes.has_default and len(plaquettes) != expected_plaquettes_number:
+        raise TQECException(
+            f"{len(plaquettes)} plaquettes have been provided, but "
+            f"{expected_plaquettes_number} were expected."
+        )
+
     # Plaquettes indices are starting at 1 in template_plaquettes. To avoid
     # offsets in the following code, we add an empty circuit at position 0.
     plaquette_circuits = {0: ScheduledCircuit(cirq.Circuit())} | {
-        i: plaquettes[i].circuit for i in _indices
+        i: plaquettes[i].circuit for i in indices
     }
 
     # Generate the ScheduledCircuit instances for each plaquette instantiation
     all_scheduled_circuits: list[ScheduledCircuit] = []
     plaquette_index: int
-    for row_index, line in enumerate(template_plaquettes):
+    for row_index, line in enumerate(plaquette_array):
         for column_index, plaquette_index in enumerate(line):
             if plaquette_index != 0:
                 scheduled_circuit = deepcopy(plaquette_circuits[plaquette_index])
diff --git a/src/tqec/circuit/detectors/flow.py b/src/tqec/circuit/detectors/flow.py
index 11b4e5e8..3bd490a9 100644
--- a/src/tqec/circuit/detectors/flow.py
+++ b/src/tqec/circuit/detectors/flow.py
@@ -223,7 +223,7 @@ class FragmentLoopFlows:
 
 
 def build_flows_from_fragments(
-    fragments: list[Fragment | FragmentLoop],
+    fragments: ty.Sequence[Fragment | FragmentLoop],
 ) -> list[FragmentFlows | FragmentLoopFlows]:
     """Compute and return the stabilizer flows of the provided fragments.
 
diff --git a/src/tqec/circuit/detectors/pauli.py b/src/tqec/circuit/detectors/pauli.py
index 4ca32bb6..eea8c0f7 100644
--- a/src/tqec/circuit/detectors/pauli.py
+++ b/src/tqec/circuit/detectors/pauli.py
@@ -13,6 +13,7 @@ import operator
 import typing as ty
 
 import stim
+
 from tqec.exceptions import TQECException
 
 PAULI_STRING_TYPE = ty.Literal["I", "X", "Y", "Z"]
@@ -154,7 +155,6 @@ class PauliString:
         """
         ret = PauliString(self._pauli_by_qubit.copy())
         for op in collapse_operators:
-            print(f"Collapsing {op} over {ret}")
             if not ret.commutes(op):
                 raise TQECException(
                     f"Cannot collapse {ret} by a non-commuting operator {op}."
diff --git a/src/tqec/circuit/operations/measurement.py b/src/tqec/circuit/operations/measurement.py
index d734cb7d..71f6e9be 100644
--- a/src/tqec/circuit/operations/measurement.py
+++ b/src/tqec/circuit/operations/measurement.py
@@ -48,9 +48,7 @@ class Measurement:
     def __lt__(self, other: object) -> bool:
         if not isinstance(other, Measurement):
             raise NotImplementedError(f"Cannot compare {type(self)} < {type(other)}.")
-        if self.offset == other.offset:
-            return self.qubit < other.qubit
-        return self.offset < other.offset
+        return (self.offset, self.qubit) < (other.offset, other.qubit)
 
     def __repr__(self) -> str:
         return f"{self.__class__.__name__}({self.qubit}, {self.offset})"
@@ -134,8 +132,8 @@ def get_measurements_from_circuit(
                 ordered_measured_qubits.append(qubit)
     measured_qubit_offset: defaultdict[cirq.GridQubit, int] = defaultdict(lambda: -1)
     measurements: list[Measurement] = []
-    for mqubit in ordered_measured_qubits:
+    for mqubit in ordered_measured_qubits[::-1]:
         index = measured_qubit_offset[mqubit]
         measured_qubit_offset[mqubit] = index - 1
         measurements.append(Measurement(mqubit, index))
-    return measurements
+    return measurements[::-1]
diff --git a/src/tqec/circuit/schedule.py b/src/tqec/circuit/schedule.py
index 1e7c8b3d..5a3d1772 100644
--- a/src/tqec/circuit/schedule.py
+++ b/src/tqec/circuit/schedule.py
@@ -546,6 +546,8 @@ def merge_scheduled_circuits(circuits: list[ScheduledCircuit]) -> cirq.Circuit:
         # is considered mergeable, and can be removed if another operation in the list
         # is considered equal (and has the mergeable tag).
         non_duplicated_operations = remove_duplicate_operations(operations)
-        all_moments.append(cirq.Moment(*non_duplicated_operations))
+        all_moments.append(
+            cirq.Moment(*sorted(non_duplicated_operations, key=lambda op: op.qubits))
+        )
 
     return cirq.Circuit(all_moments)
diff --git a/src/tqec/detectors/computation.py b/src/tqec/detectors/computation.py
new file mode 100644
index 00000000..367cf12c
--- /dev/null
+++ b/src/tqec/detectors/computation.py
@@ -0,0 +1,245 @@
+import typing as ty
+import warnings
+
+import cirq
+import numpy
+import numpy.typing as npt
+import stim
+import stimcirq
+
+from tqec.circuit.circuit import generate_circuit_from_instantiation
+from tqec.circuit.detectors.flow import build_flows_from_fragments
+from tqec.circuit.detectors.fragment import Fragment
+from tqec.circuit.detectors.match import (
+    MatchedDetector,
+    match_boundary_stabilizers,
+    match_detectors_within_fragment,
+)
+from tqec.circuit.operations.measurement import (
+    Measurement,
+    get_measurements_from_circuit,
+)
+from tqec.circuit.operations.operation import Detector
+from tqec.exceptions import TQECException, TQECWarning
+from tqec.plaquette.plaquette import Plaquettes
+from tqec.position import Displacement, Position
+from tqec.templates.base import Template
+
+
+def _get_qubit_mapping(circuit: cirq.Circuit) -> dict[cirq.GridQubit, int]:
+    """Get a mapping from qubit to indices for the provided circuit.
+
+    This function requires the whole circuit in order to avoid issues such
+    as qubits only appearing in some time steps and not the others.
+
+    Args:
+        circuit: the circuit to create the mapping for.
+
+    Raises:
+        TQECException: if any qubit of the provided circuit is not a
+            `cirq.GridQubit` instance.
+
+    Returns:
+        a mapping qubit -> index for all the qubits appearing in the circuit
+        provided. Each qubit should be mapped to a different index.
+    """
+    mapping: dict[cirq.GridQubit, int] = {}
+    for qubit in sorted(circuit.all_qubits()):
+        if not isinstance(qubit, cirq.GridQubit):
+            raise TQECException(
+                "Found a qubit that is not an instance of cirq.GridQubit."
+            )
+        mapping[qubit] = len(mapping)
+    return mapping
+
+
+def _get_measurement_offset_mapping(circuit: cirq.Circuit) -> dict[int, Measurement]:
+    """Get a mapping from measurement offsets as used and returned by the
+    package `tqec.circuit.detectors` to a `Measurement` instance.
+
+    This function returns the mapping from negative offsets that are
+    supposed to each represent a unique measurement in the circuit to
+    `Measurement` instances, that serve the same purpose but use a different
+    encoding.
+
+    Note:
+        As a sanity check, the user of this function is encouraged to check
+        if the qubits in each `Measurement` instances returned as values of
+        the mapping corresponds to the expected qubit.
+
+    Args:
+        circuit: the circuit to create the mapping for.
+
+    Raises:
+        TQECException: if any qubit of the generated circuit is not a
+            `cirq.GridQubit` instance.
+
+    Returns:
+        a mapping qubit -> index for all the qubits appearing in the circuit
+        built from the provided `instantiation` and `plaquettes_at_timestep`.
+        Each qubit should be mapped to a different index.
+    """
+    return {
+        -i - 1: m
+        for i, m in enumerate(reversed(get_measurements_from_circuit(circuit)))
+    }
+
+
+def _map_all_measurements_from_offset(
+    detectors: list[MatchedDetector],
+    measurements_by_offset: dict[int, Measurement],
+    qubits_by_indices: dict[int, cirq.GridQubit],
+) -> list[Detector]:
+    ret: list[Detector] = []
+    for d in detectors:
+        measurements: list[Measurement] = []
+        for m in d.measurements:
+            measurement = measurements_by_offset[m.offset]
+            expected_qubit = qubits_by_indices[m.qubit_index]
+            if measurement.qubit != expected_qubit:
+                raise TQECException(
+                    f"Expected qubit {expected_qubit} from index {m.qubit_index} "
+                    f"but got {measurement.qubit}."
+                )
+            measurements.append(measurement)
+        ret.append(Detector(measurements, d.coords))
+    return ret
+
+
+def _compute_detectors_within_subtemplate(
+    subtemplate: npt.NDArray[numpy.int_],
+    plaquettes_at_timestep: tuple[Plaquettes] | tuple[Plaquettes, Plaquettes],
+    increments: Displacement,
+) -> list[Detector]:
+    """Returns detectors within the subtemplate.
+
+    Args:
+        subtemplate: a square 2-dimensional array of integers with odd-length
+            sides representing the arrangement of plaquettes in a subtemplate.
+        plaquettes_at_timestep: a tuple containing either one or two collection(s)
+            of plaquettes each representing one QEC round.
+        increments: spatial increments between each `Plaquette` origin.
+
+    Returns:
+        all the detectors appearing in the subtemplate.
+    """
+    radius = subtemplate.shape[0] // 2
+    # Build subcircuit for each Plaquettes layer
+    subcircuits_cirq: list[cirq.Circuit] = []
+    complete_circuit_cirq = cirq.Circuit()
+    for plaquettes in plaquettes_at_timestep:
+        subcircuit_cirq = generate_circuit_from_instantiation(
+            subtemplate, plaquettes, increments
+        )
+        complete_circuit_cirq += subcircuit_cirq
+        subcircuits_cirq.append(subcircuit_cirq)
+
+    # Construct the different mappings we will need during the computation
+    indices_by_qubit = _get_qubit_mapping(complete_circuit_cirq)
+    qubits_by_index = {i: q for q, i in indices_by_qubit.items()}
+    coordinates_by_index: dict[int, tuple[float, ...]] = {
+        i: (q.row, q.col) for q, i in indices_by_qubit.items()
+    }
+    measurements_by_offset = _get_measurement_offset_mapping(complete_circuit_cirq)
+
+    subcircuits_stim: list[stim.Circuit] = [
+        stimcirq.cirq_circuit_to_stim_circuit(
+            subcircuit_cirq,
+            qubit_to_index_dict=ty.cast(dict[cirq.Qid, int], indices_by_qubit),
+        )
+        for subcircuit_cirq in subcircuits_cirq
+    ]
+
+    fragments = [Fragment(circ) for circ in subcircuits_stim]
+    flows = build_flows_from_fragments(fragments)
+    # Match the detectors
+    matched_detectors = match_detectors_within_fragment(flows[-1], coordinates_by_index)
+    if len(flows) == 2:
+        matched_detectors.extend(
+            match_boundary_stabilizers(flows[-2], flows[-1], coordinates_by_index)
+        )
+    detectors = _map_all_measurements_from_offset(
+        matched_detectors, measurements_by_offset, qubits_by_index
+    )
+
+    # We have two transformations to apply to it:
+    # 1. Change of coordinate system for the qubits. potential_detector is
+    #    using a coordinate system with the origin at the top-left corner of
+    #    the current sub-template, but we need to return detectors that use
+    #    the central plaquette origin as their coordinate system origin.
+    # 2. Transform the MatchedDetector instance into a Detector instance.
+    return [
+        Detector(
+            [
+                m.offset_spatially_by(-radius * increments.x, -radius * increments.y)
+                for m in detector.measurement_data
+            ],
+            detector.coordinates,
+        )
+        for detector in detectors
+    ]
+
+
+def _compute_detectors_for_fixed_radius(
+    template: Template,
+    plaquettes_at_timestep: tuple[Plaquettes] | tuple[Plaquettes, Plaquettes],
+    fixed_subtemplate_radius: int = 2,
+) -> list[Detector]:
+    unique_subtemplates = template.get_spatially_distinct_subtemplates(
+        fixed_subtemplate_radius, avoid_zero_plaquettes=True
+    )
+    increments = template.get_increments()
+    # Each detector in detectors_by_subtemplate is using a coordinate system
+    # centered on the central plaquette origin.
+    detectors_by_subtemplate: dict[int, list[Detector]] = {
+        i: _compute_detectors_within_subtemplate(
+            subtemplate, plaquettes_at_timestep, increments
+        )
+        for i, subtemplate in unique_subtemplates.subtemplates.items()
+    }
+
+    detectors_by_measurements: dict[frozenset[Measurement], Detector] = dict()
+    for i, row in enumerate(unique_subtemplates.subtemplate_indices):
+        for j, subtemplate_index in enumerate(row):
+            if subtemplate_index == 0:
+                continue
+            plaquette_origin = Position(j * increments.x, i * increments.x)
+            for d in detectors_by_subtemplate[subtemplate_index]:
+                offset_measurements = list(
+                    m.offset_spatially_by(plaquette_origin.x, plaquette_origin.y)
+                    for m in d.measurement_data
+                )
+                detectors_by_measurements[frozenset(offset_measurements)] = Detector(
+                    offset_measurements, d.coordinates
+                )
+    return list(detectors_by_measurements.values())
+
+
+def compute_detectors_in_last_timestep(
+    template: Template,
+    plaquettes_at_timestep: tuple[Plaquettes] | tuple[Plaquettes, Plaquettes],
+    subtemplate_radius_trial_range: range = range(1, 2),
+) -> list[Detector]:
+    radius_range = subtemplate_radius_trial_range
+    start, stop, step = radius_range.start, radius_range.stop, radius_range.step
+    # Check the range is monotonically increasing
+    if start < 1 or stop <= start or step < 1:
+        raise TQECException(
+            "The range of subtemplate raius must be a monotonically increasing range of positive integers."
+        )
+
+    measurements_set: set[frozenset[Measurement]] = set()
+    detectors: list[Detector]
+    for radius in radius_range:
+        detectors = _compute_detectors_for_fixed_radius(
+            template, plaquettes_at_timestep, radius
+        )
+        cur_measurements_set = {frozenset(d.measurement_data) for d in detectors}
+        if not measurements_set.issubset(cur_measurements_set):
+            raise TQECException(
+                "The detectors found for a smaller radius are not a subset of the detectors found for a larger radius."
+            )
+        if cur_measurements_set == measurements_set:
+            return detectors
+        measurements_set = cur_measurements_set
+    return detectors
diff --git a/src/tqec/plaquette/plaquette.py b/src/tqec/plaquette/plaquette.py
index a4746070..01c3771b 100644
--- a/src/tqec/plaquette/plaquette.py
+++ b/src/tqec/plaquette/plaquette.py
@@ -74,45 +74,35 @@ class Plaquette:
         return self._measurements
 
 
-def _to_plaquette_dict(
-    plaquettes: list[Plaquette] | dict[int, Plaquette] | defaultdict[int, Plaquette],
-) -> dict[int, Plaquette] | defaultdict[int, Plaquette]:
-    if isinstance(plaquettes, (dict, defaultdict)):
-        return plaquettes
-    return {i: p for i, p in enumerate(plaquettes)}
-
-
 @dataclass(frozen=True)
 class Plaquettes:
-    collection: list[Plaquette] | dict[int, Plaquette] | defaultdict[int, Plaquette]
+    collection: dict[int, Plaquette] | defaultdict[int, Plaquette]
 
     def __post_init__(self) -> None:
-        if not isinstance(self.collection, (list, dict, defaultdict)):
+        if not isinstance(self.collection, (dict, defaultdict)):
             raise TQECException(
                 f"Plaquettes initialized with {type(self.collection)} but only list, "
                 "dict and defaultdict instances are allowed."
             )
+        if 0 in self.collection:
+            raise TQECException(
+                "Found a Plaquette with index 0. This index is reserved to express "
+                '"no plaquette". Please re-number your plaquettes starting from 1.'
+            )
 
     def __getitem__(self, index: int) -> Plaquette:
         return self.collection[index]
 
     def __iter__(self) -> typing.Iterator[Plaquette]:
-        if isinstance(self.collection, list):
-            return iter(self.collection)
-        if isinstance(self.collection, defaultdict):
-            if self.collection.default_factory is not None:
-                default = self.collection.default_factory()
-                return itertools.chain(
-                    self.collection.values(), [default] if default is not None else []
-                )
-            else:
-                return iter(self.collection.values())
-        if isinstance(self.collection, dict):
-            return iter(self.collection.values())
-        else:
-            raise TQECException(
-                f"Plaquette is initialised with the wrong type: {type(self.collection)}."
+        if (
+            isinstance(self.collection, defaultdict)
+            and self.collection.default_factory is not None
+        ):
+            default = self.collection.default_factory()
+            return itertools.chain(
+                self.collection.values(), [default] if default is not None else []
             )
+        return iter(self.collection.values())
 
     @property
     def has_default(self) -> bool:
@@ -126,9 +116,5 @@ class Plaquettes:
         return len(self.collection)
 
     def __or__(self, other: Plaquettes | dict[int, Plaquette]) -> Plaquettes:
-        other_plaquettes = (
-            _to_plaquette_dict(other.collection)
-            if isinstance(other, Plaquettes)
-            else other
-        )
-        return Plaquettes(_to_plaquette_dict(self.collection) | other_plaquettes)
+        other_plaquettes = other.collection if isinstance(other, Plaquettes) else other
+        return Plaquettes(self.collection | other_plaquettes)
diff --git a/src/tqec/templates/grid.py b/src/tqec/templates/grid.py
index 4290b21b..e5d935cf 100644
--- a/src/tqec/templates/grid.py
+++ b/src/tqec/templates/grid.py
@@ -20,6 +20,31 @@ class GridTemplate(RectangularTemplate):
         k: int = 2,
         default_increments: Displacement | None = None,
     ) -> None:
+        """A grid of the same repeating template.
+
+        This class provides a way to efficiently represent a 2-dimensional grid
+        of the same `RectangularTemplate` instance tiled all over the grid.
+
+        Args:
+            index_grid: a sequence of rows representing the arrangement of the
+                provided `tiled_template` with integer indices. An entry of `0`
+                means "empty" or "no template". Then, each entry strictly above
+                `1` represent an instance of `tiled_template` on the grid.
+                Entries with the same value (e.g., in `[[1, 0], [1, 0]]`, two
+                entries have the value `1`) represent `tiled_template` instances
+                that will share the same plaquette indices.
+            tiled_template: base `RectangularTemplate` instance that will be
+                tiled over a grid.
+            k: scaling factor forwarded to the template instances composing the
+                grid. Defaults to 2.
+            default_increments: default increments between two plaquettes.
+
+        Raises:
+            TQECException: if the provided grid is empty (e.g., `[]`).
+            TQECException: if the provided grid has empty rows (e.g., `[[]]`).
+            TQECException: if the provided grid indices are not a contiguous
+                sequence of integers starting from either `0` or `1`.
+        """
         super().__init__(k, default_increments)
         if len(index_grid) == 0:
             raise TQECException("Cannot create a grid with no row.")
@@ -65,7 +90,7 @@ class GridTemplate(RectangularTemplate):
 
         Args:
             plaquette_indices: the plaquette indices that will be forwarded to
-                the underlying Shape instance's instantiate method. Defaults
+                the underlying Template instance's instantiate method. Defaults
                 to `range(1, self.expected_plaquettes_number + 1)` if `None`.
 
         Returns:
diff --git a/src/tqec/templates/subtemplates.py b/src/tqec/templates/subtemplates.py
index aebaf0ab..d1789b93 100644
--- a/src/tqec/templates/subtemplates.py
+++ b/src/tqec/templates/subtemplates.py
@@ -84,7 +84,8 @@ class UniqueSubTemplates:
     subtemplates: dict[int, npt.NDArray[numpy.int_]]
 
     def __post_init__(self) -> None:
-        indices = frozenset(numpy.unique(self.subtemplate_indices))
+        # We do not need a valid subtemplate for the 0 index.
+        indices = frozenset(numpy.unique(self.subtemplate_indices)) - {0}
         if not indices.issubset(self.subtemplates.keys()):
             raise TQECException(
                 "Found an index in subtemplate_indices that does "


Test patch:
diff --git a/src/tqec/circuit/circuit_test.py b/src/tqec/circuit/circuit_test.py
index 17037d03..e2391329 100644
--- a/src/tqec/circuit/circuit_test.py
+++ b/src/tqec/circuit/circuit_test.py
@@ -7,7 +7,7 @@ from tqec.circuit.circuit import generate_circuit
 from tqec.exceptions import TQECException
 from tqec.plaquette.enums import PlaquetteOrientation
 from tqec.plaquette.library import zz_memory_plaquette
-from tqec.plaquette.plaquette import Plaquette
+from tqec.plaquette.plaquette import Plaquette, Plaquettes
 from tqec.plaquette.qubit import PlaquetteQubits
 from tqec.templates._testing import FixedTemplate
 from tqec.templates.base import Template
@@ -33,24 +33,19 @@ def _expected_circuit(qubits: PlaquetteQubits) -> cirq.Circuit:
     )
 
 
-def test_generate_circuit_list(
-    plaquette: Plaquette, one_by_one_template: Template
-) -> None:
-    circuit = generate_circuit(one_by_one_template, [plaquette])
-    assert circuit == _expected_circuit(plaquette.qubits)
-
-
 def test_generate_circuit_dict(
     plaquette: Plaquette, one_by_one_template: Template
 ) -> None:
-    circuit = generate_circuit(one_by_one_template, {1: plaquette})
+    circuit = generate_circuit(one_by_one_template, Plaquettes({1: plaquette}))
     assert circuit == _expected_circuit(plaquette.qubits)
 
 
 def test_generate_circuit_defaultdict(
     plaquette: Plaquette, one_by_one_template: Template
 ) -> None:
-    circuit = generate_circuit(one_by_one_template, defaultdict(lambda: plaquette))
+    circuit = generate_circuit(
+        one_by_one_template, Plaquettes(defaultdict(lambda: plaquette))
+    )
     assert circuit == _expected_circuit(plaquette.qubits)
 
 
@@ -58,13 +53,13 @@ def test_generate_circuit_dict_0_indexed(
     plaquette: Plaquette, one_by_one_template: Template
 ) -> None:
     with pytest.raises(TQECException):
-        generate_circuit(one_by_one_template, {0: plaquette})
+        generate_circuit(one_by_one_template, Plaquettes({0: plaquette}))
 
 
 def test_generate_circuit_wrong_number_of_plaquettes(
     plaquette: Plaquette, one_by_one_template: Template
 ) -> None:
     with pytest.raises(TQECException):
-        generate_circuit(one_by_one_template, [plaquette, plaquette])
+        generate_circuit(one_by_one_template, Plaquettes({1: plaquette, 2: plaquette}))
     with pytest.raises(TQECException):
-        generate_circuit(one_by_one_template, [])
+        generate_circuit(one_by_one_template, Plaquettes({}))
diff --git a/src/tqec/templates/_testing.py b/src/tqec/templates/_testing.py
index eaed006f..1257c3e3 100644
--- a/src/tqec/templates/_testing.py
+++ b/src/tqec/templates/_testing.py
@@ -6,7 +6,7 @@ from typing_extensions import override
 
 from tqec.position import Displacement
 from tqec.templates.base import RectangularTemplate
-from tqec.templates.enums import TemplateOrientation, TemplateSide
+from tqec.templates.enums import TemplateSide
 from tqec.templates.scale import LinearFunction, Scalable2D