diff --git a/docs/source/channels.rst b/docs/source/channels.rst index 3bfe1a83a..cee65e930 100644 --- a/docs/source/channels.rst +++ b/docs/source/channels.rst @@ -14,10 +14,14 @@ Kraus channel .. autofunction:: depolarising_channel +.. autofunction:: amplitude_damping_channel + .. autofunction:: pauli_channel .. autofunction:: two_qubit_depolarising_channel +.. autofunction:: two_qubit_amplitude_damping_channel + .. autofunction:: two_qubit_depolarising_tensor_channel @@ -52,3 +56,14 @@ Noise model classes .. autoclass:: DepolarisingNoiseModel :members: + +.. currentmodule:: graphix.noise_models.amplitude_damping + +.. autoclass:: AmplitudeDampingNoise + :members: + +.. autoclass:: TwoQubitAmplitudeDampingNoise + :members: + +.. autoclass:: AmplitudeDampingNoiseModel + :members: 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..bfc541ab1 100644 --- a/graphix/channels.py +++ b/graphix/channels.py @@ -258,6 +258,62 @@ def two_qubit_depolarising_channel(prob: float) -> KrausChannel: ) +def amplitude_damping_channel(gamma: float) -> KrausChannel: + r"""Single-qubit amplitude damping channel. + + .. math:: + K_0 = \begin{pmatrix} 1 & 0 \\ 0 & \sqrt{1 - \gamma} \end{pmatrix}, + \quad + K_1 = \begin{pmatrix} 0 & \sqrt{\gamma} \\ 0 & 0 \end{pmatrix} + + Parameters + ---------- + gamma : float + Normalized damping parameter, between 0 and 1. + + Returns + ------- + :class:`graphix.channels.KrausChannel` + Kraus channel for amplitude damping. + """ + 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 KrausChannel( + [ + KrausData(1.0, k0), + KrausData(1.0, k1), + ] + ) + + +def two_qubit_amplitude_damping_channel(gamma: float) -> KrausChannel: + r"""Two-qubit amplitude damping channel. + + Independent amplitude damping on each qubit, with Kraus operators + :math:`K_i \otimes K_j` for :math:`i, j \in \{0, 1\}`. + + Parameters + ---------- + gamma : float + Normalized damping parameter, between 0 and 1. + + Returns + ------- + :class:`graphix.channels.KrausChannel` + Kraus channel for two-qubit amplitude damping. + """ + 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 KrausChannel( + [ + KrausData(1.0, np.kron(k0, k0)), + KrausData(1.0, np.kron(k0, k1)), + KrausData(1.0, np.kron(k1, k0)), + KrausData(1.0, np.kron(k1, k1)), + ] + ) + + def two_qubit_depolarising_tensor_channel(prob: float) -> KrausChannel: r"""Two-qubit tensor channel of single-qubit depolarising channels with same probability. 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..95eb79794 --- /dev/null +++ b/graphix/noise_models/amplitude_damping.py @@ -0,0 +1,151 @@ +"""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.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.measurements import Outcome + from graphix.noise_models.noise_model import CommandOrNoise + + +class AmplitudeDampingNoise(Noise): + """One-qubit amplitude damping noise with parameter ``gamma``.""" + + gamma = Probability() + + def __init__(self, gamma: float) -> None: + """Initialize one-qubit amplitude damping noise. + + Parameters + ---------- + gamma : float + Normalized damping parameter, 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-qubits amplitude damping noise with parameter ``gamma``.""" + + gamma = Probability() + + def __init__(self, gamma: float) -> None: + """Initialize two-qubit amplitude damping noise. + + Parameters + ---------- + gamma : float + Normalized damping parameter, 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_density_matrix.py b/tests/test_density_matrix.py index 7c67909c5..4f6884fef 100644 --- a/tests/test_density_matrix.py +++ b/tests/test_density_matrix.py @@ -12,7 +12,7 @@ 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 from graphix.fundamentals import ANGLE_PI, Plane from graphix.ops import Ops from graphix.sim.density_matrix import DensityMatrix, DensityMatrixBackend @@ -736,6 +736,56 @@ 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: + dm = DensityMatrix(randobj.rand_dm(2, fx_rng)) + rho_test = dm.rho.astype(np.complex128) + + gamma = fx_rng.uniform() + ad_channel = amplitude_damping_channel(gamma) + + assert isinstance(ad_channel, KrausChannel) + + dm.apply_channel(ad_channel, [0]) + + 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) + expected_dm = k0 @ rho_test @ k0.conj().T + k1 @ rho_test @ k1.conj().T + + assert np.allclose(expected_dm.trace(), 1.0) + assert np.allclose(dm.rho, expected_dm) + + 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) + + assert isinstance(ad_channel, KrausChannel) + + dm.apply_channel(ad_channel, [i]) + + 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) + psi_tensor = psi.reshape((2,) * nqubits) + + psi_k0 = np.tensordot(k0, psi_tensor, (1, i)) + psi_k0 = np.moveaxis(psi_k0, 0, i) + + psi_k1 = np.tensordot(k1, psi_tensor, (1, i)) + psi_k1 = np.moveaxis(psi_k1, 0, i) + + psi_k0 = np.reshape(psi_k0, (2**nqubits)) + psi_k1 = np.reshape(psi_k1, (2**nqubits)) + expected_dm = np.outer(psi_k0, psi_k0.conj()) + np.outer(psi_k1, psi_k1.conj()) + + assert np.allclose(expected_dm.trace(), 1.0) + assert np.allclose(dm.rho, expected_dm) + 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..f352dd90b 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, ) @@ -150,6 +152,51 @@ def test_2_qubit_depolarising_channel(self, fx_rng: Generator) -> None: assert np.allclose(depol_channel_2_qubit[i].coef, data[i].coef) assert np.allclose(depol_channel_2_qubit[i].operator, data[i].operator) + def test_amplitude_damping_channel(self, fx_rng: Generator) -> None: + gamma = fx_rng.uniform() + 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) + data = [ + KrausData(1.0, k0), + KrausData(1.0, k1), + ] + + 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) + + @pytest.mark.parametrize("gamma", [-0.1, 1.1]) + def test_amplitude_damping_channel_invalid_gamma(self, gamma: float) -> None: + with pytest.raises(ValueError, match="The specified channel is not normalized"): + amplitude_damping_channel(gamma) + + def test_2_qubit_amplitude_damping_channel(self, fx_rng: Generator) -> None: + gamma = fx_rng.uniform() + 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) + data = [ + KrausData(1.0, np.kron(k0, k0)), + KrausData(1.0, np.kron(k0, k1)), + KrausData(1.0, np.kron(k1, k0)), + KrausData(1.0, np.kron(k1, k1)), + ] + + 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) + def test_2_qubit_depolarising_tensor_channel(self, fx_rng: Generator) -> None: prob = fx_rng.uniform() data = [ diff --git a/tests/test_noise_model.py b/tests/test_noise_model.py index 3babd2a4c..76ba4e183 100644 --- a/tests/test_noise_model.py +++ b/tests/test_noise_model.py @@ -6,12 +6,16 @@ import pytest from graphix import Pattern -from graphix.command import CommandKind, M, N +from graphix.command import CommandKind, M, N, X, Z from graphix.noise_models import ( + AmplitudeDampingNoise, + AmplitudeDampingNoiseModel, ApplyNoise, ComposeNoiseModel, DepolarisingNoise, DepolarisingNoiseModel, + NoiseModel, + TwoQubitAmplitudeDampingNoise, TwoQubitDepolarisingNoise, ) from graphix.noise_models.noise_model import NoiselessNoiseModel @@ -100,7 +104,64 @@ 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: +def test_amplitude_damping_noise_model_transpile(fx_rng: Generator) -> None: + nqubits = 5 + depth = 5 + circuit = rand_circuit(nqubits, depth, rng=fx_rng) + pattern = circuit.transpile().pattern + noise_model = AmplitudeDampingNoiseModel( + prepare_error_prob=0.1, + x_error_prob=0.2, + z_error_prob=0.3, + entanglement_error_prob=0.4, + measure_channel_prob=0.5, + ) + noisy_pattern = noise_model.transpile(pattern, rng=fx_rng) + iterator = iter(noisy_pattern) + + def check_noise_command( + cmd: CommandOrNoise, gamma: float, two_qubits: bool, domain: set[int] | None = None + ) -> None: + assert isinstance(cmd, ApplyNoise) + if two_qubits: + assert isinstance(cmd.noise, TwoQubitAmplitudeDampingNoise) + assert cmd.noise.gamma == gamma + else: + assert isinstance(cmd.noise, AmplitudeDampingNoise) + assert cmd.noise.gamma == gamma + assert cmd.domain == domain + + for cmd in pattern: + match cmd.kind: + case CommandKind.M: + check_noise_command(next(iterator), 0.5, False) + assert next(iterator) == cmd + case CommandKind.N: + assert next(iterator) == cmd + check_noise_command(next(iterator), 0.1, False) + case CommandKind.E: + assert next(iterator) == cmd + check_noise_command(next(iterator), 0.4, True) + case CommandKind.X: + assert next(iterator) == cmd + assert isinstance(cmd, X) + check_noise_command(next(iterator), 0.2, False, cmd.domain) + case CommandKind.Z: + assert next(iterator) == cmd + assert isinstance(cmd, Z) + check_noise_command(next(iterator), 0.3, False, cmd.domain) + case _: + assert next(iterator) == cmd + + +@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() @@ -108,7 +169,6 @@ def test_confuse_result(fx_rng: Generator) -> None: 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 diff --git a/tests/test_noisy_amplitude_damping_density_matrix.py b/tests/test_noisy_amplitude_damping_density_matrix.py new file mode 100644 index 000000000..3c20d1468 --- /dev/null +++ b/tests/test_noisy_amplitude_damping_density_matrix.py @@ -0,0 +1,378 @@ +from __future__ import annotations + +from typing import TYPE_CHECKING + +import numpy as np +import numpy.typing as npt +import pytest + +from graphix.branch_selector import ConstBranchSelector, FixedBranchSelector +from graphix.channels import amplitude_damping_channel +from graphix.command import CommandKind +from graphix.fundamentals import angle_to_rad +from graphix.measurements import Measurement +from graphix.noise_models import AmplitudeDampingNoise, AmplitudeDampingNoiseModel, TwoQubitAmplitudeDampingNoise +from graphix.noise_models.noise_model import NoiselessNoiseModel +from graphix.ops import Ops +from graphix.sim.base_backend import _outcome_to_operator_matrix +from graphix.sim.density_matrix import DensityMatrix +from graphix.states import BasicStates +from graphix.transpiler import Circuit + +if TYPE_CHECKING: + from numpy.random import Generator + + from graphix.fundamentals import Angle + from graphix.measurements import Outcome + from graphix.noise_models.noise_model import Noise + from graphix.pattern import Pattern + + +def apply_amplitude_damping(rho: npt.NDArray[np.complex128], gamma: float) -> npt.NDArray[np.complex128]: + """Apply single-qubit amplitude damping to a density matrix.""" + dm = DensityMatrix(data=rho) + dm.apply_channel(amplitude_damping_channel(gamma), [0]) + return dm.rho.astype(np.complex128) + + +class _ReferenceDensityMatrixSimulator: + """Minimal density-matrix simulator for independent test references.""" + + def __init__(self) -> None: + self.dm = DensityMatrix(nqubit=0) + self.node_index: list[int] = [] + + def add_nodes(self, nodes: list[int]) -> None: + self.dm.add_nodes(len(nodes), BasicStates.PLUS) + self.node_index.extend(nodes) + + def _loc(self, node: int) -> int: + return self.node_index.index(node) + + def apply_noise(self, nodes: list[int], noise: Noise) -> None: + self.dm.apply_noise([self._loc(node) for node in nodes], noise) + + def entangle(self, edge: tuple[int, int]) -> None: + self.dm.entangle((self._loc(edge[0]), self._loc(edge[1]))) + + def measure(self, node: int, measurement: Measurement, outcome: Outcome) -> None: + bloch = measurement.to_bloch() + vec = bloch.plane.polar(bloch.angle) + operator = _outcome_to_operator_matrix(vec, outcome, symbolic=False) + loc = self._loc(node) + self.dm.evolve_single(operator, loc) + self.dm.remove_qubit(loc) + self.node_index.remove(node) + + def correct(self, node: int, operator: npt.NDArray[np.complex128]) -> None: + self.dm.evolve_single(operator, self._loc(node)) + + +def reference_hadamard_rho( + *, + prepare_gamma: float = 0.0, + entangle_gamma: float = 0.0, + measure_gamma: float = 0.0, + x_gamma: float = 0.0, + outcome: Outcome = 0, +) -> npt.NDArray[np.complex128]: + """Return the expected output density matrix for the Hadamard pattern.""" + sim = _ReferenceDensityMatrixSimulator() + sim.add_nodes([0]) + if prepare_gamma: + sim.apply_noise([0], AmplitudeDampingNoise(prepare_gamma)) + sim.add_nodes([1]) + if prepare_gamma: + sim.apply_noise([1], AmplitudeDampingNoise(prepare_gamma)) + sim.entangle((0, 1)) + if entangle_gamma: + sim.apply_noise([0, 1], TwoQubitAmplitudeDampingNoise(entangle_gamma)) + if measure_gamma: + sim.apply_noise([0], AmplitudeDampingNoise(measure_gamma)) + sim.measure(0, Measurement.X, outcome) + if outcome == 1: + sim.correct(1, Ops.X) + if x_gamma: + sim.apply_noise([1], AmplitudeDampingNoise(x_gamma)) + return sim.dm.rho.astype(np.complex128) + + +def reference_rz_rho( + alpha: Angle, + *, + prepare_gamma: float = 0.0, + entangle_gamma: float = 0.0, + measure_gamma: float = 0.0, + x_gamma: float = 0.0, + z_gamma: float = 0.0, + m0_outcome: Outcome = 0, + m1_outcome: Outcome = 0, +) -> npt.NDArray[np.complex128]: + """Return the expected output density matrix for the RZ pattern.""" + sim = _ReferenceDensityMatrixSimulator() + sim.add_nodes([0]) + if prepare_gamma: + sim.apply_noise([0], AmplitudeDampingNoise(prepare_gamma)) + sim.add_nodes([1]) + if prepare_gamma: + sim.apply_noise([1], AmplitudeDampingNoise(prepare_gamma)) + sim.add_nodes([2]) + if prepare_gamma: + sim.apply_noise([2], AmplitudeDampingNoise(prepare_gamma)) + sim.entangle((0, 1)) + if entangle_gamma: + sim.apply_noise([0, 1], TwoQubitAmplitudeDampingNoise(entangle_gamma)) + sim.entangle((1, 2)) + if entangle_gamma: + sim.apply_noise([1, 2], TwoQubitAmplitudeDampingNoise(entangle_gamma)) + if measure_gamma: + sim.apply_noise([0], AmplitudeDampingNoise(measure_gamma)) + sim.measure(0, Measurement.XY(-alpha), m0_outcome) + if measure_gamma: + sim.apply_noise([1], AmplitudeDampingNoise(measure_gamma)) + sim.measure(1, Measurement.X, m1_outcome) + if m1_outcome == 1: + sim.correct(2, Ops.X) + if x_gamma: + sim.apply_noise([2], AmplitudeDampingNoise(x_gamma)) + if m0_outcome == 1: + sim.correct(2, Ops.Z) + if z_gamma: + sim.apply_noise([2], AmplitudeDampingNoise(z_gamma)) + return sim.dm.rho.astype(np.complex128) + + +def rz_exact_res(alpha: Angle) -> npt.NDArray[np.float64]: + rad = angle_to_rad(alpha) + return 0.5 * np.array([[1, np.exp(-1j * rad)], [np.exp(1j * rad), 1]]) + + +def hpat() -> Pattern: + circ = Circuit(1) + circ.h(0) + return circ.transpile().pattern + + +def rzpat(alpha: Angle) -> Pattern: + circ = Circuit(1) + circ.rz(0, alpha) + return circ.transpile().pattern + + +class TestNoisyAmplitudeDampingDensityMatrixBackend: + """Test amplitude damping in the density-matrix backend.""" + + @pytest.mark.filterwarnings("ignore:Simulating using densitymatrix backend with no noise.") + def test_noiseless_noisy_hadamard(self, fx_rng: Generator) -> None: + hadamardpattern = hpat() + noiselessres = hadamardpattern.simulate_pattern(backend="densitymatrix", rng=fx_rng) + noisynoiselessres = hadamardpattern.simulate_pattern( + backend="densitymatrix", + noise_model=NoiselessNoiseModel(), + rng=fx_rng, + ) + assert isinstance(noiselessres, DensityMatrix) + assert np.allclose(noiselessres.rho, np.array([[1.0, 0.0], [0.0, 0.0]])) + assert isinstance(noisynoiselessres, DensityMatrix) + assert np.allclose(noisynoiselessres.rho, np.array([[1.0, 0.0], [0.0, 0.0]])) + + def test_noisy_measure_confuse_hadamard(self, fx_rng: Generator) -> None: + 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]])) + + @pytest.mark.parametrize("outcome", [0, 1]) + def test_noisy_measure_confuse_hadamard_arbitrary(self, fx_rng: Generator, outcome: Outcome) -> None: + hadamardpattern = hpat() + measure_error_pr = fx_rng.random() + res = hadamardpattern.simulate_pattern( + backend="densitymatrix", + noise_model=AmplitudeDampingNoiseModel(measure_error_prob=measure_error_pr), + branch_selector=ConstBranchSelector(outcome), + rng=fx_rng, + ) + assert isinstance(res, DensityMatrix) + assert np.allclose(res.rho, np.array([[1.0, 0.0], [0.0, 0.0]])) or np.allclose( + res.rho, + np.array([[0.0, 0.0], [0.0, 1.0]]), + ) + + def test_noisy_measure_channel_hadamard(self, fx_rng: Generator) -> None: + hadamardpattern = hpat() + measure_channel_gamma = fx_rng.random() + res = hadamardpattern.simulate_pattern( + backend="densitymatrix", + noise_model=AmplitudeDampingNoiseModel(measure_channel_prob=measure_channel_gamma), + rng=fx_rng, + ) + expected = reference_hadamard_rho(measure_gamma=measure_channel_gamma) + assert isinstance(res, DensityMatrix) + assert np.allclose(res.rho, expected) + + @pytest.mark.parametrize("outcome", [0, 1]) + def test_noisy_x_hadamard(self, fx_rng: Generator, outcome: Outcome) -> None: + hadamardpattern = hpat() + x_error_gamma = fx_rng.random() + res = hadamardpattern.simulate_pattern( + backend="densitymatrix", + noise_model=AmplitudeDampingNoiseModel(x_error_prob=x_error_gamma), + branch_selector=ConstBranchSelector(outcome), + rng=fx_rng, + ) + expected = reference_hadamard_rho(x_gamma=x_error_gamma, outcome=outcome) + assert isinstance(res, DensityMatrix) + assert np.allclose(res.rho, expected) + + @pytest.mark.parametrize("outcome", [0, 1]) + def test_noisy_entanglement_hadamard(self, fx_rng: Generator, outcome: Outcome) -> None: + hadamardpattern = hpat() + entanglement_error_gamma = fx_rng.uniform() + res = hadamardpattern.simulate_pattern( + backend="densitymatrix", + noise_model=AmplitudeDampingNoiseModel(entanglement_error_prob=entanglement_error_gamma), + branch_selector=ConstBranchSelector(outcome), + rng=fx_rng, + ) + expected = reference_hadamard_rho(entangle_gamma=entanglement_error_gamma, outcome=outcome) + assert isinstance(res, DensityMatrix) + assert np.allclose(res.rho, expected) + + @pytest.mark.parametrize("outcome", [0, 1]) + def test_noisy_preparation_hadamard(self, fx_rng: Generator, outcome: Outcome) -> None: + hadamardpattern = hpat() + prepare_error_gamma = fx_rng.random() + res = hadamardpattern.simulate_pattern( + backend="densitymatrix", + noise_model=AmplitudeDampingNoiseModel(prepare_error_prob=prepare_error_gamma), + branch_selector=ConstBranchSelector(outcome), + rng=fx_rng, + ) + expected = reference_hadamard_rho(prepare_gamma=prepare_error_gamma, outcome=outcome) + assert isinstance(res, DensityMatrix) + assert np.allclose(res.rho, expected) + + def test_amplitude_damping_channel_preserves_zero(self, fx_rng: Generator) -> None: + gamma = fx_rng.random() + rho = np.array([[1.0, 0.0], [0.0, 0.0]], dtype=np.complex128) + assert np.allclose(apply_amplitude_damping(rho, gamma), rho) + + def test_amplitude_damping_channel_on_one_state(self, fx_rng: Generator) -> None: + gamma = fx_rng.random() + rho = np.array([[0.0, 0.0], [0.0, 1.0]], dtype=np.complex128) + expected = np.array([[gamma, 0.0], [0.0, 1.0 - gamma]]) + assert np.allclose(apply_amplitude_damping(rho, gamma), expected) + + @pytest.mark.filterwarnings("ignore:Simulating using densitymatrix backend with no noise.") + def test_noiseless_noisy_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 np.allclose(noiselessres.rho, rz_exact_res(alpha)) + assert isinstance(noisynoiselessres, DensityMatrix) + assert np.allclose(noisynoiselessres.rho, rz_exact_res(alpha)) + + @pytest.mark.parametrize("outcome", [0, 1]) + def test_noisy_preparation_rz(self, fx_rng: Generator, outcome: Outcome) -> None: + alpha = fx_rng.random() + rzpattern = rzpat(alpha) + prepare_error_gamma = fx_rng.random() + res = rzpattern.simulate_pattern( + backend="densitymatrix", + noise_model=AmplitudeDampingNoiseModel(prepare_error_prob=prepare_error_gamma), + branch_selector=ConstBranchSelector(outcome), + rng=fx_rng, + ) + expected = reference_rz_rho(alpha, prepare_gamma=prepare_error_gamma, m0_outcome=outcome, m1_outcome=outcome) + assert isinstance(res, DensityMatrix) + assert np.allclose(res.rho, expected) + + @pytest.mark.parametrize("outcome", [0, 1]) + def test_noisy_entanglement_rz(self, fx_rng: Generator, outcome: Outcome) -> None: + alpha = fx_rng.random() + rzpattern = rzpat(alpha) + entanglement_error_gamma = fx_rng.uniform() + res = rzpattern.simulate_pattern( + backend="densitymatrix", + noise_model=AmplitudeDampingNoiseModel(entanglement_error_prob=entanglement_error_gamma), + branch_selector=ConstBranchSelector(outcome), + rng=fx_rng, + ) + expected = reference_rz_rho( + alpha, entangle_gamma=entanglement_error_gamma, m0_outcome=outcome, m1_outcome=outcome + ) + assert isinstance(res, DensityMatrix) + assert np.allclose(res.rho, expected) + + def test_noisy_measure_channel_rz(self, fx_rng: Generator) -> None: + alpha = fx_rng.random() + rzpattern = rzpat(alpha) + measure_channel_gamma = fx_rng.random() + res = rzpattern.simulate_pattern( + backend="densitymatrix", + noise_model=AmplitudeDampingNoiseModel(measure_channel_prob=measure_channel_gamma), + rng=fx_rng, + ) + expected = reference_rz_rho(alpha, measure_gamma=measure_channel_gamma) + assert isinstance(res, DensityMatrix) + assert np.allclose(res.rho, expected) + + @pytest.mark.parametrize("x_outcome", [0, 1]) + def test_noisy_x_rz(self, fx_rng: Generator, x_outcome: Outcome) -> None: + alpha = fx_rng.random() + rzpattern = rzpat(alpha) + x_error_gamma = fx_rng.random() + m_nodes = (cmd.node for cmd in rzpattern if cmd.kind == CommandKind.M) + results: dict[int, Outcome] = {next(m_nodes): 0, next(m_nodes): x_outcome} + res = rzpattern.simulate_pattern( + backend="densitymatrix", + noise_model=AmplitudeDampingNoiseModel(x_error_prob=x_error_gamma), + branch_selector=FixedBranchSelector(results), + rng=fx_rng, + ) + expected = reference_rz_rho(alpha, x_gamma=x_error_gamma, m1_outcome=x_outcome) + assert isinstance(res, DensityMatrix) + assert np.allclose(res.rho, expected) + + @pytest.mark.parametrize("outcome_z", [0, 1]) + @pytest.mark.parametrize("outcome_x", [0, 1]) + def test_noisy_z_rz(self, fx_rng: Generator, outcome_z: Outcome, outcome_x: Outcome) -> None: + alpha = fx_rng.random() + rzpattern = rzpat(alpha) + z_error_gamma = fx_rng.random() + results: dict[int, Outcome] = {0: outcome_z, 1: outcome_x} + res = rzpattern.simulate_pattern( + backend="densitymatrix", + noise_model=AmplitudeDampingNoiseModel(z_error_prob=z_error_gamma), + branch_selector=FixedBranchSelector(results), + rng=fx_rng, + ) + expected = reference_rz_rho(alpha, z_gamma=z_error_gamma, m0_outcome=outcome_z, m1_outcome=outcome_x) + assert isinstance(res, DensityMatrix) + assert np.allclose(res.rho, expected) + + @pytest.mark.parametrize("z_outcome", [0, 1]) + @pytest.mark.parametrize("x_outcome", [0, 1]) + def test_noisy_measure_confuse_rz(self, fx_rng: Generator, z_outcome: Outcome, x_outcome: Outcome) -> None: + alpha = fx_rng.random() + rzpattern = rzpat(alpha) + results: dict[int, Outcome] = {0: z_outcome, 1: x_outcome} + res = rzpattern.simulate_pattern( + backend="densitymatrix", + noise_model=AmplitudeDampingNoiseModel(measure_error_prob=1.0), + branch_selector=FixedBranchSelector(results), + rng=fx_rng, + ) + exact = rz_exact_res(alpha) + assert isinstance(res, DensityMatrix) + assert np.allclose(res.rho, Ops.Z @ Ops.X @ exact @ Ops.X @ Ops.Z)