Add DetectorErrorModelArrays.without_detectors

src/qldpc/decoders/dems.py

@@ -171,6 +171,15 @@ def from_arrays(
         dem_arrays.suggested_decompositions = suggested_decompositions or {}
         return dem_arrays.simplified() if simplify else dem_arrays
 
+    def copy(self) -> DetectorErrorModelArrays:
+        """Return an independent copy of this DetectorErrorModelArrays."""
+        dem_arrays = object.__new__(DetectorErrorModelArrays)
+        dem_arrays.detector_flip_matrix = self.detector_flip_matrix.copy()
+        dem_arrays.observable_flip_matrix = self.observable_flip_matrix.copy()
+        dem_arrays.error_probs = self.error_probs.copy()
+        dem_arrays.suggested_decompositions = dict(self.suggested_decompositions)
+        return dem_arrays
+
     @property
     def num_errors(self) -> int:
         """The number of distinct circuit errors."""
@@ -308,7 +317,9 @@ def to_detector_error_model(self) -> stim.DetectorErrorModel:
                 obs_targets = [stim.DemTarget.logical_observable_id(oo) for oo in observables]
                 targets = det_targets + obs_targets
 
-            dem.append("error", prob, targets)
+            # use DemInstruction (rather than dem.append("error", prob, targets)) so that error
+            # mechanisms with no targets -- valid, deterministically silent errors -- can be emitted
+            dem.append(stim.DemInstruction("error", [prob], targets))
 
         return dem
 
@@ -339,15 +350,15 @@ def simplified(self) -> DetectorErrorModelArrays:
         """Simplify this DetectorErrorModelArrays object by merging errors."""
         return DetectorErrorModelArrays(self.to_detector_error_model(), simplify=True)
 
-    def without_untriggered_detectors(self) -> DetectorErrorModelArrays:
-        """Drop all detectors that are not triggered by any error mechanism.
+    def without_detectors(self, detectors: Collection[int]) -> DetectorErrorModelArrays:
+        """Drop the given detectors.  Drop or merge error mechanisms that re empty or redundant.
 
-        Such detectors are deterministically 0, so removing them changes no sampling outcome -- it
-        only shrinks the DEM.  Error mechanisms and observables are left untouched.
+        Remaining detectors get re-indexed to range(num_remaining_detectors).
         """
-        keep = self.detector_flip_matrix.getnnz(axis=1) > 0
+        keep = np.ones(self.num_detectors, dtype=bool)
+        keep[list(detectors)] = False
         if keep.all():
-            return self
+            return self.copy()
         old_to_new = np.cumsum(keep) - 1
         suggested_decompositions = {
             error_index: _remap_decomposition_detectors(components, keep, old_to_new)
@@ -358,9 +369,17 @@ def without_untriggered_detectors(self) -> DetectorErrorModelArrays:
             self.observable_flip_matrix,
             self.error_probs,
             suggested_decompositions,
-            simplify=False,
+            simplify=True,
         )
 
+    def without_untriggered_detectors(self) -> DetectorErrorModelArrays:
+        """Drop all detectors that are not triggered by any error mechanism.
+
+        Such detectors are deterministically 0, so removing them changes no sampling outcome -- it
+        only shrinks the DEM.  Error mechanisms and observables are left untouched.
+        """
+        return self.without_detectors(np.flatnonzero(self.detector_flip_matrix.getnnz(axis=1) == 0))
+
     def with_decomposed_errors(self, *, simplify: bool = True) -> DetectorErrorModelArrays:
         """Split error mechanisms according to their suggested decompositions.
 

src/qldpc/decoders/dems_test.py

@@ -248,9 +248,11 @@ def test_post_selection() -> None:
 
 def test_without_untriggered_detectors() -> None:
     """Drop detectors that no error mechanism triggers."""
-    # with no dead detectors, the object is returned unchanged
+    # with no dead detectors, an equivalent (independent) copy is returned
     dem_arrays = decoders.DetectorErrorModelArrays(stim.DetectorErrorModel("error(0.1) D0"))
-    assert dem_arrays.without_untriggered_detectors() is dem_arrays
+    result = dem_arrays.without_untriggered_detectors()
+    assert result is not dem_arrays
+    assert result.to_dem() == dem_arrays.to_dem()
 
     # D1 is dead but appears in the decomposition D0 D1 ^ D2 D1, where it cancels; it must be
     # filtered out rather than remapped to a stale index (live detectors: D0 -> 0, D2 -> 1)