diff --git a/graphix/__init__.py b/graphix/__init__.py index 12f8da89c..d3d28a818 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 @@ -27,6 +27,7 @@ __all__ = [ "ANGLE_PI", + "AmplitudeDampingNoiseModel", "Axis", "BasicStates", "BlochMeasurement", diff --git a/graphix/channels.py b/graphix/channels.py index 8bb0ebc1b..32ff44eb5 100644 --- a/graphix/channels.py +++ b/graphix/channels.py @@ -200,6 +200,76 @@ def depolarising_channel(prob: float) -> KrausChannel: ) +def amplitude_damping_channel(gamma: 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 + ---------- + gamma : float + The damping probability, between 0 and 1. + + Returns + ------- + :class:`graphix.channels.KrausChannel` object + containing the corresponding Kraus operators + """ + if not 0 <= gamma <= 1: + raise ValueError("gamma must be between 0 and 1.") + sqrt_gamma = np.sqrt(gamma) + sqrt_one_minus_gamma = np.sqrt(1 - gamma) + return KrausChannel( + [ + KrausData(1.0, np.array([[1, 0], [0, sqrt_one_minus_gamma]], dtype=np.complex128)), + KrausData(1.0, np.array([[0, sqrt_gamma], [0, 0]], dtype=np.complex128)), + ] + ) + + +def two_qubit_amplitude_damping_channel(gamma: float) -> KrausChannel: + r"""Two-qubit amplitude damping channel (independent tensor product). + + The two-qubit channel is formed by the tensor product of two single-qubit + amplitude damping channels, yielding 4 Kraus operators: + + .. math:: + \{K_1 \otimes K_1,\; K_1 \otimes K_2,\; K_2 \otimes K_1,\; K_2 \otimes K_2\} + + Parameters + ---------- + gamma : float + The damping probability, between 0 and 1. + + Returns + ------- + :class:`graphix.channels.KrausChannel` object + containing the corresponding Kraus operators + """ + if not 0 <= gamma <= 1: + raise ValueError("gamma must be between 0 and 1.") + sqrt_gamma = np.sqrt(gamma) + sqrt_one_minus_gamma = np.sqrt(1 - gamma) + k1 = np.array([[1, 0], [0, sqrt_one_minus_gamma]], dtype=np.complex128) + k2 = np.array([[0, sqrt_gamma], [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)), + ] + ) + + def pauli_channel(px: float, py: float, pz: float) -> KrausChannel: r"""Single-qubit Pauli channel. 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..9eb6afaba --- /dev/null +++ b/graphix/noise_models/amplitude_damping.py @@ -0,0 +1,156 @@ +"""Amplitude damping noise model.""" + +from __future__ import annotations + +from typing import TYPE_CHECKING + +import typing_extensions + +from graphix.channels import amplitude_damping_channel, 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: + from collections.abc import Iterable + + from numpy.random import Generator + + from graphix.channels import KrausChannel + from graphix.measurements import Outcome + from graphix.noise_models.noise_model import CommandOrNoise + + +class AmplitudeDampingNoise(Noise): + """One-qubit amplitude damping noise with damping parameter ``gamma``.""" + + gamma = Probability() + + def __init__(self, gamma: float) -> None: + """Initialize one-qubit amplitude damping noise. + + Parameters + ---------- + gamma : float + Damping parameter of the noise, between 0 and 1. + """ + self.gamma = gamma + + @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.gamma) + + +class TwoQubitAmplitudeDampingNoise(Noise): + """Two-qubit amplitude damping noise with damping parameter ``gamma``.""" + + gamma = Probability() + + def __init__(self, gamma: float) -> None: + """Initialize two-qubit amplitude damping noise. + + Parameters + ---------- + gamma : float + Damping parameter of the noise, between 0 and 1. + """ + self.gamma = gamma + + @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.gamma) + + +class AmplitudeDampingNoiseModel(NoiseModel): + """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, + 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_error_prob = measure_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: + """Assign wrong measurement result cmd = "M".""" + 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_kraus.py b/tests/test_kraus.py index f7e4983f0..f4e0fa57d 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,49 @@ 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(0, 0.5) + sqrt_gamma = np.sqrt(gamma) + sqrt_one_minus_gamma = np.sqrt(1 - gamma) + data_k1 = np.array([[1, 0], [0, sqrt_one_minus_gamma]], dtype=np.complex128) + data_k2 = np.array([[0, sqrt_gamma], [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 + assert np.allclose(ad_channel[0].operator, data_k1) + assert np.allclose(ad_channel[0].coef, 1.0) + assert np.allclose(ad_channel[1].operator, data_k2) + assert np.allclose(ad_channel[1].coef, 1.0) + + def test_two_qubit_amplitude_damping_channel(self, fx_rng: Generator) -> None: + gamma = fx_rng.uniform(0, 0.5) + sqrt_gamma = np.sqrt(gamma) + sqrt_one_minus_gamma = np.sqrt(1 - gamma) + k1 = np.array([[1, 0], [0, sqrt_one_minus_gamma]], dtype=np.complex128) + k2 = np.array([[0, sqrt_gamma], [0, 0]], dtype=np.complex128) + ad2_channel = two_qubit_amplitude_damping_channel(gamma) + assert isinstance(ad2_channel, KrausChannel) + assert ad2_channel.nqubit == 2 + assert len(ad2_channel) == 4 + expected = [np.kron(k1, k1), np.kron(k1, k2), np.kron(k2, k1), np.kron(k2, k2)] + for i, exp in enumerate(expected): + assert np.allclose(ad2_channel[i].operator, exp) + assert np.allclose(ad2_channel[i].coef, 1.0) + + def test_amplitude_damping_channel_edge_cases(self) -> None: + # gamma=0 → only K1 (identity) + c0 = amplitude_damping_channel(0.0) + assert np.allclose(c0[0].operator, np.eye(2)) + assert np.allclose(c0[0].coef, 1.0) + assert len(c0) == 2 + # gamma=1 → K1 = |0><0|, K2 = |0><1| + c1 = amplitude_damping_channel(1.0) + assert np.allclose(c1[0].operator, np.array([[1, 0], [0, 0]])) + assert np.allclose(c1[1].operator, np.array([[0, 1], [0, 0]])) + # invalid gamma + with pytest.raises(ValueError): + amplitude_damping_channel(-0.1) + with pytest.raises(ValueError): + amplitude_damping_channel(1.1) diff --git a/tests/test_noise_model.py b/tests/test_noise_model.py index 3babd2a4c..3c831aa71 100644 --- a/tests/test_noise_model.py +++ b/tests/test_noise_model.py @@ -8,10 +8,13 @@ from graphix import Pattern from graphix.command import CommandKind, M, N from graphix.noise_models import ( + AmplitudeDampingNoise, + AmplitudeDampingNoiseModel, ApplyNoise, ComposeNoiseModel, DepolarisingNoise, DepolarisingNoiseModel, + TwoQubitAmplitudeDampingNoise, TwoQubitDepolarisingNoise, ) from graphix.noise_models.noise_model import NoiselessNoiseModel @@ -100,6 +103,83 @@ 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_amplitude_damping_noise_nqubits() -> None: + ad_noise = AmplitudeDampingNoise(0.3) + assert ad_noise.nqubits == 1 + + +def test_two_qubit_amplitude_damping_noise_nqubits() -> None: + ad2_noise = TwoQubitAmplitudeDampingNoise(0.3) + assert ad2_noise.nqubits == 2 + + +def test_amplitude_damping_noise_to_kraus_channel(fx_rng: Generator) -> None: + gamma = fx_rng.uniform(0, 0.5) + ad_noise = AmplitudeDampingNoise(gamma) + channel = ad_noise.to_kraus_channel() + assert channel.nqubit == 1 + assert len(channel) == 2 + + +def test_amplitude_damping_noise_model_default() -> None: + model = AmplitudeDampingNoiseModel() + assert model.prepare_error_prob == 0.0 + assert model.measure_error_prob == 0.0 + assert not model.entanglement_error_prob + assert not model.measure_channel_prob + + +def test_amplitude_damping_noise_model_properties(fx_rng: Generator) -> None: + gamma = fx_rng.uniform(0, 0.5) + model = AmplitudeDampingNoiseModel( + prepare_error_prob=gamma, + x_error_prob=gamma, + z_error_prob=gamma, + entanglement_error_prob=gamma, + measure_channel_prob=gamma, + measure_error_prob=gamma, + ) + assert model.prepare_error_prob == gamma + assert model.x_error_prob == gamma + assert model.z_error_prob == gamma + assert model.entanglement_error_prob == gamma + assert model.measure_channel_prob == gamma + assert model.measure_error_prob == gamma + + +def test_amplitude_damping_noise_model_input_nodes(fx_rng: Generator) -> None: + gamma = fx_rng.uniform(0, 0.5) + model = AmplitudeDampingNoiseModel(prepare_error_prob=gamma) + nodes = model.input_nodes([0, 1]) + assert len(nodes) == 2 + for node_cmd in nodes: + assert isinstance(node_cmd, ApplyNoise) + assert isinstance(node_cmd.noise, AmplitudeDampingNoise) + assert node_cmd.noise.gamma == gamma + + +def test_amplitude_damping_noise_model_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=AmplitudeDampingNoiseModel(), + rng=fx_rng, + measure_method=measure_method, + ) + assert measure_method.results[0] == 0 + noise_model = AmplitudeDampingNoiseModel(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_confuse_result(fx_rng: Generator) -> None: # Pattern that measures 0 on qubit 0 with probability 1. pattern = Pattern(cmds=[N(0), M(0)]) diff --git a/tests/test_noisy_density_matrix.py b/tests/test_noisy_density_matrix.py index 104923265..af41008a9 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, DepolarisingNoiseModel from graphix.noise_models.noise_model import NoiselessNoiseModel from graphix.ops import Ops from graphix.sim.density_matrix import DensityMatrix @@ -524,3 +524,281 @@ 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) ) + + # --- Amplitude damping noise model tests --- + + def test_noiseless_amplitude_damping_hadamard(self, fx_rng: Generator) -> None: + """AmplitudeDampingNoiseModel with all-zero defaults behaves like noiseless.""" + hadamardpattern = hpat() + res = hadamardpattern.simulate_pattern( + backend="densitymatrix", + noise_model=AmplitudeDampingNoiseModel(), + rng=fx_rng, + ) + assert isinstance(res, DensityMatrix) + assert np.allclose(res.rho, np.array([[1.0, 0.0], [0.0, 0.0]])) + + def test_amplitude_damping_noise_changes_hadamard(self, fx_rng: Generator) -> None: + """Amplitude damping prepare noise changes the state from the noiseless result.""" + hadamardpattern = hpat() + gamma = fx_rng.uniform(0.1, 0.5) + res = hadamardpattern.simulate_pattern( + backend="densitymatrix", + noise_model=AmplitudeDampingNoiseModel(prepare_error_prob=gamma), + rng=fx_rng, + ) + assert isinstance(res, DensityMatrix) + # State should be different from noiseless |0> + assert not np.allclose(res.rho, np.array([[1.0, 0.0], [0.0, 0.0]])) + # Density matrix should be valid: trace 1, Hermitian + assert np.isclose(np.trace(res.rho), 1.0) + assert np.allclose(res.rho, res.rho.conj().T) + + def test_amplitude_damping_confuse_measurement(self, fx_rng: Generator) -> None: + """measure_error_prob=1 flips measurement outcome.""" + hadamardpattern = hpat() + res = hadamardpattern.simulate_pattern( + backend="densitymatrix", + noise_model=AmplitudeDampingNoiseModel(measure_error_prob=1.0), + rng=fx_rng, + ) + assert isinstance(res, DensityMatrix) + assert np.allclose(res.rho, np.array([[0.0, 0.0], [0.0, 1.0]])) + + # ----- Analytical comparison tests for amplitude damping ----- + + @staticmethod + def _ad_kraus(gamma: float) -> tuple[npt.NDArray[np.complex128], npt.NDArray[np.complex128]]: + """Return the two Kraus operators for the single-qubit amplitude damping channel.""" + k0 = np.array([[1.0, 0.0], [0.0, np.sqrt(1.0 - gamma)]], dtype=np.complex128) + k1 = np.array([[0.0, np.sqrt(gamma)], [0.0, 0.0]], dtype=np.complex128) + return k0, k1 + + @staticmethod + def _apply_ad_2q(rho: npt.NDArray[np.complex128], gamma: float, qubit: int) -> npt.NDArray[np.complex128]: + """Apply amplitude damping with parameter gamma to one qubit of a 2-qubit density matrix.""" + k0, k1 = TestNoisyDensityMatrixBackend._ad_kraus(gamma) + if qubit == 0: + op0 = np.kron(k0, np.eye(2, dtype=np.complex128)) + op1 = np.kron(k1, np.eye(2, dtype=np.complex128)) + else: + op0 = np.kron(np.eye(2, dtype=np.complex128), k0) + op1 = np.kron(np.eye(2, dtype=np.complex128), k1) + return op0 @ rho @ op0.conj().T + op1 @ rho @ op1.conj().T + + @staticmethod + def _partial_trace_qubit0( + rho: npt.NDArray[np.complex128], + ) -> npt.NDArray[np.complex128]: + """Trace out qubit 0 of a 2-qubit density matrix.""" + res = np.zeros((2, 2), dtype=np.complex128) + for k in range(2): + e = np.zeros(2, dtype=np.complex128) + e[k] = 1.0 + proj = np.kron(e.reshape(-1, 1), np.eye(2, dtype=np.complex128)) + res += proj.conj().T @ rho @ proj + return res + + @staticmethod + def _amplitude_damping_measure_analytical(gamma: float, outcome: Outcome) -> npt.NDArray[np.complex128]: + """Compute expected hpat output for AmplitudeDampingNoiseModel(measure_channel_prob=gamma). + + Applies AD on node 0 after CZ, then measures in X basis, then + applies X(1,{0}) correction when outcome is odd. + """ + rho_plus = 0.5 * np.ones((2, 2), dtype=np.complex128) + rho_plus[0, 1] = 0.5 + rho_plus[1, 0] = 0.5 + rho_initial = np.kron(rho_plus, rho_plus) + + cz = Ops.CZ + rho_cz = cz @ rho_initial @ cz.conj().T + + rho_noisy = TestNoisyDensityMatrixBackend._apply_ad_2q(rho_cz, gamma, qubit=0) + + if outcome == 0: + proj = np.array([[0.5, 0.5], [0.5, 0.5]], dtype=np.complex128) + else: + proj = np.array([[0.5, -0.5], [-0.5, 0.5]], dtype=np.complex128) + proj_full = np.kron(proj, np.eye(2, dtype=np.complex128)) + rho_meas = proj_full @ rho_noisy @ proj_full.conj().T + prob = np.trace(rho_meas).real + rho_out = TestNoisyDensityMatrixBackend._partial_trace_qubit0(rho_meas) / prob + + if outcome % 2 == 1: + x_mat = Ops.X + rho_out = x_mat @ rho_out @ x_mat.conj().T + + return rho_out # type: ignore[no-any-return] + + @staticmethod + def _amplitude_damping_prepare_analytical(gamma: float, outcome: Outcome) -> npt.NDArray[np.complex128]: + """Compute expected hpat output for AmplitudeDampingNoiseModel(prepare_error_prob=gamma). + + Applies AD to |+⟩ on each node independently, then CZ, measure, correct. + """ + k0, k1 = TestNoisyDensityMatrixBackend._ad_kraus(gamma) + rho_plus = 0.5 * np.ones((2, 2), dtype=np.complex128) + rho_plus[0, 1] = 0.5 + rho_plus[1, 0] = 0.5 + + rho_ad = k0 @ rho_plus @ k0.conj().T + k1 @ rho_plus @ k1.conj().T + rho_prep = np.kron(rho_ad, rho_ad) + + cz = Ops.CZ + rho_cz = cz @ rho_prep @ cz.conj().T + + if outcome == 0: + proj = np.array([[0.5, 0.5], [0.5, 0.5]], dtype=np.complex128) + else: + proj = np.array([[0.5, -0.5], [-0.5, 0.5]], dtype=np.complex128) + proj_full = np.kron(proj, np.eye(2, dtype=np.complex128)) + rho_meas = proj_full @ rho_cz @ proj_full.conj().T + prob = np.trace(rho_meas).real + rho_out = TestNoisyDensityMatrixBackend._partial_trace_qubit0(rho_meas) / prob + + if outcome % 2 == 1: + x_mat = Ops.X + rho_out = x_mat @ rho_out @ x_mat.conj().T + + return rho_out # type: ignore[no-any-return] + + @staticmethod + def _amplitude_damping_entanglement_analytical(gamma: float, outcome: Outcome) -> npt.NDArray[np.complex128]: + """Compute expected hpat output for AmplitudeDampingNoiseModel(entanglement_error_prob=gamma). + + Applies independent AD to each qubit after CZ (tensor Kraus product), + then measures and corrects. + """ + rho_plus = 0.5 * np.ones((2, 2), dtype=np.complex128) + rho_plus[0, 1] = 0.5 + rho_plus[1, 0] = 0.5 + rho_initial = np.kron(rho_plus, rho_plus) + + cz = Ops.CZ + rho_cz = cz @ rho_initial @ cz.conj().T + + k0, k1 = TestNoisyDensityMatrixBackend._ad_kraus(gamma) + kraus_2q = [ + np.kron(k0, k0), + np.kron(k0, k1), + np.kron(k1, k0), + np.kron(k1, k1), + ] + rho_noisy = np.zeros_like(rho_cz, dtype=np.complex128) + for op in kraus_2q: + rho_noisy += op @ rho_cz @ op.conj().T + + if outcome == 0: + proj = np.array([[0.5, 0.5], [0.5, 0.5]], dtype=np.complex128) + else: + proj = np.array([[0.5, -0.5], [-0.5, 0.5]], dtype=np.complex128) + proj_full = np.kron(proj, np.eye(2, dtype=np.complex128)) + rho_meas = proj_full @ rho_noisy @ proj_full.conj().T + prob = np.trace(rho_meas).real + rho_out = TestNoisyDensityMatrixBackend._partial_trace_qubit0(rho_meas) / prob + + if outcome % 2 == 1: + x_mat = Ops.X + rho_out = x_mat @ rho_out @ x_mat.conj().T + + return rho_out # type: ignore[no-any-return] + + @staticmethod + def _amplitude_damping_x_error_analytical(gamma: float, outcome: Outcome) -> npt.NDArray[np.complex128]: + """Compute expected hpat output for AmplitudeDampingNoiseModel(x_error_prob=gamma). + + AD is applied to the output node AFTER the X correction (per the + noise model's cmds -> ApplyNoise ordering). For hpat the post-correction + output is |0>, and AD is the identity on |0>. + """ + rho_plus = 0.5 * np.ones((2, 2), dtype=np.complex128) + rho_plus[0, 1] = 0.5 + rho_plus[1, 0] = 0.5 + rho_initial = np.kron(rho_plus, rho_plus) + + cz = Ops.CZ + rho_cz = cz @ rho_initial @ cz.conj().T + + if outcome == 0: + proj = np.array([[0.5, 0.5], [0.5, 0.5]], dtype=np.complex128) + else: + proj = np.array([[0.5, -0.5], [-0.5, 0.5]], dtype=np.complex128) + proj_full = np.kron(proj, np.eye(2, dtype=np.complex128)) + rho_meas = proj_full @ rho_cz @ proj_full.conj().T + prob = np.trace(rho_meas).real + rho_out = TestNoisyDensityMatrixBackend._partial_trace_qubit0(rho_meas) / prob + + if outcome % 2 == 1: + x_mat = Ops.X + rho_out = x_mat @ rho_out @ x_mat.conj().T + + k0, k1 = TestNoisyDensityMatrixBackend._ad_kraus(gamma) + return k0 @ rho_out @ k0.conj().T + k1 @ rho_out @ k1.conj().T # type: ignore[no-any-return] + + # --- Tests using analytical formulas --- + + @pytest.mark.parametrize("outcome", [0, 1]) + def test_amplitude_damping_measure_analytical(self, fx_rng: Generator, outcome: Outcome) -> None: + """Compare measure_channel simulation against analytical Kraus computation.""" + gamma = fx_rng.uniform(0.0, 0.9) + hadamardpattern = hpat() + expected = TestNoisyDensityMatrixBackend._amplitude_damping_measure_analytical(gamma, outcome) + res = hadamardpattern.simulate_pattern( + backend="densitymatrix", + noise_model=AmplitudeDampingNoiseModel(measure_channel_prob=gamma), + branch_selector=ConstBranchSelector(outcome), + rng=fx_rng, + ) + assert isinstance(res, DensityMatrix) + assert np.allclose(res.rho, expected, atol=1e-10) + + @pytest.mark.parametrize("outcome", [0, 1]) + def test_amplitude_damping_prepare_analytical(self, fx_rng: Generator, outcome: Outcome) -> None: + """Compare prepare_error simulation against analytical Kraus computation.""" + gamma = fx_rng.uniform(0.0, 0.9) + hadamardpattern = hpat() + expected = TestNoisyDensityMatrixBackend._amplitude_damping_prepare_analytical(gamma, outcome) + res = hadamardpattern.simulate_pattern( + backend="densitymatrix", + noise_model=AmplitudeDampingNoiseModel(prepare_error_prob=gamma), + branch_selector=ConstBranchSelector(outcome), + rng=fx_rng, + ) + assert isinstance(res, DensityMatrix) + assert np.allclose(res.rho, expected, atol=1e-10) + + @pytest.mark.parametrize("outcome", [0, 1]) + def test_amplitude_damping_entanglement_analytical(self, fx_rng: Generator, outcome: Outcome) -> None: + """Compare entanglement_error simulation against analytical Kraus computation.""" + gamma = fx_rng.uniform(0.0, 0.9) + hadamardpattern = hpat() + expected = TestNoisyDensityMatrixBackend._amplitude_damping_entanglement_analytical(gamma, outcome) + res = hadamardpattern.simulate_pattern( + backend="densitymatrix", + noise_model=AmplitudeDampingNoiseModel(entanglement_error_prob=gamma), + branch_selector=ConstBranchSelector(outcome), + rng=fx_rng, + ) + assert isinstance(res, DensityMatrix) + assert np.allclose(res.rho, expected, atol=1e-10) + + @pytest.mark.parametrize("outcome", [0, 1]) + def test_amplitude_damping_x_error_analytical(self, fx_rng: Generator, outcome: Outcome) -> None: + """Compare x_error simulation against analytical Kraus computation. + + Note: for the hpat pattern, the post-correction state of the output + node is always |0⟩, and AD does not affect |0⟩. So the output is + always |0⟩ regardless of gamma. + """ + gamma = fx_rng.uniform(0.0, 0.9) + hadamardpattern = hpat() + expected = TestNoisyDensityMatrixBackend._amplitude_damping_x_error_analytical(gamma, outcome) + res = hadamardpattern.simulate_pattern( + backend="densitymatrix", + noise_model=AmplitudeDampingNoiseModel(x_error_prob=gamma), + branch_selector=ConstBranchSelector(outcome), + rng=fx_rng, + ) + assert isinstance(res, DensityMatrix) + assert np.allclose(res.rho, expected, atol=1e-10)