diff --git a/CHANGELOG.md b/CHANGELOG.md index c7c280f1d..6233cf9cb 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 diff --git a/graphix/channels.py b/graphix/channels.py index 8bb0ebc1b..5cefb914e 100644 --- a/graphix/channels.py +++ b/graphix/channels.py @@ -296,3 +296,50 @@ 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(left, right)) for left in (k1, k2) for right in (k1, 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..ea23f9311 --- /dev/null +++ b/graphix/noise_models/amplitude_damping.py @@ -0,0 +1,155 @@ +"""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 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, + 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( + 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 _: # 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: + """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_density_matrix.py b/tests/test_density_matrix.py index 7c67909c5..25bd94d95 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 @@ -569,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() @@ -647,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() @@ -736,6 +742,95 @@ 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 = np.asarray(dm.rho, dtype=np.complex128) + + 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) + + 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) + 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..c3f3f9b58 100644 --- a/tests/test_noise_model.py +++ b/tests/test_noise_model.py @@ -6,15 +6,19 @@ import pytest from graphix import Pattern -from graphix.command import CommandKind, M, N +from graphix.clifford import Clifford +from graphix.command import C, CommandKind, E, M, N, S, T, X, Z from graphix.noise_models import ( + AmplitudeDampingNoise, + AmplitudeDampingNoiseModel, ApplyNoise, ComposeNoiseModel, DepolarisingNoise, DepolarisingNoiseModel, + 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 @@ -100,7 +104,61 @@ 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 _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 out[0].kind == CommandKind.N + _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_noise(out[1], TwoQubitAmplitudeDampingNoise, 0.3, [0, 1]) + + # M: noise applied BEFORE measurement + out = model.command(M(node=0)) + _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_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( + "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,9 +166,103 @@ 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 ) assert measure_method.results[0] == 1 + + +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) + + +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) diff --git a/tests/test_noisy_density_matrix.py b/tests/test_noisy_density_matrix.py index 104923265..6f38919c6 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 @@ -9,10 +10,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 DepolarisingNoiseModel -from graphix.noise_models.noise_model import NoiselessNoiseModel +from graphix.noise_models import AmplitudeDampingNoiseModel, DepolarisingNoiseModel +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 from graphix.transpiler import Circuit if TYPE_CHECKING: @@ -524,3 +526,233 @@ 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: + """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.""" + + @staticmethod + 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), + ] + + @staticmethod + def _kron_all(ops: list[npt.NDArray[np.complex128]]) -> npt.NDArray[np.complex128]: + """Kronecker product of a list of operators, left to right.""" + 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.""" + return sum( + (k @ rho @ k.conj().T for k in self._kraus_operators(gamma)), + np.zeros((2, 2), dtype=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.""" + dim = 2**nqubit + out = np.zeros((dim, dim), dtype=np.complex128) + for k in self._kraus_operators(gamma): + ops = [Ops.I] * nqubit + ops[qubit] = k + full = self._kron_all(ops) + out += full @ rho @ full.conj().T + return out + + 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``.""" + 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 = [Ops.I] * nqubit + ops[qubit_a] = left + ops[qubit_b] = right + full = self._kron_all(ops) + out += full @ rho @ full.conj().T + return out + + @staticmethod + 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 _hadamard_expected(self, step: str, gamma: float, outcome: Outcome) -> npt.NDArray[np.complex128]: + """Output density matrix of the Hadamard pattern with damping at ``step``.""" + 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 = Ops.CZ @ rho @ Ops.CZ + if step == "entangle": + rho = self._damp_two_qubits_in_register(rho, 0, 1, 2, gamma) + if step == "measure": + 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), 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)) + if outcome == 1: + reduced = Ops.X @ reduced @ Ops.X + if step == "xcorr": + reduced = self._damp_one_qubit(reduced, gamma) + return reduced + + @pytest.mark.parametrize("outcome", [0, 1]) + @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_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}), + branch_selector=ConstBranchSelector(outcome), + rng=fx_rng, + ) + assert isinstance(res, DensityMatrix) + 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. + """ + 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(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(Ops.I, Ops.CZ) + rho = cz12 @ rho @ cz12 + if step == "entangle": + rho = self._damp_two_qubits_in_register(rho, 1, 2, 3, gamma) + + if step == "measure": + rho = self._damp_qubit_in_register(rho, 0, 3, gamma) + 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([Ops.I, self._xy_projector(outcome_x, 0.0), Ops.I]) + 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 = Ops.X @ reduced @ Ops.X + if step == "xcorr": + reduced = self._damp_one_qubit(reduced, gamma) + if outcome_z == 1: + 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]) + @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))