diff --git a/source/compiler/qsc_eval/src/backend/noise_tests.rs b/source/compiler/qsc_eval/src/backend/noise_tests.rs index 13449c07af..48f0d4d1a0 100644 --- a/source/compiler/qsc_eval/src/backend/noise_tests.rs +++ b/source/compiler/qsc_eval/src/backend/noise_tests.rs @@ -10,7 +10,7 @@ use crate::{ use expect_test::{Expect, expect}; use num_bigint::BigUint; use num_complex::Complex; -use qdk_simulators::noise_config::{NoiseConfig, NoiseTable, encode_pauli}; +use qdk_simulators::noise_config::{LossPolicy, NoiseConfig, NoiseTable, encode_pauli}; use std::fmt::Write; #[test] @@ -258,6 +258,7 @@ fn noise_config_with_single_qubit_fault( qubits: 1, pauli_strings: vec![encode_pauli(pauli)], probabilities: vec![1.0], + on_loss: LossPolicy::Skip, }; set_gate(&mut config, table); config @@ -274,6 +275,7 @@ fn noise_config_with_two_qubit_fault( qubits: 2, pauli_strings: vec![encode_pauli(pauli)], probabilities: vec![1.0], + on_loss: LossPolicy::Skip, }; set_gate(&mut config, table); config @@ -525,6 +527,7 @@ fn noise_config_mz_with_loss() { let mut config = NoiseConfig::NOISELESS; config.mz = NoiseTable { qubits: 1, + on_loss: LossPolicy::Skip, pauli_strings: vec![encode_pauli("L")], probabilities: vec![1.0], }; @@ -545,6 +548,7 @@ fn noise_config_gate_loss_causes_measurement_loss() { let mut config = NoiseConfig::NOISELESS; config.x = NoiseTable { qubits: 1, + on_loss: LossPolicy::Skip, pauli_strings: vec![encode_pauli("L")], probabilities: vec![1.0], }; diff --git a/source/qdk_package/qdk/_native.pyi b/source/qdk_package/qdk/_native.pyi index d4d0636e25..d6c10360ad 100644 --- a/source/qdk_package/qdk/_native.pyi +++ b/source/qdk_package/qdk/_native.pyi @@ -868,6 +868,25 @@ class QirInstruction: ... class IdleNoiseParams: s_probability: float +class LossPolicy(Enum): + """ + Specifies the behavior of a multi-qubit gate when at least one of its + qubit operands is lost. + """ + + # If any operand of a gate is lost, skip the gate entirely. + SKIP: int + # If any operand of a gate is lost, propagate the loss to the other operands. + PROPAGATE: int + # For multi-qubit rotations, degrade the unitary to its single-qubit version + # on the surviving operand (e.g. rxx -> rx). Falls back to SKIP for gates with + # no single-qubit reduction (cx, cy, cz, swap, and single-qubit gates). + DEGRADE: int + # Skip the gate and instead apply an S adjoint to each surviving operand. + RESIDUAL_S_DAGGER: int + # Apply the unitary anyway, ignoring the loss. + APPLY_ANYWAY: int + class NoiseTable: # Deprecated. Setting `loss` distributes the per-qubit loss probability # across the correlated loss fault strings ('L' for a single-qubit @@ -876,6 +895,7 @@ class NoiseTable: # `loss` reconstructs that per-qubit probability. Prefer setting the loss # fault strings directly via `set_pauli_noise`. loss: float + on_loss: LossPolicy def __init__(self, num_qubits: int): """ diff --git a/source/qdk_package/qdk/simulation/__init__.py b/source/qdk_package/qdk/simulation/__init__.py index 8e5701332c..8e93ba3f2e 100644 --- a/source/qdk_package/qdk/simulation/__init__.py +++ b/source/qdk_package/qdk/simulation/__init__.py @@ -15,6 +15,10 @@ to individual gate intrinsics to model depolarizing, bit-flip, phase-flip, or correlated noise channels. +- :class:`~qdk.simulation.LossPolicy` — selects how a gate behaves when one of + its qubit operands is lost. Assign it to a noise table's ``on_loss`` attribute + (e.g. ``noise.cx.on_loss = LossPolicy.SKIP``). + - :func:`~qdk.simulation.run_qir` — simulates QIR as given in one of three backend simulators: clifford, gpu or cpu. @@ -26,7 +30,7 @@ """ from .._device._atom import NeutralAtomDevice -from ._simulation import NoiseConfig, run_qir +from ._simulation import NoiseConfig, LossPolicy, run_qir from ._noisy_simulator import ( NoisySimulatorError, DensityMatrixSimulator, @@ -40,6 +44,7 @@ __all__ = [ "NeutralAtomDevice", "NoiseConfig", + "LossPolicy", "run_qir", "NoisySimulatorError", "Operation", diff --git a/source/qdk_package/qdk/simulation/_simulation.py b/source/qdk_package/qdk/simulation/_simulation.py index 4aa6bf367d..68847400ac 100644 --- a/source/qdk_package/qdk/simulation/_simulation.py +++ b/source/qdk_package/qdk/simulation/_simulation.py @@ -15,6 +15,7 @@ run_cpu_adaptive, run_cpu_full_state, NoiseConfig, + LossPolicy, GpuContext, try_create_gpu_adapter, Result, diff --git a/source/qdk_package/src/interpreter.rs b/source/qdk_package/src/interpreter.rs index a77fc2064c..ce90c14e50 100644 --- a/source/qdk_package/src/interpreter.rs +++ b/source/qdk_package/src/interpreter.rs @@ -23,7 +23,7 @@ use crate::{ }, noisy_simulator::register_noisy_simulator_submodule, qir_simulation::{ - IdleNoiseParams, NoiseConfig, NoiseTable, QirInstruction, QirInstructionId, + IdleNoiseParams, LossPolicy, NoiseConfig, NoiseTable, QirInstruction, QirInstructionId, cpu_simulators::{ run_clifford, run_clifford_adaptive, run_cpu_adaptive, run_cpu_full_state, }, @@ -104,6 +104,7 @@ fn verify_classes_are_sendable() { is_send::(); is_send::(); is_send::(); + is_send::(); } #[pymodule] @@ -133,6 +134,7 @@ fn _native<'a>(py: Python<'a>, m: &Bound<'a, PyModule>) -> PyResult<()> { m.add_class::()?; m.add_class::()?; m.add_class::()?; + m.add_class::()?; m.add_function(wrap_pyfunction!(physical_estimates, m)?)?; m.add_function(wrap_pyfunction!(run_clifford, m)?)?; m.add_function(wrap_pyfunction!(try_create_gpu_adapter, m)?)?; diff --git a/source/qdk_package/src/qir_simulation.rs b/source/qdk_package/src/qir_simulation.rs index 7b819b63ca..ada487cfd8 100644 --- a/source/qdk_package/src/qir_simulation.rs +++ b/source/qdk_package/src/qir_simulation.rs @@ -9,7 +9,7 @@ use crate::qir_simulation::correlated_noise::parse_noise_table; use num_traits::{Float, Unsigned}; use pyo3::{ - Bound, FromPyObject, Py, PyRef, PyResult, Python, + Bound, FromPyObject, Py, PyAny, PyRef, PyResult, Python, exceptions::{PyAttributeError, PyKeyError, PyTypeError, PyValueError}, pybacked::PyBackedStr, pyclass, pymethods, @@ -88,6 +88,47 @@ pub enum QirInstruction { ), } +/// Specifies the behavior of a multi-qubit gate when at least one of its +/// qubit operands is lost. Mirrors [`qdk_simulators::noise_config::LossPolicy`]. +#[derive(Clone, Copy, Debug, PartialEq, Eq)] +#[pyclass(eq, eq_int, from_py_object, module = "qdk._native")] +pub enum LossPolicy { + #[pyo3(name = "SKIP")] + Skip = 0, + #[pyo3(name = "PROPAGATE")] + Propagate = 1, + #[pyo3(name = "DEGRADE")] + Degrade = 2, + #[pyo3(name = "RESIDUAL_S_DAGGER")] + ResidualSDagger = 3, + #[pyo3(name = "APPLY_ANYWAY")] + ApplyAnyway = 4, +} + +impl From for qdk_simulators::noise_config::LossPolicy { + fn from(value: LossPolicy) -> Self { + match value { + LossPolicy::Skip => Self::Skip, + LossPolicy::Propagate => Self::Propagate, + LossPolicy::Degrade => Self::Degrade, + LossPolicy::ResidualSDagger => Self::ResidualSDagger, + LossPolicy::ApplyAnyway => Self::ApplyAnyway, + } + } +} + +impl From for LossPolicy { + fn from(value: qdk_simulators::noise_config::LossPolicy) -> Self { + match value { + qdk_simulators::noise_config::LossPolicy::Skip => Self::Skip, + qdk_simulators::noise_config::LossPolicy::Propagate => Self::Propagate, + qdk_simulators::noise_config::LossPolicy::Degrade => Self::Degrade, + qdk_simulators::noise_config::LossPolicy::ResidualSDagger => Self::ResidualSDagger, + qdk_simulators::noise_config::LossPolicy::ApplyAnyway => Self::ApplyAnyway, + } + } +} + #[derive(Debug)] #[pyclass(module = "qdk._native")] pub struct NoiseConfig { @@ -357,6 +398,8 @@ impl From for IdleNoiseParams { #[pyclass(from_py_object, module = "qdk._native")] pub struct NoiseTable { qubits: u32, + #[pyo3(get)] + pub on_loss: LossPolicy, pauli_noise: FxHashMap, } @@ -522,6 +565,7 @@ impl NoiseTable { NoiseTable { qubits: num_qubits, pauli_noise: FxHashMap::default(), + on_loss: LossPolicy::Skip, } } @@ -553,11 +597,22 @@ E.g.: `noise_config.cz.IL`" /// /// for arbitrary pauli fields. Setting an element that was /// previously set overrides that entry with the new value. - fn __setattr__(&mut self, name: &str, value: Probability) -> PyResult<()> { - if name == "loss" { - self.set_loss(value) - } else { - self.set_pauli_noise_elt(&name.to_uppercase(), value) + /// + /// The `on_loss` attribute is special-cased to accept a `LossPolicy`. + fn __setattr__(&mut self, name: &str, value: &Bound<'_, PyAny>) -> PyResult<()> { + match name { + "on_loss" => { + if self.qubits < 2 { + Err(PyAttributeError::new_err( + "Loss policies only apply to multi-qubit gates.".to_string(), + )) + } else { + self.on_loss = value.extract::()?; + Ok(()) + } + } + "loss" => self.set_loss(value.extract::()?), + _ => self.set_pauli_noise_elt(&name.to_uppercase(), value.extract::()?), } } @@ -649,6 +704,7 @@ impl From for qdk_simulators::noise_config::NoiseTable qubits: value.qubits, pauli_strings, probabilities, + on_loss: value.on_loss.into(), } } } @@ -667,6 +723,7 @@ fn from_noise_table_ref( qubits: value.qubits, pauli_strings, probabilities, + on_loss: value.on_loss.into(), } } @@ -686,6 +743,7 @@ impl From> for NoiseTable NoiseTable { qubits: value.qubits, pauli_noise, + on_loss: value.on_loss.into(), } } } diff --git a/source/qdk_package/src/qir_simulation/correlated_noise.rs b/source/qdk_package/src/qir_simulation/correlated_noise.rs index faff78f93b..bc1add013e 100644 --- a/source/qdk_package/src/qir_simulation/correlated_noise.rs +++ b/source/qdk_package/src/qir_simulation/correlated_noise.rs @@ -8,7 +8,7 @@ use rustc_hash::FxHashMap; use std::fmt; use std::str::FromStr; -use crate::qir_simulation::NoiseTable; +use crate::qir_simulation::{LossPolicy, NoiseTable}; /// Errors that can occur while parsing a noise-table CSV. #[derive(Debug)] @@ -84,6 +84,7 @@ pub fn parse_noise_table(contents: &str) -> Result { return Ok(NoiseTable { qubits: qubits.unwrap_or(0), pauli_noise, + on_loss: LossPolicy::Skip, }); } @@ -156,6 +157,7 @@ pub fn parse_noise_table(contents: &str) -> Result { Ok(NoiseTable { qubits, pauli_noise, + on_loss: LossPolicy::Skip, }) } diff --git a/source/qdk_package/tests/test_clifford_simulator.py b/source/qdk_package/tests/test_clifford_simulator.py index 04417c9c31..80d8801ccd 100644 --- a/source/qdk_package/tests/test_clifford_simulator.py +++ b/source/qdk_package/tests/test_clifford_simulator.py @@ -303,147 +303,69 @@ def test_clifford_run_no_noise(): assert output == [[Result.Zero] * 16], "Expected result of 0s with pi/2 angles." +QSHARP_OP_25_QUBITS = """ +operation Test() : Result[] { + use qs = Qubit[25]; X(qs[0]); CZ(qs[23], qs[24]); MResetEachZ(qs) +}""" + + def test_clifford_run_bitflip_noise(): """Bitflip noise for Clifford simulator.""" qsharp.init(target_profile=TargetProfile.Base) - qsharp.eval(read_file_relative("CliffordIsing.qs")) + qsharp.eval(QSHARP_OP_25_QUBITS) - p_noise = 0.005 + p_noise = 0.2 noise = NoiseConfig() - noise.rx.set_bitflip(p_noise) - noise.rzz.set_pauli_noise("XX", p_noise) - noise.mresetz.set_bitflip(p_noise) - - output = qsharp.run( - "IsingModel2DEvolution(4, 4, PI() / 2.0, PI() / 2.0, 10.0, 10)", - shots=10_000, - noise=noise, - seed=17, - type="clifford", - ) - result = [result_array_to_string(cast(Sequence[Result], x)) for x in output] - # Reasonable results obtained from manual run + noise.x.set_bitflip(p_noise) + noise.cz.set_pauli_noise("XX", p_noise) expect = { - "0000000001000000": 0.0084, - "0000010000000000": 0.0079, - "0001000000000000": 0.0087, - "0000100000000000": 0.0096, - "0000000000000000": 0.1412, - "0011000000000000": 0.0066, - "0000000001100000": 0.0082, - "1100000000000000": 0.0072, - "0000000000000011": 0.0083, - "0000000000100010": 0.0091, - "0000000000000010": 0.0074, - "0100010000000000": 0.0058, - "0000000001000100": 0.0078, - "0000011000000000": 0.0067, - "0010000000000000": 0.0089, - "0000000000000110": 0.0085, - "0000000000100000": 0.0091, - "0000000010001000": 0.007, - "0000001000000000": 0.008, - "0000100010000000": 0.0078, - "1000100000000000": 0.0086, - "1000000000000000": 0.0067, - "0000000000010001": 0.0073, - "0001000100000000": 0.0085, - "0000110000000000": 0.0075, - "0000000000000001": 0.0076, - "0110000000000000": 0.0073, - "0010001000000000": 0.0068, - "0100000000000000": 0.0087, - "0000000100000000": 0.0066, - "0000010001000000": 0.007, - "0000000000000100": 0.0068, - "0000001000100000": 0.0067, - "0000000011000000": 0.0102, - "0000000000010000": 0.0087, - "0000000000110000": 0.0081, - "0000000010000000": 0.008, - "0000000100010000": 0.008, - "0000001100000000": 0.0075, - "0000000000001000": 0.0088, - "0000000000001100": 0.0066, + "1000000000000000000000000": (1 - p_noise) ** 2, # No noise + "0000000000000000000000000": p_noise * (1 - p_noise), # X bitflip + "1000000000000000000000011": (1 - p_noise) * p_noise, # CZ bitflip + "0000000000000000000000011": p_noise**2, # X & CZ bitflip } + + output = qsharp.run("Test()", shots=500, noise=noise, seed=17, type="clifford") + result = [result_array_to_string(cast(Sequence[Result], x)) for x in output] expect_distribution( result, expect, - tolerance=0.005, + tolerance=0.02, ) # Same execution should work with the operation itself. - output = qsharp.run( - qdk.code.IsingModel2DEvolution, - 10_000, - 4, - 4, - math.pi / 2, - math.pi / 2, - 10.0, - 10, - noise=noise, - seed=17, - type="clifford", - ) + output = qsharp.run(qdk.code.Test, 500, noise=noise, seed=17, type="clifford") result = [result_array_to_string(cast(Sequence[Result], x)) for x in output] expect_distribution( result, expect, - tolerance=0.005, + tolerance=0.02, ) def test_clifford_run_mixed_noise(): qsharp.init(target_profile=TargetProfile.Base) - qsharp.eval(read_file_relative("CliffordIsing.qs")) + qsharp.eval(QSHARP_OP_25_QUBITS) + p_noise = 0.2 noise = NoiseConfig() - noise.rx.set_bitflip(0.008) - noise.rx.loss = 0.005 - noise.rzz.set_depolarizing(0.008) - noise.rzz.loss = 0.005 + noise.x.set_bitflip(p_noise) + noise.cz.XI = p_noise + noise.cz.IL = p_noise - output = qsharp.run( - "IsingModel2DEvolution(4, 4, PI() / 2.0, PI() / 2.0, 4.0, 4)", - shots=10_000, - noise=noise, - seed=228, - type="clifford", - ) + output = qsharp.run("Test()", shots=500, noise=noise, seed=17, type="clifford") result = [result_array_to_string(cast(Sequence[Result], x)) for x in output] expect_distribution( result, - # Reasonable results obtained from manual run { - "0000000000-00000": 0.01, - "0000000000001000": 0.0055, - "000000000-000000": 0.0104, - "0000000000000000": 0.0854, - "0100000000000000": 0.0062, - "0000-00000000000": 0.0098, - "-000000000000000": 0.0066, - "0-00000000000000": 0.0084, - "00000-0000000000": 0.0116, - "00000000000-0000": 0.0069, - "00-0000000000000": 0.0104, - "0000000001000000": 0.0057, - "00000000-0000000": 0.0108, - "0010000000000000": 0.0054, - "000000-000000000": 0.0113, - "0000000000010000": 0.0067, - "00000000000000-0": 0.0092, - "000000000000-000": 0.0072, - "0000000000100000": 0.0074, - "0000000010000000": 0.0056, - "0000010000000000": 0.0065, - "0001000000000000": 0.0052, - "0000000000000-00": 0.0087, - "0000000-00000000": 0.0081, - "000000000000000-": 0.0052, - "0000000100000000": 0.0052, + "1000000000000000000000000": (1 - p_noise) * (1 - 2 * p_noise), # No noise + "0000000000000000000000000": p_noise * (1 - 2 * p_noise), # X bitflip + "1000000000000000000000010": (1 - p_noise) * p_noise, # CZ bitflip + "100000000000000000000000-": (1 - p_noise) * p_noise, # CZ loss + "0000000000000000000000010": p_noise**2, # X bitflip + CZ bitflip + "000000000000000000000000-": p_noise**2, # X bitflip + CZ loss }, - tolerance=0.005, + tolerance=0.02, ) diff --git a/source/qdk_package/tests/test_simulators_gates_noisy.py b/source/qdk_package/tests/test_simulators_gates_noisy.py index 24030677aa..7d2964489e 100644 --- a/source/qdk_package/tests/test_simulators_gates_noisy.py +++ b/source/qdk_package/tests/test_simulators_gates_noisy.py @@ -7,7 +7,7 @@ from qdk import qsharp from qdk._interpreter import compile from qdk import Result, TargetProfile -from qdk.simulation import run_qir as _run_qir, NoiseConfig +from qdk.simulation import run_qir as _run_qir, NoiseConfig, LossPolicy from qdk.simulation._simulation import try_create_gpu_adapter from typing import Literal, List, Optional, TypeAlias @@ -256,6 +256,239 @@ def test_two_qubit_loss(sim_type): ) +# =========================================================================== +# Loss-policy (on_loss) tests +# =========================================================================== +# +# These exercise the per-gate `NoiseConfig..on_loss` behavior. +# +# A qubit is lost deterministically by giving a single-qubit gate a loss +# probability of 1.0 and then applying that gate. The gate under test then sees +# a lost operand and applies its configured policy. All outcomes are +# deterministic, so a single shot is sufficient. +# +# Need to test: CX, CY, CZ, RXX, RYY, RZZ, SWAP +# with: SKIP, PROPAGATE, DEGRADE, RESIDUAL_S_DAGGER +# +# TEST 0: C*, R**, and SWAP default loss policies +# TEST 1: C*, R**, and SWAP with SKIP do not apply unitary +# TEST 2: C*, R**, and SWAP with PROPAGATE lose first qubit +# TEST 3: C*, R**, and SWAP with PROPAGATE lose second qubit +# TEST 4: C* and SWAP with DEGRADE behave like skip +# TEST 5: R** with DEGRADE behave like R* +# TEST 6: C*, R**, and SWAP with RESIDUAL_S_DAGGER do not apply unitary +# TEST 7: SWAP always exchanges loss flags and qubit states + +# Two-qubit gate-call fragments, grouped by how they reduce when exactly one +# operand is lost. Each entry is (NoiseConfig attribute, Q# gate call). +CONTROLLED_GATES = [ + ("cx", "CNOT(qs[0], qs[1])"), + ("cy", "CY(qs[0], qs[1])"), + ("cz", "CZ(qs[0], qs[1])"), +] +ROTATION_GATES = [ + ("rxx", "Rxx(Std.Math.PI(), qs[0], qs[1])"), + ("ryy", "Ryy(Std.Math.PI(), qs[0], qs[1])"), + ("rzz", "Rzz(Std.Math.PI(), qs[0], qs[1])"), +] +SWAP_GATE = ("swap", "SWAP(qs[0], qs[1])") +ALL_GATES = CONTROLLED_GATES + ROTATION_GATES + [SWAP_GATE] + +CONTROLLED_IDS = [attr for attr, _ in CONTROLLED_GATES] +ROTATION_IDS = [attr for attr, _ in ROTATION_GATES] +SWAP_ID = SWAP_GATE[0] +ALL_IDS = [attr for attr, _ in ALL_GATES] + +# Rotation gates under DEGRADE reduce to their single-qubit version on the +# survivor. With theta = PI the degraded rotation flips the survivor's measured +# bit to 1, but Rz only adds phase, so Rzz is prepared/measured in the X basis. +ROTATION_DEGRADE_RECIPES = [ + ("rxx", "Rxx(Std.Math.PI(), qs[0], qs[1])", "", ""), + ("ryy", "Ryy(Std.Math.PI(), qs[0], qs[1])", "", ""), + ("rzz", "Rzz(Std.Math.PI(), qs[0], qs[1])", "H(qs[1]);", "H(qs[1]);"), +] + + +def run_loss_policy_scenario( + gate: str, + sim_type: SimType, + *, + attr: str = "", + on_loss=None, + prep: str = "", + post: str = "", + lose: int = 0, +) -> str: + """ + Lose one operand of a two-qubit gate deterministically, apply the gate, and + measure both qubits. + + The qubit at index *lose* is taken out via a Y gate configured with + ``loss = 1.0``; the survivor can therefore be prepared with any non-Y gate + through *prep* (and post-processed through *post*). Returns the single + deterministic shot as a two-character string for + ``[MResetZ(qs[0]), MResetZ(qs[1])]``. + """ + noise = NoiseConfig() + noise.y.loss = 1.0 + if on_loss is not None: + setattr(getattr(noise, attr), "on_loss", on_loss) + source = ( + f"{{use qs = Qubit[2]; {prep} Y(qs[{lose}]); {gate}; {post} " + f"[MResetZ(qs[0]), MResetZ(qs[1])]}}" + ) + return compile_and_run(source, shots=1, seed=SEED, noise=noise, sim_type=sim_type)[ + 0 + ] + + +# TEST 0: C*, R**, and SWAP default loss policies +def test_on_loss_defaults(): + noise = NoiseConfig() + assert noise.cx.on_loss == LossPolicy.SKIP + assert noise.cy.on_loss == LossPolicy.SKIP + assert noise.cz.on_loss == LossPolicy.SKIP + assert noise.rxx.on_loss == LossPolicy.DEGRADE + assert noise.ryy.on_loss == LossPolicy.DEGRADE + assert noise.rzz.on_loss == LossPolicy.DEGRADE + assert noise.swap.on_loss == LossPolicy.APPLY_ANYWAY + + +# TEST 1: SKIP never applies the unitary. C*/R** leave the survivor (|1>) +# untouched ("-1"); SWAP also skips the unitary but still exchanges the loss +# flag, so the survivor's |1> ends up on the now-lost qubit ("0-"). +@pytest.mark.parametrize("sim_type", SIM_TYPES) +@pytest.mark.parametrize( + "attr,gate,expected", + [(*elt, "-1") for elt in ALL_GATES], + ids=ALL_IDS, +) +def test_on_loss_skip_does_not_apply_unitary(attr, gate, expected, sim_type): + res = run_loss_policy_scenario( + gate, sim_type, attr=attr, on_loss=LossPolicy.SKIP, prep="X(qs[1]);" + ) + assert res == expected + + +# TEST 2: PROPAGATE loses the surviving operand too. Here the first operand +# (qs[0]) is the one originally lost, and the loss propagates to qs[1], so both +# operands measure as loss. +@pytest.mark.parametrize("sim_type", SIM_TYPES) +@pytest.mark.parametrize("attr,gate", ALL_GATES, ids=ALL_IDS) +def test_on_loss_propagate_lose_first(attr, gate, sim_type): + res = run_loss_policy_scenario( + gate, sim_type, attr=attr, on_loss=LossPolicy.PROPAGATE, lose=0 + ) + assert res == "--" + + +# TEST 3: same as TEST 2 but the second operand (qs[1]) is the one originally +# lost; the loss still propagates to the survivor, so both measure as loss. +@pytest.mark.parametrize("sim_type", SIM_TYPES) +@pytest.mark.parametrize("attr,gate", ALL_GATES, ids=ALL_IDS) +def test_on_loss_propagate_lose_second(attr, gate, sim_type): + res = run_loss_policy_scenario( + gate, sim_type, attr=attr, on_loss=LossPolicy.PROPAGATE, lose=1 + ) + assert res == "--" + + +# TEST 4: DEGRADE has no single-qubit reduction for controlled gates or SWAP, +# so it falls back to SKIP. C* leave the survivor (|1>) untouched ("-1"); SWAP +# skips the unitary but exchanges the loss flag ("0-"). +@pytest.mark.parametrize("sim_type", SIM_TYPES) +@pytest.mark.parametrize( + "attr,gate,expected", + [(*elt, "-1") for elt in CONTROLLED_GATES + [SWAP_GATE]], + ids=CONTROLLED_IDS + [SWAP_ID], +) +def test_on_loss_degrade_behaves_like_skip(attr, gate, expected, sim_type): + res = run_loss_policy_scenario( + gate, sim_type, attr=attr, on_loss=LossPolicy.DEGRADE, prep="X(qs[1]);" + ) + assert res == expected + + +# TEST 5: DEGRADE reduces a two-qubit rotation to its single-qubit version on +# the survivor (Rxx->Rx, Ryy->Ry, Rzz->Rz). With theta = PI the reduced +# rotation flips the survivor's measured bit; Rz only adds phase, so Rzz is +# measured in the X basis (H before and after). +@pytest.mark.parametrize("sim_type", SIM_TYPES) +@pytest.mark.parametrize( + "attr,gate,prep,post", ROTATION_DEGRADE_RECIPES, ids=ROTATION_IDS +) +def test_on_loss_degrade_reduces_rotation(attr, gate, prep, post, sim_type): + res = run_loss_policy_scenario( + gate, sim_type, attr=attr, on_loss=LossPolicy.DEGRADE, prep=prep, post=post + ) + assert res == "-1" + + +# TEST 6a: RESIDUAL_S_DAGGER skips the unitary but applies an S-dagger to the +# survivor in place. The survivor is prepared in |+i> = S H |0>; S-dagger maps it +# back to |+>, and the final H rotates it to |0>, so the survivor (still at +# qs[1] for C*/R**) measures as 0 ("-0"). +@pytest.mark.parametrize("sim_type", SIM_TYPES) +@pytest.mark.parametrize( + "attr,gate", + CONTROLLED_GATES + ROTATION_GATES, + ids=CONTROLLED_IDS + ROTATION_IDS, +) +def test_on_loss_residual_s_dagger_applies_s_adjoint(attr, gate, sim_type): + res = run_loss_policy_scenario( + gate, + sim_type, + attr=attr, + on_loss=LossPolicy.RESIDUAL_S_DAGGER, + prep="H(qs[1]); S(qs[1]);", + post="H(qs[1]);", + ) + assert res == "-0" + + +# TEST 6b: SWAP under RESIDUAL_S_DAGGER applies the full swap, then the S-dagger +# to the survivor, then exchanges the loss flag. Unlike C*/R**, the swap +# relocates the survivor to qs[0] (the originally-lost operand) and moves the +# loss flag to qs[1]. The survivor |+i> = S H |0> maps under S-dagger to |+>, so +# the post-processing H is applied to qs[0] (the survivor's new location) to +# rotate it to |0>: qs[0] measures 0 and the now-lost qs[1] measures as loss +# ("0-"). Measuring the survivor in place (its X-basis phase) is what verifies +# the S-dagger was actually applied. +@pytest.mark.parametrize("sim_type", SIM_TYPES) +def test_on_loss_swap_residual_s_dagger_applies_s_adjoint(sim_type): + res = run_loss_policy_scenario( + SWAP_GATE[1], + sim_type, + attr=SWAP_GATE[0], + on_loss=LossPolicy.RESIDUAL_S_DAGGER, + prep="H(qs[1]); S(qs[1]);", + post="H(qs[0]);", + ) + assert res == "0-" + + +# TEST 7: SWAP always exchanges the loss flag between its operands, whatever the +# policy. The survivor (|1>) becomes the lost qubit and the originally-lost +# qubit returns as a reset |0> ("0-"); under PROPAGATE both end up lost ("--"). +@pytest.mark.parametrize("sim_type", SIM_TYPES) +@pytest.mark.parametrize( + "on_loss,expected", + [ + (LossPolicy.SKIP, "-1"), + (LossPolicy.PROPAGATE, "--"), + (LossPolicy.DEGRADE, "-1"), + (LossPolicy.RESIDUAL_S_DAGGER, "1-"), + (LossPolicy.APPLY_ANYWAY, "1-"), + ], + ids=["skip", "propagate", "residual_s_dagger", "degrade", "apply_anyway"], +) +def test_on_loss_swap_swaps_loss_flag(on_loss, expected, sim_type): + res = run_loss_policy_scenario( + SWAP_GATE[1], sim_type, attr="swap", on_loss=on_loss, prep="X(qs[1]);" + ) + assert res == expected + + # =========================================================================== # Correlated loss tests ('L' in a noise string) # =========================================================================== diff --git a/source/qdk_package/tests/test_sparse_simulator.py b/source/qdk_package/tests/test_sparse_simulator.py index b2e36c4a29..0105b392bf 100644 --- a/source/qdk_package/tests/test_sparse_simulator.py +++ b/source/qdk_package/tests/test_sparse_simulator.py @@ -3,7 +3,7 @@ from collections import Counter from pathlib import Path -from typing import Sequence, cast +from typing import Dict, Sequence, cast import math import random @@ -43,6 +43,52 @@ def result_array_to_string(results: Sequence[Result]) -> str: return "".join(chars) +def format_expectation(actual: Dict[str, float], expect: Dict[str, float]): + return f"Expected distribution:\n {expect}\n\nActual distribution:\n {actual}" + + +def assert_err(msg: str, actual: Dict[str, float], expect: Dict[str, float]): + return msg + "\n\n" + format_expectation(actual, expect) + + +def assert_distributions_eq( + actual: Dict[str, float], expect: Dict[str, float], tolerance: float +): + # Prune values that are smaller than the tolerance. + actual = {key: val for key, val in actual.items() if val > tolerance} + expect = {key: val for key, val in expect.items() if val > tolerance} + + for key in actual: + assert key in expect, assert_err( + f"Unexpected measurement string: '{key}'.", actual, expect + ) + + for key in expect: + assert key in actual, assert_err( + f"Missing measurement string: '{key}'", actual, expect + ) + + tolerance_percent = int(tolerance * 100) + for key in actual: + assert abs(actual[key] - expect[key]) < tolerance, assert_err( + f"Probability for {key} outside {tolerance_percent}% tolerance.", + actual, + expect, + ) + + +def expect_distribution( + results, + expected: Dict[str, float], + *, + tolerance: float = 0.01, +): + histogram = Counter(results) + total = sum(histogram.values()) + actual = {key: val / total for key, val in histogram.items()} + assert_distributions_eq(actual, expected, tolerance) + + def test_sparse_no_noise(): """Simple test that sparse simulator works without noise.""" qsharp.init(target_profile=TargetProfile.Base) @@ -54,49 +100,60 @@ def test_sparse_no_noise(): assert output == [[Result.Zero] * 16], "Expected result of 0s with pi/2 angles." +QSHARP_OP_25_QUBITS = """ +operation Test() : Result[] { + use qs = Qubit[25]; X(qs[0]); CZ(qs[23], qs[24]); MResetEachZ(qs) +}""" + + def test_sparse_bitflip_noise(): - """Bitflip noise for sparse simulator.""" + """Bitflip noise for Clifford simulator.""" qsharp.init(target_profile=TargetProfile.Base) - qsharp.eval(read_file_relative("CliffordIsing.qs")) + qsharp.eval(QSHARP_OP_25_QUBITS) - p_noise = 0.005 + p_noise = 0.2 noise = NoiseConfig() - noise.rx.set_bitflip(p_noise) - noise.rzz.set_pauli_noise("XX", p_noise) - noise.mresetz.set_bitflip(p_noise) - - output = run( - "IsingModel2DEvolution(4, 4, PI() / 2.0, PI() / 2.0, 10.0, 10)", - shots=1, - noise=noise, - seed=17, - ) + noise.x.set_bitflip(p_noise) + noise.cz.set_pauli_noise("XX", p_noise) + + output = qsharp.run("Test()", shots=1000, noise=noise, seed=17) result = [result_array_to_string(cast(Sequence[Result], x)) for x in output] - print(result) - # Reasonable results obtained from manual run - assert result == ["1000110000000000"] + expect_distribution( + result, + { + "1000000000000000000000000": (1 - p_noise) ** 2, # No noise + "0000000000000000000000000": p_noise * (1 - p_noise), # X bitflip + "1000000000000000000000011": (1 - p_noise) * p_noise, # CZ bitflip + "0000000000000000000000011": p_noise**2, # X & CZ bitflip + }, + tolerance=0.02, + ) def test_sparse_mixed_noise(): qsharp.init(target_profile=TargetProfile.Base) - qsharp.eval(read_file_relative("CliffordIsing.qs")) + qsharp.eval(QSHARP_OP_25_QUBITS) + p_noise = 0.2 noise = NoiseConfig() - noise.rz.set_bitflip(0.008) - noise.rz.loss = 0.005 - noise.rzz.set_depolarizing(0.008) - noise.rzz.loss = 0.005 - - output = run( - "IsingModel2DEvolution(4, 4, PI() / 2.0, PI() / 2.0, 4.0, 4)", - shots=1, - noise=noise, - seed=23, - ) + noise.x.set_bitflip(p_noise) + noise.cz.XI = p_noise + noise.cz.IL = p_noise + + output = qsharp.run("Test()", shots=1000, noise=noise, seed=17) result = [result_array_to_string(cast(Sequence[Result], x)) for x in output] - print(result) - # Reasonable results obtained from manual run - assert result == ["0000000000-00011"] + expect_distribution( + result, + { + "1000000000000000000000000": (1 - p_noise) * (1 - 2 * p_noise), # No noise + "0000000000000000000000000": p_noise * (1 - 2 * p_noise), # X bitflip + "1000000000000000000000010": (1 - p_noise) * p_noise, # CZ bitflip + "100000000000000000000000-": (1 - p_noise) * p_noise, # CZ loss + "0000000000000000000000010": p_noise**2, # X bitflip + CZ bitflip + "000000000000000000000000-": p_noise**2, # X bitflip + CZ loss + }, + tolerance=0.02, + ) def test_sparse_isolated_loss(): diff --git a/source/simulators/src/cpu_full_state_simulator.rs b/source/simulators/src/cpu_full_state_simulator.rs index 1ca4f5260c..72f8f8917e 100644 --- a/source/simulators/src/cpu_full_state_simulator.rs +++ b/source/simulators/src/cpu_full_state_simulator.rs @@ -3,7 +3,7 @@ use crate::{ MeasurementResult, QubitID, Simulator, - noise_config::{CumulativeNoiseConfig, Fault, FaultTerm, IntrinsicID}, + noise_config::{CumulativeNoiseConfig, Fault, FaultTerm, IntrinsicID, LossPolicy}, }; use core::f64; use nalgebra::Complex; @@ -495,6 +495,15 @@ impl NoisySimulator { } } + /// Applies an `S` adjoint to the given target + /// Used by the [`LossPolicy::ResidualSDagger`] behavior. + fn residual_s_dagger(&mut self, target: QubitID) { + self.apply_idle_noise(target); + self.state + .apply_operation(&S_ADJ, &[target]) + .expect("apply_operation should succeed"); + } + /// Records a z-measurement on the given `target`. fn record_mz(&mut self, target: QubitID, result_id: QubitID) { let measurement = self.mz_impl(target); @@ -550,6 +559,13 @@ impl NoisySimulator { MeasurementResult::Zero } } + + fn loss_impl(&mut self, target: QubitID) { + if !self.loss[target] { + self.mresetz_impl(target); + self.loss[target] = true; + } + } } /// Design decision: Why is this a macro? @@ -760,36 +776,84 @@ impl Simulator for NoisySimulator { } fn cx(&mut self, control: QubitID, target: QubitID) { - if !self.loss[control] && !self.loss[target] { - self.apply_idle_noise(control); - self.apply_idle_noise(target); - self.state - .apply_operation(&CX, &[control, target]) - .expect("apply_operation should succeed"); + match (self.loss[control], self.loss[target]) { + (true, true) => (), + (true, false) | (false, true) => { + let remaining_qubit = if self.loss[control] { target } else { control }; + self.apply_idle_noise(remaining_qubit); + match self.noise_config.cx.on_loss { + LossPolicy::Skip | LossPolicy::Degrade => (), + LossPolicy::Propagate => self.loss_impl(remaining_qubit), + LossPolicy::ResidualSDagger => self.residual_s_dagger(remaining_qubit), + LossPolicy::ApplyAnyway => self + .state + .apply_operation(&CX, &[control, target]) + .expect("apply_operation should succeed"), + } + } + (false, false) => { + self.apply_idle_noise(control); + self.apply_idle_noise(target); + self.state + .apply_operation(&CX, &[control, target]) + .expect("apply_operation should succeed"); + } } // We still apply operation faults to non-lost qubits. apply_noise!(self, cx, &[control, target]); } fn cy(&mut self, control: QubitID, target: QubitID) { - if !self.loss[control] && !self.loss[target] { - self.apply_idle_noise(control); - self.apply_idle_noise(target); - self.state - .apply_operation(&CY, &[control, target]) - .expect("apply_operation should succeed"); + match (self.loss[control], self.loss[target]) { + (true, true) => (), + (true, false) | (false, true) => { + let remaining_qubit = if self.loss[control] { target } else { control }; + self.apply_idle_noise(remaining_qubit); + match self.noise_config.cy.on_loss { + LossPolicy::Skip | LossPolicy::Degrade => (), + LossPolicy::Propagate => self.loss_impl(remaining_qubit), + LossPolicy::ResidualSDagger => self.residual_s_dagger(remaining_qubit), + LossPolicy::ApplyAnyway => self + .state + .apply_operation(&CY, &[control, target]) + .expect("apply_operation should succeed"), + } + } + (false, false) => { + self.apply_idle_noise(control); + self.apply_idle_noise(target); + self.state + .apply_operation(&CY, &[control, target]) + .expect("apply_operation should succeed"); + } } // We still apply operation faults to non-lost qubits. apply_noise!(self, cy, &[control, target]); } fn cz(&mut self, control: QubitID, target: QubitID) { - if !self.loss[control] && !self.loss[target] { - self.apply_idle_noise(control); - self.apply_idle_noise(target); - self.state - .apply_operation(&CZ, &[control, target]) - .expect("apply_operation should succeed"); + match (self.loss[control], self.loss[target]) { + (true, true) => (), + (true, false) | (false, true) => { + let remaining_qubit = if self.loss[control] { target } else { control }; + self.apply_idle_noise(remaining_qubit); + match self.noise_config.cz.on_loss { + LossPolicy::Skip | LossPolicy::Degrade => (), + LossPolicy::Propagate => self.loss_impl(remaining_qubit), + LossPolicy::ResidualSDagger => self.residual_s_dagger(remaining_qubit), + LossPolicy::ApplyAnyway => self + .state + .apply_operation(&CZ, &[control, target]) + .expect("apply_operation should succeed"), + } + } + (false, false) => { + self.apply_idle_noise(control); + self.apply_idle_noise(target); + self.state + .apply_operation(&CZ, &[control, target]) + .expect("apply_operation should succeed"); + } } // We still apply operation faults to non-lost qubits. apply_noise!(self, cz, &[control, target]); @@ -798,65 +862,120 @@ impl Simulator for NoisySimulator { fn rxx(&mut self, angle: f64, q1: QubitID, q2: QubitID) { match (self.loss[q1], self.loss[q2]) { (true, true) => (), - (true, false) => self.rx(angle, q2), - (false, true) => self.rx(angle, q1), + (true, false) | (false, true) => { + let remaining_qubit = if self.loss[q1] { q2 } else { q1 }; + self.apply_idle_noise(remaining_qubit); + match self.noise_config.rxx.on_loss { + LossPolicy::Skip => (), + LossPolicy::Degrade => return self.rx(angle, remaining_qubit), + LossPolicy::Propagate => self.loss_impl(remaining_qubit), + LossPolicy::ResidualSDagger => self.residual_s_dagger(remaining_qubit), + LossPolicy::ApplyAnyway => self + .state + .apply_operation(&rxx(angle), &[q1, q2]) + .expect("apply_operation should succeed"), + } + } (false, false) => { self.apply_idle_noise(q1); self.apply_idle_noise(q2); self.state .apply_operation(&rxx(angle), &[q1, q2]) .expect("apply_operation should succeed"); - apply_noise!(self, rxx, &[q1, q2]); } } + apply_noise!(self, rxx, &[q1, q2]); } fn ryy(&mut self, angle: f64, q1: QubitID, q2: QubitID) { match (self.loss[q1], self.loss[q2]) { (true, true) => (), - (true, false) => self.ry(angle, q2), - (false, true) => self.ry(angle, q1), + (true, false) | (false, true) => { + let remaining_qubit = if self.loss[q1] { q2 } else { q1 }; + self.apply_idle_noise(remaining_qubit); + match self.noise_config.ryy.on_loss { + LossPolicy::Skip => (), + LossPolicy::Degrade => return self.ry(angle, remaining_qubit), + LossPolicy::Propagate => self.loss_impl(remaining_qubit), + LossPolicy::ResidualSDagger => self.residual_s_dagger(remaining_qubit), + LossPolicy::ApplyAnyway => self + .state + .apply_operation(&ryy(angle), &[q1, q2]) + .expect("apply_operation should succeed"), + } + } (false, false) => { self.apply_idle_noise(q1); self.apply_idle_noise(q2); self.state .apply_operation(&ryy(angle), &[q1, q2]) .expect("apply_operation should succeed"); - apply_noise!(self, ryy, &[q1, q2]); } } + apply_noise!(self, ryy, &[q1, q2]); } fn rzz(&mut self, angle: f64, q1: QubitID, q2: QubitID) { match (self.loss[q1], self.loss[q2]) { (true, true) => (), - (true, false) => self.rz(angle, q2), - (false, true) => self.rz(angle, q1), + (true, false) | (false, true) => { + let remaining_qubit = if self.loss[q1] { q2 } else { q1 }; + self.apply_idle_noise(remaining_qubit); + match self.noise_config.rzz.on_loss { + LossPolicy::Skip => (), + LossPolicy::Degrade => return self.rz(angle, remaining_qubit), + LossPolicy::Propagate => self.loss_impl(remaining_qubit), + LossPolicy::ResidualSDagger => self.residual_s_dagger(remaining_qubit), + LossPolicy::ApplyAnyway => self + .state + .apply_operation(&rzz(angle), &[q1, q2]) + .expect("apply_operation should succeed"), + } + } (false, false) => { self.apply_idle_noise(q1); self.apply_idle_noise(q2); self.state .apply_operation(&rzz(angle), &[q1, q2]) .expect("apply_operation should succeed"); - apply_noise!(self, rzz, &[q1, q2]); } } + apply_noise!(self, rzz, &[q1, q2]); } fn swap(&mut self, q1: QubitID, q2: QubitID) { + // There are three kinds of swaps: + // 1. A logical swap, also called a relabel. + // 2. A swap by physically exchanging the location of the qubits. + // 3. An exchange of information by doing three CX. + // + // This method is concerned with the kinds (1) and (2), since (3) + // gets decomposed into other instructions before making it to the simulator. + // In both (1) and (2), the loss state of the qubits gets exchanged. + match (self.loss[q1], self.loss[q2]) { (true, true) => (), - (true, false) => { - self.apply_idle_noise(q2); - self.state - .apply_operation(&SWAP, &[q1, q2]) - .expect("apply_operation should succeed"); - } - (false, true) => { - self.apply_idle_noise(q1); - self.state - .apply_operation(&SWAP, &[q1, q2]) - .expect("apply_operation should succeed"); + (true, false) | (false, true) => { + let lost_qubit = if self.loss[q1] { q1 } else { q2 }; + let remaining_qubit = if self.loss[q1] { q2 } else { q1 }; + self.apply_idle_noise(remaining_qubit); + match self.noise_config.swap.on_loss { + LossPolicy::Skip | LossPolicy::Degrade => (), + LossPolicy::Propagate => self.loss_impl(remaining_qubit), + LossPolicy::ResidualSDagger => { + self.state + .apply_operation(&SWAP, &[q1, q2]) + .expect("apply_operation should succeed"); + self.residual_s_dagger(lost_qubit); + self.loss.swap(q1, q2); + } + LossPolicy::ApplyAnyway => { + self.state + .apply_operation(&SWAP, &[q1, q2]) + .expect("apply_operation should succeed"); + self.loss.swap(q1, q2); + } + } } (false, false) => { self.apply_idle_noise(q1); @@ -864,17 +983,9 @@ impl Simulator for NoisySimulator { self.state .apply_operation(&SWAP, &[q1, q2]) .expect("apply_operation should succeed"); + self.loss.swap(q1, q2); } } - // There are three kinds of swaps: - // 1. A logical swap, also called a relabel. - // 2. A swap by physically exchanging the location of the qubits. - // 3. An exchange of information by doing three CX. - // - // This method is concerned with the kinds (1) and (2), since (3) - // gets decomposed into other instructions before making it to the simulator. - // In both (1) and (2), the loss state of the qubits gets exchanged. - self.loss.swap(q1, q2); // Is up to the user if swap is a virtual operation or not. // If they don't specify noise/loss probability for swap, then it is virtual. diff --git a/source/simulators/src/gpu_full_state_simulator/common.wgsl b/source/simulators/src/gpu_full_state_simulator/common.wgsl index fb3eccacbf..cbd62f0d75 100644 --- a/source/simulators/src/gpu_full_state_simulator/common.wgsl +++ b/source/simulators/src/gpu_full_state_simulator/common.wgsl @@ -331,6 +331,265 @@ fn get_pauli_noise_idx(op_idx: u32) -> u32 { return 0u; } +// From the starting index given, return the next index if loss noise, else 0 +fn get_loss_idx(op_idx: u32) -> u32 { + if (arrayLength(&ops) > (op_idx + 1)) { + let op = &ops[op_idx + 1]; + if (op.id == OPID_LOSS_NOISE) { + return op_idx + 1u; + } + } + return 0u; +} + +// Loss policy values. These are stamped onto a gate op's `q3` field by the host +// (see `LossPolicy::as_u32` on the Rust side) and tell the shader how to handle +// the gate when one of its operands is lost. `0` means "no policy stamped", +// which the shader treats the same as SKIP. +const LOSS_POLICY_SKIP = 0u; +const LOSS_POLICY_PROPAGATE = 1u; +const LOSS_POLICY_DEGRADE = 2u; +const LOSS_POLICY_RESIDUAL_S_DAGGER = 3u; +const LOSS_POLICY_APPLY_ANYWAY = 4u; + +// Returns true if the gate at `op_idx` touches at least one lost qubit. +// `q1`/`q2` are the (resolved) operands of the gate. +fn gate_has_lost_operand(shot_idx: u32, op_idx: u32, q1: u32, q2: u32) -> bool { + let shot = &shots[shot_idx]; + let op = &ops[op_idx]; + if (shot.qubit_state[q1].heat == -1.0) { + return true; + } + let is_2q = !is_1q_op(op.id); + return is_2q && (shot.qubit_state[q2].heat == -1.0); +} + +// Builds a 4x4 (in shot.unitary) that applies the 1-qubit matrix `m` (given as +// m00,m01,m10,m11) to `target_is_q2 ? q2 : q1` and identity to the other qubit +// of the pair. The lost qubit is in the |0> state, so the identity factor keeps +// it there. The 2-qubit basis is |q1 q2>, so the row/col index is +// (2 * q1_bit + q2_bit). +fn set_1q_on_pair_unitary(shot_idx: u32, target_is_q2: bool, + m00: vec2f, m01: vec2f, m10: vec2f, m11: vec2f) { + let shot = &shots[shot_idx]; + // Zero the whole 4x4 first. + for (var i = 0u; i < 16u; i++) { + shot.unitary[i] = vec2f(0.0, 0.0); + } + if target_is_q2 { + // Acts on q2 (low bit): block-diagonal diag(M, M). + // Top-left block (q1 = 0): + shot.unitary[0] = m00; shot.unitary[1] = m01; + shot.unitary[4] = m10; shot.unitary[5] = m11; + // Bottom-right block (q1 = 1): + shot.unitary[10] = m00; shot.unitary[11] = m01; + shot.unitary[14] = m10; shot.unitary[15] = m11; + } else { + // Acts on q1 (high bit): M (x) I. + shot.unitary[0] = m00; shot.unitary[2] = m01; + shot.unitary[8] = m10; shot.unitary[10] = m11; + shot.unitary[5] = m00; shot.unitary[7] = m01; + shot.unitary[13] = m10; shot.unitary[15] = m11; + } +} + +// Multiplies one row of the 4x4 pair unitary (in shot.unitary) by -i, in place. +// Folding a diag(1, -i) = S-dagger factor on one qubit into a 2-qubit matrix +// scales the rows whose target-qubit bit is 1 by -i. For a complex entry +// (x + y i), (x + y i) * -i = y - x i. +fn scale_pair_unitary_row_by_neg_i(shot_idx: u32, row: u32) { + let shot = &shots[shot_idx]; + for (var c = 0u; c < 4u; c++) { + let e = shot.unitary[row * 4u + c]; + shot.unitary[row * 4u + c] = vec2f(e.y, -e.x); + } +} + +// Sets up the shot to execute a 2-qubit shot-buffer op on the gate's operands. +fn finish_2q_shot_buffer(shot_idx: u32, op_idx: u32, q1: u32, q2: u32) { + let shot = &shots[shot_idx]; + shot.op_idx = op_idx; + shot.op_type = OPID_SHOT_BUFF_2Q; + shot.qubits_updated_last_op_mask = (1u << q1) | (1u << q2); +} + +// Loses a single surviving `qubit` for the PROPAGATE policy: samples a +// measurement outcome, collapses the qubit to that outcome and resets it to +// |0>, and marks it lost (heat = -1.0). The collapse is expressed as a 2-qubit +// tensor on the gate's operands (reset on `qubit`, identity on the lost +// partner, which is already in |0>), reusing the standard shot-buffer execute +// path. `qubit` must be one of the gate's two operands `q1`/`q2`. +fn propagate_loss_to_qubit(shot_idx: u32, op_idx: u32, q1: u32, q2: u32, qubit: u32) { + let shot = &shots[shot_idx]; + + let result = select(1u, 0u, shot.rand_measure < shot.qubit_state[qubit].zero_probability); + + // Reset instrument (project + move |1> into |0> slot), same as MResetZ: + // result==0: [[1,0],[0,0]] + // result==1: [[0,1],[0,0]] + let m00 = select(vec2f(1.0, 0.0), vec2f(0.0, 0.0), result == 1u); + let m01 = select(vec2f(0.0, 0.0), vec2f(1.0, 0.0), result == 1u); + let m10 = vec2f(0.0, 0.0); + let m11 = vec2f(0.0, 0.0); + + let target_is_q2 = (qubit == q2); + set_1q_on_pair_unitary(shot_idx, target_is_q2, m00, m01, m10, m11); + + // Renormalize by the measured branch probability. + shot.renormalize = select( + 1.0 / sqrt(shot.qubit_state[qubit].zero_probability), + 1.0 / sqrt(shot.qubit_state[qubit].one_probability), + result == 1u); + + // Mark the qubit lost and clear its definite-state bits so the probability + // pass recomputes it. + shot.qubit_state[qubit].heat = -1.0; + shot.qubit_is_0_mask = shot.qubit_is_0_mask & ~(1u << qubit); + shot.qubit_is_1_mask = shot.qubit_is_1_mask & ~(1u << qubit); + + finish_2q_shot_buffer(shot_idx, op_idx, q1, q2); +} + +// Handles a gate whose operand(s) include at least one lost qubit, according to +// the loss policy stamped on the op's `q3` field. `q1`/`q2` are the (resolved) +// operands. Returns true if the gate was fully handled here (the caller should +// return), or false if normal processing should continue (only for +// APPLY_ANYWAY, which runs the gate as usual). +fn handle_lost_operand_policy(shot_idx: u32, op_idx: u32, q1: u32, q2: u32) -> bool { + let shot = &shots[shot_idx]; + let op = &ops[op_idx]; + let is_1q = is_1q_op(op.id); + let is_2q = !is_1q; + let policy = op.policy; + + // Loss policies only make sense for multi-qubit gates. + // If this is a single-qubit gate, skip it entirely. + if (is_1q) { + shot.op_type = OPID_ID; + shot.op_idx = op_idx; + return true; + } + + let q1_lost = shot.qubit_state[q1].heat == -1.0; + let q2_lost = is_2q && (shot.qubit_state[q2].heat == -1.0); + let has_survivor = is_2q && !(q1_lost && q2_lost); + // The surviving operand (only meaningful when has_survivor is true). + let survivor = select(q1, q2, q1_lost); + let survivor_is_q2 = q1_lost; + + // SWAP is special: it physically relocates the two qubits, so their loss + // state is always exchanged regardless of the policy (the policy only + // governs whether the unitary runs). Handle it explicitly here. + if (op.id == OPID_SWAP) { + switch policy { + case LOSS_POLICY_PROPAGATE { + propagate_loss_to_qubit(shot_idx, op_idx, q1, q2, survivor); + return true; + } + case LOSS_POLICY_RESIDUAL_S_DAGGER { + // Match the CPU/stabilizer SWAP + residual S-dagger semantics: + // 1. Apply the full SWAP (shot.unitary already holds it). + // 2. Apply S-dagger = diag(1, -i) to the (originally) lost + // operand's position, which after the SWAP holds the + // survivor's amplitudes. + // 3. Exchange the per-qubit loss flag (heat) of the operands. + + // Fold the S-dagger into the SWAP matrix by scaling, by -i, the + // two rows of the |q1 q2> pair matrix whose lost-qubit bit is 1. + // q1 is the high bit (rows 2, 3); q2 is the low bit (rows 1, 3). + let lost_row = select(1u, 2u, q1_lost); + scale_pair_unitary_row_by_neg_i(shot_idx, lost_row); + scale_pair_unitary_row_by_neg_i(shot_idx, 3u); + // Exchange the per-qubit loss flag (heat) of the two operands. + let heat1 = shot.qubit_state[q1].heat; + shot.qubit_state[q1].heat = shot.qubit_state[q2].heat; + shot.qubit_state[q2].heat = heat1; + // The 2-qubit execute path skips amplitudes for qubits known to be + // in a definite state, which would skip the amplitudes SWAP needs to move. + // Clear those bits for both operands so the swap is actually applied. + shot.qubit_is_0_mask = shot.qubit_is_0_mask & ~((1u << q1) | (1u << q2)); + shot.qubit_is_1_mask = shot.qubit_is_1_mask & ~((1u << q1) | (1u << q2)); + // shot.unitary now holds (S-dagger on lost) * SWAP. + finish_2q_shot_buffer(shot_idx, op_idx, q1, q2); + return true; + } + case LOSS_POLICY_APPLY_ANYWAY { + // Exchange the per-qubit loss flag (heat) of the two operands. + let heat1 = shot.qubit_state[q1].heat; + shot.qubit_state[q1].heat = shot.qubit_state[q2].heat; + shot.qubit_state[q2].heat = heat1; + // The 2-qubit execute path skips amplitudes for qubits known to be + // in a definite state, which would skip the amplitudes SWAP needs to move. + // Clear those bits for both operands so the swap is actually applied. + shot.qubit_is_0_mask = shot.qubit_is_0_mask & ~((1u << q1) | (1u << q2)); + shot.qubit_is_1_mask = shot.qubit_is_1_mask & ~((1u << q1) | (1u << q2)); + // shot.unitary already holds the SWAP matrix (set by the caller). + finish_2q_shot_buffer(shot_idx, op_idx, q1, q2); + return true; + } + default { + shot.op_type = OPID_ID; + shot.op_idx = op_idx; + return true; + } + } + } + + // APPLY_ANYWAY: run the gate as if nothing was lost. + if (policy == LOSS_POLICY_APPLY_ANYWAY) { + return false; + } + + if (policy == LOSS_POLICY_PROPAGATE && has_survivor) { + propagate_loss_to_qubit(shot_idx, op_idx, q1, q2, survivor); + return true; + } + + if (policy == LOSS_POLICY_RESIDUAL_S_DAGGER && has_survivor) { + // Apply S-dagger = diag(1, -i) to the surviving operand. + set_1q_on_pair_unitary(shot_idx, survivor_is_q2, + vec2f(1.0, 0.0), vec2f(0.0, 0.0), + vec2f(0.0, 0.0), vec2f(0.0, -1.0)); + finish_2q_shot_buffer(shot_idx, op_idx, q1, q2); + return true; + } + + if (policy == LOSS_POLICY_DEGRADE && has_survivor && + (op.id == OPID_RXX || op.id == OPID_RYY || op.id == OPID_RZZ)) { + // Degrade the two-qubit rotation to its single-qubit version on the + // survivor. The op's unitary[0] holds cos(θ/2) for Rxx/Ryy; we recover + // the angle to build the 1-qubit rotation matrix. + let cos_half = op.unitary[0].x; + if (op.id == OPID_RXX) { + // Rx(θ) = [[c, -i s], [-i s, c]], where s = sin(θ/2). + let s = op.unitary[3].y * -1.0; // unitary[3] = (0, -sin(θ/2)) + set_1q_on_pair_unitary(shot_idx, survivor_is_q2, + vec2f(cos_half, 0.0), vec2f(0.0, -s), + vec2f(0.0, -s), vec2f(cos_half, 0.0)); + } else if (op.id == OPID_RYY) { + // Ry(θ) = [[c, -s], [s, c]], where s = sin(θ/2). + let s = op.unitary[3].y; // unitary[3] = (0, sin(θ/2)) for Ryy + set_1q_on_pair_unitary(shot_idx, survivor_is_q2, + vec2f(cos_half, 0.0), vec2f(-s, 0.0), + vec2f(s, 0.0), vec2f(cos_half, 0.0)); + } else { + // Rzz -> Rz(θ). The GPU Rz convention is [[1, 0], [0, e^{iθ}]], + // and unitary[5] = e^{iθ} holds the full-angle phase. + let phase = op.unitary[5]; + set_1q_on_pair_unitary(shot_idx, survivor_is_q2, + vec2f(1.0, 0.0), vec2f(0.0, 0.0), + vec2f(0.0, 0.0), phase); + } + finish_2q_shot_buffer(shot_idx, op_idx, q1, q2); + return true; + } + + // SKIP (and any policy with no applicable survivor, e.g. DEGRADE on a + // controlled gate, or a single-qubit gate): skip the gate entirely. + shot.op_type = OPID_ID; + shot.op_idx = op_idx; + return true; +} fn apply_1q_pauli_noise(shot_idx: u32, op_idx: u32, noise_idx: u32) { // NOTE: Assumes that whatever prepared the program ensured that noise_op.q1 matches op.q1 and diff --git a/source/simulators/src/gpu_full_state_simulator/gpu_context.rs b/source/simulators/src/gpu_full_state_simulator/gpu_context.rs index 56df1b87e7..c87500504f 100644 --- a/source/simulators/src/gpu_full_state_simulator/gpu_context.rs +++ b/source/simulators/src/gpu_full_state_simulator/gpu_context.rs @@ -10,7 +10,7 @@ use crate::bytecode::AdaptiveProgram; use crate::correlated_noise::NoiseTables; use crate::gpu_resources::GpuResources; use crate::noise_config::NoiseConfig; -use crate::noise_mapping::{expand_correlated_loss_commits, get_noise_ops}; +use crate::noise_mapping::{expand_correlated_loss_commits, get_noise_ops, loss_policy_u32}; use crate::shader_types::{ self, DiagnosticsData, InterpreterState, MAX_ALLOCA_SIZE, MAX_BUFFER_SIZE, MAX_QUBIT_COUNT, MAX_QUBITS_PER_WORKGROUP, MAX_REGISTERS, MAX_SHOTS_PER_BATCH, MIN_QUBIT_COUNT, MIN_REGISTERS, @@ -931,8 +931,14 @@ fn add_noise_config_to_ops(ops: &[Op], noise: &NoiseConfig) -> Vec let mut noisy_ops: Vec = Vec::with_capacity(ops.len() + 1); for op in ops { - let mut add_ops: Vec = vec![*op]; - // If there's a NoiseConfig, and we get noise for this op, append it. + // Stamp the configured loss policy onto the gate op so the shader can + // decide how to handle the gate when one of its operands is lost. + let mut gate_op = *op; + if let Some(policy) = loss_policy_u32(op, noise) { + gate_op.policy = policy; + } + let mut add_ops: Vec = vec![gate_op]; + // If there's a NoiseConfig, and we get noise for this op, append it // The base path dispatches ops linearly, so it needs explicit // loss-commit ops to perform any deferred qubit loss. if let Some(noise_ops) = get_noise_ops(op, noise, true) { @@ -989,7 +995,13 @@ fn add_noise_to_adaptive_ops( let new_idx = noisy_ops.len() as u32; index_map.push(new_idx); - noisy_ops.push(*op); + // Stamp the configured loss policy onto the gate op so the shader can + // decide how to handle the gate when one of its operands is lost. + let mut gate_op = *op; + if let Some(policy) = loss_policy_u32(op, noise) { + gate_op.policy = policy; + } + noisy_ops.push(gate_op); // Append the Pauli/loss sampler op (no loss-commit ops; see above). if let Some(noise_ops) = get_noise_ops(op, noise, false) { diff --git a/source/simulators/src/gpu_full_state_simulator/noise_mapping.rs b/source/simulators/src/gpu_full_state_simulator/noise_mapping.rs index 03ddff440e..8d595cf445 100644 --- a/source/simulators/src/gpu_full_state_simulator/noise_mapping.rs +++ b/source/simulators/src/gpu_full_state_simulator/noise_mapping.rs @@ -60,22 +60,12 @@ fn get_noise_op(op: &Op, noise_table: &NoiseTable) -> Op { noise_op } -/// Builds the noise ops to insert after `op` for the given config, or `None` -/// if the gate is noiseless. -/// -/// `emit_loss_commits` controls whether loss-commit ops are appended after the -/// categorical sampler op. The base (non-adaptive) path dispatches ops linearly -/// and needs an explicit loss-commit op per qubit to perform the deferred -/// measure + reset, so it passes `true`. The adaptive path instead drains -/// `pending_loss_mask` inside the interpreter loop, so it passes `false` to -/// avoid emitting loss-commit ops that would never be dispatched (which would -/// otherwise roughly double the op pool for circuits with loss on every gate). -#[must_use] -pub fn get_noise_ops( +/// Returns the [`NoiseTable`] in `noise_config` that applies to the given op, +/// or `None` if the op has no associated noise table (e.g. a noise op itself). +fn noise_table_for<'a>( op: &Op, - noise_config: &NoiseConfig, - emit_loss_commits: bool, -) -> Option> { + noise_config: &'a NoiseConfig, +) -> Option<&'a NoiseTable> { let noise_table = match op.id { ops::ID => &noise_config.i, ops::X => &noise_config.x, @@ -103,6 +93,38 @@ pub fn get_noise_ops( ops::MRESETZ | ops::RESETZ => &noise_config.mresetz, _ => return None, }; + Some(noise_table) +} + +/// Returns the [`LossPolicy`] configured for the given op's gate, encoded as a +/// `u32` for the GPU shader (see [`LossPolicy::as_u32`]). Returns `None` for +/// ops that have no associated gate noise table. +/// +/// The shader reads this from the gate op's `q3` field to decide how to handle +/// the gate when one of its operands is lost. +#[must_use] +pub fn loss_policy_u32(op: &Op, noise_config: &NoiseConfig) -> Option { + noise_table_for(op, noise_config).map(|table| table.on_loss.as_u32()) +} + +/// Builds the noise ops to insert after `op` for the given config, or `None` +/// if the gate is noiseless. +/// +/// `emit_loss_commits` controls whether loss-commit ops are appended after the +/// categorical sampler op. The base (non-adaptive) path dispatches ops linearly +/// and needs an explicit loss-commit op per qubit to perform the deferred +/// measure + reset, so it passes `true`. The adaptive path instead drains +/// `pending_loss_mask` inside the interpreter loop, so it passes `false` to +/// avoid emitting loss-commit ops that would never be dispatched (which would +/// otherwise roughly double the op pool for circuits with loss on every gate). +#[must_use] +pub fn get_noise_ops( + op: &Op, + noise_config: &NoiseConfig, + emit_loss_commits: bool, +) -> Option> { + let noise_table = noise_table_for(op, noise_config)?; + if noise_table.is_noiseless() { return None; } diff --git a/source/simulators/src/gpu_full_state_simulator/shader_types.rs b/source/simulators/src/gpu_full_state_simulator/shader_types.rs index 97a33cdc4f..74aaf9f528 100644 --- a/source/simulators/src/gpu_full_state_simulator/shader_types.rs +++ b/source/simulators/src/gpu_full_state_simulator/shader_types.rs @@ -298,6 +298,10 @@ pub struct Op { pub q1: u32, pub q2: u32, pub q3: u32, // For ccx + pub policy: u32, + pub pad0: u32, // for 16-byte alignment + pub pad1: u32, // for 16-byte alignment + pub pad2: u32, // for 16-byte alignment pub r00: f32, pub i00: f32, pub r01: f32, @@ -333,7 +337,7 @@ pub struct Op { } // safety check to make sure Op is the correct size with padding at compile time -const _: () = assert!(std::mem::size_of::() == 144); +const _: () = assert!(std::mem::size_of::() == 160); impl Default for Op { fn default() -> Self { @@ -342,6 +346,10 @@ impl Default for Op { q1: 0, q2: 0, q3: 0, + policy: 0, + pad0: 0, + pad1: 0, + pad2: 0, r00: 0.0, i00: 0.0, r01: 0.0, @@ -693,10 +701,14 @@ impl Op { /// in the op's matrix-storage floats. The host and shader agree that flat /// slot `k` maps to WGSL `unitary[k / 2][k % 2]`. pub fn set_noise_prob_slot(&mut self, slot: usize, prob: f32) { - // Op is 4 leading u32 fields followed by 32 matrix floats - // (r00, i00, r01, i01, ...), so matrix flat slot `k` lives at index `4 + k`. - let floats: &mut [f32; 36] = bytemuck::cast_mut(self); - floats[4 + slot] = prob; + // The matrix floats (r00, i00, r01, i01, ...) follow the leading u32 + // header fields. Derive the lengths from the struct layout so this stays + // correct if header fields are added or removed, and so the `cast_mut` + // below never hits a size mismatch. + const OP_FLOATS: usize = std::mem::size_of::() / std::mem::size_of::(); + const MATRIX_OFFSET: usize = std::mem::offset_of!(Op, r00) / std::mem::size_of::(); + let floats: &mut [f32; OP_FLOATS] = bytemuck::cast_mut(self); + floats[MATRIX_OFFSET + slot] = prob; } #[must_use] @@ -928,16 +940,22 @@ impl Op { /// in matrix flat slot `i`). #[must_use] pub fn correlated_noise_qubit(&self, index: u32) -> u32 { - let floats: &[f32; 36] = bytemuck::cast_ref(self); - // The 4 leading u32 fields precede the 32 matrix floats. Qubit ids are - // stored as exact f32 values (range limited to 32), mirroring how the - // shader reads them back as u32. + // The matrix floats (r00, i00, r01, i01, ...) follow the leading u32 + // header fields. Derive the lengths from the struct layout so this stays + // correct if header fields are added or removed, and so the `cast_mut` + // below never hits a size mismatch. + const OP_FLOATS: usize = std::mem::size_of::() / std::mem::size_of::(); + const MATRIX_OFFSET: usize = std::mem::offset_of!(Op, r00) / std::mem::size_of::(); + let floats: &[f32; OP_FLOATS] = bytemuck::cast_ref(self); + + // Qubit ids are stored as exact f32 values (range limited to 32), + // mirroring how the shader reads them back as u32. #[allow( clippy::cast_possible_truncation, clippy::cast_sign_loss, reason = "qubit ids are small non-negative integers stored exactly as f32" )] - let qubit = floats[4 + index as usize] as u32; + let qubit = floats[MATRIX_OFFSET + index as usize] as u32; qubit } diff --git a/source/simulators/src/gpu_full_state_simulator/simulator_adaptive.wgsl b/source/simulators/src/gpu_full_state_simulator/simulator_adaptive.wgsl index 1ac9f8e2d4..7d5e862e03 100644 --- a/source/simulators/src/gpu_full_state_simulator/simulator_adaptive.wgsl +++ b/source/simulators/src/gpu_full_state_simulator/simulator_adaptive.wgsl @@ -80,10 +80,14 @@ struct Op { q1: u32, q2: u32, q3: u32, + policy: u32, + pad0: u32, + pad1: u32, + pad2: u32, // Entries in the unitary are: 00, 01, 02, 03, 10, 11, 12, 13, 20, ..., 32, 33 // 1q matrix elements are stored in: 00, 01, 10, 11 (i.e., indices 0, 1, 4, and 5) unitary: array, -} // Struct size: 4 * 4 + 16 * 8 = 144 bytes (which is aligned to 16 bytes) +} // Struct size: 4 * 8 + 16 * 8 = 160 bytes (which is aligned to 16 bytes) @group(0) @binding(2) var ops: array; @@ -1617,6 +1621,16 @@ fn prepare_op(@builtin(global_invocation_id) globalId: vec3) { shot.op_idx = op_idx; shot.op_type = op.id; + // If any operand is lost, dispatch the gate's configured loss + // policy (stamped on op.policy). For most policies this fully handles + // the op; APPLY_ANYWAY returns false so the gate runs as usual. + if gate_has_lost_operand(shot_idx, op_idx, q1, q2) { + if handle_lost_operand_policy(shot_idx, op_idx, q1, q2) { + shots[shot_idx].interp.status = STATUS_RUNNING; + return; + } + } + // Check for noise ops after this gate in the ops pool let pauli_op_idx = get_pauli_noise_idx(op_idx); diff --git a/source/simulators/src/gpu_full_state_simulator/simulator_base.wgsl b/source/simulators/src/gpu_full_state_simulator/simulator_base.wgsl index e1f966630d..27fc874679 100644 --- a/source/simulators/src/gpu_full_state_simulator/simulator_base.wgsl +++ b/source/simulators/src/gpu_full_state_simulator/simulator_base.wgsl @@ -72,10 +72,14 @@ struct Op { q1: u32, q2: u32, q3: u32, + policy: u32, + pad0: u32, + pad1: u32, + pad2: u32, // Entries in the unitary are: 00, 01, 02, 03, 10, 11, 12, 13, 20, ..., 32, 33 // 1q matrix elements are stored in: 00, 01, 10, 11 (i.e., indices 0, 1, 4, and 5) unitary: array, -} // Struct size: 4 * 4 + 16 * 8 = 144 bytes (which is aligned to 16 bytes) +} // Struct size: 4 * 4 + 16 * 8 = 160 bytes (which is aligned to 16 bytes) @group(0) @binding(2) var ops: array; @@ -310,13 +314,14 @@ fn prepare_op(@builtin(global_invocation_id) globalId: vec3) { return; } - // Before doing further work, if any qubit for the gate is lost, just skip by marking the op as ID - if (shot.qubit_state[op.q1].heat == -1.0) || - (op.id == OPID_CX || op.id == OPID_CY || op.id == OPID_CZ || op.id == OPID_SWAP || op.id == OPID_RXX || op.id == OPID_RYY || op.id == OPID_RZZ || op.id == OPID_MAT2Q) && - (shot.qubit_state[op.q2].heat == -1.0) { - shot.op_type = OPID_ID; - shot.op_idx = op_idx; - return; + // Before doing further work, if any qubit for the gate is lost, dispatch + // the gate's configured loss policy (stamped on op.policy). For most policies + // this fully handles the op; APPLY_ANYWAY returns false so the gate runs as + // usual below. + if (gate_has_lost_operand(shot_idx, op_idx, op.q1, op.q2)) { + if (handle_lost_operand_policy(shot_idx, op_idx, op.q1, op.q2)) { + return; + } } if pauli_op_idx != 0 { diff --git a/source/simulators/src/noise_config.rs b/source/simulators/src/noise_config.rs index 7841909b45..22b97300fe 100644 --- a/source/simulators/src/noise_config.rs +++ b/source/simulators/src/noise_config.rs @@ -63,6 +63,51 @@ impl CumulativeNoiseConfig { } } +/// Specifies the behavior of a multi-qubit gate when at least one of its qubit +/// operands is lost. +/// +/// This lets users experiment with different lost-qubit gate behaviors +/// from Python (via the per-gate `on_loss` field of the noise config) +/// without modifying and recompiling the simulator. +#[derive(Clone, Copy, Debug, Default, PartialEq, Eq)] +pub enum LossPolicy { + /// If any of the qubit operands of a gate is lost, skip the gate entirely. + #[default] + Skip, + /// If any of the qubit operands of a gate is lost, propagate the loss to + /// the other operands (the surviving operands are measured, reset, and + /// flagged as lost) and skip the gate. + Propagate, + /// For multi-qubit rotations, degrade the unitary to its single-qubit + /// version applied to the surviving operand (e.g. `rxx` -> `rx`). + /// For gates with no single-qubit reduction (`cx`, `cy`, `cz`, `swap`, + /// and single-qubit gates) this falls back to [`LossPolicy::Skip`]. + Degrade, + /// Skip the gate and instead apply an `S` adjoint to each surviving operand. + ResidualSDagger, + /// Apply the unitary anyway, ignoring the loss. Lost operands (already + /// measured and reset to the zero state) remain flagged as lost. + ApplyAnyway, +} + +impl LossPolicy { + /// Encodes the policy as a `u32` for transport to the GPU shader. + /// + /// The values match the Python `LossPolicy` enum (`SKIP = 0` .. + /// `APPLY_ANYWAY = 4`). The value `0` is reserved by the shader to mean + /// "no policy stamped" and is never produced here. + #[must_use] + pub fn as_u32(self) -> u32 { + match self { + Self::Skip => 0, + Self::Propagate => 1, + Self::Degrade => 2, + Self::ResidualSDagger => 3, + Self::ApplyAnyway => 4, + } + } +} + /// Noise description for each operation. /// /// This is the format in which the user config files are @@ -126,10 +171,10 @@ impl NoiseConfig { cx: NoiseTable::::noiseless(2), cy: NoiseTable::::noiseless(2), cz: NoiseTable::::noiseless(2), - rxx: NoiseTable::::noiseless(2), - ryy: NoiseTable::::noiseless(2), - rzz: NoiseTable::::noiseless(2), - swap: NoiseTable::::noiseless(2), + rxx: NoiseTable::::noiseless_with_loss_policy(2, LossPolicy::Degrade), + ryy: NoiseTable::::noiseless_with_loss_policy(2, LossPolicy::Degrade), + rzz: NoiseTable::::noiseless_with_loss_policy(2, LossPolicy::Degrade), + swap: NoiseTable::::noiseless_with_loss_policy(2, LossPolicy::ApplyAnyway), ccx: NoiseTable::::noiseless(3), mov: NoiseTable::::noiseless(1), mz: NoiseTable::::noiseless(1), @@ -180,15 +225,23 @@ pub struct NoiseTable { pub qubits: u32, pub pauli_strings: Vec, pub probabilities: Vec, + /// The behavior of this gate when at least one of its operands is lost. + pub on_loss: LossPolicy, } impl NoiseTable { #[must_use] pub const fn noiseless(qubits: u32) -> Self { + Self::noiseless_with_loss_policy(qubits, LossPolicy::Skip) + } + + #[must_use] + pub const fn noiseless_with_loss_policy(qubits: u32, on_loss: LossPolicy) -> Self { Self { qubits, pauli_strings: Vec::new(), probabilities: Vec::new(), + on_loss, } } @@ -285,6 +338,8 @@ impl From> for CumulativeNoiseConfig { /// This is the internal format used by the simulator. #[derive(Default)] pub struct CumulativeNoiseTable { + /// The behavior of this gate when at least one of its operands is lost. + pub on_loss: LossPolicy, pub sampler: Sampler, } @@ -297,6 +352,7 @@ impl From> for CumulativeNoiseTable { .map(|p| Fault::from((p, qubits))); let probs = value.probabilities.into_iter().map(uq1_63::from_prob); Self { + on_loss: value.on_loss, sampler: Sampler::new(choices, probs), } } diff --git a/source/simulators/src/stabilizer_simulator.rs b/source/simulators/src/stabilizer_simulator.rs index de5be2a7c9..bb6ce5a1a6 100644 --- a/source/simulators/src/stabilizer_simulator.rs +++ b/source/simulators/src/stabilizer_simulator.rs @@ -7,7 +7,7 @@ pub mod operation; use crate::{ MeasurementResult, NearlyZero, QubitID, Simulator, - noise_config::{CumulativeNoiseConfig, Fault, FaultTerm, IntrinsicID}, + noise_config::{CumulativeNoiseConfig, Fault, FaultTerm, IntrinsicID, LossPolicy}, }; use operation::Operation; use paulimer::{ @@ -191,6 +191,13 @@ impl StabilizerSimulator { self.state.pauli(&observable); } + /// Applies an `S` adjoint to the given target + /// Used by the [`LossPolicy::ResidualSDagger`] behavior. + fn residual_s_dagger(&mut self, target: QubitID) { + self.apply_idle_noise(target); + self.state.apply_unitary(UnitaryOp::SqrtZInv, &[target]); + } + /// Records a z-measurement on the given `target`. fn record_mz(&mut self, target: QubitID, result_id: QubitID) { let measurement = self.mz_impl(target); @@ -246,6 +253,13 @@ impl StabilizerSimulator { MeasurementResult::Zero } } + + fn loss_impl(&mut self, target: QubitID) { + if !self.loss[target] { + self.mresetz_impl(target); + self.loss[target] = true; + } + } } impl Simulator for StabilizerSimulator { @@ -329,35 +343,83 @@ impl Simulator for StabilizerSimulator { } fn cx(&mut self, control: QubitID, target: QubitID) { - if !self.loss[control] && !self.loss[target] { - self.apply_idle_noise(control); - self.apply_idle_noise(target); - self.state - .apply_unitary(UnitaryOp::ControlledX, &[control, target]); + match (self.loss[control], self.loss[target]) { + (true, true) => (), + (true, false) | (false, true) => { + let remaining_qubit = if self.loss[control] { target } else { control }; + self.apply_idle_noise(remaining_qubit); + match self.noise_config.cx.on_loss { + LossPolicy::Skip | LossPolicy::Degrade => (), + LossPolicy::Propagate => self.loss_impl(remaining_qubit), + LossPolicy::ResidualSDagger => self.residual_s_dagger(remaining_qubit), + LossPolicy::ApplyAnyway => self + .state + .apply_unitary(UnitaryOp::ControlledX, &[control, target]), + } + } + (false, false) => { + self.apply_idle_noise(control); + self.apply_idle_noise(target); + self.state + .apply_unitary(UnitaryOp::ControlledX, &[control, target]); + } } // We still apply operation faults to non-lost qubits. apply_noise!(self, cx, &[control, target]); } fn cy(&mut self, control: QubitID, target: QubitID) { - if !self.loss[control] && !self.loss[target] { - self.apply_idle_noise(control); - self.apply_idle_noise(target); - self.state.apply_unitary(UnitaryOp::SqrtZInv, &[target]); - self.state - .apply_unitary(UnitaryOp::ControlledX, &[control, target]); - self.state.apply_unitary(UnitaryOp::SqrtZ, &[target]); + match (self.loss[control], self.loss[target]) { + (true, true) => (), + (true, false) | (false, true) => { + let remaining_qubit = if self.loss[control] { target } else { control }; + self.apply_idle_noise(remaining_qubit); + match self.noise_config.cy.on_loss { + LossPolicy::Skip | LossPolicy::Degrade => (), + LossPolicy::Propagate => self.loss_impl(remaining_qubit), + LossPolicy::ResidualSDagger => self.residual_s_dagger(remaining_qubit), + LossPolicy::ApplyAnyway => { + self.state.apply_unitary(UnitaryOp::SqrtZInv, &[target]); + self.state + .apply_unitary(UnitaryOp::ControlledX, &[control, target]); + self.state.apply_unitary(UnitaryOp::SqrtZ, &[target]); + } + } + } + (false, false) => { + self.apply_idle_noise(control); + self.apply_idle_noise(target); + self.state.apply_unitary(UnitaryOp::SqrtZInv, &[target]); + self.state + .apply_unitary(UnitaryOp::ControlledX, &[control, target]); + self.state.apply_unitary(UnitaryOp::SqrtZ, &[target]); + } } // We still apply operation faults to non-lost qubits. apply_noise!(self, cy, &[control, target]); } fn cz(&mut self, control: QubitID, target: QubitID) { - if !self.loss[control] && !self.loss[target] { - self.apply_idle_noise(control); - self.apply_idle_noise(target); - self.state - .apply_unitary(UnitaryOp::ControlledZ, &[control, target]); + match (self.loss[control], self.loss[target]) { + (true, true) => (), + (true, false) | (false, true) => { + let remaining_qubit = if self.loss[control] { target } else { control }; + self.apply_idle_noise(remaining_qubit); + match self.noise_config.cz.on_loss { + LossPolicy::Skip | LossPolicy::Degrade => (), + LossPolicy::Propagate => self.loss_impl(remaining_qubit), + LossPolicy::ResidualSDagger => self.residual_s_dagger(remaining_qubit), + LossPolicy::ApplyAnyway => self + .state + .apply_unitary(UnitaryOp::ControlledZ, &[control, target]), + } + } + (false, false) => { + self.apply_idle_noise(control); + self.apply_idle_noise(target); + self.state + .apply_unitary(UnitaryOp::ControlledZ, &[control, target]); + } } // We still apply operation faults to non-lost qubits. apply_noise!(self, cz, &[control, target]); @@ -420,8 +482,34 @@ impl Simulator for StabilizerSimulator { fn rxx(&mut self, angle: f64, q1: QubitID, q2: QubitID) { match (self.loss[q1], self.loss[q2]) { (true, true) => (), - (true, false) => self.rx(angle, q2), - (false, true) => self.rx(angle, q1), + (true, false) | (false, true) => { + let remaining_qubit = if self.loss[q1] { q2 } else { q1 }; + self.apply_idle_noise(remaining_qubit); + match self.noise_config.rxx.on_loss { + LossPolicy::Skip => (), + LossPolicy::Degrade => return self.rx(angle, remaining_qubit), + LossPolicy::Propagate => self.loss_impl(remaining_qubit), + LossPolicy::ResidualSDagger => self.residual_s_dagger(remaining_qubit), + LossPolicy::ApplyAnyway => { + // We can only perform rotations by multiples of PI / 2 in the stabilizer, so normalize the angle + // and check to see if it is supported. + let unitary = unitary_from_normalized_angle( + angle, + UnitaryOp::Z, + UnitaryOp::SqrtZ, + UnitaryOp::SqrtZInv, + ); + // NOTE: We perform the Rxx gate by changing basis to Y and performing the decomposition of Rzz. + self.state.apply_unitary(UnitaryOp::Hadamard, &[q1]); + self.state.apply_unitary(UnitaryOp::Hadamard, &[q2]); + self.state.apply_unitary(UnitaryOp::ControlledX, &[q2, q1]); + self.state.apply_unitary(unitary, &[q1]); + self.state.apply_unitary(UnitaryOp::ControlledX, &[q2, q1]); + self.state.apply_unitary(UnitaryOp::Hadamard, &[q1]); + self.state.apply_unitary(UnitaryOp::Hadamard, &[q2]); + } + } + } (false, false) => { self.apply_idle_noise(q1); self.apply_idle_noise(q2); @@ -442,17 +530,42 @@ impl Simulator for StabilizerSimulator { self.state.apply_unitary(UnitaryOp::ControlledX, &[q2, q1]); self.state.apply_unitary(UnitaryOp::Hadamard, &[q1]); self.state.apply_unitary(UnitaryOp::Hadamard, &[q2]); - - apply_noise!(self, rxx, &[q1, q2]); } } + apply_noise!(self, rxx, &[q1, q2]); } fn ryy(&mut self, angle: f64, q1: QubitID, q2: QubitID) { match (self.loss[q1], self.loss[q2]) { (true, true) => (), - (true, false) => self.ry(angle, q2), - (false, true) => self.ry(angle, q1), + (true, false) | (false, true) => { + let remaining_qubit = if self.loss[q1] { q2 } else { q1 }; + self.apply_idle_noise(remaining_qubit); + match self.noise_config.ryy.on_loss { + LossPolicy::Skip => (), + LossPolicy::Degrade => return self.ry(angle, remaining_qubit), + LossPolicy::Propagate => self.loss_impl(remaining_qubit), + LossPolicy::ResidualSDagger => self.residual_s_dagger(remaining_qubit), + LossPolicy::ApplyAnyway => { + // We can only perform rotations by multiples of PI / 2 in the stabilizer, so normalize the angle + // and check to see if it is supported. + let unitary = unitary_from_normalized_angle( + angle, + UnitaryOp::Z, + UnitaryOp::SqrtZ, + UnitaryOp::SqrtZInv, + ); + // NOTE: We perform the Ryy gate by changing basis to Z and performing the decomposition of Rzz. + self.state.apply_unitary(UnitaryOp::SqrtX, &[q1]); + self.state.apply_unitary(UnitaryOp::SqrtX, &[q2]); + self.state.apply_unitary(UnitaryOp::ControlledX, &[q2, q1]); + self.state.apply_unitary(unitary, &[q1]); + self.state.apply_unitary(UnitaryOp::ControlledX, &[q2, q1]); + self.state.apply_unitary(UnitaryOp::SqrtXInv, &[q1]); + self.state.apply_unitary(UnitaryOp::SqrtXInv, &[q2]); + } + } + } (false, false) => { self.apply_idle_noise(q1); self.apply_idle_noise(q2); @@ -473,17 +586,37 @@ impl Simulator for StabilizerSimulator { self.state.apply_unitary(UnitaryOp::ControlledX, &[q2, q1]); self.state.apply_unitary(UnitaryOp::SqrtXInv, &[q1]); self.state.apply_unitary(UnitaryOp::SqrtXInv, &[q2]); - - apply_noise!(self, ryy, &[q1, q2]); } } + apply_noise!(self, ryy, &[q1, q2]); } fn rzz(&mut self, angle: f64, q1: QubitID, q2: QubitID) { match (self.loss[q1], self.loss[q2]) { (true, true) => (), - (true, false) => self.rz(angle, q2), - (false, true) => self.rz(angle, q1), + (true, false) | (false, true) => { + let remaining_qubit = if self.loss[q1] { q2 } else { q1 }; + self.apply_idle_noise(remaining_qubit); + match self.noise_config.rzz.on_loss { + LossPolicy::Skip => (), + LossPolicy::Degrade => return self.rz(angle, remaining_qubit), + LossPolicy::Propagate => self.loss_impl(remaining_qubit), + LossPolicy::ResidualSDagger => self.residual_s_dagger(remaining_qubit), + LossPolicy::ApplyAnyway => { + // We can only perform rotations by multiples of PI / 2 in the stabilizer, so normalize the angle + // and check to see if it is supported. + let unitary = unitary_from_normalized_angle( + angle, + UnitaryOp::Z, + UnitaryOp::SqrtZ, + UnitaryOp::SqrtZInv, + ); + self.state.apply_unitary(UnitaryOp::ControlledX, &[q2, q1]); + self.state.apply_unitary(unitary, &[q1]); + self.state.apply_unitary(UnitaryOp::ControlledX, &[q2, q1]); + } + } + } (false, false) => { self.apply_idle_noise(q1); self.apply_idle_noise(q2); @@ -499,38 +632,48 @@ impl Simulator for StabilizerSimulator { self.state.apply_unitary(UnitaryOp::ControlledX, &[q2, q1]); self.state.apply_unitary(unitary, &[q1]); self.state.apply_unitary(UnitaryOp::ControlledX, &[q2, q1]); - - apply_noise!(self, rzz, &[q1, q2]); } } + apply_noise!(self, rzz, &[q1, q2]); } fn swap(&mut self, q1: QubitID, q2: QubitID) { + // There are three kinds of swaps: + // 1. A logical swap, also called a relabel. + // 2. A swap by physically exchanging the location of the qubits. + // 3. An exchange of information by doing three CX. + // + // This method is concerned with the kinds (1) and (2), since (3) + // gets decomposed into other instructions before making it to the simulator. + // In both (1) and (2), the loss state of the qubits gets exchanged. + match (self.loss[q1], self.loss[q2]) { (true, true) => (), - (true, false) => { - self.apply_idle_noise(q2); - self.state.apply_permutation(&[1, 0], &[q1, q2]); - } - (false, true) => { - self.apply_idle_noise(q1); - self.state.apply_permutation(&[1, 0], &[q1, q2]); + (true, false) | (false, true) => { + let lost_qubit = if self.loss[q1] { q1 } else { q2 }; + let remaining_qubit = if self.loss[q1] { q2 } else { q1 }; + self.apply_idle_noise(remaining_qubit); + match self.noise_config.swap.on_loss { + LossPolicy::Skip | LossPolicy::Degrade => (), + LossPolicy::Propagate => self.loss_impl(remaining_qubit), + LossPolicy::ResidualSDagger => { + self.state.apply_permutation(&[1, 0], &[q1, q2]); + self.residual_s_dagger(lost_qubit); + self.loss.swap(q1, q2); + } + LossPolicy::ApplyAnyway => { + self.state.apply_permutation(&[1, 0], &[q1, q2]); + self.loss.swap(q1, q2); + } + } } (false, false) => { self.apply_idle_noise(q1); self.apply_idle_noise(q2); self.state.apply_permutation(&[1, 0], &[q1, q2]); + self.loss.swap(q1, q2); } } - // There are three kinds of swaps: - // 1. A logical swap, also called a relabel. - // 2. A swap by physically exchanging the location of the qubits. - // 3. An exchange of information by doing three CX. - // - // This method is concerned with the kinds (1) and (2), since (3) - // gets decomposed into other instructions before making it to the simulator. - // In both (1) and (2), the loss state of the qubits gets exchanged. - self.loss.swap(q1, q2); // Is up to the user if swap is a virtual operation or not. // If they don't specify noise/loss probability for swap, then it is virtual.