From a9b30c94505cb9228f97f7f617fe80998de332dd Mon Sep 17 00:00:00 2001
From: acasta-yhliu <yhliu2000@outlook.com>
Date: Fri, 27 Mar 2026 14:00:12 -0400
Subject: [PATCH 01/13] [update] start implementation of the TargetedNoiseRule
 #426

TargetedNoiseRule class is appended to the code, and implmented is_targeted_noisy_op
---
 src/qldpc/circuits/noise_model.py | 79 +++++++++++++++++++++++++++++++
 1 file changed, 79 insertions(+)

diff --git a/src/qldpc/circuits/noise_model.py b/src/qldpc/circuits/noise_model.py
index 70062257a..b8ae22324 100644
--- a/src/qldpc/circuits/noise_model.py
+++ b/src/qldpc/circuits/noise_model.py
@@ -292,6 +292,85 @@ def noisy_operation(
         return noisy_op, noise_after
 
 
+class TargetedNoiseRule(NoiseRule):
+    """Describes how to add noise to a specific circuit instruction on specific qubits.
+
+    Unlike NoiseRule, which applies to all operations of a given type, this rule matches only an
+    exact gate-and-target combination, allowing fine-grained per-operation noise overrides.
+    """
+
+    def __init__(
+        self,
+        *,
+        noisy_op: stim.CircuitInstruction,
+        tags: Collection[str] | None = None,
+        noise: stim.Circuit = stim.Circuit(),
+        readout_error: float = 0,
+        reset_error: float = 0,
+    ):
+        """Initializes a targeted noise rule for a specific circuit instruction.
+
+        Args:
+            noisy_op: The circuit instruction that this rule targets.  Defines the gate name and
+                qubit targets to match against.  Gate args on this instruction are ignored during
+                matching.
+            tags: If not None, only match operations whose tag exactly matches one of the given
+                strings.  If None, match operations regardless of their tag.
+            noise: An explicit noise circuit to append after the matched operation.  Defaults to an
+                empty circuit (no noise).
+            readout_error: The probability that a measurement result is reported incorrectly.  Only
+                allowed when noisy_op is a measurement.
+            reset_error: The probability that a qubit is reset to the wrong state.  Only allowed
+                when noisy_op is a reset.
+
+        Raises:
+            ValueError: If readout_error or reset_error is not between 0 and 1 (inclusive).
+        """
+        super().__init__(readout_error=readout_error, reset_error=reset_error)
+        self.noisy_op = noisy_op
+        self.tags: frozenset[str] | None = frozenset(tags) if tags is not None else None
+        self.noise = noise
+
+    def is_targeted_noisy_op(self, op: stim.CircuitInstruction) -> bool:
+        """Returns whether the given operation matches this rule's target instruction.
+
+        Two operations match if they have the same gate name and the same qubit target values in the
+        same order.  Gate args are ignored.  If self.tags is not None, op.tag must exactly match one
+        of the given tags.
+
+        Args:
+            op: The circuit instruction to check.
+
+        Returns:
+            True if op matches this rule's target instruction.  False otherwise.
+        """
+        if op.name != self.noisy_op.name:
+            return False
+        if self.tags is not None and op.tag not in self.tags:
+            return False
+        return op.targets_copy() == self.noisy_op.targets_copy()
+
+    def noisy_operation(
+        self,
+        op: stim.CircuitInstruction,
+        *,
+        immune_qubits: set[int] = set(),
+    ) -> tuple[stim.CircuitInstruction, stim.Circuit]:
+        """Apply this targeted noise rule to the given operation.
+
+        Args:
+            op: The operation to add noise to.
+            immune_qubits: Set of qubit indices that should not have noise applied to them.  If any
+                target qubit of the matched operation is immune, no noise is emitted.
+
+        Returns:
+            stim.CircuitInstruction: The given operation, possibly modified to account for readout
+                or reset errors.
+            stim.Circuit: Noise operations that should follow the given operation.
+        """
+        ...
+
+
 class NoiseModel:
     """A model that defines how to add noise to quantum circuits.
 

From b7fd15ce0b60c0e14b9571ee243a6c8ba8824ff1 Mon Sep 17 00:00:00 2001
From: acasta-yhliu <yhliu2000@outlook.com>
Date: Fri, 27 Mar 2026 14:07:10 -0400
Subject: [PATCH 02/13] [update] partially implement
 `TargetedNoiseRule.noisy_operation` #426

The implementation follows `NoiseRule.noisy_operation` for now
---
 src/qldpc/circuits/noise_model.py | 39 ++++++++++++++++++++++++++++++-
 1 file changed, 38 insertions(+), 1 deletion(-)

diff --git a/src/qldpc/circuits/noise_model.py b/src/qldpc/circuits/noise_model.py
index b8ae22324..e4149b0fb 100644
--- a/src/qldpc/circuits/noise_model.py
+++ b/src/qldpc/circuits/noise_model.py
@@ -368,7 +368,44 @@ def noisy_operation(
                 or reset errors.
             stim.Circuit: Noise operations that should follow the given operation.
         """
-        ...
+        if not self.is_targeted_noisy_op(op):
+            return op, stim.Circuit()
+
+        targets = op.targets_copy()
+        if immune_qubits and any(
+            (
+                target.is_qubit_target
+                or target.is_x_target
+                or target.is_y_target
+                or target.is_z_target
+            )
+            and target.value in immune_qubits
+            for target in targets
+        ):
+            return op, stim.Circuit()
+
+        args = op.gate_args_copy()
+        if self.readout_error:
+            assert op.name in JUST_MEASURE_OPS or op.name in MEASURE_AND_RESET_OPS
+            if not args:
+                args = [self.readout_error]
+            else:
+                assert len(args) == 1
+                # combine bit-flip probabilities
+                args = [1 - (1 - self.readout_error) * (1 - args[0])]
+
+        noisy_op = stim.CircuitInstruction(op.name, targets, args, tag=op.tag)
+        noise_after = stim.Circuit()
+
+        if self.reset_error:
+            assert op.name in JUST_RESET_OPS or op.name in MEASURE_AND_RESET_OPS
+            qubit_targets = [target.value for target in targets if not target.is_combiner]
+            error_name = ("X" if _get_standardized_name(op)[-1] != "X" else "Z") + "_ERROR"
+            noise_after.append(stim.CircuitInstruction(error_name, qubit_targets, [self.reset_error]))
+
+        noise_after += self.noise
+
+        return noisy_op, noise_after
 
 
 class NoiseModel:

From 192e88158f714179482dc7598f15889cec6c3c72 Mon Sep 17 00:00:00 2001
From: "Michael A. Perlin" <michael.perlin@jpmchase.com>
Date: Sat, 28 Mar 2026 22:50:59 -0400
Subject: [PATCH 03/13] read-through pass

---
 src/qldpc/circuits/noise_model.py | 42 +++++++++++++++++--------------
 1 file changed, 23 insertions(+), 19 deletions(-)

diff --git a/src/qldpc/circuits/noise_model.py b/src/qldpc/circuits/noise_model.py
index e4149b0fb..90cece01e 100644
--- a/src/qldpc/circuits/noise_model.py
+++ b/src/qldpc/circuits/noise_model.py
@@ -198,7 +198,7 @@ def __init__(
         after: dict[str, float | Iterable[float]] = {},
         readout_error: float = 0,
         reset_error: float = 0,
-    ):
+    ) -> None:
         """Initializes a noise rule with specified error channels.
 
         Args:
@@ -303,40 +303,45 @@ def __init__(
         self,
         *,
         noisy_op: stim.CircuitInstruction,
-        tags: Collection[str] | None = None,
-        noise: stim.Circuit = stim.Circuit(),
+        noise: stim.Circuit,
         readout_error: float = 0,
         reset_error: float = 0,
-    ):
+        tags: Collection[str] | None = None,
+    ) -> None:
         """Initializes a targeted noise rule for a specific circuit instruction.
 
         Args:
             noisy_op: The circuit instruction that this rule targets.  Defines the gate name and
                 qubit targets to match against.  Gate args on this instruction are ignored during
                 matching.
+            noise: An explicit noise circuit to append after the matched operation.
+            readout_error: The probability that a measurement result is reported incorrectly.  Only
+                allowed for operations that produce measurement results.
+            reset_error: The probability that a qubit is reset to the wrong state.  Only allowed for
+                operations that reset qubits.
             tags: If not None, only match operations whose tag exactly matches one of the given
                 strings.  If None, match operations regardless of their tag.
-            noise: An explicit noise circuit to append after the matched operation.  Defaults to an
-                empty circuit (no noise).
-            readout_error: The probability that a measurement result is reported incorrectly.  Only
-                allowed when noisy_op is a measurement.
-            reset_error: The probability that a qubit is reset to the wrong state.  Only allowed
-                when noisy_op is a reset.
 
         Raises:
             ValueError: If readout_error or reset_error is not between 0 and 1 (inclusive).
         """
         super().__init__(readout_error=readout_error, reset_error=reset_error)
         self.noisy_op = noisy_op
-        self.tags: frozenset[str] | None = frozenset(tags) if tags is not None else None
         self.noise = noise
+        self.tags: frozenset[str] | None = frozenset(tags) if tags is not None else None
+
+    def __bool__(self) -> bool:
+        """Is this noise rule nontrivial?"""
+        nontrivial_noise = bool(self.noise) or bool(self.readout_error) or bool(self.reset_error)
+        nontrivial_targets = bool(self.noisy_op.targets_copy())
+        return nontrivial_noise and nontrivial_targets
 
     def is_targeted_noisy_op(self, op: stim.CircuitInstruction) -> bool:
-        """Returns whether the given operation matches this rule's target instruction.
+        """Determine whether the given operation matches this rule's target instruction.
 
         Two operations match if they have the same gate name and the same qubit target values in the
-        same order.  Gate args are ignored.  If self.tags is not None, op.tag must exactly match one
-        of the given tags.
+        same order.  Gate args are ignored.  If self.tags is not None, op.tag must exactly match an
+        element of self.tags.
 
         Args:
             op: The circuit instruction to check.
@@ -344,9 +349,7 @@ def is_targeted_noisy_op(self, op: stim.CircuitInstruction) -> bool:
         Returns:
             True if op matches this rule's target instruction.  False otherwise.
         """
-        if op.name != self.noisy_op.name:
-            return False
-        if self.tags is not None and op.tag not in self.tags:
+        if op.name != self.noisy_op.name or (self.tags is not None and op.tag not in self.tags):
             return False
         return op.targets_copy() == self.noisy_op.targets_copy()
 
@@ -401,7 +404,8 @@ def noisy_operation(
             assert op.name in JUST_RESET_OPS or op.name in MEASURE_AND_RESET_OPS
             qubit_targets = [target.value for target in targets if not target.is_combiner]
             error_name = ("X" if _get_standardized_name(op)[-1] != "X" else "Z") + "_ERROR"
-            noise_after.append(stim.CircuitInstruction(error_name, qubit_targets, [self.reset_error]))
+            error_op = stim.CircuitInstruction(error_name, qubit_targets, [self.reset_error])
+            noise_after.append(error_op)
 
         noise_after += self.noise
 
@@ -426,7 +430,7 @@ def __init__(
         idle_error: float | None = None,
         additional_error_waiting_for_m_or_r: float | None = None,
         rules: dict[str, NoiseRule] | None = None,
-    ):
+    ) -> None:
         """Initializes a noise model with specified parameters.
 
         Args:

From 1b8a5248ba8ce9f59d143a000c203c6fc88549bc Mon Sep 17 00:00:00 2001
From: "Michael A. Perlin" <michael.perlin@jpmchase.com>
Date: Mon, 30 Mar 2026 09:36:05 -0400
Subject: [PATCH 04/13] use NoiseRule.noisy_operation inside
 TargetedNoiseRule.noisy_operation

---
 src/qldpc/circuits/noise_model.py | 71 ++++++++-----------------------
 1 file changed, 17 insertions(+), 54 deletions(-)

diff --git a/src/qldpc/circuits/noise_model.py b/src/qldpc/circuits/noise_model.py
index 90cece01e..9148fd382 100644
--- a/src/qldpc/circuits/noise_model.py
+++ b/src/qldpc/circuits/noise_model.py
@@ -48,6 +48,7 @@
 from __future__ import annotations
 
 import collections
+import dataclasses
 from collections.abc import Collection, Iterable, Iterator
 
 import stim
@@ -192,6 +193,10 @@ class NoiseRule:
     applied to a particular type of quantum operation.
     """
 
+    after: dict[str, float | Iterable[float]] = dataclasses.field(default_factory=dict)
+    readout_error: float
+    reset_error: float
+
     def __init__(
         self,
         *,
@@ -241,7 +246,7 @@ def __bool__(self) -> bool:
         return bool(self.after) or bool(self.readout_error) or bool(self.reset_error)
 
     def noisy_operation(
-        self, op: stim.CircuitInstruction, *, immune_qubits: set[int] = set()
+        self, op: stim.CircuitInstruction
     ) -> tuple[stim.CircuitInstruction, stim.Circuit]:
         """Apply this noise rule to the given operation.
 
@@ -253,26 +258,15 @@ def noisy_operation(
             stim.Circuit: Noise operations that should follow the given operation.
         """
         targets = op.targets_copy()
-        if immune_qubits and any(
-            (
-                target.is_qubit_target
-                or target.is_x_target
-                or target.is_y_target
-                or target.is_z_target
-            )
-            and target.value in immune_qubits
-            for target in targets
-        ):
-            return op, stim.Circuit()
-
         args = op.gate_args_copy()
+
         if self.readout_error:
             assert op.name in JUST_MEASURE_OPS or op.name in MEASURE_AND_RESET_OPS
             if not args:
                 args = [self.readout_error]
             else:
                 assert len(args) == 1
-                # combine bit-flip probabilities
+                # combine the bit-flip probabilities self.readout_error and args[0]
                 args = [1 - (1 - self.readout_error) * (1 - args[0])]
 
         noisy_op = stim.CircuitInstruction(op.name, targets, args, tag=op.tag)
@@ -299,6 +293,10 @@ class TargetedNoiseRule(NoiseRule):
     exact gate-and-target combination, allowing fine-grained per-operation noise overrides.
     """
 
+    noisy_op: stim.CircuitInstruction
+    noise: stim.Circuit
+    tags: frozenset[str] | None
+
     def __init__(
         self,
         *,
@@ -354,10 +352,7 @@ def is_targeted_noisy_op(self, op: stim.CircuitInstruction) -> bool:
         return op.targets_copy() == self.noisy_op.targets_copy()
 
     def noisy_operation(
-        self,
-        op: stim.CircuitInstruction,
-        *,
-        immune_qubits: set[int] = set(),
+        self, op: stim.CircuitInstruction
     ) -> tuple[stim.CircuitInstruction, stim.Circuit]:
         """Apply this targeted noise rule to the given operation.
 
@@ -373,42 +368,8 @@ def noisy_operation(
         """
         if not self.is_targeted_noisy_op(op):
             return op, stim.Circuit()
-
-        targets = op.targets_copy()
-        if immune_qubits and any(
-            (
-                target.is_qubit_target
-                or target.is_x_target
-                or target.is_y_target
-                or target.is_z_target
-            )
-            and target.value in immune_qubits
-            for target in targets
-        ):
-            return op, stim.Circuit()
-
-        args = op.gate_args_copy()
-        if self.readout_error:
-            assert op.name in JUST_MEASURE_OPS or op.name in MEASURE_AND_RESET_OPS
-            if not args:
-                args = [self.readout_error]
-            else:
-                assert len(args) == 1
-                # combine bit-flip probabilities
-                args = [1 - (1 - self.readout_error) * (1 - args[0])]
-
-        noisy_op = stim.CircuitInstruction(op.name, targets, args, tag=op.tag)
-        noise_after = stim.Circuit()
-
-        if self.reset_error:
-            assert op.name in JUST_RESET_OPS or op.name in MEASURE_AND_RESET_OPS
-            qubit_targets = [target.value for target in targets if not target.is_combiner]
-            error_name = ("X" if _get_standardized_name(op)[-1] != "X" else "Z") + "_ERROR"
-            error_op = stim.CircuitInstruction(error_name, qubit_targets, [self.reset_error])
-            noise_after.append(error_op)
-
+        noisy_op, noise_after = super().noisy_operation(op)
         noise_after += self.noise
-
         return noisy_op, noise_after
 
 
@@ -609,10 +570,12 @@ def _inplace_append_noisy_moment(
             if immune_op_tag in op.tag or (rule := self.get_noise_rule(op)) is None:
                 circuit.append(op)
             else:
-                noisy_op, after = rule.noisy_operation(op, immune_qubits=immune_qubits)
+                noisy_op, after = rule.noisy_operation(op)
                 circuit.append(noisy_op)
                 noise_after_moment += after
 
+        # TODO: deal with immmune_qubits
+
         circuit += noise_after_moment
 
         moment_was_noisy = any(immune_op_tag not in op.tag for op in moment)

From c33ca429c193bc2bac4c691065287ae2b03cffe7 Mon Sep 17 00:00:00 2001
From: "Michael A. Perlin" <michael.perlin@jpmchase.com>
Date: Mon, 30 Mar 2026 09:36:33 -0400
Subject: [PATCH 05/13] remove class variables

---
 src/qldpc/circuits/noise_model.py | 8 --------
 1 file changed, 8 deletions(-)

diff --git a/src/qldpc/circuits/noise_model.py b/src/qldpc/circuits/noise_model.py
index 9148fd382..a97244d2d 100644
--- a/src/qldpc/circuits/noise_model.py
+++ b/src/qldpc/circuits/noise_model.py
@@ -193,10 +193,6 @@ class NoiseRule:
     applied to a particular type of quantum operation.
     """
 
-    after: dict[str, float | Iterable[float]] = dataclasses.field(default_factory=dict)
-    readout_error: float
-    reset_error: float
-
     def __init__(
         self,
         *,
@@ -293,10 +289,6 @@ class TargetedNoiseRule(NoiseRule):
     exact gate-and-target combination, allowing fine-grained per-operation noise overrides.
     """
 
-    noisy_op: stim.CircuitInstruction
-    noise: stim.Circuit
-    tags: frozenset[str] | None
-
     def __init__(
         self,
         *,

From d29fe994472dd4788a2f8af37c5bf1ada028b917 Mon Sep 17 00:00:00 2001
From: "Michael A. Perlin" <michael.perlin@jpmchase.com>
Date: Mon, 30 Mar 2026 09:41:07 -0400
Subject: [PATCH 06/13] remove unused import

---
 src/qldpc/circuits/noise_model.py | 1 -
 1 file changed, 1 deletion(-)

diff --git a/src/qldpc/circuits/noise_model.py b/src/qldpc/circuits/noise_model.py
index a97244d2d..5bd95b35a 100644
--- a/src/qldpc/circuits/noise_model.py
+++ b/src/qldpc/circuits/noise_model.py
@@ -48,7 +48,6 @@
 from __future__ import annotations
 
 import collections
-import dataclasses
 from collections.abc import Collection, Iterable, Iterator
 
 import stim

From b99929ace025fc58d0e4534127260dc4f019f733 Mon Sep 17 00:00:00 2001
From: "Michael A. Perlin" <michael.perlin@jpmchase.com>
Date: Mon, 30 Mar 2026 09:43:56 -0400
Subject: [PATCH 07/13] update TODO

---
 src/qldpc/circuits/noise_model.py | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/src/qldpc/circuits/noise_model.py b/src/qldpc/circuits/noise_model.py
index 5bd95b35a..caf487fe7 100644
--- a/src/qldpc/circuits/noise_model.py
+++ b/src/qldpc/circuits/noise_model.py
@@ -565,7 +565,10 @@ def _inplace_append_noisy_moment(
                 circuit.append(noisy_op)
                 noise_after_moment += after
 
-        # TODO: deal with immmune_qubits
+        # TODO: post-process noise_after_moment using immune_qubits
+        # see comments:
+        # - https://github.com/qLDPCOrg/qLDPC/issues/426#issuecomment-4144385907
+        # - https://github.com/qLDPCOrg/qLDPC/issues/426#issuecomment-4149324799
 
         circuit += noise_after_moment
 

From d0f6f207fb37b6088417b214f24f485b0cf932cd Mon Sep 17 00:00:00 2001
From: acasta-yhliu <yhliu2000@outlook.com>
Date: Mon, 30 Mar 2026 10:37:16 -0400
Subject: [PATCH 08/13] [update] remove argument "immune_qubits" in both
 `NoiseRule` and `TargetedNoiseRule`

An `immunize_qubits` function is appended and called in `_inplace_append_noisy_moment`. Detailed implementation is remain incomplete.
---
 src/qldpc/circuits/noise_model.py | 93 ++++++++++++++++++++++++++-----
 1 file changed, 79 insertions(+), 14 deletions(-)

diff --git a/src/qldpc/circuits/noise_model.py b/src/qldpc/circuits/noise_model.py
index caf487fe7..13e730628 100644
--- a/src/qldpc/circuits/noise_model.py
+++ b/src/qldpc/circuits/noise_model.py
@@ -185,6 +185,39 @@ def as_noiseless_circuit(circuit: stim.Circuit) -> stim.Circuit:
     return noiseless_circuit
 
 
+def immunize_noise(noise: stim.Circuit, immune_qubits: set[int]) -> stim.Circuit:
+    """Return a copy of a flat noise circuit with instructions targeting immune qubits removed.
+
+    An instruction is removed if any of its qubit targets belongs to immune_qubits.
+
+    Args:
+        noise: A flat noise circuit (no repeat blocks) to filter.
+        immune_qubits: Set of qubit indices that should not have noise applied to them.
+
+    Returns:
+        stim.Circuit: A filtered copy of the input circuit.
+    """
+    if not immune_qubits:
+        return noise
+    result = stim.Circuit()
+    for instruction in noise:
+        assert isinstance(instruction, stim.CircuitInstruction)
+        if not any(
+            (
+                target.is_qubit_target
+                or target.is_x_target
+                or target.is_y_target
+                or target.is_z_target
+            )
+            and target.value in immune_qubits
+            for target in instruction.targets_copy()
+            if not target.is_combiner
+        ):
+            result.append(instruction)
+    return result
+
+
+
 class NoiseRule:
     """Describes how to add noise to an operation.
 
@@ -242,9 +275,13 @@ def __bool__(self) -> bool:
 
     def noisy_operation(
         self, op: stim.CircuitInstruction
+        self, op: stim.CircuitInstruction
     ) -> tuple[stim.CircuitInstruction, stim.Circuit]:
         """Apply this noise rule to the given operation.
 
+        Immunity to noise is not handled here.  The caller is responsible for filtering the
+        returned noise circuit using immunize_noise if needed.
+
         Args:
             op: The operation to add noise to.
 
@@ -267,19 +304,34 @@ def noisy_operation(
         noisy_op = stim.CircuitInstruction(op.name, targets, args, tag=op.tag)
         noise_after = stim.Circuit()
 
-        qubit_targets = [target.value for target in targets if not target.is_combiner]
         if self.reset_error:
             assert op.name in JUST_RESET_OPS or op.name in MEASURE_AND_RESET_OPS
+            qubit_targets = [target.value for target in targets if not target.is_combiner]
             error_name = ("X" if _get_standardized_name(op)[-1] != "X" else "Z") + "_ERROR"
-            error_op = stim.CircuitInstruction(error_name, qubit_targets, [self.reset_error])
-            noise_after.append(error_op)
+            noise_after.append(stim.CircuitInstruction(error_name, qubit_targets, [self.reset_error]))
 
-        for op_name, args in self.after.items():
-            error_op = stim.CircuitInstruction(op_name, qubit_targets, args)
-            noise_after.append(error_op)
+        noise_after += self._build_noise_after(op)
 
         return noisy_op, noise_after
 
+    def _build_noise_after(self, op: stim.CircuitInstruction) -> stim.Circuit:
+        """Build the extra noise circuit to append after the given operation.
+
+        Subclasses may override this to customize the noise that follows an operation.  Reset errors
+        are excluded here and handled separately in noisy_operation.
+
+        Args:
+            op: The operation being applied.
+
+        Returns:
+            stim.Circuit: Additional noise to append after the operation.
+        """
+        qubit_targets = [target.value for target in op.targets_copy() if not target.is_combiner]
+        noise = stim.Circuit()
+        for op_name, args in self.after.items():
+            noise.append(stim.CircuitInstruction(op_name, qubit_targets, args))
+        return noise
+
 
 class TargetedNoiseRule(NoiseRule):
     """Describes how to add noise to a specific circuit instruction on specific qubits.
@@ -343,14 +395,16 @@ def is_targeted_noisy_op(self, op: stim.CircuitInstruction) -> bool:
         return op.targets_copy() == self.noisy_op.targets_copy()
 
     def noisy_operation(
-        self, op: stim.CircuitInstruction
+        self,
+        op: stim.CircuitInstruction,
     ) -> tuple[stim.CircuitInstruction, stim.Circuit]:
         """Apply this targeted noise rule to the given operation.
 
+        Immunity to noise is not handled here.  The caller is responsible for filtering the
+        returned noise circuit using immunize_noise if needed.
+
         Args:
             op: The operation to add noise to.
-            immune_qubits: Set of qubit indices that should not have noise applied to them.  If any
-                target qubit of the matched operation is immune, no noise is emitted.
 
         Returns:
             stim.CircuitInstruction: The given operation, possibly modified to account for readout
@@ -359,9 +413,10 @@ def noisy_operation(
         """
         if not self.is_targeted_noisy_op(op):
             return op, stim.Circuit()
-        noisy_op, noise_after = super().noisy_operation(op)
-        noise_after += self.noise
-        return noisy_op, noise_after
+        return super().noisy_operation(op)
+
+    def _build_noise_after(self, op: stim.CircuitInstruction) -> stim.Circuit:
+        return self.noise
 
 
 class NoiseModel:
@@ -562,8 +617,18 @@ def _inplace_append_noisy_moment(
                 circuit.append(op)
             else:
                 noisy_op, after = rule.noisy_operation(op)
-                circuit.append(noisy_op)
-                noise_after_moment += after
+                op_targets_immune = immune_qubits and any(
+                    (
+                        target.is_qubit_target
+                        or target.is_x_target
+                        or target.is_y_target
+                        or target.is_z_target
+                    )
+                    and target.value in immune_qubits
+                    for target in op.targets_copy()
+                )
+                circuit.append(op if op_targets_immune else noisy_op)
+                noise_after_moment += immunize_noise(after, immune_qubits)
 
         # TODO: post-process noise_after_moment using immune_qubits
         # see comments:

From af6a194e0cae42c2496edde8926efe3422bbeb59 Mon Sep 17 00:00:00 2001
From: acasta-yhliu <yhliu2000@outlook.com>
Date: Mon, 30 Mar 2026 11:27:06 -0400
Subject: [PATCH 09/13] [fix] fixed coverage and format issue

tests are added to noise_model_test.py for TargetedNoiseRule, the behavior of immunize_qubits is still remain incomplete.
---
 src/qldpc/circuits/__init__.py         |   2 +
 src/qldpc/circuits/noise_model.py      |   5 +-
 src/qldpc/circuits/noise_model_test.py | 117 +++++++++++++++++++++++++
 3 files changed, 122 insertions(+), 2 deletions(-)

diff --git a/src/qldpc/circuits/__init__.py b/src/qldpc/circuits/__init__.py
index a4e82fd3b..4b470342e 100644
--- a/src/qldpc/circuits/__init__.py
+++ b/src/qldpc/circuits/__init__.py
@@ -25,6 +25,7 @@
     NoiseModel,
     NoiseRule,
     SI1000NoiseModel,
+    TargetedNoiseRule,
     as_noiseless_circuit,
 )
 from .syndrome_measurement import (
@@ -59,6 +60,7 @@
     "NoiseModel",
     "NoiseRule",
     "SI1000NoiseModel",
+    "TargetedNoiseRule",
     "as_noiseless_circuit",
     "EdgeColoring",
     "EdgeColoringXZ",
diff --git a/src/qldpc/circuits/noise_model.py b/src/qldpc/circuits/noise_model.py
index 13e730628..2b8e14a94 100644
--- a/src/qldpc/circuits/noise_model.py
+++ b/src/qldpc/circuits/noise_model.py
@@ -217,7 +217,6 @@ def immunize_noise(noise: stim.Circuit, immune_qubits: set[int]) -> stim.Circuit
     return result
 
 
-
 class NoiseRule:
     """Describes how to add noise to an operation.
 
@@ -308,7 +307,9 @@ def noisy_operation(
             assert op.name in JUST_RESET_OPS or op.name in MEASURE_AND_RESET_OPS
             qubit_targets = [target.value for target in targets if not target.is_combiner]
             error_name = ("X" if _get_standardized_name(op)[-1] != "X" else "Z") + "_ERROR"
-            noise_after.append(stim.CircuitInstruction(error_name, qubit_targets, [self.reset_error]))
+            noise_after.append(
+                stim.CircuitInstruction(error_name, qubit_targets, [self.reset_error])
+            )
 
         noise_after += self._build_noise_after(op)
 
diff --git a/src/qldpc/circuits/noise_model_test.py b/src/qldpc/circuits/noise_model_test.py
index a430a3f16..68bcad775 100644
--- a/src/qldpc/circuits/noise_model_test.py
+++ b/src/qldpc/circuits/noise_model_test.py
@@ -103,6 +103,123 @@ def test_gate_errors() -> None:
         noise_model.noisy_circuit(circuit, insert_ticks=False)
 
 
+def test_targeted_gate_errors() -> None:
+    """TargetedNoiseRule applies noise only to exact gate-and-qubit matches."""
+
+    # only CX 0 1 gets depolarizing noise; CX 1 2 is unaffected
+    circuit = stim.Circuit("""
+        CX 0 1 1 2
+        TICK
+        M 0
+        RX 1
+        MR 2
+    """)
+
+    targeted_cx = circuits.TargetedNoiseRule(
+        noisy_op=stim.CircuitInstruction("CX", [0, 1]),
+        noise=stim.Circuit("DEPOLARIZE2(0.2) 0 1"),
+    )
+    noise_model = circuits.NoiseModel(clifford_2q_error=targeted_cx)
+    assert _circuits_are_equivalent(
+        stim.Circuit("""
+        CX 0 1
+        DEPOLARIZE2(0.2) 0 1
+        CX 1 2
+        TICK
+        M 0
+        RX 1
+        MR 2
+    """),
+        noise_model.noisy_circuit(circuit),
+    )
+
+    # only CX 0 1 gets depolarizing noise but on 1 2
+    targeted_cx = circuits.TargetedNoiseRule(
+        noisy_op=stim.CircuitInstruction("CX", [0, 1]),
+        noise=stim.Circuit("DEPOLARIZE2(0.2) 1 2"),
+    )
+    noise_model = circuits.NoiseModel(clifford_2q_error=targeted_cx)
+    assert _circuits_are_equivalent(
+        stim.Circuit("""
+        CX 0 1
+        DEPOLARIZE2(0.2) 1 2
+        CX 1 2
+        TICK
+        M 0
+        RX 1
+        MR 2
+    """),
+        noise_model.noisy_circuit(circuit),
+    )
+
+    # targeted readout and reset errors on specific qubits
+    circuit = stim.Circuit("""
+        H 0
+        CX 0 1 1 2
+        TICK
+        M 0
+        RX 1
+        MR 2
+    """)
+    noise_model = circuits.NoiseModel(
+        clifford_1q_error=0.1,
+        clifford_2q_error=0.2,
+        rules={
+            "MZ": circuits.TargetedNoiseRule(
+                noisy_op=stim.CircuitInstruction("MZ", [0]), readout_error=0.3
+            ),
+            "RX": circuits.TargetedNoiseRule(
+                noisy_op=stim.CircuitInstruction("RX", [1]), reset_error=0.4
+            ),
+        },
+    )
+    assert _circuits_are_equivalent(
+        stim.Circuit("""
+        H 0
+        DEPOLARIZE1(0.1) 0
+        CX 0 1
+        DEPOLARIZE2(0.2) 0 1
+        CX 1 2
+        DEPOLARIZE2(0.2) 1 2
+        TICK
+        MZ(0.3) 0
+        RX 1
+        MR 2
+        Z_ERROR(0.4) 1
+    """),
+        noise_model.noisy_circuit(circuit),
+    )
+
+    # # tags: only apply noise to operations with the matching tag
+    circuit = stim.Circuit("""
+        H[ancilla] 0
+        H 1
+        S 0
+        TICK
+        CX 0 1
+    """)
+    noise_model = circuits.NoiseModel(
+        clifford_2q_error=0.2,
+        clifford_1q_error=circuits.TargetedNoiseRule(
+            noisy_op=stim.CircuitInstruction("H", [0]),
+            tags={"ancilla"},
+            noise=stim.Circuit("DEPOLARIZE1(0.1) 0"),
+        ),
+    )
+    noisy_circuit = stim.Circuit("""
+        H[ancilla] 0
+        H 1
+        DEPOLARIZE1(0.1) 0
+        S 0
+        TICK
+        CX 0 1
+        DEPOLARIZE2(0.2) 0 1
+    """)
+    print("=====Expect", noisy_circuit)
+    print("=====But got", noise_model.noisy_circuit(circuit))
+    assert _circuits_are_equivalent(noisy_circuit, noise_model.noisy_circuit(circuit))
+
+
 def test_idle_errors() -> None:
     """Add idling errors to a circuit."""
 

From 383b7955183c82ff8b9441edd55a48331e0d533f Mon Sep 17 00:00:00 2001
From: acasta-yhliu <yhliu2000@outlook.com>
Date: Mon, 30 Mar 2026 11:36:20 -0400
Subject: [PATCH 10/13] [fix] duplicate argument and invalid syntax during
 merging

---
 src/qldpc/circuits/noise_model.py      | 1 -
 src/qldpc/circuits/noise_model_test.py | 4 ++--
 2 files changed, 2 insertions(+), 3 deletions(-)

diff --git a/src/qldpc/circuits/noise_model.py b/src/qldpc/circuits/noise_model.py
index 2b8e14a94..26cd58139 100644
--- a/src/qldpc/circuits/noise_model.py
+++ b/src/qldpc/circuits/noise_model.py
@@ -274,7 +274,6 @@ def __bool__(self) -> bool:
 
     def noisy_operation(
         self, op: stim.CircuitInstruction
-        self, op: stim.CircuitInstruction
     ) -> tuple[stim.CircuitInstruction, stim.Circuit]:
         """Apply this noise rule to the given operation.
 
diff --git a/src/qldpc/circuits/noise_model_test.py b/src/qldpc/circuits/noise_model_test.py
index 68bcad775..0a96655e6 100644
--- a/src/qldpc/circuits/noise_model_test.py
+++ b/src/qldpc/circuits/noise_model_test.py
@@ -166,10 +166,10 @@ def test_targeted_gate_errors() -> None:
         clifford_2q_error=0.2,
         rules={
             "MZ": circuits.TargetedNoiseRule(
-                noisy_op=stim.CircuitInstruction("MZ", [0]), readout_error=0.3
+                noisy_op=stim.CircuitInstruction("MZ", [0]), noise=stim.Circuit(), readout_error=0.3
             ),
             "RX": circuits.TargetedNoiseRule(
-                noisy_op=stim.CircuitInstruction("RX", [1]), reset_error=0.4
+                noisy_op=stim.CircuitInstruction("RX", [1]), noise=stim.Circuit(), reset_error=0.4
             ),
         },
     )

From 110ae058d266b167ff2a40918a6c2136452984a4 Mon Sep 17 00:00:00 2001
From: acasta-yhliu <yhliu2000@outlook.com>
Date: Tue, 31 Mar 2026 10:22:31 -0400
Subject: [PATCH 11/13] [update] implemented immunize_noise to support both
 DEPOLARIZE2 and PAULI_CHANNEL_2

tests are modified to cover new written code, need further refactor to make the code cleaner
---
 src/qldpc/circuits/noise_model.py      | 129 ++++++++++++++++++-------
 src/qldpc/circuits/noise_model_test.py |  69 +++++++++++++
 2 files changed, 164 insertions(+), 34 deletions(-)

diff --git a/src/qldpc/circuits/noise_model.py b/src/qldpc/circuits/noise_model.py
index 26cd58139..89a4623f7 100644
--- a/src/qldpc/circuits/noise_model.py
+++ b/src/qldpc/circuits/noise_model.py
@@ -48,6 +48,7 @@
 from __future__ import annotations
 
 import collections
+import warnings
 from collections.abc import Collection, Iterable, Iterator
 
 import stim
@@ -185,14 +186,85 @@ def as_noiseless_circuit(circuit: stim.Circuit) -> stim.Circuit:
     return noiseless_circuit
 
 
-def immunize_noise(noise: stim.Circuit, immune_qubits: set[int]) -> stim.Circuit:
+
+# PAULI_CHANNEL_2 arg order: IX(0), IY(1), IZ(2), XI(3), XX(4), XY(5), XZ(6),
+#                             YI(7), YX(8), YY(9), YZ(10), ZI(11), ZX(12), ZY(13), ZZ(14)
+# Marginal indices when the second qubit is immune (surviving: first qubit, non-cross terms: XI, YI, ZI):
+_PC2_SECOND_IMMUNE_INDICES = (3, 7, 11)
+# Marginal indices when the first qubit is immune (surviving: second qubit, non-cross terms: IX, IY, IZ):
+_PC2_FIRST_IMMUNE_INDICES = (0, 1, 2)
+
+
+def _marginalize_2q_noise(
+    noise_op: stim.CircuitInstruction, immune_qubits: set[int], *, marginalize: bool
+) -> stim.Circuit:
+    """Filter or marginalize a 2-qubit noise instruction over immune qubits.
+
+    Processes each pair of targets independently.  Pairs with no immune qubits are kept as-is.
+    Pairs where both qubits are immune are dropped.  For partially-immune pairs, if marginalize is
+    True the surviving qubit receives a 1-qubit marginal (ignoring cross-Pauli terms); if False the
+    pair is dropped.  Only DEPOLARIZE2 and PAULI_CHANNEL_2 support marginalization; other 2-qubit
+    channels emit a warning and are dropped for partially-immune pairs.
+
+    Args:
+        noise_op: A 2-qubit noise instruction.
+        immune_qubits: Set of qubit indices that should not receive noise.
+        marginalize: If True, emit a 1-qubit marginal for partially-immune pairs.
+
+    Returns:
+        stim.Circuit: The filtered/marginalized circuit for this instruction.
+    """
+    result = stim.Circuit()
+    name = noise_op.name
+    args = noise_op.gate_args_copy()
+    targets = noise_op.targets_copy()
+
+    for i in range(0, len(targets), 2):
+        q1, q2 = targets[i], targets[i + 1]
+        q1_immune = q1.value in immune_qubits
+        q2_immune = q2.value in immune_qubits
+
+        if not q1_immune and not q2_immune:
+            result.append(stim.CircuitInstruction(name, [q1, q2], args))
+        elif q1_immune and q2_immune:
+            pass  # both immune: skip
+        elif not marginalize:
+            pass  # partially immune, no marginalization: drop the pair
+        elif name == "DEPOLARIZE2":
+            p = args[0]
+            surviving = q2.value if q1_immune else q1.value
+            result.append(stim.CircuitInstruction("DEPOLARIZE1", [surviving], [p / 5]))
+        elif name == "PAULI_CHANNEL_2":
+            if q2_immune:
+                marginal = [args[idx] for idx in _PC2_SECOND_IMMUNE_INDICES]
+                surviving = q1.value
+            else:
+                marginal = [args[idx] for idx in _PC2_FIRST_IMMUNE_INDICES]
+                surviving = q2.value
+            result.append(stim.CircuitInstruction("PAULI_CHANNEL_1", [surviving], marginal))
+        else:  # pragma: no cover
+            warnings.warn(
+                f"Cannot marginalize {name} over immune qubits; noise is dropped.",
+                stacklevel=2,
+            )
+
+    return result
+
+
+def immunize_noise(
+    noise: stim.Circuit, immune_qubits: set[int], *, marginalize: bool = False
+) -> stim.Circuit:
     """Return a copy of a flat noise circuit with instructions targeting immune qubits removed.
 
-    An instruction is removed if any of its qubit targets belongs to immune_qubits.
+    An instruction is removed if any of its qubit targets belongs to immune_qubits.  If marginalize
+    is True, DEPOLARIZE2 and PAULI_CHANNEL_2 instructions where only one qubit in a pair is immune
+    are replaced by the marginal 1-qubit noise channel on the surviving qubit (ignoring cross-Pauli
+    terms).
 
     Args:
         noise: A flat noise circuit (no repeat blocks) to filter.
         immune_qubits: Set of qubit indices that should not have noise applied to them.
+        marginalize: If True, marginalize 2-qubit noise onto surviving qubits instead of removing.
 
     Returns:
         stim.Circuit: A filtered copy of the input circuit.
@@ -200,20 +272,17 @@ def immunize_noise(noise: stim.Circuit, immune_qubits: set[int]) -> stim.Circuit
     if not immune_qubits:
         return noise
     result = stim.Circuit()
-    for instruction in noise:
-        assert isinstance(instruction, stim.CircuitInstruction)
-        if not any(
-            (
-                target.is_qubit_target
-                or target.is_x_target
-                or target.is_y_target
-                or target.is_z_target
-            )
-            and target.value in immune_qubits
-            for target in instruction.targets_copy()
-            if not target.is_combiner
-        ):
-            result.append(instruction)
+    for noise_op in noise:
+        assert isinstance(noise_op, stim.CircuitInstruction)
+        if all(t.value not in immune_qubits for t in noise_op.targets_copy() if not t.is_combiner):
+            result.append(noise_op)
+        elif stim.gate_data(noise_op.name).is_two_qubit_gate:
+            result += _marginalize_2q_noise(noise_op, immune_qubits, marginalize=marginalize)
+        else:
+            # 1-qubit noise with multiple targets: keep only non-immune targets
+            surviving = [t for t in noise_op.targets_copy() if t.value not in immune_qubits]
+            if surviving:
+                result.append(stim.CircuitInstruction(noise_op.name, surviving, noise_op.gate_args_copy()))
     return result
 
 
@@ -520,6 +589,7 @@ def noisy_circuit(
         immune_qubits: Collection[int] | None = None,
         immune_op_tag: str = DEFAULT_IMMUNE_OP_TAG,
         insert_ticks: bool = True,
+        marginalize: bool = False
     ) -> stim.Circuit:
         f"""Construct a noisy version of the given circuit.
 
@@ -537,6 +607,7 @@ def noisy_circuit(
                 noiseless.  Default: "{DEFAULT_IMMUNE_OP_TAG}".
             insert_ticks: If True, automatically inserts TICK operations to prevent qubit reuse
                 conflicts.  If False, assumes that this preprocessing is not necessary.
+            marginalize: If True, marginalize 2-qubit noise onto surviving qubits instead of removing.
 
         Returns:
             stim.Circuit: A noisy version of the input circuit.
@@ -569,6 +640,7 @@ def noisy_circuit(
                         moment_or_repeat_block.body_copy(),
                         system_qubits=system_qubits,
                         immune_qubits=immune_qubits,
+                        marginalize=marginalize,
                     )
                     noisy_body.append("TICK")
                     noisy_circuit.append(
@@ -585,6 +657,7 @@ def noisy_circuit(
                     system_qubits=system_qubits,
                     immune_qubits=immune_qubits,
                     immune_op_tag=immune_op_tag,
+                    marginalize=marginalize,
                 )
 
         return noisy_circuit
@@ -597,6 +670,7 @@ def _inplace_append_noisy_moment(
         system_qubits: set[int],
         immune_qubits: set[int],
         immune_op_tag: str,
+        marginalize: bool,
     ) -> None:
         """Apps noise to a moment and appends it to a circuit (in-place).
 
@@ -610,6 +684,8 @@ def _inplace_append_noisy_moment(
             immune_qubits: Set of all qubits that should not have noise applied to them.
             immune_op_tag: If an operation contains this string in its tag, that operation is
                 noiseless.
+            marginalize: If True, marginalize 2-qubit noise onto surviving qubits instead of
+                removing.
         """
         noise_after_moment = stim.Circuit()
         for op in moment:
@@ -617,25 +693,10 @@ def _inplace_append_noisy_moment(
                 circuit.append(op)
             else:
                 noisy_op, after = rule.noisy_operation(op)
-                op_targets_immune = immune_qubits and any(
-                    (
-                        target.is_qubit_target
-                        or target.is_x_target
-                        or target.is_y_target
-                        or target.is_z_target
-                    )
-                    and target.value in immune_qubits
-                    for target in op.targets_copy()
-                )
-                circuit.append(op if op_targets_immune else noisy_op)
-                noise_after_moment += immunize_noise(after, immune_qubits)
-
-        # TODO: post-process noise_after_moment using immune_qubits
-        # see comments:
-        # - https://github.com/qLDPCOrg/qLDPC/issues/426#issuecomment-4144385907
-        # - https://github.com/qLDPCOrg/qLDPC/issues/426#issuecomment-4149324799
+                circuit.append(noisy_op)
+                noise_after_moment += after
 
-        circuit += noise_after_moment
+        circuit += immunize_noise(noise_after_moment, immune_qubits, marginalize=marginalize)
 
         moment_was_noisy = any(immune_op_tag not in op.tag for op in moment)
         if moment_was_noisy and self.idle_error or self.additional_error_waiting_for_m_or_r:
diff --git a/src/qldpc/circuits/noise_model_test.py b/src/qldpc/circuits/noise_model_test.py
index 0a96655e6..0c5181e3d 100644
--- a/src/qldpc/circuits/noise_model_test.py
+++ b/src/qldpc/circuits/noise_model_test.py
@@ -249,11 +249,14 @@ def test_immunity() -> None:
     # qubits can be immune to errors
     circuit = stim.Circuit("""
         H 0 1
+        CNOT 1 2
     """)
     noise_model = circuits.DepolarizingNoiseModel(0.1, include_idling_error=False)
     noisy_circuit = stim.Circuit("""
         H 0 1
+        CNOT 1 2
         DEPOLARIZE1(0.1) 1
+        DEPOLARIZE2(0.1) 1 2
     """)
     assert _circuits_are_equivalent(
         noisy_circuit, noise_model.noisy_circuit(circuit, immune_qubits=[0], insert_ticks=False)
@@ -291,6 +294,72 @@ def test_immunity() -> None:
     assert noise_model.noisy_circuit(noiseless_circuit).to_tableau() == tableau
 
 
+def test_immune_marginalize() -> None:
+    circuit = stim.Circuit("""
+        CNOT 0 1
+        CNOT 1 2
+        CNOT 3 4
+    """)
+    noise_model = circuits.DepolarizingNoiseModel(0.1, include_idling_error=False)
+
+    # error marginalize on part of DEPOLARIZE2 on qubit 0 1
+    noisy_circuit_marginalize = stim.Circuit("""
+        CNOT 0 1
+        CNOT 1 2
+        CNOT 3 4
+        DEPOLARIZE1(0.02) 2
+        DEPOLARIZE2(0.1) 3 4
+    """)
+    assert _circuits_are_equivalent(
+        noisy_circuit_marginalize,
+        noise_model.noisy_circuit(
+            circuit, immune_qubits=[0, 1], insert_ticks=False, marginalize=True
+        ),
+    )
+
+    # turn off marginalization
+    noisy_circuit_marginalize = stim.Circuit("""
+        CNOT 0 1
+        CNOT 1 2
+        CNOT 3 4
+        DEPOLARIZE2(0.1) 3 4
+    """)
+    assert _circuits_are_equivalent(
+        noisy_circuit_marginalize,
+        noise_model.noisy_circuit(
+            circuit, immune_qubits=[0, 1], insert_ticks=False, marginalize=False
+        ),
+    )
+
+    # test PAULI_CHANNEL_2
+    circuit = stim.Circuit("""
+        CNOT 0 1
+    """)
+    noise_model = circuits.NoiseModel(
+        clifford_2q_error=circuits.NoiseRule(
+            after={
+                "PAULI_CHANNEL_2": [0.01, 0.02, 0.03, 0.04, 0, 0, 0, 0.05, 0, 0, 0, 0.06, 0, 0, 0]
+            }
+        )
+    )
+
+    assert _circuits_are_equivalent(
+        stim.Circuit("""
+            CNOT 0 1
+            PAULI_CHANNEL_1(0.01, 0.02, 0.03) 1
+        """),
+        noise_model.noisy_circuit(circuit, immune_qubits=[0], insert_ticks=False, marginalize=True),
+    )
+
+    assert _circuits_are_equivalent(
+        stim.Circuit("""
+            CNOT 0 1
+            PAULI_CHANNEL_1(0.04, 0.05, 0.06) 0
+        """),
+        noise_model.noisy_circuit(circuit, immune_qubits=[1], insert_ticks=False, marginalize=True),
+    )
+
+
 def test_classical_controls() -> None:
     """Classically controled gates get special treatment."""
     noise_model: circuits.NoiseModel

From 3748d29031d9deb38e6780828d31bfde31f674ca Mon Sep 17 00:00:00 2001
From: "Michael A. Perlin" <michael.perlin@jpmchase.com>
Date: Wed, 24 Jun 2026 14:12:04 -0400
Subject: [PATCH 12/13] nits

---
 src/qldpc/circuits/noise_model.py | 22 +++++++++-------------
 1 file changed, 9 insertions(+), 13 deletions(-)

diff --git a/src/qldpc/circuits/noise_model.py b/src/qldpc/circuits/noise_model.py
index 89a4623f7..68dbc8fe0 100644
--- a/src/qldpc/circuits/noise_model.py
+++ b/src/qldpc/circuits/noise_model.py
@@ -186,7 +186,6 @@ def as_noiseless_circuit(circuit: stim.Circuit) -> stim.Circuit:
     return noiseless_circuit
 
 
-
 # PAULI_CHANNEL_2 arg order: IX(0), IY(1), IZ(2), XI(3), XX(4), XY(5), XZ(6),
 #                             YI(7), YX(8), YY(9), YZ(10), ZI(11), ZX(12), ZY(13), ZZ(14)
 # Marginal indices when the second qubit is immune (surviving: first qubit, non-cross terms: XI, YI, ZI):
@@ -282,7 +281,9 @@ def immunize_noise(
             # 1-qubit noise with multiple targets: keep only non-immune targets
             surviving = [t for t in noise_op.targets_copy() if t.value not in immune_qubits]
             if surviving:
-                result.append(stim.CircuitInstruction(noise_op.name, surviving, noise_op.gate_args_copy()))
+                result.append(
+                    stim.CircuitInstruction(noise_op.name, surviving, noise_op.gate_args_copy())
+                )
     return result
 
 
@@ -346,9 +347,6 @@ def noisy_operation(
     ) -> tuple[stim.CircuitInstruction, stim.Circuit]:
         """Apply this noise rule to the given operation.
 
-        Immunity to noise is not handled here.  The caller is responsible for filtering the
-        returned noise circuit using immunize_noise if needed.
-
         Args:
             op: The operation to add noise to.
 
@@ -387,7 +385,7 @@ def _build_noise_after(self, op: stim.CircuitInstruction) -> stim.Circuit:
         """Build the extra noise circuit to append after the given operation.
 
         Subclasses may override this to customize the noise that follows an operation.  Reset errors
-        are excluded here and handled separately in noisy_operation.
+        are excluded here and handled separately in NoiseRule.noisy_operation.
 
         Args:
             op: The operation being applied.
@@ -464,14 +462,10 @@ def is_targeted_noisy_op(self, op: stim.CircuitInstruction) -> bool:
         return op.targets_copy() == self.noisy_op.targets_copy()
 
     def noisy_operation(
-        self,
-        op: stim.CircuitInstruction,
+        self, op: stim.CircuitInstruction
     ) -> tuple[stim.CircuitInstruction, stim.Circuit]:
         """Apply this targeted noise rule to the given operation.
 
-        Immunity to noise is not handled here.  The caller is responsible for filtering the
-        returned noise circuit using immunize_noise if needed.
-
         Args:
             op: The operation to add noise to.
 
@@ -588,8 +582,8 @@ def noisy_circuit(
         system_qubits: Collection[int] | None = None,
         immune_qubits: Collection[int] | None = None,
         immune_op_tag: str = DEFAULT_IMMUNE_OP_TAG,
+        marginalize: bool = False,
         insert_ticks: bool = True,
-        marginalize: bool = False
     ) -> stim.Circuit:
         f"""Construct a noisy version of the given circuit.
 
@@ -607,7 +601,9 @@ def noisy_circuit(
                 noiseless.  Default: "{DEFAULT_IMMUNE_OP_TAG}".
             insert_ticks: If True, automatically inserts TICK operations to prevent qubit reuse
                 conflicts.  If False, assumes that this preprocessing is not necessary.
-            marginalize: If True, marginalize 2-qubit noise onto surviving qubits instead of removing.
+            marginalize: If True, marginalize 2-qubit noise on the absense of errors on noise-immune
+                qubits.  If False, gates that address noise-immune qubits are also noiseless.
+                Default: False.
 
         Returns:
             stim.Circuit: A noisy version of the input circuit.

From a388b0f578846502d3a5321595b410cbf7a6e0ae Mon Sep 17 00:00:00 2001
From: "Michael A. Perlin" <michael.perlin@jpmchase.com>
Date: Wed, 24 Jun 2026 14:18:46 -0400
Subject: [PATCH 13/13] some rearranging

---
 src/qldpc/circuits/noise_model.py | 766 +++++++++++++++---------------
 1 file changed, 384 insertions(+), 382 deletions(-)

diff --git a/src/qldpc/circuits/noise_model.py b/src/qldpc/circuits/noise_model.py
index 68dbc8fe0..0f28dcd5e 100644
--- a/src/qldpc/circuits/noise_model.py
+++ b/src/qldpc/circuits/noise_model.py
@@ -1,8 +1,10 @@
 """Implementation of noise models for Stim circuits
 
 The main components of this module are:
-- NoiseRule: Defines how to add noise to individual operations.
 - NoiseModel: Defines how noise is added to circuits.
+- NoiseRule: Defines how to add noise to operations according to their gate type.
+- TargetedNoiseRule: Defines how to add noise to individual operations according to gate type and
+    targeted qubits.
 - Built-in noise models: DepolarizingNoiseModel and the superconducting-inspired SI1000NoiseModel.
 
 Examples of basic usage with a predefined noise model:
@@ -186,410 +188,114 @@ def as_noiseless_circuit(circuit: stim.Circuit) -> stim.Circuit:
     return noiseless_circuit
 
 
-# PAULI_CHANNEL_2 arg order: IX(0), IY(1), IZ(2), XI(3), XX(4), XY(5), XZ(6),
-#                             YI(7), YX(8), YY(9), YZ(10), ZI(11), ZX(12), ZY(13), ZZ(14)
-# Marginal indices when the second qubit is immune (surviving: first qubit, non-cross terms: XI, YI, ZI):
-_PC2_SECOND_IMMUNE_INDICES = (3, 7, 11)
-# Marginal indices when the first qubit is immune (surviving: second qubit, non-cross terms: IX, IY, IZ):
-_PC2_FIRST_IMMUNE_INDICES = (0, 1, 2)
-
-
-def _marginalize_2q_noise(
-    noise_op: stim.CircuitInstruction, immune_qubits: set[int], *, marginalize: bool
-) -> stim.Circuit:
-    """Filter or marginalize a 2-qubit noise instruction over immune qubits.
-
-    Processes each pair of targets independently.  Pairs with no immune qubits are kept as-is.
-    Pairs where both qubits are immune are dropped.  For partially-immune pairs, if marginalize is
-    True the surviving qubit receives a 1-qubit marginal (ignoring cross-Pauli terms); if False the
-    pair is dropped.  Only DEPOLARIZE2 and PAULI_CHANNEL_2 support marginalization; other 2-qubit
-    channels emit a warning and are dropped for partially-immune pairs.
-
-    Args:
-        noise_op: A 2-qubit noise instruction.
-        immune_qubits: Set of qubit indices that should not receive noise.
-        marginalize: If True, emit a 1-qubit marginal for partially-immune pairs.
-
-    Returns:
-        stim.Circuit: The filtered/marginalized circuit for this instruction.
-    """
-    result = stim.Circuit()
-    name = noise_op.name
-    args = noise_op.gate_args_copy()
-    targets = noise_op.targets_copy()
-
-    for i in range(0, len(targets), 2):
-        q1, q2 = targets[i], targets[i + 1]
-        q1_immune = q1.value in immune_qubits
-        q2_immune = q2.value in immune_qubits
-
-        if not q1_immune and not q2_immune:
-            result.append(stim.CircuitInstruction(name, [q1, q2], args))
-        elif q1_immune and q2_immune:
-            pass  # both immune: skip
-        elif not marginalize:
-            pass  # partially immune, no marginalization: drop the pair
-        elif name == "DEPOLARIZE2":
-            p = args[0]
-            surviving = q2.value if q1_immune else q1.value
-            result.append(stim.CircuitInstruction("DEPOLARIZE1", [surviving], [p / 5]))
-        elif name == "PAULI_CHANNEL_2":
-            if q2_immune:
-                marginal = [args[idx] for idx in _PC2_SECOND_IMMUNE_INDICES]
-                surviving = q1.value
-            else:
-                marginal = [args[idx] for idx in _PC2_FIRST_IMMUNE_INDICES]
-                surviving = q2.value
-            result.append(stim.CircuitInstruction("PAULI_CHANNEL_1", [surviving], marginal))
-        else:  # pragma: no cover
-            warnings.warn(
-                f"Cannot marginalize {name} over immune qubits; noise is dropped.",
-                stacklevel=2,
-            )
-
-    return result
-
-
-def immunize_noise(
-    noise: stim.Circuit, immune_qubits: set[int], *, marginalize: bool = False
-) -> stim.Circuit:
-    """Return a copy of a flat noise circuit with instructions targeting immune qubits removed.
-
-    An instruction is removed if any of its qubit targets belongs to immune_qubits.  If marginalize
-    is True, DEPOLARIZE2 and PAULI_CHANNEL_2 instructions where only one qubit in a pair is immune
-    are replaced by the marginal 1-qubit noise channel on the surviving qubit (ignoring cross-Pauli
-    terms).
-
-    Args:
-        noise: A flat noise circuit (no repeat blocks) to filter.
-        immune_qubits: Set of qubit indices that should not have noise applied to them.
-        marginalize: If True, marginalize 2-qubit noise onto surviving qubits instead of removing.
-
-    Returns:
-        stim.Circuit: A filtered copy of the input circuit.
-    """
-    if not immune_qubits:
-        return noise
-    result = stim.Circuit()
-    for noise_op in noise:
-        assert isinstance(noise_op, stim.CircuitInstruction)
-        if all(t.value not in immune_qubits for t in noise_op.targets_copy() if not t.is_combiner):
-            result.append(noise_op)
-        elif stim.gate_data(noise_op.name).is_two_qubit_gate:
-            result += _marginalize_2q_noise(noise_op, immune_qubits, marginalize=marginalize)
-        else:
-            # 1-qubit noise with multiple targets: keep only non-immune targets
-            surviving = [t for t in noise_op.targets_copy() if t.value not in immune_qubits]
-            if surviving:
-                result.append(
-                    stim.CircuitInstruction(noise_op.name, surviving, noise_op.gate_args_copy())
-                )
-    return result
-
-
-class NoiseRule:
-    """Describes how to add noise to an operation.
+class NoiseModel:
+    """A model that defines how to add noise to quantum circuits.
 
-    This class encapsulates the noise channels and measurement error probabilities that should be
-    applied to a particular type of quantum operation.
+    This class provides a framework for adding various types of noise to quantum circuits, including
+    gate errors, readout errors, reset errors, and idling errors.  Classically controlled operations
+    are assumed to NOT occur, so the corresponding qubits pick up idling errors, if applicable.
     """
 
     def __init__(
         self,
+        clifford_1q_error: NoiseRule | float | None = None,
+        clifford_2q_error: NoiseRule | float | None = None,
+        readout_error: float | None = None,
+        reset_error: float | None = None,
         *,
-        after: dict[str, float | Iterable[float]] = {},
-        readout_error: float = 0,
-        reset_error: float = 0,
+        idle_error: float | None = None,
+        additional_error_waiting_for_m_or_r: float | None = None,
+        rules: dict[str, NoiseRule] | None = None,
     ) -> None:
-        """Initializes a noise rule with specified error channels.
+        """Initializes a noise model with specified parameters.
 
         Args:
-            after: A dictionary mapping noise channel names to their probability arguments.  For
-                example, {"DEPOLARIZE2": 0.01, "PAULI_CHANNEL_1": [0.02, 0, 0]} will add two-qubit
-                depolarization with parameter 0.01, followed by 2% bit-flip noise.  These noise
-                channels occur after all other operations in the moment and are applied to the same
-                targets as the relevant operation.
-            readout_error: The probability that a measurement result is reported incorrectly.  Only
-                allowed for operations that produce measurement results.
-            reset_error: The probability that a qubit is reset to the wrong state.  Only allowed for
-                operations that reset qubits.
-
-        Raises:
-            ValueError: If any noise channel name is not recognized or if any net probability of an
-                error is not between 0 and 1 (inclusive).
+            clifford_1q_error: Default noise rule or depolarization probability for one-qubit
+                unitary Clifford gates.
+            clifford_2q_error: Default noise rule or depolarization probability for two-qubit
+                unitary Clifford gates.
+            readout_error: Default probability of flipping measurement results.
+            reset_error: Default probability of resetting qubits to the wrong state.
+            idle_error: Probability of depolarization for each idling qubit in any given moment.
+            additional_error_waiting_for_m_or_r: Additional depolarization probability applied to
+                qubits that are waiting while other qubits undergo measurement or reset operations.
+            rules: Dictionary mapping specific gate names to their noise rules.  Overrides all other
+                rules for unitary, measurement, and reset gates.
         """
-        self.readout_error = readout_error
-        if not (0 <= readout_error <= 1):
-            raise ValueError(f"{readout_error=} is not between 0 and 1")
-
-        self.reset_error = reset_error
-        if not (0 <= reset_error <= 1):
-            raise ValueError(f"{reset_error=} is not between 0 and 1")
+        if not (isinstance(clifford_1q_error, NoiseRule) or clifford_1q_error is None):
+            clifford_1q_error = NoiseRule(after={"DEPOLARIZE1": clifford_1q_error})
+        if not (isinstance(clifford_2q_error, NoiseRule) or clifford_2q_error is None):
+            clifford_2q_error = NoiseRule(after={"DEPOLARIZE2": clifford_2q_error})
 
-        self.after = {
-            op: tuple(prob_or_probs) if isinstance(prob_or_probs, Iterable) else (prob_or_probs,)
-            for op, prob_or_probs in after.items()
-        }
-        for op, probs in self.after.items():
-            if OP_TYPES[op] != NOISE:
-                raise ValueError(f"Invalid or unrecognized noise channel {op} in {after=}")
-            if not (0 <= sum(probs) <= 1):
-                raise ValueError(
-                    f"The net probability of an error is not between 0 and 1 in {after=}"
-                )
+        self.rules = rules
+        self.clifford_1q_error = clifford_1q_error
+        self.clifford_2q_error = clifford_2q_error
+        self.readout_error = readout_error or 0
+        self.reset_error = reset_error or 0
+        self.idle_error = idle_error
+        self.additional_error_waiting_for_m_or_r = additional_error_waiting_for_m_or_r
 
     def __bool__(self) -> bool:
-        """Is this noise rule nontrivial?"""
-        return bool(self.after) or bool(self.readout_error) or bool(self.reset_error)
+        """Is this noise model nontrivial?"""
+        return (
+            bool(self.rules)
+            or bool(self.clifford_1q_error)
+            or bool(self.clifford_2q_error)
+            or bool(self.readout_error)
+            or bool(self.reset_error)
+            or bool(self.idle_error)
+            or bool(self.additional_error_waiting_for_m_or_r)
+        )
 
-    def noisy_operation(
-        self, op: stim.CircuitInstruction
-    ) -> tuple[stim.CircuitInstruction, stim.Circuit]:
-        """Apply this noise rule to the given operation.
+    def get_noise_rule(self, op: stim.CircuitInstruction) -> NoiseRule | None:
+        """Determines the noise rule to apply to a specific operation.
 
         Args:
-            op: The operation to add noise to.
+            op: The circuit instruction to find a noise rule for.
 
         Returns:
-            stim.CircuitInstruction: The given operation possibly modified to account for noise.
-            stim.Circuit: Noise operations that should follow the given operation.
+            The NoiseRule to apply for the given operation, or None for no noise.
         """
-        targets = op.targets_copy()
-        args = op.gate_args_copy()
-
-        if self.readout_error:
-            assert op.name in JUST_MEASURE_OPS or op.name in MEASURE_AND_RESET_OPS
-            if not args:
-                args = [self.readout_error]
-            else:
-                assert len(args) == 1
-                # combine the bit-flip probabilities self.readout_error and args[0]
-                args = [1 - (1 - self.readout_error) * (1 - args[0])]
+        if OP_TYPES[op.name] == ANNOTATION or _involves_classical_bits(op):
+            return None
 
-        noisy_op = stim.CircuitInstruction(op.name, targets, args, tag=op.tag)
-        noise_after = stim.Circuit()
+        if self.rules is not None:
+            rule = self.rules.get(_get_standardized_name(op)) or self.rules.get(
+                op.name
+            )  # allows for an MPP rule, but first checks for rules such as MXY
+            if rule is not None:
+                return rule
 
-        if self.reset_error:
-            assert op.name in JUST_RESET_OPS or op.name in MEASURE_AND_RESET_OPS
-            qubit_targets = [target.value for target in targets if not target.is_combiner]
-            error_name = ("X" if _get_standardized_name(op)[-1] != "X" else "Z") + "_ERROR"
-            noise_after.append(
-                stim.CircuitInstruction(error_name, qubit_targets, [self.reset_error])
-            )
+        op_type = OP_TYPES[op.name]
+        if self.clifford_1q_error is not None and op_type == CLIFFORD_1Q:
+            return self.clifford_1q_error
+        if self.clifford_2q_error is not None and op_type == CLIFFORD_2Q:
+            return self.clifford_2q_error
 
-        noise_after += self._build_noise_after(op)
+        if self.readout_error and op.name in JUST_MEASURE_OPS:
+            return NoiseRule(readout_error=self.readout_error)
+        if self.reset_error and op.name in JUST_RESET_OPS:
+            return NoiseRule(reset_error=self.reset_error)
+        if (self.readout_error or self.reset_error) and op.name in MEASURE_AND_RESET_OPS:
+            return NoiseRule(readout_error=self.readout_error, reset_error=self.reset_error)
 
-        return noisy_op, noise_after
+        return None
 
-    def _build_noise_after(self, op: stim.CircuitInstruction) -> stim.Circuit:
-        """Build the extra noise circuit to append after the given operation.
+    def noisy_circuit(
+        self,
+        circuit: stim.Circuit,
+        *,
+        system_qubits: Collection[int] | None = None,
+        immune_qubits: Collection[int] | None = None,
+        immune_op_tag: str = DEFAULT_IMMUNE_OP_TAG,
+        marginalize: bool = False,
+        insert_ticks: bool = True,
+    ) -> stim.Circuit:
+        f"""Construct a noisy version of the given circuit.
 
-        Subclasses may override this to customize the noise that follows an operation.  Reset errors
-        are excluded here and handled separately in NoiseRule.noisy_operation.
-
-        Args:
-            op: The operation being applied.
-
-        Returns:
-            stim.Circuit: Additional noise to append after the operation.
-        """
-        qubit_targets = [target.value for target in op.targets_copy() if not target.is_combiner]
-        noise = stim.Circuit()
-        for op_name, args in self.after.items():
-            noise.append(stim.CircuitInstruction(op_name, qubit_targets, args))
-        return noise
-
-
-class TargetedNoiseRule(NoiseRule):
-    """Describes how to add noise to a specific circuit instruction on specific qubits.
-
-    Unlike NoiseRule, which applies to all operations of a given type, this rule matches only an
-    exact gate-and-target combination, allowing fine-grained per-operation noise overrides.
-    """
-
-    def __init__(
-        self,
-        *,
-        noisy_op: stim.CircuitInstruction,
-        noise: stim.Circuit,
-        readout_error: float = 0,
-        reset_error: float = 0,
-        tags: Collection[str] | None = None,
-    ) -> None:
-        """Initializes a targeted noise rule for a specific circuit instruction.
-
-        Args:
-            noisy_op: The circuit instruction that this rule targets.  Defines the gate name and
-                qubit targets to match against.  Gate args on this instruction are ignored during
-                matching.
-            noise: An explicit noise circuit to append after the matched operation.
-            readout_error: The probability that a measurement result is reported incorrectly.  Only
-                allowed for operations that produce measurement results.
-            reset_error: The probability that a qubit is reset to the wrong state.  Only allowed for
-                operations that reset qubits.
-            tags: If not None, only match operations whose tag exactly matches one of the given
-                strings.  If None, match operations regardless of their tag.
-
-        Raises:
-            ValueError: If readout_error or reset_error is not between 0 and 1 (inclusive).
-        """
-        super().__init__(readout_error=readout_error, reset_error=reset_error)
-        self.noisy_op = noisy_op
-        self.noise = noise
-        self.tags: frozenset[str] | None = frozenset(tags) if tags is not None else None
-
-    def __bool__(self) -> bool:
-        """Is this noise rule nontrivial?"""
-        nontrivial_noise = bool(self.noise) or bool(self.readout_error) or bool(self.reset_error)
-        nontrivial_targets = bool(self.noisy_op.targets_copy())
-        return nontrivial_noise and nontrivial_targets
-
-    def is_targeted_noisy_op(self, op: stim.CircuitInstruction) -> bool:
-        """Determine whether the given operation matches this rule's target instruction.
-
-        Two operations match if they have the same gate name and the same qubit target values in the
-        same order.  Gate args are ignored.  If self.tags is not None, op.tag must exactly match an
-        element of self.tags.
-
-        Args:
-            op: The circuit instruction to check.
-
-        Returns:
-            True if op matches this rule's target instruction.  False otherwise.
-        """
-        if op.name != self.noisy_op.name or (self.tags is not None and op.tag not in self.tags):
-            return False
-        return op.targets_copy() == self.noisy_op.targets_copy()
-
-    def noisy_operation(
-        self, op: stim.CircuitInstruction
-    ) -> tuple[stim.CircuitInstruction, stim.Circuit]:
-        """Apply this targeted noise rule to the given operation.
-
-        Args:
-            op: The operation to add noise to.
-
-        Returns:
-            stim.CircuitInstruction: The given operation, possibly modified to account for readout
-                or reset errors.
-            stim.Circuit: Noise operations that should follow the given operation.
-        """
-        if not self.is_targeted_noisy_op(op):
-            return op, stim.Circuit()
-        return super().noisy_operation(op)
-
-    def _build_noise_after(self, op: stim.CircuitInstruction) -> stim.Circuit:
-        return self.noise
-
-
-class NoiseModel:
-    """A model that defines how to add noise to quantum circuits.
-
-    This class provides a framework for adding various types of noise to quantum circuits, including
-    gate errors, readout errors, reset errors, and idling errors.  Classically controlled operations
-    are assumed to NOT occur, so the corresponding qubits pick up idling errors, if applicable.
-    """
-
-    def __init__(
-        self,
-        clifford_1q_error: NoiseRule | float | None = None,
-        clifford_2q_error: NoiseRule | float | None = None,
-        readout_error: float | None = None,
-        reset_error: float | None = None,
-        *,
-        idle_error: float | None = None,
-        additional_error_waiting_for_m_or_r: float | None = None,
-        rules: dict[str, NoiseRule] | None = None,
-    ) -> None:
-        """Initializes a noise model with specified parameters.
-
-        Args:
-            clifford_1q_error: Default noise rule or depolarization probability for one-qubit unitary
-                Clifford gates.
-            clifford_2q_error: Default noise rule or depolarization probability for two-qubit unitary
-                Clifford gates.
-            readout_error: Default probability of flipping measurement results.
-            reset_error: Default probability of resetting qubits to the wrong state.
-            idle_error: Probability of depolarization for each idling qubit in any given moment.
-            additional_error_waiting_for_m_or_r: Additional depolarization probability applied to
-                qubits that are waiting while other qubits undergo measurement or reset operations.
-            rules: Dictionary mapping specific gate names to their noise rules.  Overrides all other
-                rules for unitary, measurement, and reset gates.
-        """
-        if not (isinstance(clifford_1q_error, NoiseRule) or clifford_1q_error is None):
-            clifford_1q_error = NoiseRule(after={"DEPOLARIZE1": clifford_1q_error})
-        if not (isinstance(clifford_2q_error, NoiseRule) or clifford_2q_error is None):
-            clifford_2q_error = NoiseRule(after={"DEPOLARIZE2": clifford_2q_error})
-
-        self.rules = rules
-        self.clifford_1q_error = clifford_1q_error
-        self.clifford_2q_error = clifford_2q_error
-        self.readout_error = readout_error or 0
-        self.reset_error = reset_error or 0
-        self.idle_error = idle_error
-        self.additional_error_waiting_for_m_or_r = additional_error_waiting_for_m_or_r
-
-    def __bool__(self) -> bool:
-        """Is this noise model nontrivial?"""
-        return (
-            bool(self.rules)
-            or bool(self.clifford_1q_error)
-            or bool(self.clifford_2q_error)
-            or bool(self.readout_error)
-            or bool(self.reset_error)
-            or bool(self.idle_error)
-            or bool(self.additional_error_waiting_for_m_or_r)
-        )
-
-    def get_noise_rule(self, op: stim.CircuitInstruction) -> NoiseRule | None:
-        """Determines the noise rule to apply to a specific operation.
-
-        Args:
-            op: The circuit instruction to find a noise rule for.
-
-        Returns:
-            The NoiseRule to apply for the given operation, or None for no noise.
-        """
-        if OP_TYPES[op.name] == ANNOTATION or _involves_classical_bits(op):
-            return None
-
-        if self.rules is not None:
-            rule = self.rules.get(_get_standardized_name(op)) or self.rules.get(
-                op.name
-            )  # allows for an MPP rule, but first checks for rules such as MXY
-            if rule is not None:
-                return rule
-
-        op_type = OP_TYPES[op.name]
-        if self.clifford_1q_error is not None and op_type == CLIFFORD_1Q:
-            return self.clifford_1q_error
-        if self.clifford_2q_error is not None and op_type == CLIFFORD_2Q:
-            return self.clifford_2q_error
-
-        if self.readout_error and op.name in JUST_MEASURE_OPS:
-            return NoiseRule(readout_error=self.readout_error)
-        if self.reset_error and op.name in JUST_RESET_OPS:
-            return NoiseRule(reset_error=self.reset_error)
-        if (self.readout_error or self.reset_error) and op.name in MEASURE_AND_RESET_OPS:
-            return NoiseRule(readout_error=self.readout_error, reset_error=self.reset_error)
-
-        return None
-
-    def noisy_circuit(
-        self,
-        circuit: stim.Circuit,
-        *,
-        system_qubits: Collection[int] | None = None,
-        immune_qubits: Collection[int] | None = None,
-        immune_op_tag: str = DEFAULT_IMMUNE_OP_TAG,
-        marginalize: bool = False,
-        insert_ticks: bool = True,
-    ) -> stim.Circuit:
-        f"""Construct a noisy version of the given circuit.
-
-        This method first uses TICKs to split the input circuit into moments of operations that can
-        be applied in parallel, thereby preventing qubit reuse conflicts.  Noise is then applied to
-        each operation according to the rules of this NoiseModel.
+        This method first uses TICKs to split the input circuit into moments of operations that can
+        be applied in parallel, thereby preventing qubit reuse conflicts.  Noise is then applied to
+        each operation according to the rules of this NoiseModel.
 
         Args:
             circuit: The circuit to apply noise to.
@@ -713,8 +419,8 @@ def _inplace_append_idle_errors(
     ) -> None:
         """Append idling errors from the given moment to the given circuit.
 
-        This method identifies which qubits are idle during a moment and applies depolarization noise
-        to them according to the noise model parameters.
+        This method identifies which qubits are idle during a moment and applies depolarization
+        noise to them according to the noise model parameters.
 
         Args:
             circuit: The circuit to append idle error operations to.
@@ -812,6 +518,302 @@ def __init__(self, p: float) -> None:
         )
 
 
+class NoiseRule:
+    """Describes how to add noise to an operation.
+
+    This class encapsulates the noise channels and measurement error probabilities that should be
+    applied to a particular type of quantum operation.
+    """
+
+    def __init__(
+        self,
+        *,
+        after: dict[str, float | Iterable[float]] = {},
+        readout_error: float = 0,
+        reset_error: float = 0,
+    ) -> None:
+        """Initializes a noise rule with specified error channels.
+
+        Args:
+            after: A dictionary mapping noise channel names to their probability arguments.  For
+                example, {"DEPOLARIZE2": 0.01, "PAULI_CHANNEL_1": [0.02, 0, 0]} will add two-qubit
+                depolarization with parameter 0.01, followed by 2% bit-flip noise.  These noise
+                channels occur after all other operations in the moment and are applied to the same
+                targets as the relevant operation.
+            readout_error: The probability that a measurement result is reported incorrectly.  Only
+                allowed for operations that produce measurement results.
+            reset_error: The probability that a qubit is reset to the wrong state.  Only allowed for
+                operations that reset qubits.
+
+        Raises:
+            ValueError: If any noise channel name is not recognized or if any net probability of an
+                error is not between 0 and 1 (inclusive).
+        """
+        self.readout_error = readout_error
+        if not (0 <= readout_error <= 1):
+            raise ValueError(f"{readout_error=} is not between 0 and 1")
+
+        self.reset_error = reset_error
+        if not (0 <= reset_error <= 1):
+            raise ValueError(f"{reset_error=} is not between 0 and 1")
+
+        self.after = {
+            op: tuple(prob_or_probs) if isinstance(prob_or_probs, Iterable) else (prob_or_probs,)
+            for op, prob_or_probs in after.items()
+        }
+        for op, probs in self.after.items():
+            if OP_TYPES[op] != NOISE:
+                raise ValueError(f"Invalid or unrecognized noise channel {op} in {after=}")
+            if not (0 <= sum(probs) <= 1):
+                raise ValueError(
+                    f"The net probability of an error is not between 0 and 1 in {after=}"
+                )
+
+    def __bool__(self) -> bool:
+        """Is this noise rule nontrivial?"""
+        return bool(self.after) or bool(self.readout_error) or bool(self.reset_error)
+
+    def noisy_operation(
+        self, op: stim.CircuitInstruction
+    ) -> tuple[stim.CircuitInstruction, stim.Circuit]:
+        """Apply this noise rule to the given operation.
+
+        Args:
+            op: The operation to add noise to.
+
+        Returns:
+            stim.CircuitInstruction: The given operation possibly modified to account for noise.
+            stim.Circuit: Noise operations that should follow the given operation.
+        """
+        targets = op.targets_copy()
+        args = op.gate_args_copy()
+
+        if self.readout_error:
+            assert op.name in JUST_MEASURE_OPS or op.name in MEASURE_AND_RESET_OPS
+            if not args:
+                args = [self.readout_error]
+            else:
+                assert len(args) == 1
+                # combine the bit-flip probabilities self.readout_error and args[0]
+                args = [1 - (1 - self.readout_error) * (1 - args[0])]
+
+        noisy_op = stim.CircuitInstruction(op.name, targets, args, tag=op.tag)
+        noise_after = stim.Circuit()
+
+        if self.reset_error:
+            assert op.name in JUST_RESET_OPS or op.name in MEASURE_AND_RESET_OPS
+            qubit_targets = [target.value for target in targets if not target.is_combiner]
+            error_name = ("X" if _get_standardized_name(op)[-1] != "X" else "Z") + "_ERROR"
+            noise_after.append(
+                stim.CircuitInstruction(error_name, qubit_targets, [self.reset_error])
+            )
+
+        noise_after += self._build_noise_after(op)
+
+        return noisy_op, noise_after
+
+    def _build_noise_after(self, op: stim.CircuitInstruction) -> stim.Circuit:
+        """Build the extra noise circuit to append after the given operation.
+
+        Subclasses may override this to customize the noise that follows an operation.  Reset errors
+        are excluded here and handled separately in NoiseRule.noisy_operation.
+
+        Args:
+            op: The operation being applied.
+
+        Returns:
+            stim.Circuit: Additional noise to append after the operation.
+        """
+        qubit_targets = [target.value for target in op.targets_copy() if not target.is_combiner]
+        noise = stim.Circuit()
+        for op_name, args in self.after.items():
+            noise.append(stim.CircuitInstruction(op_name, qubit_targets, args))
+        return noise
+
+
+class TargetedNoiseRule(NoiseRule):
+    """Describes how to add noise to a specific circuit instruction on specific qubits.
+
+    Unlike NoiseRule, which applies to all operations of a given type, this rule matches only an
+    exact gate-and-target combination, allowing fine-grained per-operation noise overrides.
+    """
+
+    def __init__(
+        self,
+        *,
+        noisy_op: stim.CircuitInstruction,
+        noise: stim.Circuit,
+        readout_error: float = 0,
+        reset_error: float = 0,
+        tags: Collection[str] | None = None,
+    ) -> None:
+        """Initializes a targeted noise rule for a specific circuit instruction.
+
+        Args:
+            noisy_op: The circuit instruction that this rule targets.  Defines the gate name and
+                qubit targets to match against.  Gate args on this instruction are ignored during
+                matching.
+            noise: An explicit noise circuit to append after the matched operation.
+            readout_error: The probability that a measurement result is reported incorrectly.  Only
+                allowed for operations that produce measurement results.
+            reset_error: The probability that a qubit is reset to the wrong state.  Only allowed for
+                operations that reset qubits.
+            tags: If not None, only match operations whose tag exactly matches one of the given
+                strings.  If None, match operations regardless of their tag.
+
+        Raises:
+            ValueError: If readout_error or reset_error is not between 0 and 1 (inclusive).
+        """
+        super().__init__(readout_error=readout_error, reset_error=reset_error)
+        self.noisy_op = noisy_op
+        self.noise = noise
+        self.tags: frozenset[str] | None = frozenset(tags) if tags is not None else None
+
+    def __bool__(self) -> bool:
+        """Is this noise rule nontrivial?"""
+        nontrivial_noise = bool(self.noise) or bool(self.readout_error) or bool(self.reset_error)
+        nontrivial_targets = bool(self.noisy_op.targets_copy())
+        return nontrivial_noise and nontrivial_targets
+
+    def is_targeted_noisy_op(self, op: stim.CircuitInstruction) -> bool:
+        """Determine whether the given operation matches this rule's target instruction.
+
+        Two operations match if they have the same gate name and the same qubit target values in the
+        same order.  Gate args are ignored.  If self.tags is not None, op.tag must exactly match an
+        element of self.tags.
+
+        Args:
+            op: The circuit instruction to check.
+
+        Returns:
+            True if op matches this rule's target instruction.  False otherwise.
+        """
+        if op.name != self.noisy_op.name or (self.tags is not None and op.tag not in self.tags):
+            return False
+        return op.targets_copy() == self.noisy_op.targets_copy()
+
+    def noisy_operation(
+        self, op: stim.CircuitInstruction
+    ) -> tuple[stim.CircuitInstruction, stim.Circuit]:
+        """Apply this targeted noise rule to the given operation.
+
+        Args:
+            op: The operation to add noise to.
+
+        Returns:
+            stim.CircuitInstruction: The given operation, possibly modified to account for readout
+                or reset errors.
+            stim.Circuit: Noise operations that should follow the given operation.
+        """
+        if not self.is_targeted_noisy_op(op):
+            return op, stim.Circuit()
+        return super().noisy_operation(op)
+
+    def _build_noise_after(self, op: stim.CircuitInstruction) -> stim.Circuit:
+        return self.noise
+
+
+# PAULI_CHANNEL_2 arg order:
+# IX(0), IY(1), IZ(2), XI(3), XX(4), XY(5), XZ(6), YI(7), YX(8), YY(9), YZ(10), ZI(11), ZX(12), ZY(13), ZZ(14)
+# Marginal indices when the second qubit is immune (surviving: first qubit, non-cross terms: XI, YI, ZI):
+_PC2_SECOND_IMMUNE_INDICES = (3, 7, 11)
+# Marginal indices when the first qubit is immune (surviving: second qubit, non-cross terms: IX, IY, IZ):
+_PC2_FIRST_IMMUNE_INDICES = (0, 1, 2)
+
+
+def _marginalize_2q_noise(
+    noise_op: stim.CircuitInstruction, immune_qubits: set[int], *, marginalize: bool
+) -> stim.Circuit:
+    """Filter or marginalize a 2-qubit noise instruction over immune qubits.
+
+    Processes each pair of targets independently.  Pairs with no immune qubits are kept as-is.
+    Pairs where both qubits are immune are dropped.  For partially-immune pairs, if marginalize is
+    True the surviving qubit receives a 1-qubit marginal (ignoring cross-Pauli terms); if False the
+    pair is dropped.  Only DEPOLARIZE2 and PAULI_CHANNEL_2 support marginalization; other 2-qubit
+    channels emit a warning and are dropped for partially-immune pairs.
+
+    Args:
+        noise_op: A 2-qubit noise instruction.
+        immune_qubits: Set of qubit indices that should not receive noise.
+        marginalize: If True, emit a 1-qubit marginal for partially-immune pairs.
+
+    Returns:
+        stim.Circuit: The filtered/marginalized circuit for this instruction.
+    """
+    result = stim.Circuit()
+    name = noise_op.name
+    args = noise_op.gate_args_copy()
+    targets = noise_op.targets_copy()
+
+    for i in range(0, len(targets), 2):
+        q1, q2 = targets[i], targets[i + 1]
+        q1_immune = q1.value in immune_qubits
+        q2_immune = q2.value in immune_qubits
+
+        if not q1_immune and not q2_immune:
+            result.append(stim.CircuitInstruction(name, [q1, q2], args))
+        elif q1_immune and q2_immune:
+            pass  # both immune: skip
+        elif not marginalize:
+            pass  # partially immune, no marginalization: drop the pair
+        elif name == "DEPOLARIZE2":
+            p = args[0]
+            surviving = q2.value if q1_immune else q1.value
+            result.append(stim.CircuitInstruction("DEPOLARIZE1", [surviving], [p / 5]))
+        elif name == "PAULI_CHANNEL_2":
+            if q2_immune:
+                marginal = [args[idx] for idx in _PC2_SECOND_IMMUNE_INDICES]
+                surviving = q1.value
+            else:
+                marginal = [args[idx] for idx in _PC2_FIRST_IMMUNE_INDICES]
+                surviving = q2.value
+            result.append(stim.CircuitInstruction("PAULI_CHANNEL_1", [surviving], marginal))
+        else:  # pragma: no cover
+            warnings.warn(
+                f"Cannot marginalize {name} over immune qubits; noise is dropped.",
+                stacklevel=2,
+            )
+
+    return result
+
+
+def immunize_noise(
+    noise: stim.Circuit, immune_qubits: set[int], *, marginalize: bool = False
+) -> stim.Circuit:
+    """Return a copy of a flat noise circuit with instructions targeting immune qubits removed.
+
+    An instruction is removed if any of its qubit targets belongs to immune_qubits.  If marginalize
+    is True, DEPOLARIZE2 and PAULI_CHANNEL_2 instructions where only one qubit in a pair is immune
+    are replaced by the marginal 1-qubit noise channel on the surviving qubit (ignoring cross-Pauli
+    terms).
+
+    Args:
+        noise: A flat noise circuit (no repeat blocks) to filter.
+        immune_qubits: Set of qubit indices that should not have noise applied to them.
+        marginalize: If True, marginalize 2-qubit noise onto surviving qubits instead of removing.
+
+    Returns:
+        stim.Circuit: A filtered copy of the input circuit.
+    """
+    if not immune_qubits:
+        return noise
+    result = stim.Circuit()
+    for noise_op in noise:
+        assert isinstance(noise_op, stim.CircuitInstruction)
+        if all(t.value not in immune_qubits for t in noise_op.targets_copy() if not t.is_combiner):
+            result.append(noise_op)
+        elif stim.gate_data(noise_op.name).is_two_qubit_gate:
+            result += _marginalize_2q_noise(noise_op, immune_qubits, marginalize=marginalize)
+        else:
+            # 1-qubit noise with multiple targets: keep only non-immune targets
+            surviving = [t for t in noise_op.targets_copy() if t.value not in immune_qubits]
+            if surviving:
+                result.append(
+                    stim.CircuitInstruction(noise_op.name, surviving, noise_op.gate_args_copy())
+                )
+    return result
+
+
 def _get_standardized_name(op: stim.CircuitInstruction) -> str:
     """Stardardized name of a circuit instruction.