From c33e65b2ca3b1ce988d33687be7c775fb2a5e4c6 Mon Sep 17 00:00:00 2001 From: Nandan Patel Date: Fri, 12 Jun 2026 16:36:56 +0800 Subject: [PATCH 01/11] Add amplitude damping noise model (#497) --- graphix/channels.py | 52 ++++++++ graphix/noise_models/__init__.py | 8 ++ graphix/noise_models/amplitude_damping.py | 153 ++++++++++++++++++++++ tests/test_density_matrix.py | 94 ++++++++++++- tests/test_kraus.py | 40 ++++++ tests/test_noise_model.py | 116 +++++++++++++++- 6 files changed, 461 insertions(+), 2 deletions(-) create mode 100644 graphix/noise_models/amplitude_damping.py diff --git a/graphix/channels.py b/graphix/channels.py index 8bb0ebc1b..9339a2dbf 100644 --- a/graphix/channels.py +++ b/graphix/channels.py @@ -296,3 +296,55 @@ def two_qubit_depolarising_tensor_channel(prob: float) -> KrausChannel: KrausData(prob / 3.0, np.kron(Ops.Z, Ops.Y)), ] ) +def amplitude_damping_channel(prob: float) -> KrausChannel: + r"""Single-qubit amplitude damping channel. + + .. math:: + K_1 = \begin{pmatrix} 1 & 0 \\ 0 & \sqrt{1-\gamma} \end{pmatrix}, \quad + K_2 = \begin{pmatrix} 0 & \sqrt{\gamma} \\ 0 & 0 \end{pmatrix} + + Parameters + ---------- + prob : float + The damping parameter :math:`\gamma` associated to the channel. + + Returns + ------- + :class:`graphix.channels.KrausChannel` object + containing the corresponding Kraus operators + """ + return KrausChannel( + [ + KrausData(1.0, np.array([[1.0, 0.0], [0.0, np.sqrt(1 - prob)]], dtype=np.complex128)), + KrausData(1.0, np.array([[0.0, np.sqrt(prob)], [0.0, 0.0]], dtype=np.complex128)), + ] + ) + + +def two_qubit_amplitude_damping_channel(prob: float) -> KrausChannel: + r"""Two-qubit amplitude damping channel. + + Tensor product of two independent single-qubit amplitude damping channels + with the same damping parameter :math:`\gamma`, giving the four Kraus + operators :math:`\{K_i \otimes K_j\}` for :math:`i, j \in \{1, 2\}`. + + Parameters + ---------- + prob : float + The damping parameter :math:`\gamma` associated to the channel. + + Returns + ------- + :class:`graphix.channels.KrausChannel` object + containing the corresponding Kraus operators + """ + k1 = np.array([[1.0, 0.0], [0.0, np.sqrt(1 - prob)]], dtype=np.complex128) + k2 = np.array([[0.0, np.sqrt(prob)], [0.0, 0.0]], dtype=np.complex128) + return KrausChannel( + [ + KrausData(1.0, np.kron(k1, k1)), + KrausData(1.0, np.kron(k1, k2)), + KrausData(1.0, np.kron(k2, k1)), + KrausData(1.0, np.kron(k2, k2)), + ] + ) diff --git a/graphix/noise_models/__init__.py b/graphix/noise_models/__init__.py index 1d74beaac..fb7a375e3 100644 --- a/graphix/noise_models/__init__.py +++ b/graphix/noise_models/__init__.py @@ -4,6 +4,11 @@ from typing import TYPE_CHECKING +from graphix.noise_models.amplitude_damping import ( + AmplitudeDampingNoise, + AmplitudeDampingNoiseModel, + TwoQubitAmplitudeDampingNoise, +) from graphix.noise_models.depolarising import DepolarisingNoise, DepolarisingNoiseModel, TwoQubitDepolarisingNoise from graphix.noise_models.noise_model import ( ApplyNoise, @@ -16,11 +21,14 @@ from graphix.noise_models.noise_model import CommandOrNoise as CommandOrNoise __all__ = [ + "AmplitudeDampingNoise", + "AmplitudeDampingNoiseModel", "ApplyNoise", "ComposeNoiseModel", "DepolarisingNoise", "DepolarisingNoiseModel", "Noise", "NoiseModel", + "TwoQubitAmplitudeDampingNoise", "TwoQubitDepolarisingNoise", ] diff --git a/graphix/noise_models/amplitude_damping.py b/graphix/noise_models/amplitude_damping.py new file mode 100644 index 000000000..b280019b3 --- /dev/null +++ b/graphix/noise_models/amplitude_damping.py @@ -0,0 +1,153 @@ +"""Amplitude damping noise model.""" + +from __future__ import annotations + +from typing import TYPE_CHECKING + +import typing_extensions + +from graphix.channels import ( + KrausChannel, + amplitude_damping_channel, + two_qubit_amplitude_damping_channel, +) +from graphix.command import BaseM, CommandKind +from graphix.noise_models.noise_model import ApplyNoise, Noise, NoiseModel +from graphix.utils import Probability + +if TYPE_CHECKING: + from collections.abc import Iterable + + from numpy.random import Generator + + from graphix.measurements import Outcome + from graphix.noise_models.noise_model import CommandOrNoise + + +class AmplitudeDampingNoise(Noise): + """One-qubit amplitude damping noise with damping parameter ``prob``.""" + + prob = Probability() + + def __init__(self, prob: float) -> None: + r"""Initialize one-qubit amplitude damping noise. + + Parameters + ---------- + prob : float + Damping parameter :math:`\\gamma` of the noise, between 0 and 1. + """ + self.prob = prob + + @property + @typing_extensions.override + def nqubits(self) -> int: + """Return the number of qubits targetted by the noise element.""" + return 1 + + @typing_extensions.override + def to_kraus_channel(self) -> KrausChannel: + """Return the Kraus channel describing the noise element.""" + return amplitude_damping_channel(self.prob) + + +class TwoQubitAmplitudeDampingNoise(Noise): + """Two-qubit amplitude damping noise with damping parameter ``prob``.""" + + prob = Probability() + + def __init__(self, prob: float) -> None: + r"""Initialize two-qubit amplitude damping noise. + + Parameters + ---------- + prob : float + Damping parameter :math:`\\gamma` of the noise, between 0 and 1. + """ + self.prob = prob + + @property + @typing_extensions.override + def nqubits(self) -> int: + """Return the number of qubits targetted by the noise element.""" + return 2 + + @typing_extensions.override + def to_kraus_channel(self) -> KrausChannel: + """Return the Kraus channel describing the noise element.""" + return two_qubit_amplitude_damping_channel(self.prob) + + +class AmplitudeDampingNoiseModel(NoiseModel): + r"""Amplitude damping noise model. + + :param NoiseModel: Parent abstract class class:`NoiseModel` + :type NoiseModel: class + """ + + def __init__( + self, + prepare_error_prob: float = 0.0, + x_error_prob: float = 0.0, + z_error_prob: float = 0.0, + entanglement_error_prob: float = 0.0, + measure_channel_prob: float = 0.0, + ) -> None: + self.prepare_error_prob = prepare_error_prob + self.x_error_prob = x_error_prob + self.z_error_prob = z_error_prob + self.entanglement_error_prob = entanglement_error_prob + self.measure_channel_prob = measure_channel_prob + + @typing_extensions.override + def input_nodes( + self, nodes: Iterable[int], rng: Generator | None = None, *, stacklevel: int = 1 + ) -> list[CommandOrNoise]: + """Return the noise to apply to input nodes.""" + return [ApplyNoise(noise=AmplitudeDampingNoise(self.prepare_error_prob), nodes=[node]) for node in nodes] + + @typing_extensions.override + def command( + self, cmd: CommandOrNoise, rng: Generator | None = None, *, stacklevel: int = 1 + ) -> list[CommandOrNoise]: + """Return the noise to apply to the command ``cmd``.""" + match cmd.kind: + case CommandKind.N: + return [cmd, ApplyNoise(noise=AmplitudeDampingNoise(self.prepare_error_prob), nodes=[cmd.node])] + case CommandKind.E: + return [ + cmd, + ApplyNoise( + noise=TwoQubitAmplitudeDampingNoise(self.entanglement_error_prob), nodes=list(cmd.nodes) + ), + ] + case CommandKind.M: + return [ApplyNoise(noise=AmplitudeDampingNoise(self.measure_channel_prob), nodes=[cmd.node]), cmd] + case CommandKind.X: + return [ + cmd, + ApplyNoise(noise=AmplitudeDampingNoise(self.x_error_prob), nodes=[cmd.node], domain=cmd.domain), + ] + case CommandKind.Z: + return [ + cmd, + ApplyNoise(noise=AmplitudeDampingNoise(self.z_error_prob), nodes=[cmd.node], domain=cmd.domain), + ] + case CommandKind.C | CommandKind.T | CommandKind.ApplyNoise: + return [cmd] + case CommandKind.S: + raise ValueError("Unexpected signal!") + case _: + typing_extensions.assert_never(cmd.kind) + + @typing_extensions.override + def confuse_result( + self, cmd: BaseM, result: Outcome, rng: Generator | None = None, *, stacklevel: int = 1 + ) -> Outcome: + r"""Return the measurement result unchanged. + + Amplitude damping is a purely quantum channel and introduces no + classical readout error. Compose this model with another noise + model (e.g. via :class:`ComposeNoiseModel`) to add readout error. + """ + return result diff --git a/tests/test_density_matrix.py b/tests/test_density_matrix.py index 7c67909c5..df0d5b3be 100644 --- a/tests/test_density_matrix.py +++ b/tests/test_density_matrix.py @@ -12,7 +12,13 @@ import graphix.random_objects as randobj from graphix import command from graphix.branch_selector import ConstBranchSelector -from graphix.channels import KrausChannel, dephasing_channel, depolarising_channel +from graphix.channels import ( + KrausChannel, + amplitude_damping_channel, + dephasing_channel, + depolarising_channel, + two_qubit_amplitude_damping_channel, +) from graphix.fundamentals import ANGLE_PI, Plane from graphix.ops import Ops from graphix.sim.density_matrix import DensityMatrix, DensityMatrixBackend @@ -736,6 +742,92 @@ def test_apply_depolarising_channel(self, fx_rng: Generator) -> None: assert np.allclose(expected_dm.trace(), 1.0) assert np.allclose(dm.rho, expected_dm) + def test_apply_amplitude_damping_channel(self, fx_rng: Generator) -> None: + # check on single qubit first, against the by-hand Kraus sum + dm = DensityMatrix(randobj.rand_dm(2, fx_rng)) + rho_test = dm.rho + + gamma = fx_rng.uniform() + ad_channel = amplitude_damping_channel(gamma) + + assert isinstance(ad_channel, KrausChannel) + + dm.apply_channel(ad_channel, [0]) + + k1 = np.array([[1.0, 0.0], [0.0, np.sqrt(1 - gamma)]], dtype=np.complex128) + k2 = np.array([[0.0, np.sqrt(gamma)], [0.0, 0.0]], dtype=np.complex128) + expected_dm = k1 @ rho_test @ k1.conj().T + k2 @ rho_test @ k2.conj().T + + assert np.allclose(expected_dm.trace(), 1.0) + assert np.allclose(dm.rho, expected_dm) + + # check embedded in a larger random register + nqubits = int(fx_rng.integers(2, 5)) + i = int(fx_rng.integers(0, nqubits)) + + psi = _randstate_raw(nqubits, fx_rng) + psi /= np.sqrt(np.sum(np.abs(psi) ** 2)) + dm = DensityMatrix(data=np.outer(psi, psi.conj())) + + gamma = fx_rng.uniform() + ad_channel = amplitude_damping_channel(gamma) + dm.apply_channel(ad_channel, [i]) + + expected_dm = np.zeros((2**nqubits, 2**nqubits), dtype=np.complex128) + for elem in ad_channel: + psi_evolved = np.tensordot(elem.operator, psi.reshape((2,) * nqubits), (1, i)) + psi_evolved = np.moveaxis(psi_evolved, 0, i).reshape(2**nqubits) + expected_dm += elem.coef * np.conj(elem.coef) * np.outer(psi_evolved, psi_evolved.conj()) + + assert np.allclose(expected_dm.trace(), 1.0) + assert np.allclose(dm.rho, expected_dm) + + @pytest.mark.parametrize("gamma", [0.0, 0.2, 0.5, 0.9, 1.0]) + def test_amplitude_damping_ground_state_fixed(self, gamma: float) -> None: + #``|0><0|`` is a fixed point of amplitude damping for any gamma. + ket0 = np.array([[1.0, 0.0], [0.0, 0.0]], dtype=np.complex128) + dm = DensityMatrix(data=BasicStates.ZERO) + dm.apply_channel(amplitude_damping_channel(gamma), [0]) + assert np.allclose(dm.rho, ket0) + + @pytest.mark.parametrize("gamma", [0.0, 0.2, 0.5, 0.9, 1.0]) + def test_amplitude_damping_excited_state_decays(self, gamma: float) -> None: + #``|1><1| -> (1 - gamma)|1><1| + gamma|0><0|`` (directional T1 decay). + ket0 = np.array([[1.0, 0.0], [0.0, 0.0]], dtype=np.complex128) + ket1 = np.array([[0.0, 0.0], [0.0, 1.0]], dtype=np.complex128) + dm = DensityMatrix(data=BasicStates.ONE) + dm.apply_channel(amplitude_damping_channel(gamma), [0]) + assert np.allclose(dm.rho, (1 - gamma) * ket1 + gamma * ket0) + + @pytest.mark.parametrize("gamma", [0.1, 0.4, 0.8]) + def test_amplitude_damping_coherence_decay(self, gamma: float) -> None: + #Off-diagonal coherences scale by ``sqrt(1 - gamma)`` (distinct from dephasing). + dm = DensityMatrix(data=BasicStates.PLUS) + dm.apply_channel(amplitude_damping_channel(gamma), [0]) + assert np.isclose(dm.rho[0, 1], np.sqrt(1 - gamma) / 2) + assert np.isclose(dm.rho[1, 0], np.sqrt(1 - gamma) / 2) + + @pytest.mark.parametrize("gamma", [0.0, 0.25, 0.6, 1.0]) + def test_apply_two_qubit_amplitude_damping_channel(self, gamma: float, fx_rng: Generator) -> None: + #The two-qubit channel equals independent damping on each factor. + a = _randstate_raw(1, fx_rng) + a /= np.sqrt(np.sum(np.abs(a) ** 2)) + b = _randstate_raw(1, fx_rng) + b /= np.sqrt(np.sum(np.abs(b) ** 2)) + rho_a = np.outer(a, a.conj()) + rho_b = np.outer(b, b.conj()) + + dm = DensityMatrix(data=np.kron(rho_a, rho_b)) + dm.apply_channel(two_qubit_amplitude_damping_channel(gamma), [0, 1]) + + k1 = np.array([[1.0, 0.0], [0.0, np.sqrt(1 - gamma)]], dtype=np.complex128) + k2 = np.array([[0.0, np.sqrt(gamma)], [0.0, 0.0]], dtype=np.complex128) + single = lambda rho: k1 @ rho @ k1.conj().T + k2 @ rho @ k2.conj().T # noqa: E731 + expected = np.kron(single(rho_a), single(rho_b)) + + assert np.allclose(expected.trace(), 1.0) + assert np.allclose(dm.rho, expected) + def test_apply_random_channel_one_qubit(self, fx_rng: Generator) -> None: """Test using complex parameters.""" # check against statevector backend by hand for now. diff --git a/tests/test_kraus.py b/tests/test_kraus.py index f7e4983f0..8a6b9fd0f 100644 --- a/tests/test_kraus.py +++ b/tests/test_kraus.py @@ -9,8 +9,10 @@ from graphix.channels import ( KrausChannel, KrausData, + amplitude_damping_channel, dephasing_channel, depolarising_channel, + two_qubit_amplitude_damping_channel, two_qubit_depolarising_channel, two_qubit_depolarising_tensor_channel, ) @@ -180,3 +182,41 @@ def test_2_qubit_depolarising_tensor_channel(self, fx_rng: Generator) -> None: for i in range(len(depol_tensor_channel_2_qubit)): assert np.allclose(depol_tensor_channel_2_qubit[i].coef, data[i].coef) assert np.allclose(depol_tensor_channel_2_qubit[i].operator, data[i].operator) + + def test_amplitude_damping_channel(self, fx_rng: Generator) -> None: + gamma = fx_rng.uniform() + data = [ + KrausData(1.0, np.array([[1.0, 0.0], [0.0, np.sqrt(1 - gamma)]], dtype=np.complex128)), + KrausData(1.0, np.array([[0.0, np.sqrt(gamma)], [0.0, 0.0]], dtype=np.complex128)), + ] + + ad_channel = amplitude_damping_channel(gamma) + + assert isinstance(ad_channel, KrausChannel) + assert ad_channel.nqubit == 1 + assert len(ad_channel) == 2 + + for i in range(len(ad_channel)): + assert np.allclose(ad_channel[i].coef, data[i].coef) + assert np.allclose(ad_channel[i].operator, data[i].operator) + + def test_2_qubit_amplitude_damping_channel(self, fx_rng: Generator) -> None: + gamma = fx_rng.uniform() + k1 = np.array([[1.0, 0.0], [0.0, np.sqrt(1 - gamma)]], dtype=np.complex128) + k2 = np.array([[0.0, np.sqrt(gamma)], [0.0, 0.0]], dtype=np.complex128) + data = [ + KrausData(1.0, np.kron(k1, k1)), + KrausData(1.0, np.kron(k1, k2)), + KrausData(1.0, np.kron(k2, k1)), + KrausData(1.0, np.kron(k2, k2)), + ] + + ad_channel_2_qubit = two_qubit_amplitude_damping_channel(gamma) + + assert isinstance(ad_channel_2_qubit, KrausChannel) + assert ad_channel_2_qubit.nqubit == 2 + assert len(ad_channel_2_qubit) == 4 + + for i in range(len(ad_channel_2_qubit)): + assert np.allclose(ad_channel_2_qubit[i].coef, data[i].coef) + assert np.allclose(ad_channel_2_qubit[i].operator, data[i].operator) diff --git a/tests/test_noise_model.py b/tests/test_noise_model.py index 3babd2a4c..bccbad640 100644 --- a/tests/test_noise_model.py +++ b/tests/test_noise_model.py @@ -6,12 +6,15 @@ import pytest from graphix import Pattern -from graphix.command import CommandKind, M, N +from graphix.command import CommandKind, E, M, N, X from graphix.noise_models import ( + AmplitudeDampingNoise, + AmplitudeDampingNoiseModel, ApplyNoise, ComposeNoiseModel, DepolarisingNoise, DepolarisingNoiseModel, + TwoQubitAmplitudeDampingNoise, TwoQubitDepolarisingNoise, ) from graphix.noise_models.noise_model import NoiselessNoiseModel @@ -114,3 +117,114 @@ def test_confuse_result(fx_rng: Generator) -> None: backend="densitymatrix", noise_model=noise_model, rng=fx_rng, measure_method=measure_method ) assert measure_method.results[0] == 1 + + +def test_amplitude_damping_command_injection() -> None: + """Amplitude damping noise is injected at the correct command positions.""" + model = AmplitudeDampingNoiseModel( + prepare_error_prob=0.1, + x_error_prob=0.2, + entanglement_error_prob=0.3, + measure_channel_prob=0.4, + ) + + # N: noise applied AFTER preparation + out = model.command(N(node=0)) + assert len(out) == 2 + assert out[0].kind == CommandKind.N + assert isinstance(out[1], ApplyNoise) + assert isinstance(out[1].noise, AmplitudeDampingNoise) + assert out[1].nodes == [0] + + # E: two-qubit noise applied AFTER entanglement + out = model.command(E(nodes=(0, 1))) + assert out[0].kind == CommandKind.E + assert isinstance(out[1], ApplyNoise) + assert isinstance(out[1].noise, TwoQubitAmplitudeDampingNoise) + assert out[1].noise.nqubits == 2 + + # M: noise applied BEFORE measurement + out = model.command(M(node=0)) + assert isinstance(out[0], ApplyNoise) + assert isinstance(out[0].noise, AmplitudeDampingNoise) + assert out[1].kind == CommandKind.M + + # X: correction kept, noise conditioned on the same domain + out = model.command(X(node=0, domain={1, 2})) + assert out[0].kind == CommandKind.X + assert isinstance(out[1], ApplyNoise) + assert out[1].domain == {1, 2} + + +@pytest.mark.parametrize("outcome", [0, 1]) +def test_amplitude_damping_confuse_result_is_identity(outcome: int) -> None: + """Amplitude damping introduces no classical readout error.""" + model = AmplitudeDampingNoiseModel() + assert model.confuse_result(M(node=0), outcome) == outcome + + +def test_compose_amplitude_damping_depolarising_transpile(fx_rng: Generator) -> None: + #Compose an amplitude damping and a depolarising model, and verify that each composed command injects the two models' noise in order. + nqubits = 5 + depth = 5 + circuit = rand_circuit(nqubits, depth, rng=fx_rng) + pattern = circuit.transpile().pattern + noise_model = ComposeNoiseModel( + [AmplitudeDampingNoiseModel(x_error_prob=0.5), DepolarisingNoiseModel(z_error_prob=0.5)] + ) + noisy_pattern = noise_model.transpile(pattern, rng=fx_rng) + iterator = iter(noisy_pattern) + + def check_ad(cmd: CommandOrNoise, prob: float, two_qubits: bool) -> None: + assert isinstance(cmd, ApplyNoise) + if two_qubits: + assert isinstance(cmd.noise, TwoQubitAmplitudeDampingNoise) + else: + assert isinstance(cmd.noise, AmplitudeDampingNoise) + assert cmd.noise.prob == prob + + def check_depol(cmd: CommandOrNoise, prob: float, two_qubits: bool) -> None: + assert isinstance(cmd, ApplyNoise) + if two_qubits: + assert isinstance(cmd.noise, TwoQubitDepolarisingNoise) + else: + assert isinstance(cmd.noise, DepolarisingNoise) + assert cmd.noise.prob == prob + + for cmd in pattern: + if cmd.kind == CommandKind.M: + # measurement-channel noise injected BEFORE M; prepend reverses the + # order relative to appended slots, so AD precedes depol here + check_ad(next(iterator), 0, False) + check_depol(next(iterator), 0, False) + assert next(iterator) == cmd + match cmd.kind: + case CommandKind.N: + check_depol(next(iterator), 0, False) + check_ad(next(iterator), 0, False) + case CommandKind.E: + check_depol(next(iterator), 0, True) + check_ad(next(iterator), 0, True) + case CommandKind.X: + # depol carries 0 on X, AD carries x_error_prob=0.5 + check_depol(next(iterator), 0, False) + check_ad(next(iterator), 0.5, False) + case CommandKind.Z: + # depol carries z_error_prob=0.5, AD carries 0 on Z + check_depol(next(iterator), 0.5, False) + check_ad(next(iterator), 0, False) + + +def test_compose_amplitude_damping_depolarising_simulation(fx_rng: Generator) -> None: + """A composed model with both noiseless-configured models reproduces the ideal state.""" + nqubits = 5 + depth = 5 + circuit = rand_circuit(nqubits, depth, rng=fx_rng) + state = circuit.simulate_statevector().statevec + pattern = circuit.transpile().pattern + pattern.standardize() + pattern.minimize_space() + # both models default to noiseless (all probs 0) + noise_model = ComposeNoiseModel([AmplitudeDampingNoiseModel(), DepolarisingNoiseModel()]) + state_mbqc = pattern.simulate_pattern(backend="densitymatrix", noise_model=noise_model, rng=fx_rng) + assert np.abs(np.dot(state_mbqc.flatten().conjugate(), DensityMatrix(state).rho.flatten())) == pytest.approx(1) From bec57de45f2456d36476e793a39221a8315c1090 Mon Sep 17 00:00:00 2001 From: Nandan Patel Date: Mon, 15 Jun 2026 16:21:41 +0800 Subject: [PATCH 02/11] Apply ruff format --- graphix/channels.py | 2 ++ tests/test_density_matrix.py | 16 ++++++++++------ tests/test_noise_model.py | 5 +++-- 3 files changed, 15 insertions(+), 8 deletions(-) diff --git a/graphix/channels.py b/graphix/channels.py index 9339a2dbf..8b3204044 100644 --- a/graphix/channels.py +++ b/graphix/channels.py @@ -296,6 +296,8 @@ def two_qubit_depolarising_tensor_channel(prob: float) -> KrausChannel: KrausData(prob / 3.0, np.kron(Ops.Z, Ops.Y)), ] ) + + def amplitude_damping_channel(prob: float) -> KrausChannel: r"""Single-qubit amplitude damping channel. diff --git a/tests/test_density_matrix.py b/tests/test_density_matrix.py index df0d5b3be..ba6f2eb58 100644 --- a/tests/test_density_matrix.py +++ b/tests/test_density_matrix.py @@ -8,6 +8,7 @@ import numpy as np import numpy.typing as npt import pytest +from scipy.special import k1 import graphix.random_objects as randobj from graphix import command @@ -745,7 +746,7 @@ def test_apply_depolarising_channel(self, fx_rng: Generator) -> None: def test_apply_amplitude_damping_channel(self, fx_rng: Generator) -> None: # check on single qubit first, against the by-hand Kraus sum dm = DensityMatrix(randobj.rand_dm(2, fx_rng)) - rho_test = dm.rho + rho_test = np.asarray(dm.rho, dtype=np.complex128) gamma = fx_rng.uniform() ad_channel = amplitude_damping_channel(gamma) @@ -784,7 +785,7 @@ def test_apply_amplitude_damping_channel(self, fx_rng: Generator) -> None: @pytest.mark.parametrize("gamma", [0.0, 0.2, 0.5, 0.9, 1.0]) def test_amplitude_damping_ground_state_fixed(self, gamma: float) -> None: - #``|0><0|`` is a fixed point of amplitude damping for any gamma. + # ``|0><0|`` is a fixed point of amplitude damping for any gamma. ket0 = np.array([[1.0, 0.0], [0.0, 0.0]], dtype=np.complex128) dm = DensityMatrix(data=BasicStates.ZERO) dm.apply_channel(amplitude_damping_channel(gamma), [0]) @@ -792,7 +793,7 @@ def test_amplitude_damping_ground_state_fixed(self, gamma: float) -> None: @pytest.mark.parametrize("gamma", [0.0, 0.2, 0.5, 0.9, 1.0]) def test_amplitude_damping_excited_state_decays(self, gamma: float) -> None: - #``|1><1| -> (1 - gamma)|1><1| + gamma|0><0|`` (directional T1 decay). + # ``|1><1| -> (1 - gamma)|1><1| + gamma|0><0|`` (directional T1 decay). ket0 = np.array([[1.0, 0.0], [0.0, 0.0]], dtype=np.complex128) ket1 = np.array([[0.0, 0.0], [0.0, 1.0]], dtype=np.complex128) dm = DensityMatrix(data=BasicStates.ONE) @@ -801,7 +802,7 @@ def test_amplitude_damping_excited_state_decays(self, gamma: float) -> None: @pytest.mark.parametrize("gamma", [0.1, 0.4, 0.8]) def test_amplitude_damping_coherence_decay(self, gamma: float) -> None: - #Off-diagonal coherences scale by ``sqrt(1 - gamma)`` (distinct from dephasing). + # Off-diagonal coherences scale by ``sqrt(1 - gamma)`` (distinct from dephasing). dm = DensityMatrix(data=BasicStates.PLUS) dm.apply_channel(amplitude_damping_channel(gamma), [0]) assert np.isclose(dm.rho[0, 1], np.sqrt(1 - gamma) / 2) @@ -809,7 +810,7 @@ def test_amplitude_damping_coherence_decay(self, gamma: float) -> None: @pytest.mark.parametrize("gamma", [0.0, 0.25, 0.6, 1.0]) def test_apply_two_qubit_amplitude_damping_channel(self, gamma: float, fx_rng: Generator) -> None: - #The two-qubit channel equals independent damping on each factor. + # The two-qubit channel equals independent damping on each factor. a = _randstate_raw(1, fx_rng) a /= np.sqrt(np.sum(np.abs(a) ** 2)) b = _randstate_raw(1, fx_rng) @@ -822,7 +823,10 @@ def test_apply_two_qubit_amplitude_damping_channel(self, gamma: float, fx_rng: G k1 = np.array([[1.0, 0.0], [0.0, np.sqrt(1 - gamma)]], dtype=np.complex128) k2 = np.array([[0.0, np.sqrt(gamma)], [0.0, 0.0]], dtype=np.complex128) - single = lambda rho: k1 @ rho @ k1.conj().T + k2 @ rho @ k2.conj().T # noqa: E731 + + def single(rho: npt.NDArray[np.complex128]) -> npt.NDArray[np.complex128]: + return k1 @ rho @ k1.conj().T + k2 @ rho @ k2.conj().T + expected = np.kron(single(rho_a), single(rho_b)) assert np.allclose(expected.trace(), 1.0) diff --git a/tests/test_noise_model.py b/tests/test_noise_model.py index bccbad640..384f9b547 100644 --- a/tests/test_noise_model.py +++ b/tests/test_noise_model.py @@ -18,6 +18,7 @@ TwoQubitDepolarisingNoise, ) from graphix.noise_models.noise_model import NoiselessNoiseModel +from graphix.measurements import Outcome from graphix.random_objects import rand_circuit from graphix.sim.density_matrix import DensityMatrix from graphix.simulator import DefaultMeasureMethod @@ -157,14 +158,14 @@ def test_amplitude_damping_command_injection() -> None: @pytest.mark.parametrize("outcome", [0, 1]) -def test_amplitude_damping_confuse_result_is_identity(outcome: int) -> None: +def test_amplitude_damping_confuse_result_is_identity(outcome: Outcome) -> None: """Amplitude damping introduces no classical readout error.""" model = AmplitudeDampingNoiseModel() assert model.confuse_result(M(node=0), outcome) == outcome def test_compose_amplitude_damping_depolarising_transpile(fx_rng: Generator) -> None: - #Compose an amplitude damping and a depolarising model, and verify that each composed command injects the two models' noise in order. + # Compose an amplitude damping and a depolarising model, and verify that each composed command injects the two models' noise in order. nqubits = 5 depth = 5 circuit = rand_circuit(nqubits, depth, rng=fx_rng) From 2dc824d781896ce158ecba3fd862c5bb11df9f52 Mon Sep 17 00:00:00 2001 From: Nandan Patel Date: Mon, 15 Jun 2026 17:56:30 +0800 Subject: [PATCH 03/11] Add per-step analytic tests and fix docstring raw-string lint --- tests/test_noisy_density_matrix.py | 145 ++++++++++++++++++++++++++++- 1 file changed, 144 insertions(+), 1 deletion(-) diff --git a/tests/test_noisy_density_matrix.py b/tests/test_noisy_density_matrix.py index 104923265..24ce8b069 100644 --- a/tests/test_noisy_density_matrix.py +++ b/tests/test_noisy_density_matrix.py @@ -9,7 +9,7 @@ from graphix.branch_selector import ConstBranchSelector, FixedBranchSelector from graphix.command import CommandKind from graphix.fundamentals import angle_to_rad -from graphix.noise_models import DepolarisingNoiseModel +from graphix.noise_models import AmplitudeDampingNoiseModel, ComposeNoiseModel, DepolarisingNoiseModel from graphix.noise_models.noise_model import NoiselessNoiseModel from graphix.ops import Ops from graphix.sim.density_matrix import DensityMatrix @@ -524,3 +524,146 @@ def test_noisy_measure_confuse_rz_arbitrary( or np.allclose(res.rho, Ops.Z @ exact @ Ops.Z) or np.allclose(res.rho, Ops.Z @ Ops.X @ exact @ Ops.X @ Ops.Z) ) + + +class TestNoisyDensityMatrixAmplitudeDamping: + """Amplitude damping noise model on the density matrix backend.""" + + @pytest.mark.filterwarnings("ignore:Simulating using densitymatrix backend with no noise.") + def test_noiseless_amplitude_damping_hadamard(self, fx_rng: Generator) -> None: + """A zero-parameter amplitude damping model reproduces the noiseless result.""" + + hadamardpattern = hpat() + noiselessres = hadamardpattern.simulate_pattern(backend="densitymatrix", rng=fx_rng) + noisynoiselessres = hadamardpattern.simulate_pattern( + backend="densitymatrix", + noise_model=AmplitudeDampingNoiseModel(), + rng=fx_rng, + ) + assert isinstance(noiselessres, DensityMatrix) + assert isinstance(noisynoiselessres, DensityMatrix) + assert np.allclose(noiselessres.rho, np.array([[1.0, 0.0], [0.0, 0.0]])) + assert np.allclose(noisynoiselessres.rho, np.array([[1.0, 0.0], [0.0, 0.0]])) + + @pytest.mark.filterwarnings("ignore:Simulating using densitymatrix backend with no noise.") + def test_noiseless_amplitude_damping_rz(self, fx_rng: Generator) -> None: + + alpha = fx_rng.random() + rzpattern = rzpat(alpha) + noiselessres = rzpattern.simulate_pattern(backend="densitymatrix", rng=fx_rng) + noisynoiselessres = rzpattern.simulate_pattern( + backend="densitymatrix", + noise_model=AmplitudeDampingNoiseModel(), + rng=fx_rng, + ) + assert isinstance(noiselessres, DensityMatrix) + assert isinstance(noisynoiselessres, DensityMatrix) + assert np.allclose(noiselessres.rho, rz_exact_res(alpha)) + assert np.allclose(noisynoiselessres.rho, rz_exact_res(alpha)) + + def test_compose_amplitude_damping_with_depolarising_runs(self, fx_rng: Generator) -> None: + """ComposeNoiseModel of depolarising + amplitude damping runs and stays normalized.""" + + hadamardpattern = hpat() + depol = DepolarisingNoiseModel(prepare_error_prob=0.2) + ad = AmplitudeDampingNoiseModel(prepare_error_prob=0.2) + composed = ComposeNoiseModel([depol, ad]) + + r_depol = hadamardpattern.simulate_pattern(backend="densitymatrix", noise_model=depol, rng=fx_rng) + r_both = hadamardpattern.simulate_pattern(backend="densitymatrix", noise_model=composed, rng=fx_rng) + assert isinstance(r_depol, DensityMatrix) + assert isinstance(r_both, DensityMatrix) + assert np.isclose(r_both.rho.trace(), 1.0) + assert not np.allclose(r_depol.rho, r_both.rho) + + +class TestAmplitudeDampingAnalytic: + """Analytic per-step verification of amplitude damping on the Hadamard pattern.""" + + @staticmethod + def _kraus(gamma: float) -> list[npt.NDArray[np.complex128]]: + return [ + np.array([[1.0, 0.0], [0.0, np.sqrt(1 - gamma)]], dtype=np.complex128), + np.array([[0.0, np.sqrt(gamma)], [0.0, 0.0]], dtype=np.complex128), + ] + + def _ad_single(self, rho: npt.NDArray[np.complex128], gamma: float) -> npt.NDArray[np.complex128]: + return sum((k @ rho @ k.conj().T for k in self._kraus(gamma)), np.zeros((2, 2), dtype=np.complex128)) + + def _ad_on(self, rho: npt.NDArray[np.complex128], qubit: int, gamma: float) -> npt.NDArray[np.complex128]: + eye = np.eye(2, dtype=np.complex128) + out = np.zeros((4, 4), dtype=np.complex128) + for k in self._kraus(gamma): + full = np.kron(k, eye) if qubit == 0 else np.kron(eye, k) + out += full @ rho @ full.conj().T + return out + + def _ad_two(self, rho: npt.NDArray[np.complex128], gamma: float) -> npt.NDArray[np.complex128]: + out = np.zeros((4, 4), dtype=np.complex128) + for left in self._kraus(gamma): + for right in self._kraus(gamma): + full = np.kron(left, right) + out += full @ rho @ full.conj().T + return out + + @staticmethod + def _proj_x(outcome: Outcome) -> npt.NDArray[np.complex128]: + vec = np.array([1.0, 1.0 if outcome == 0 else -1.0], dtype=np.complex128) / np.sqrt(2) + return np.outer(vec, vec.conj()) + + def _expected(self, step: str, gamma: float, outcome: Outcome) -> npt.NDArray[np.complex128]: + eye = np.eye(2, dtype=np.complex128) + pauli_x = np.array([[0.0, 1.0], [1.0, 0.0]], dtype=np.complex128) + cz = np.diag([1.0, 1.0, 1.0, -1.0]).astype(np.complex128) + plus = np.array([1.0, 1.0], dtype=np.complex128) / np.sqrt(2) + + rho = np.kron(np.outer(plus, plus.conj()), np.outer(plus, plus.conj())) + if step == "prep": + rho = self._ad_on(rho, 0, gamma) + rho = self._ad_on(rho, 1, gamma) + rho = cz @ rho @ cz.conj().T + if step == "entangle": + rho = self._ad_two(rho, gamma) + if step == "measure": + rho = self._ad_on(rho, 0, gamma) + proj = np.kron(self._proj_x(outcome), eye) + rho = proj @ rho @ proj.conj().T + rho = rho / np.trace(rho) + reduced: npt.NDArray[np.complex128] = np.einsum("ijik->jk", rho.reshape(2, 2, 2, 2)) + if outcome == 1: + reduced = pauli_x @ reduced @ pauli_x.conj().T + if step == "xcorr" and outcome == 1: + reduced = self._ad_single(reduced, gamma) + return reduced + + @pytest.mark.parametrize("outcome", [0, 1]) + @pytest.mark.parametrize("gamma", [0.0, 0.3, 0.7, 1.0]) + @pytest.mark.parametrize( + ("step", "param"), + [ + ("prep", "prepare_error_prob"), + ("entangle", "entanglement_error_prob"), + ("measure", "measure_channel_prob"), + ("xcorr", "x_error_prob"), + ], + ) + def test_amplitude_damping_step_matches_analytic( + self, step: str, param: str, gamma: float, outcome: Outcome, fx_rng: Generator + ) -> None: + res = hpat().simulate_pattern( + backend="densitymatrix", + noise_model=AmplitudeDampingNoiseModel(**{param: gamma}), + branch_selector=ConstBranchSelector(outcome), + rng=fx_rng, + ) + assert isinstance(res, DensityMatrix) + assert np.allclose(res.rho, self._expected(step, gamma, outcome)) + + # For the measurement step the result has a clean branch-independent + # closed form: damping the qubit before the X-basis measurement leaves + # the state diagonal with the coherence factor ``sqrt(1 - gamma)`` + # showing through. Assert it explicitly as a human-readable cross-check. + if step == "measure": + sqrt_term = np.sqrt(1 - gamma) + closed_form = np.array([[(1 + sqrt_term) / 2, 0.0], [0.0, (1 - sqrt_term) / 2]], dtype=np.complex128) + assert np.allclose(res.rho, closed_form) From af099714521810c7cb52a955274d0e26ef5c30c9 Mon Sep 17 00:00:00 2001 From: Nandan Patel Date: Mon, 15 Jun 2026 23:18:59 +0800 Subject: [PATCH 04/11] Fix lint and type errors in amplitude damping tests --- tests/test_density_matrix.py | 5 ++--- tests/test_noise_model.py | 2 +- tests/test_noisy_density_matrix.py | 7 ++++--- 3 files changed, 7 insertions(+), 7 deletions(-) diff --git a/tests/test_density_matrix.py b/tests/test_density_matrix.py index ba6f2eb58..25bd94d95 100644 --- a/tests/test_density_matrix.py +++ b/tests/test_density_matrix.py @@ -8,7 +8,6 @@ import numpy as np import numpy.typing as npt import pytest -from scipy.special import k1 import graphix.random_objects as randobj from graphix import command @@ -576,7 +575,7 @@ def test_apply_dephasing_channel(self, fx_rng: Generator) -> None: dm = DensityMatrix(randobj.rand_dm(2, fx_rng)) # copy of initial dm - rho_test = dm.rho + rho_test = np.asarray(dm.rho, dtype=np.complex128) # create dephasing channel prob = fx_rng.uniform() @@ -654,7 +653,7 @@ def test_apply_depolarising_channel(self, fx_rng: Generator) -> None: dm = DensityMatrix(randobj.rand_dm(2, fx_rng)) # copy of initial dm - rho_test = dm.rho + rho_test = np.asarray(dm.rho, dtype=np.complex128) # create dephasing channel prob = fx_rng.uniform() diff --git a/tests/test_noise_model.py b/tests/test_noise_model.py index 384f9b547..d86dd6d67 100644 --- a/tests/test_noise_model.py +++ b/tests/test_noise_model.py @@ -18,7 +18,6 @@ TwoQubitDepolarisingNoise, ) from graphix.noise_models.noise_model import NoiselessNoiseModel -from graphix.measurements import Outcome from graphix.random_objects import rand_circuit from graphix.sim.density_matrix import DensityMatrix from graphix.simulator import DefaultMeasureMethod @@ -26,6 +25,7 @@ if TYPE_CHECKING: from numpy.random import Generator + from graphix.measurements import Outcome from graphix.noise_models import CommandOrNoise diff --git a/tests/test_noisy_density_matrix.py b/tests/test_noisy_density_matrix.py index 24ce8b069..e2c4dd799 100644 --- a/tests/test_noisy_density_matrix.py +++ b/tests/test_noisy_density_matrix.py @@ -532,7 +532,6 @@ class TestNoisyDensityMatrixAmplitudeDamping: @pytest.mark.filterwarnings("ignore:Simulating using densitymatrix backend with no noise.") def test_noiseless_amplitude_damping_hadamard(self, fx_rng: Generator) -> None: """A zero-parameter amplitude damping model reproduces the noiseless result.""" - hadamardpattern = hpat() noiselessres = hadamardpattern.simulate_pattern(backend="densitymatrix", rng=fx_rng) noisynoiselessres = hadamardpattern.simulate_pattern( @@ -563,7 +562,6 @@ def test_noiseless_amplitude_damping_rz(self, fx_rng: Generator) -> None: def test_compose_amplitude_damping_with_depolarising_runs(self, fx_rng: Generator) -> None: """ComposeNoiseModel of depolarising + amplitude damping runs and stays normalized.""" - hadamardpattern = hpat() depol = DepolarisingNoiseModel(prepare_error_prob=0.2) ad = AmplitudeDampingNoiseModel(prepare_error_prob=0.2) @@ -617,7 +615,10 @@ def _expected(self, step: str, gamma: float, outcome: Outcome) -> npt.NDArray[np cz = np.diag([1.0, 1.0, 1.0, -1.0]).astype(np.complex128) plus = np.array([1.0, 1.0], dtype=np.complex128) / np.sqrt(2) - rho = np.kron(np.outer(plus, plus.conj()), np.outer(plus, plus.conj())) + rho: npt.NDArray[np.complex128] = np.kron(np.outer(plus, plus.conj()), np.outer(plus, plus.conj())).astype( + np.complex128 + ) + if step == "prep": rho = self._ad_on(rho, 0, gamma) rho = self._ad_on(rho, 1, gamma) From dda1110f11c20bf6ca1c5799f652c17da9495fd5 Mon Sep 17 00:00:00 2001 From: Nandan Patel Date: Mon, 15 Jun 2026 23:43:59 +0800 Subject: [PATCH 05/11] Add tests for amplitude damping noise model coverage --- tests/test_noise_model.py | 30 +++++++++++++++++++++++++++++- 1 file changed, 29 insertions(+), 1 deletion(-) diff --git a/tests/test_noise_model.py b/tests/test_noise_model.py index d86dd6d67..aa2573136 100644 --- a/tests/test_noise_model.py +++ b/tests/test_noise_model.py @@ -6,7 +6,8 @@ import pytest from graphix import Pattern -from graphix.command import CommandKind, E, M, N, X +from graphix.command import C, CommandKind, E, M, N, S, T, X +from graphix.clifford import Clifford from graphix.noise_models import ( AmplitudeDampingNoise, AmplitudeDampingNoiseModel, @@ -229,3 +230,30 @@ def test_compose_amplitude_damping_depolarising_simulation(fx_rng: Generator) -> noise_model = ComposeNoiseModel([AmplitudeDampingNoiseModel(), DepolarisingNoiseModel()]) state_mbqc = pattern.simulate_pattern(backend="densitymatrix", noise_model=noise_model, rng=fx_rng) assert np.abs(np.dot(state_mbqc.flatten().conjugate(), DensityMatrix(state).rho.flatten())) == pytest.approx(1) + +def test_amplitude_damping_nqubits() -> None: + # covers line 46: AmplitudeDampingNoise.nqubits returning 1 + assert AmplitudeDampingNoise(0.1).nqubits == 1 + assert TwoQubitAmplitudeDampingNoise(0.1).nqubits == 2 + + +def test_amplitude_damping_command_passthrough() -> None: + # covers line 136: the C / T / ApplyNoise case returns [cmd] unchanged + model = AmplitudeDampingNoiseModel() + + c_cmd = C(node=0, clifford=Clifford.H) + assert model.command(c_cmd) == [c_cmd] + + t_cmd = T() + assert model.command(t_cmd) == [t_cmd] + + apply_noise = ApplyNoise(noise=AmplitudeDampingNoise(0.1), nodes=[0]) + assert model.command(apply_noise) == [apply_noise] + + +def test_amplitude_damping_command_rejects_signal() -> None: + # covers lines 138-139: the CommandKind.S case raises ValueError + model = AmplitudeDampingNoiseModel() + s_cmd = S(node=0) + with pytest.raises(ValueError, match="Unexpected signal"): + model.command(s_cmd) From 54cfebe111d5c9b8be50eded58faea90ea372993 Mon Sep 17 00:00:00 2001 From: Nandan Patel Date: Wed, 17 Jun 2026 10:07:04 +0800 Subject: [PATCH 06/11] Add measure_error_prob, RZ analytic tests, and branch coverage for amplitude damping --- graphix/noise_models/amplitude_damping.py | 16 +- tests/test_noise_model.py | 44 ++--- tests/test_noisy_density_matrix.py | 231 +++++++++++++--------- 3 files changed, 170 insertions(+), 121 deletions(-) diff --git a/graphix/noise_models/amplitude_damping.py b/graphix/noise_models/amplitude_damping.py index b280019b3..ea23f9311 100644 --- a/graphix/noise_models/amplitude_damping.py +++ b/graphix/noise_models/amplitude_damping.py @@ -12,7 +12,9 @@ two_qubit_amplitude_damping_channel, ) from graphix.command import BaseM, CommandKind +from graphix.measurements import toggle_outcome from graphix.noise_models.noise_model import ApplyNoise, Noise, NoiseModel +from graphix.rng import ensure_rng from graphix.utils import Probability if TYPE_CHECKING: @@ -92,12 +94,14 @@ def __init__( z_error_prob: float = 0.0, entanglement_error_prob: float = 0.0, measure_channel_prob: float = 0.0, + measure_error_prob: float = 0.0, ) -> None: self.prepare_error_prob = prepare_error_prob self.x_error_prob = x_error_prob self.z_error_prob = z_error_prob self.entanglement_error_prob = entanglement_error_prob self.measure_channel_prob = measure_channel_prob + self.measure_error_prob = measure_error_prob @typing_extensions.override def input_nodes( @@ -137,17 +141,15 @@ def command( return [cmd] case CommandKind.S: raise ValueError("Unexpected signal!") - case _: + case _: # pragma: no cover typing_extensions.assert_never(cmd.kind) @typing_extensions.override def confuse_result( self, cmd: BaseM, result: Outcome, rng: Generator | None = None, *, stacklevel: int = 1 ) -> Outcome: - r"""Return the measurement result unchanged. - - Amplitude damping is a purely quantum channel and introduces no - classical readout error. Compose this model with another noise - model (e.g. via :class:`ComposeNoiseModel`) to add readout error. - """ + """Assign wrong measurement result with probability ``measure_error_prob``.""" + rng = ensure_rng(rng, stacklevel=stacklevel + 1) + if rng.uniform() < self.measure_error_prob: + return toggle_outcome(result) return result diff --git a/tests/test_noise_model.py b/tests/test_noise_model.py index aa2573136..bd7f2ba3c 100644 --- a/tests/test_noise_model.py +++ b/tests/test_noise_model.py @@ -6,8 +6,8 @@ import pytest from graphix import Pattern -from graphix.command import C, CommandKind, E, M, N, S, T, X from graphix.clifford import Clifford +from graphix.command import C, CommandKind, E, M, N, S, T, X from graphix.noise_models import ( AmplitudeDampingNoise, AmplitudeDampingNoiseModel, @@ -18,7 +18,7 @@ TwoQubitAmplitudeDampingNoise, TwoQubitDepolarisingNoise, ) -from graphix.noise_models.noise_model import NoiselessNoiseModel +from graphix.noise_models.noise_model import NoiselessNoiseModel, NoiseModel from graphix.random_objects import rand_circuit from graphix.sim.density_matrix import DensityMatrix from graphix.simulator import DefaultMeasureMethod @@ -26,7 +26,6 @@ if TYPE_CHECKING: from numpy.random import Generator - from graphix.measurements import Outcome from graphix.noise_models import CommandOrNoise @@ -105,22 +104,6 @@ def test_compose_noise_model_simulation(fx_rng: Generator) -> None: assert np.abs(np.dot(state_mbqc.flatten().conjugate(), DensityMatrix(state).rho.flatten())) == pytest.approx(1) -def test_confuse_result(fx_rng: Generator) -> None: - # Pattern that measures 0 on qubit 0 with probability 1. - pattern = Pattern(cmds=[N(0), M(0)]) - measure_method = DefaultMeasureMethod() - pattern.simulate_pattern( - backend="densitymatrix", noise_model=NoiselessNoiseModel(), rng=fx_rng, measure_method=measure_method - ) - assert measure_method.results[0] == 0 - noise_model = DepolarisingNoiseModel(measure_error_prob=1) - measure_method = DefaultMeasureMethod() - pattern.simulate_pattern( - backend="densitymatrix", noise_model=noise_model, rng=fx_rng, measure_method=measure_method - ) - assert measure_method.results[0] == 1 - - def test_amplitude_damping_command_injection() -> None: """Amplitude damping noise is injected at the correct command positions.""" model = AmplitudeDampingNoiseModel( @@ -158,11 +141,23 @@ def test_amplitude_damping_command_injection() -> None: assert out[1].domain == {1, 2} -@pytest.mark.parametrize("outcome", [0, 1]) -def test_amplitude_damping_confuse_result_is_identity(outcome: Outcome) -> None: - """Amplitude damping introduces no classical readout error.""" - model = AmplitudeDampingNoiseModel() - assert model.confuse_result(M(node=0), outcome) == outcome +@pytest.mark.parametrize( + "noise_model", + [DepolarisingNoiseModel(measure_error_prob=1), AmplitudeDampingNoiseModel(measure_error_prob=1)], +) +def test_confuse_result(fx_rng: Generator, noise_model: NoiseModel) -> None: + # Pattern that measures 0 on qubit 0 with probability 1. + pattern = Pattern(cmds=[N(0), M(0)]) + measure_method = DefaultMeasureMethod() + pattern.simulate_pattern( + backend="densitymatrix", noise_model=NoiselessNoiseModel(), rng=fx_rng, measure_method=measure_method + ) + assert measure_method.results[0] == 0 + measure_method = DefaultMeasureMethod() + pattern.simulate_pattern( + backend="densitymatrix", noise_model=noise_model, rng=fx_rng, measure_method=measure_method + ) + assert measure_method.results[0] == 1 def test_compose_amplitude_damping_depolarising_transpile(fx_rng: Generator) -> None: @@ -231,6 +226,7 @@ def test_compose_amplitude_damping_depolarising_simulation(fx_rng: Generator) -> state_mbqc = pattern.simulate_pattern(backend="densitymatrix", noise_model=noise_model, rng=fx_rng) assert np.abs(np.dot(state_mbqc.flatten().conjugate(), DensityMatrix(state).rho.flatten())) == pytest.approx(1) + def test_amplitude_damping_nqubits() -> None: # covers line 46: AmplitudeDampingNoise.nqubits returning 1 assert AmplitudeDampingNoise(0.1).nqubits == 1 diff --git a/tests/test_noisy_density_matrix.py b/tests/test_noisy_density_matrix.py index e2c4dd799..d9de81f2f 100644 --- a/tests/test_noisy_density_matrix.py +++ b/tests/test_noisy_density_matrix.py @@ -9,10 +9,11 @@ from graphix.branch_selector import ConstBranchSelector, FixedBranchSelector from graphix.command import CommandKind from graphix.fundamentals import angle_to_rad -from graphix.noise_models import AmplitudeDampingNoiseModel, ComposeNoiseModel, DepolarisingNoiseModel +from graphix.noise_models import AmplitudeDampingNoiseModel, DepolarisingNoiseModel from graphix.noise_models.noise_model import NoiselessNoiseModel from graphix.ops import Ops from graphix.sim.density_matrix import DensityMatrix +from graphix.states import BasicStates from graphix.transpiler import Circuit if TYPE_CHECKING: @@ -526,119 +527,97 @@ def test_noisy_measure_confuse_rz_arbitrary( ) -class TestNoisyDensityMatrixAmplitudeDamping: - """Amplitude damping noise model on the density matrix backend.""" - - @pytest.mark.filterwarnings("ignore:Simulating using densitymatrix backend with no noise.") - def test_noiseless_amplitude_damping_hadamard(self, fx_rng: Generator) -> None: - """A zero-parameter amplitude damping model reproduces the noiseless result.""" - hadamardpattern = hpat() - noiselessres = hadamardpattern.simulate_pattern(backend="densitymatrix", rng=fx_rng) - noisynoiselessres = hadamardpattern.simulate_pattern( - backend="densitymatrix", - noise_model=AmplitudeDampingNoiseModel(), - rng=fx_rng, - ) - assert isinstance(noiselessres, DensityMatrix) - assert isinstance(noisynoiselessres, DensityMatrix) - assert np.allclose(noiselessres.rho, np.array([[1.0, 0.0], [0.0, 0.0]])) - assert np.allclose(noisynoiselessres.rho, np.array([[1.0, 0.0], [0.0, 0.0]])) - - @pytest.mark.filterwarnings("ignore:Simulating using densitymatrix backend with no noise.") - def test_noiseless_amplitude_damping_rz(self, fx_rng: Generator) -> None: - - alpha = fx_rng.random() - rzpattern = rzpat(alpha) - noiselessres = rzpattern.simulate_pattern(backend="densitymatrix", rng=fx_rng) - noisynoiselessres = rzpattern.simulate_pattern( - backend="densitymatrix", - noise_model=AmplitudeDampingNoiseModel(), - rng=fx_rng, - ) - assert isinstance(noiselessres, DensityMatrix) - assert isinstance(noisynoiselessres, DensityMatrix) - assert np.allclose(noiselessres.rho, rz_exact_res(alpha)) - assert np.allclose(noisynoiselessres.rho, rz_exact_res(alpha)) - - def test_compose_amplitude_damping_with_depolarising_runs(self, fx_rng: Generator) -> None: - """ComposeNoiseModel of depolarising + amplitude damping runs and stays normalized.""" - hadamardpattern = hpat() - depol = DepolarisingNoiseModel(prepare_error_prob=0.2) - ad = AmplitudeDampingNoiseModel(prepare_error_prob=0.2) - composed = ComposeNoiseModel([depol, ad]) - - r_depol = hadamardpattern.simulate_pattern(backend="densitymatrix", noise_model=depol, rng=fx_rng) - r_both = hadamardpattern.simulate_pattern(backend="densitymatrix", noise_model=composed, rng=fx_rng) - assert isinstance(r_depol, DensityMatrix) - assert isinstance(r_both, DensityMatrix) - assert np.isclose(r_both.rho.trace(), 1.0) - assert not np.allclose(r_depol.rho, r_both.rho) - - class TestAmplitudeDampingAnalytic: - """Analytic per-step verification of amplitude damping on the Hadamard pattern.""" + """Analytic per-step verification of amplitude damping on the Hadamard and RZ patterns.""" @staticmethod - def _kraus(gamma: float) -> list[npt.NDArray[np.complex128]]: + def _kraus_operators(gamma: float) -> list[npt.NDArray[np.complex128]]: return [ np.array([[1.0, 0.0], [0.0, np.sqrt(1 - gamma)]], dtype=np.complex128), np.array([[0.0, np.sqrt(gamma)], [0.0, 0.0]], dtype=np.complex128), ] - def _ad_single(self, rho: npt.NDArray[np.complex128], gamma: float) -> npt.NDArray[np.complex128]: - return sum((k @ rho @ k.conj().T for k in self._kraus(gamma)), np.zeros((2, 2), dtype=np.complex128)) + @staticmethod + def _kron_all(ops: list[npt.NDArray[np.complex128]]) -> npt.NDArray[np.complex128]: + """Kronecker product of a list of operators, left to right.""" + result = ops[0] + for op in ops[1:]: + result = np.kron(result, op).astype(np.complex128) + return result + + def _damp_one_qubit(self, rho: npt.NDArray[np.complex128], gamma: float) -> npt.NDArray[np.complex128]: + """Apply single-qubit amplitude damping to a 2x2 density matrix.""" + return sum( + (k @ rho @ k.conj().T for k in self._kraus_operators(gamma)), + np.zeros((2, 2), dtype=np.complex128), + ) - def _ad_on(self, rho: npt.NDArray[np.complex128], qubit: int, gamma: float) -> npt.NDArray[np.complex128]: + def _damp_qubit_in_register( + self, rho: npt.NDArray[np.complex128], qubit: int, nqubit: int, gamma: float + ) -> npt.NDArray[np.complex128]: + """Apply amplitude damping to a single ``qubit`` of an ``nqubit`` register.""" eye = np.eye(2, dtype=np.complex128) - out = np.zeros((4, 4), dtype=np.complex128) - for k in self._kraus(gamma): - full = np.kron(k, eye) if qubit == 0 else np.kron(eye, k) + dim = 2**nqubit + out = np.zeros((dim, dim), dtype=np.complex128) + for k in self._kraus_operators(gamma): + ops = [eye] * nqubit + ops[qubit] = k + full = self._kron_all(ops) out += full @ rho @ full.conj().T return out - def _ad_two(self, rho: npt.NDArray[np.complex128], gamma: float) -> npt.NDArray[np.complex128]: - out = np.zeros((4, 4), dtype=np.complex128) - for left in self._kraus(gamma): - for right in self._kraus(gamma): - full = np.kron(left, right) + def _damp_two_qubits_in_register( + self, rho: npt.NDArray[np.complex128], qubit_a: int, qubit_b: int, nqubit: int, gamma: float + ) -> npt.NDArray[np.complex128]: + """Apply the two-qubit amplitude damping channel to ``qubit_a`` and ``qubit_b``.""" + eye = np.eye(2, dtype=np.complex128) + dim = 2**nqubit + out = np.zeros((dim, dim), dtype=np.complex128) + for left in self._kraus_operators(gamma): + for right in self._kraus_operators(gamma): + ops = [eye] * nqubit + ops[qubit_a] = left + ops[qubit_b] = right + full = self._kron_all(ops) out += full @ rho @ full.conj().T return out @staticmethod - def _proj_x(outcome: Outcome) -> npt.NDArray[np.complex128]: - vec = np.array([1.0, 1.0 if outcome == 0 else -1.0], dtype=np.complex128) / np.sqrt(2) + def _xy_projector(outcome: Outcome, angle_rad: float) -> npt.NDArray[np.complex128]: + """Projector onto the XY-plane measurement basis state for the given outcome.""" + sign = 1.0 if outcome == 0 else -1.0 + vec = np.array([1.0, sign * np.exp(1j * angle_rad)], dtype=np.complex128) / np.sqrt(2) return np.outer(vec, vec.conj()) - def _expected(self, step: str, gamma: float, outcome: Outcome) -> npt.NDArray[np.complex128]: + def _hadamard_expected(self, step: str, gamma: float, outcome: Outcome) -> npt.NDArray[np.complex128]: + """Output density matrix of the Hadamard pattern with damping at ``step``.""" eye = np.eye(2, dtype=np.complex128) - pauli_x = np.array([[0.0, 1.0], [1.0, 0.0]], dtype=np.complex128) - cz = np.diag([1.0, 1.0, 1.0, -1.0]).astype(np.complex128) - plus = np.array([1.0, 1.0], dtype=np.complex128) / np.sqrt(2) - - rho: npt.NDArray[np.complex128] = np.kron(np.outer(plus, plus.conj()), np.outer(plus, plus.conj())).astype( - np.complex128 - ) + pauli_x = np.asarray(Ops.X, dtype=np.complex128) + cz = np.asarray(Ops.CZ, dtype=np.complex128) + rho = np.asarray(DensityMatrix(data=[BasicStates.PLUS, BasicStates.PLUS]).rho, dtype=np.complex128) if step == "prep": - rho = self._ad_on(rho, 0, gamma) - rho = self._ad_on(rho, 1, gamma) + rho = self._damp_qubit_in_register(rho, 0, 2, gamma) + rho = self._damp_qubit_in_register(rho, 1, 2, gamma) rho = cz @ rho @ cz.conj().T if step == "entangle": - rho = self._ad_two(rho, gamma) + rho = self._damp_two_qubits_in_register(rho, 0, 1, 2, gamma) if step == "measure": - rho = self._ad_on(rho, 0, gamma) - proj = np.kron(self._proj_x(outcome), eye) + rho = self._damp_qubit_in_register(rho, 0, 2, gamma) + + # measured qubit 0 in the X basis (XY plane, angle 0) + proj = np.kron(self._xy_projector(outcome, 0.0), eye) rho = proj @ rho @ proj.conj().T rho = rho / np.trace(rho) + reduced: npt.NDArray[np.complex128] = np.einsum("ijik->jk", rho.reshape(2, 2, 2, 2)) if outcome == 1: reduced = pauli_x @ reduced @ pauli_x.conj().T if step == "xcorr" and outcome == 1: - reduced = self._ad_single(reduced, gamma) + reduced = self._damp_one_qubit(reduced, gamma) return reduced @pytest.mark.parametrize("outcome", [0, 1]) - @pytest.mark.parametrize("gamma", [0.0, 0.3, 0.7, 1.0]) @pytest.mark.parametrize( ("step", "param"), [ @@ -648,9 +627,10 @@ def _expected(self, step: str, gamma: float, outcome: Outcome) -> npt.NDArray[np ("xcorr", "x_error_prob"), ], ) - def test_amplitude_damping_step_matches_analytic( - self, step: str, param: str, gamma: float, outcome: Outcome, fx_rng: Generator + def test_amplitude_damping_step_matches_analytic_hadamard( + self, step: str, param: str, outcome: Outcome, fx_rng: Generator ) -> None: + gamma = fx_rng.random() res = hpat().simulate_pattern( backend="densitymatrix", noise_model=AmplitudeDampingNoiseModel(**{param: gamma}), @@ -658,13 +638,84 @@ def test_amplitude_damping_step_matches_analytic( rng=fx_rng, ) assert isinstance(res, DensityMatrix) - assert np.allclose(res.rho, self._expected(step, gamma, outcome)) + assert np.allclose(res.rho, self._hadamard_expected(step, gamma, outcome)) + + def _rz_expected( + self, step: str, gamma: float, alpha: Angle, outcome_z: Outcome, outcome_x: Outcome + ) -> npt.NDArray[np.complex128]: + """Output density matrix of the RZ pattern with damping at ``step``. + + The RZ pattern measures node 0 (determining the Z correction) and node 1 + (determining the X correction) before the output lands on node 2. + """ + eye = np.eye(2, dtype=np.complex128) + pauli_x = np.asarray(Ops.X, dtype=np.complex128) + pauli_z = np.asarray(Ops.Z, dtype=np.complex128) + cz = np.asarray(Ops.CZ, dtype=np.complex128) + rad = angle_to_rad(alpha) + + rho = np.asarray( + DensityMatrix(data=[BasicStates.PLUS, BasicStates.PLUS, BasicStates.PLUS]).rho, + dtype=np.complex128, + ) + if step == "prep": + for qubit in (0, 1, 2): + rho = self._damp_qubit_in_register(rho, qubit, 3, gamma) + + cz01 = np.kron(cz, eye) + rho = cz01 @ rho @ cz01.conj().T + if step == "entangle": + rho = self._damp_two_qubits_in_register(rho, 0, 1, 3, gamma) + cz12 = np.kron(eye, cz) + rho = cz12 @ rho @ cz12.conj().T + if step == "entangle": + rho = self._damp_two_qubits_in_register(rho, 1, 2, 3, gamma) - # For the measurement step the result has a clean branch-independent - # closed form: damping the qubit before the X-basis measurement leaves - # the state diagonal with the coherence factor ``sqrt(1 - gamma)`` - # showing through. Assert it explicitly as a human-readable cross-check. if step == "measure": - sqrt_term = np.sqrt(1 - gamma) - closed_form = np.array([[(1 + sqrt_term) / 2, 0.0], [0.0, (1 - sqrt_term) / 2]], dtype=np.complex128) - assert np.allclose(res.rho, closed_form) + rho = self._damp_qubit_in_register(rho, 0, 3, gamma) + proj0 = self._kron_all([self._xy_projector(outcome_z, -rad), eye, eye]) + rho = proj0 @ rho @ proj0.conj().T + if step == "measure": + rho = self._damp_qubit_in_register(rho, 1, 3, gamma) + proj1 = self._kron_all([eye, self._xy_projector(outcome_x, 0.0), eye]) + rho = proj1 @ rho @ proj1.conj().T + rho = rho / np.trace(rho) + + reduced: npt.NDArray[np.complex128] = np.einsum("ijkijl->kl", rho.reshape(2, 2, 2, 2, 2, 2)) + if outcome_x == 1: + reduced = pauli_x @ reduced @ pauli_x.conj().T + if step == "xcorr" and outcome_x == 1: + reduced = self._damp_one_qubit(reduced, gamma) + if outcome_z == 1: + reduced = pauli_z @ reduced @ pauli_z.conj().T + if step == "zcorr" and outcome_z == 1: + reduced = self._damp_one_qubit(reduced, gamma) + return reduced + + @pytest.mark.parametrize("outcome_x", [0, 1]) + @pytest.mark.parametrize("outcome_z", [0, 1]) + @pytest.mark.parametrize( + ("step", "param"), + [ + ("prep", "prepare_error_prob"), + ("entangle", "entanglement_error_prob"), + ("measure", "measure_channel_prob"), + ("xcorr", "x_error_prob"), + ("zcorr", "z_error_prob"), + ], + ) + def test_amplitude_damping_step_matches_analytic_rz( + self, step: str, param: str, outcome_z: Outcome, outcome_x: Outcome, fx_rng: Generator + ) -> None: + gamma = fx_rng.random() + alpha = fx_rng.random() + rzpattern = rzpat(alpha) + results: dict[int, Outcome] = {0: outcome_z, 1: outcome_x} + res = rzpattern.simulate_pattern( + backend="densitymatrix", + noise_model=AmplitudeDampingNoiseModel(**{param: gamma}), + branch_selector=FixedBranchSelector(results), + rng=fx_rng, + ) + assert isinstance(res, DensityMatrix) + assert np.allclose(res.rho, self._rz_expected(step, gamma, alpha, outcome_z, outcome_x)) From c3f8c7c37184d9954ade255b447a40dec9938198 Mon Sep 17 00:00:00 2001 From: Nandan Patel Date: Wed, 17 Jun 2026 22:19:21 +0800 Subject: [PATCH 07/11] Used the Ops functionality, ensured all function calls are standardized --- graphix/channels.py | 9 +-- tests/test_noisy_density_matrix.py | 99 +++++++++++++++++++++--------- 2 files changed, 71 insertions(+), 37 deletions(-) diff --git a/graphix/channels.py b/graphix/channels.py index 8b3204044..5cefb914e 100644 --- a/graphix/channels.py +++ b/graphix/channels.py @@ -342,11 +342,4 @@ def two_qubit_amplitude_damping_channel(prob: float) -> KrausChannel: """ k1 = np.array([[1.0, 0.0], [0.0, np.sqrt(1 - prob)]], dtype=np.complex128) k2 = np.array([[0.0, np.sqrt(prob)], [0.0, 0.0]], dtype=np.complex128) - return KrausChannel( - [ - KrausData(1.0, np.kron(k1, k1)), - KrausData(1.0, np.kron(k1, k2)), - KrausData(1.0, np.kron(k2, k1)), - KrausData(1.0, np.kron(k2, k2)), - ] - ) + return KrausChannel([KrausData(1.0, np.kron(left, right)) for left in (k1, k2) for right in (k1, k2)]) diff --git a/tests/test_noisy_density_matrix.py b/tests/test_noisy_density_matrix.py index d9de81f2f..ff7512f91 100644 --- a/tests/test_noisy_density_matrix.py +++ b/tests/test_noisy_density_matrix.py @@ -1,5 +1,6 @@ from __future__ import annotations +import functools from typing import TYPE_CHECKING import numpy as np @@ -10,7 +11,7 @@ from graphix.command import CommandKind from graphix.fundamentals import angle_to_rad from graphix.noise_models import AmplitudeDampingNoiseModel, DepolarisingNoiseModel -from graphix.noise_models.noise_model import NoiselessNoiseModel +from graphix.noise_models.noise_model import ComposeNoiseModel, NoiselessNoiseModel from graphix.ops import Ops from graphix.sim.density_matrix import DensityMatrix from graphix.states import BasicStates @@ -527,6 +528,55 @@ def test_noisy_measure_confuse_rz_arbitrary( ) +class TestNoisyDensityMatrixAmplitudeDamping: + """Amplitude damping noise model on the density matrix backend.""" + + @pytest.mark.filterwarnings("ignore:Simulating using densitymatrix backend with no noise.") + def test_noiseless_amplitude_damping_hadamard(self, fx_rng: Generator) -> None: + """A zero-parameter amplitude damping model reproduces the noiseless result.""" + hadamardpattern = hpat() + noiselessres = hadamardpattern.simulate_pattern(backend="densitymatrix", rng=fx_rng) + noisynoiselessres = hadamardpattern.simulate_pattern( + backend="densitymatrix", + noise_model=AmplitudeDampingNoiseModel(), + rng=fx_rng, + ) + assert isinstance(noiselessres, DensityMatrix) + assert isinstance(noisynoiselessres, DensityMatrix) + assert np.allclose(noiselessres.rho, np.array([[1.0, 0.0], [0.0, 0.0]])) + assert np.allclose(noisynoiselessres.rho, np.array([[1.0, 0.0], [0.0, 0.0]])) + + @pytest.mark.filterwarnings("ignore:Simulating using densitymatrix backend with no noise.") + def test_noiseless_amplitude_damping_rz(self, fx_rng: Generator) -> None: + """A zero-parameter amplitude damping model reproduces the noiseless RZ result.""" + alpha = fx_rng.random() + rzpattern = rzpat(alpha) + noiselessres = rzpattern.simulate_pattern(backend="densitymatrix", rng=fx_rng) + noisynoiselessres = rzpattern.simulate_pattern( + backend="densitymatrix", + noise_model=AmplitudeDampingNoiseModel(), + rng=fx_rng, + ) + assert isinstance(noiselessres, DensityMatrix) + assert isinstance(noisynoiselessres, DensityMatrix) + assert np.allclose(noiselessres.rho, rz_exact_res(alpha)) + assert np.allclose(noisynoiselessres.rho, rz_exact_res(alpha)) + + def test_compose_amplitude_damping_with_depolarising_runs(self, fx_rng: Generator) -> None: + """ComposeNoiseModel of depolarising + amplitude damping runs and stays normalized.""" + hadamardpattern = hpat() + depol = DepolarisingNoiseModel(prepare_error_prob=0.2) + ad = AmplitudeDampingNoiseModel(prepare_error_prob=0.2) + composed = ComposeNoiseModel([depol, ad]) + + r_depol = hadamardpattern.simulate_pattern(backend="densitymatrix", noise_model=depol, rng=fx_rng) + r_both = hadamardpattern.simulate_pattern(backend="densitymatrix", noise_model=composed, rng=fx_rng) + assert isinstance(r_depol, DensityMatrix) + assert isinstance(r_both, DensityMatrix) + assert np.isclose(r_both.rho.trace(), 1.0) + assert not np.allclose(r_depol.rho, r_both.rho) + + class TestAmplitudeDampingAnalytic: """Analytic per-step verification of amplitude damping on the Hadamard and RZ patterns.""" @@ -540,10 +590,7 @@ def _kraus_operators(gamma: float) -> list[npt.NDArray[np.complex128]]: @staticmethod def _kron_all(ops: list[npt.NDArray[np.complex128]]) -> npt.NDArray[np.complex128]: """Kronecker product of a list of operators, left to right.""" - result = ops[0] - for op in ops[1:]: - result = np.kron(result, op).astype(np.complex128) - return result + return functools.reduce(lambda a, b: np.kron(a, b).astype(np.complex128), ops) def _damp_one_qubit(self, rho: npt.NDArray[np.complex128], gamma: float) -> npt.NDArray[np.complex128]: """Apply single-qubit amplitude damping to a 2x2 density matrix.""" @@ -556,7 +603,7 @@ def _damp_qubit_in_register( self, rho: npt.NDArray[np.complex128], qubit: int, nqubit: int, gamma: float ) -> npt.NDArray[np.complex128]: """Apply amplitude damping to a single ``qubit`` of an ``nqubit`` register.""" - eye = np.eye(2, dtype=np.complex128) + eye = np.asarray(Ops.I, dtype=np.complex128) dim = 2**nqubit out = np.zeros((dim, dim), dtype=np.complex128) for k in self._kraus_operators(gamma): @@ -570,7 +617,7 @@ def _damp_two_qubits_in_register( self, rho: npt.NDArray[np.complex128], qubit_a: int, qubit_b: int, nqubit: int, gamma: float ) -> npt.NDArray[np.complex128]: """Apply the two-qubit amplitude damping channel to ``qubit_a`` and ``qubit_b``.""" - eye = np.eye(2, dtype=np.complex128) + eye = np.asarray(Ops.I, dtype=np.complex128) dim = 2**nqubit out = np.zeros((dim, dim), dtype=np.complex128) for left in self._kraus_operators(gamma): @@ -591,15 +638,12 @@ def _xy_projector(outcome: Outcome, angle_rad: float) -> npt.NDArray[np.complex1 def _hadamard_expected(self, step: str, gamma: float, outcome: Outcome) -> npt.NDArray[np.complex128]: """Output density matrix of the Hadamard pattern with damping at ``step``.""" - eye = np.eye(2, dtype=np.complex128) - pauli_x = np.asarray(Ops.X, dtype=np.complex128) - cz = np.asarray(Ops.CZ, dtype=np.complex128) - + eye = np.asarray(Ops.I, dtype=np.complex128) rho = np.asarray(DensityMatrix(data=[BasicStates.PLUS, BasicStates.PLUS]).rho, dtype=np.complex128) if step == "prep": rho = self._damp_qubit_in_register(rho, 0, 2, gamma) rho = self._damp_qubit_in_register(rho, 1, 2, gamma) - rho = cz @ rho @ cz.conj().T + rho = Ops.CZ @ rho @ Ops.CZ if step == "entangle": rho = self._damp_two_qubits_in_register(rho, 0, 1, 2, gamma) if step == "measure": @@ -612,9 +656,9 @@ def _hadamard_expected(self, step: str, gamma: float, outcome: Outcome) -> npt.N reduced: npt.NDArray[np.complex128] = np.einsum("ijik->jk", rho.reshape(2, 2, 2, 2)) if outcome == 1: - reduced = pauli_x @ reduced @ pauli_x.conj().T - if step == "xcorr" and outcome == 1: - reduced = self._damp_one_qubit(reduced, gamma) + reduced = Ops.X @ reduced @ Ops.X + if step == "xcorr": + reduced = self._damp_one_qubit(reduced, gamma) return reduced @pytest.mark.parametrize("outcome", [0, 1]) @@ -648,10 +692,7 @@ def _rz_expected( The RZ pattern measures node 0 (determining the Z correction) and node 1 (determining the X correction) before the output lands on node 2. """ - eye = np.eye(2, dtype=np.complex128) - pauli_x = np.asarray(Ops.X, dtype=np.complex128) - pauli_z = np.asarray(Ops.Z, dtype=np.complex128) - cz = np.asarray(Ops.CZ, dtype=np.complex128) + eye = np.asarray(Ops.I, dtype=np.complex128) rad = angle_to_rad(alpha) rho = np.asarray( @@ -662,12 +703,12 @@ def _rz_expected( for qubit in (0, 1, 2): rho = self._damp_qubit_in_register(rho, qubit, 3, gamma) - cz01 = np.kron(cz, eye) - rho = cz01 @ rho @ cz01.conj().T + cz01 = np.kron(Ops.CZ, Ops.I) + rho = cz01 @ rho @ cz01 if step == "entangle": rho = self._damp_two_qubits_in_register(rho, 0, 1, 3, gamma) - cz12 = np.kron(eye, cz) - rho = cz12 @ rho @ cz12.conj().T + cz12 = np.kron(Ops.I, Ops.CZ) + rho = cz12 @ rho @ cz12 if step == "entangle": rho = self._damp_two_qubits_in_register(rho, 1, 2, 3, gamma) @@ -683,13 +724,13 @@ def _rz_expected( reduced: npt.NDArray[np.complex128] = np.einsum("ijkijl->kl", rho.reshape(2, 2, 2, 2, 2, 2)) if outcome_x == 1: - reduced = pauli_x @ reduced @ pauli_x.conj().T - if step == "xcorr" and outcome_x == 1: - reduced = self._damp_one_qubit(reduced, gamma) + reduced = Ops.X @ reduced @ Ops.X + if step == "xcorr": + reduced = self._damp_one_qubit(reduced, gamma) if outcome_z == 1: - reduced = pauli_z @ reduced @ pauli_z.conj().T - if step == "zcorr" and outcome_z == 1: - reduced = self._damp_one_qubit(reduced, gamma) + reduced = Ops.Z @ reduced @ Ops.Z + if step == "zcorr": + reduced = self._damp_one_qubit(reduced, gamma) return reduced @pytest.mark.parametrize("outcome_x", [0, 1]) From ebad847e60c40ec4e073786d6aad0a00a8de2b6f Mon Sep 17 00:00:00 2001 From: Nandan Patel Date: Thu, 18 Jun 2026 00:01:07 +0800 Subject: [PATCH 08/11] Made reviewer changes --- tests/test_noise_model.py | 37 ++++++++++++++++++++---------- tests/test_noisy_density_matrix.py | 16 +++++-------- 2 files changed, 31 insertions(+), 22 deletions(-) diff --git a/tests/test_noise_model.py b/tests/test_noise_model.py index bd7f2ba3c..c3f3f9b58 100644 --- a/tests/test_noise_model.py +++ b/tests/test_noise_model.py @@ -7,7 +7,7 @@ from graphix import Pattern from graphix.clifford import Clifford -from graphix.command import C, CommandKind, E, M, N, S, T, X +from graphix.command import C, CommandKind, E, M, N, S, T, X, Z from graphix.noise_models import ( AmplitudeDampingNoise, AmplitudeDampingNoiseModel, @@ -104,41 +104,54 @@ def test_compose_noise_model_simulation(fx_rng: Generator) -> None: assert np.abs(np.dot(state_mbqc.flatten().conjugate(), DensityMatrix(state).rho.flatten())) == pytest.approx(1) +def _assert_noise( + cmd: CommandOrNoise, + noise_type: type[AmplitudeDampingNoise | TwoQubitAmplitudeDampingNoise], + gamma: float, + nodes: list[int], + domain: set[int] | None = None, +) -> None: + assert isinstance(cmd, ApplyNoise) + assert isinstance(cmd.noise, noise_type) + assert cmd.noise.prob == gamma + assert cmd.nodes == nodes + assert cmd.domain == domain + + def test_amplitude_damping_command_injection() -> None: """Amplitude damping noise is injected at the correct command positions.""" model = AmplitudeDampingNoiseModel( prepare_error_prob=0.1, x_error_prob=0.2, + z_error_prob=0.5, entanglement_error_prob=0.3, measure_channel_prob=0.4, ) # N: noise applied AFTER preparation out = model.command(N(node=0)) - assert len(out) == 2 assert out[0].kind == CommandKind.N - assert isinstance(out[1], ApplyNoise) - assert isinstance(out[1].noise, AmplitudeDampingNoise) - assert out[1].nodes == [0] + _assert_noise(out[1], AmplitudeDampingNoise, 0.1, [0]) # E: two-qubit noise applied AFTER entanglement out = model.command(E(nodes=(0, 1))) assert out[0].kind == CommandKind.E - assert isinstance(out[1], ApplyNoise) - assert isinstance(out[1].noise, TwoQubitAmplitudeDampingNoise) - assert out[1].noise.nqubits == 2 + _assert_noise(out[1], TwoQubitAmplitudeDampingNoise, 0.3, [0, 1]) # M: noise applied BEFORE measurement out = model.command(M(node=0)) - assert isinstance(out[0], ApplyNoise) - assert isinstance(out[0].noise, AmplitudeDampingNoise) + _assert_noise(out[0], AmplitudeDampingNoise, 0.4, [0]) assert out[1].kind == CommandKind.M # X: correction kept, noise conditioned on the same domain out = model.command(X(node=0, domain={1, 2})) assert out[0].kind == CommandKind.X - assert isinstance(out[1], ApplyNoise) - assert out[1].domain == {1, 2} + _assert_noise(out[1], AmplitudeDampingNoise, 0.2, [0], {1, 2}) + + # Z: correction kept, noise conditioned on the same domain + out = model.command(Z(node=0, domain={3})) + assert out[0].kind == CommandKind.Z + _assert_noise(out[1], AmplitudeDampingNoise, 0.5, [0], {3}) @pytest.mark.parametrize( diff --git a/tests/test_noisy_density_matrix.py b/tests/test_noisy_density_matrix.py index ff7512f91..6f38919c6 100644 --- a/tests/test_noisy_density_matrix.py +++ b/tests/test_noisy_density_matrix.py @@ -603,11 +603,10 @@ def _damp_qubit_in_register( self, rho: npt.NDArray[np.complex128], qubit: int, nqubit: int, gamma: float ) -> npt.NDArray[np.complex128]: """Apply amplitude damping to a single ``qubit`` of an ``nqubit`` register.""" - eye = np.asarray(Ops.I, dtype=np.complex128) dim = 2**nqubit out = np.zeros((dim, dim), dtype=np.complex128) for k in self._kraus_operators(gamma): - ops = [eye] * nqubit + ops = [Ops.I] * nqubit ops[qubit] = k full = self._kron_all(ops) out += full @ rho @ full.conj().T @@ -617,12 +616,11 @@ def _damp_two_qubits_in_register( self, rho: npt.NDArray[np.complex128], qubit_a: int, qubit_b: int, nqubit: int, gamma: float ) -> npt.NDArray[np.complex128]: """Apply the two-qubit amplitude damping channel to ``qubit_a`` and ``qubit_b``.""" - eye = np.asarray(Ops.I, dtype=np.complex128) dim = 2**nqubit out = np.zeros((dim, dim), dtype=np.complex128) for left in self._kraus_operators(gamma): for right in self._kraus_operators(gamma): - ops = [eye] * nqubit + ops = [Ops.I] * nqubit ops[qubit_a] = left ops[qubit_b] = right full = self._kron_all(ops) @@ -638,7 +636,6 @@ def _xy_projector(outcome: Outcome, angle_rad: float) -> npt.NDArray[np.complex1 def _hadamard_expected(self, step: str, gamma: float, outcome: Outcome) -> npt.NDArray[np.complex128]: """Output density matrix of the Hadamard pattern with damping at ``step``.""" - eye = np.asarray(Ops.I, dtype=np.complex128) rho = np.asarray(DensityMatrix(data=[BasicStates.PLUS, BasicStates.PLUS]).rho, dtype=np.complex128) if step == "prep": rho = self._damp_qubit_in_register(rho, 0, 2, gamma) @@ -650,8 +647,8 @@ def _hadamard_expected(self, step: str, gamma: float, outcome: Outcome) -> npt.N rho = self._damp_qubit_in_register(rho, 0, 2, gamma) # measured qubit 0 in the X basis (XY plane, angle 0) - proj = np.kron(self._xy_projector(outcome, 0.0), eye) - rho = proj @ rho @ proj.conj().T + proj = np.kron(self._xy_projector(outcome, 0.0), Ops.I) + rho = proj @ rho @ proj rho = rho / np.trace(rho) reduced: npt.NDArray[np.complex128] = np.einsum("ijik->jk", rho.reshape(2, 2, 2, 2)) @@ -692,7 +689,6 @@ def _rz_expected( The RZ pattern measures node 0 (determining the Z correction) and node 1 (determining the X correction) before the output lands on node 2. """ - eye = np.asarray(Ops.I, dtype=np.complex128) rad = angle_to_rad(alpha) rho = np.asarray( @@ -714,11 +710,11 @@ def _rz_expected( if step == "measure": rho = self._damp_qubit_in_register(rho, 0, 3, gamma) - proj0 = self._kron_all([self._xy_projector(outcome_z, -rad), eye, eye]) + proj0 = self._kron_all([self._xy_projector(outcome_z, -rad), Ops.I, Ops.I]) rho = proj0 @ rho @ proj0.conj().T if step == "measure": rho = self._damp_qubit_in_register(rho, 1, 3, gamma) - proj1 = self._kron_all([eye, self._xy_projector(outcome_x, 0.0), eye]) + proj1 = self._kron_all([Ops.I, self._xy_projector(outcome_x, 0.0), Ops.I]) rho = proj1 @ rho @ proj1.conj().T rho = rho / np.trace(rho) From ff8a83b3d80829add0da4f4a7d5904c7a91c1082 Mon Sep 17 00:00:00 2001 From: Nandan Patel Date: Thu, 18 Jun 2026 10:03:07 +0800 Subject: [PATCH 09/11] Updated CHANGELOG.md --- CHANGELOG.md | 2 ++ 1 file changed, 2 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index ca9380e82..689852461 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -11,6 +11,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - #476: Introduced new methods `OpenGraph.extract_circuit`, `CliffordMap.to_tableau` and new function `graphix.circ_ext.compilation.cm_berg_pass`. Circuit extraction can be done natively in Graphix. +- #545: Added an amplitude damping noise model. Introduces `amplitude_damping_channel` / `two_qubit_amplitude_damping_channel`, the `AmplitudeDampingNoise` / `TwoQubitAmplitudeDampingNoise` noise elements, and `AmplitudeDampingNoiseModel`. + - #490: Introduced new `Instruction` and `Command` namespace classes for instruction and command instantiation. - #505 From 1d5b095d1144f1383c03a4eba0ece0faf2ab6f7e Mon Sep 17 00:00:00 2001 From: Nandan Patel Date: Thu, 18 Jun 2026 17:23:30 +0800 Subject: [PATCH 10/11] Added noise model to init file --- graphix/__init__.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/graphix/__init__.py b/graphix/__init__.py index 12f8da89c..8ca30e8ff 100644 --- a/graphix/__init__.py +++ b/graphix/__init__.py @@ -13,7 +13,7 @@ from graphix.graphsim import GraphState from graphix.instruction import Instruction from graphix.measurements import BlochMeasurement, Measurement, PauliMeasurement -from graphix.noise_models import DepolarisingNoiseModel, NoiseModel +from graphix.noise_models import AmplitudeDampingNoiseModel, DepolarisingNoiseModel, NoiseModel from graphix.opengraph import OpenGraph from graphix.optimization import StandardizedPattern from graphix.parameter import Placeholder @@ -26,6 +26,7 @@ from graphix.transpiler import Circuit __all__ = [ + "AmplitudeDampingNoiseModel", "ANGLE_PI", "Axis", "BasicStates", From 9ff69763494023acee2ce31872be45cc6a7085fe Mon Sep 17 00:00:00 2001 From: Nandan Patel Date: Thu, 18 Jun 2026 20:09:44 +0800 Subject: [PATCH 11/11] Fix sorting in __init__.py --- graphix/__init__.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/graphix/__init__.py b/graphix/__init__.py index 8ca30e8ff..d3d28a818 100644 --- a/graphix/__init__.py +++ b/graphix/__init__.py @@ -26,8 +26,8 @@ from graphix.transpiler import Circuit __all__ = [ - "AmplitudeDampingNoiseModel", "ANGLE_PI", + "AmplitudeDampingNoiseModel", "Axis", "BasicStates", "BlochMeasurement",