Noise

A real quantum computer is never perfect. dense_evolution.noise has three ways to put that imperfection into a simulation, depending on what you need from it.

Step 1. Build a circuit

import numpy as np
import dense_evolution as de

qasm = 'OPENQASM 2.0; include "qelib1.inc"; qreg q[2]; creg c[2]; h q[0]; cx q[0],q[1]; measure q -> c;'
circuit = de.QASMParser().parse(qasm)
sim = de.DenseSVSimulator(2)
sim.run_circuit_jit(circuit.to_tuples())
sv = np.asarray(sim.get_statevector())

sv is the exact, noiseless result. Every step below starts from it.

Step 2. Add noise to it

from dense_evolution.noise import NoiseModel

rng = np.random.default_rng(0)
sv_noisy = NoiseModel.apply_to_sv(sv.copy(), 2, "depolarizing", 0.1, rng=rng)
round(float(np.vdot(sv_noisy, sv_noisy).real), 4)

1.0

NoiseModel.apply_to_sv takes the statevector, the qubit count, which error model to simulate, and how strong it is. The result is still a valid, correctly normalised state — noise redistributes probability, it never breaks it.

Step 3. See every available model

NoiseModel.MODELS

['ideal', 'depolarizing', 'bitflip', 'phaseflip', 'amplitude_damping', 'combined']

Swap "depolarizing" above for any other name to use it instead. p is that model's error probability, except for "amplitude_damping", where p is the decay rate.

NoiseModel.kraus_description("bitflip")["physical"]

'Bit flip σ_x with probability p'

Step 4. Use a real device's own noise instead of a made-up number

from qiskit_ibm_runtime.fake_provider import FakeSherbrooke
from dense_evolution.interop import noise_model_from_qiskit_backend

backend = FakeSherbrooke()
spec = noise_model_from_qiskit_backend(backend)
len(spec)

652

Each entry in spec is a real, measured error rate for one gate on one qubit of that device — pass any one of them as model/p/qubits to NoiseModel.apply_to_sv above, instead of guessing a number.

Step 5. Make the noise strength part of a gradient

import jax
from dense_evolution.noise import NoiseSpec

key = jax.random.PRNGKey(0)
spec = NoiseSpec(model="depolarizing", p=0.05, jax_key=key, qubits=[0, 1])
spec

NoiseSpec(model='depolarizing', p=0.05, qubits=(0, 1))

Pass a NoiseSpec as circuit_to_energy_fn's noise= argument to trace noise strength through jax.grad along with every other parameter, instead of applying it as a separate step outside the gradient.

Step 6. Search for a worst-case coherent error against a QEC code

The three steps above are all stochastic, single-qubit errors. This one is different: a coherent error spread across every qubit at once, searched for by gradient ascent instead of sampled at random.

import jax.numpy as jnp
from dense_evolution.noise.coherent_attack import craft_adversarial_delta_constrained

x_stabilizers = ["IIIXXXX", "IXXIIXX", "XIXIXIX"]
sv_code = jnp.ones(2 ** 7, dtype=jnp.complex128) / jnp.sqrt(2 ** 7)
delta, leakage, _ = craft_adversarial_delta_constrained(
    sv_code, x_stabilizers, epsilon=0.5, linf_cap=0.1, n_steps=50
)
delta.shape

(7,)

delta[q] is the coherent rz angle the search found for qubit q. The linf_cap argument matters: without it (craft_adversarial_delta, no per-qubit cap), the search degenerately concentrates the whole error on one qubit — which a real decoder corrects exactly every time, making that "worst case" actually harmless. Capping each qubit's share forces the search into directions that spread across multiple qubits, which is where real decoder failures happen.

Step 7. Simulate a cosmic-ray burst hitting the chip

A cosmic ray or gamma ray striking the chip briefly raises every qubit's decay rate, then lets it fall back to normal — this profile reproduces a real measured event from a published 26-qubit device.

import numpy as np
from dense_evolution.mitigation.zne import cosmic_ray_burst_profile

baseline = 0.001
times_us = np.array([0.0, 10.0, 1000.0, 50000.0])
profile = cosmic_ray_burst_profile(times_us, baseline_gamma=baseline)
[round(float(x), 6) for x in profile]

[0.001, 0.002487, 0.003599, 0.001372]

At t=0 (the impact instant) the decay rate is still just baseline. It climbs to over 3x baseline within about a millisecond, then relaxes back down over tens of milliseconds. Feed the result into continuous_dissipative_evolve alongside amplitude_damping_channel to inject a real burst shape into a circuit or QEC study, instead of a constant noise rate.

Step 8. Apply noise to a density matrix instead of a statevector

Steps 2-7 above all work on a statevector. Density-matrix ZNE (see Density-matrix ZNE healing) needs noise applied directly to a density matrix instead.

rho = np.outer(sv, sv.conj())

from dense_evolution.noise import global_depolarizing_channel

rho_noisy = global_depolarizing_channel(rho, 0.1)
round(float(np.trace(rho_noisy).real), 6)

1.0

global_depolarizing_channel mixes the whole register toward the fully mixed state as one unit — a different physical process from NoiseModel's "depolarizing", which acts independently per qubit. Use this one for a state-prep/measurement (SPAM) error reported as a single number over the whole register, not per-qubit gate noise. amplitude_damping_channel(rho, gamma) is the density-matrix equivalent of Step 3's "amplitude_damping", single-qubit only, and is what Step 7's burst profile is meant to drive.

Step 9. Make the noise strength oscillate instead of scaling smoothly

from dense_evolution.noise import oscillating_p_eff

[round(float(oscillating_p_eff(base_p=0.1, factor=f, freq=2.0, amp=0.5)), 4) for f in [0.0, 1.0, 2.0, 3.0]]

[0.1, 0.15, 0.1, 0.05]

Most mitigation techniques (Richardson extrapolation, for one) assume noise grows smoothly as you scale it up. oscillating_p_eff builds a deliberately non-smooth noise-vs-scale relationship instead, to check whether a technique still works when that assumption doesn't hold.

---

Details

Kraus formulas

Model	Kraus operators
depolarizing	{√(1-p)I, √(p/3)X, √(p/3)Y, √(p/3)Z}
bitflip	{√(1-p)I, √p·X}
phaseflip	{√(1-p)I, √p·Z}
amplitude_damping	K0=diag(1,√(1-γ)), K1=[[0,√γ],[0,0]]
combined	depolarizing(p/2) then amplitude_damping(p/3)
ideal	identity

Every channel draws one fire/no-fire decision per qubit per shot (plus one Pauli choice for depolarizing/combined's depolarizing sub-step) — the same single-Pauli-per-qubit-per-shot convention STIM's DEPOLARIZE1(p) uses. Prior to v8.1.57, each channel instead drew one independent decision per computational-basis amplitude pair, which over-decohered entangled states (a measured value dropped from 1.0 to 0.31 at p=0.15 on one test case) — fixed by drawing one decision per qubit and applying it uniformly.

Photon loss is not a separate channel

Photon loss on a dual-rail-encoded qubit is exactly this library's amplitude_damping channel (Step 3 above) — there is no separate photon noise model, and none is needed.

Moved here from mitigation.zne

global_depolarizing_channel, amplitude_damping_channel, cosmic_ray_burst_profile, and oscillating_p_eff used to live in dense_evolution.mitigation.zne — they generate noise, they don't mitigate it, so that was the wrong home. dense_evolution.mitigation.zne still re-exports all four for backward compatibility.

Coherent adversarial noise: an honest negative result

craft_adversarial_delta (no L-infinity cap) was tested against a real decoder for the Steane [[7,1,3]] code and found to converge to a direction with zero real decoder-failure rate — a coherent error concentrated on a single qubit always collapses to something a distance-3 code corrects exactly, so the search's own "worst case" is actually safe. Random noise directions of the same L2 budget, spread across multiple qubits, failed the decoder readily by comparison. craft_adversarial_delta_constrained's L-infinity cap fixes this by forbidding that degenerate solution. Promoted from Dense-Evolution-Discovery's Steane investigation, generalized from a fixed 7-qubit table to any stabilizer list.

noise

Every way to put noise into a Dense-Evolution simulation, in one place.

• NoiseModel (.kraus_channels, dispatching to one file per channel under .kraus: ideal, depolarizing, bitflip, phaseflip, amplitude_damping, combined) -- 6 stochastic single-qubit Kraus channels, applied directly to a statevector via apply_to_sv.
• NoiseSpec (.differentiable) -- the native JAX-differentiable representation of a noise configuration, so noise strength itself can be a traced/differentiable value inside circuit_to_energy_fn.
• Coherent adversarial noise (.coherent_attack) -- a genuinely continuous, multi-qubit coherent error channel (apply_rz_all) and a JAX-differentiable search (craft_adversarial_delta, craft_adversarial_delta_constrained) for a worst-case direction against a stabilizer code's syndrome -- promoted from Dense-Evolution-Discovery's Steane [[7,1,3]] investigation, including its honest negative result (see coherent_attack's module docstring).
• Density-matrix channels (.density_matrix_channels) -- global_depolarizing_channel, amplitude_damping_channel: noise applied directly to a density matrix instead of a statevector, for density-matrix ZNE's noise ensemble.
• cosmic_ray_burst_profile (.cosmic_ray) -- a real, time-dependent noise-strength profile for a cosmic-ray-induced quasiparticle burst.
• oscillating_p_eff (.oscillating) -- a noise strength that oscillates instead of scaling smoothly, for stress-testing mitigation techniques that assume smoothness.

Real device noise from a Qiskit backend's own calibration data (noise_model_from_qiskit_backend) lives in dense_evolution.interop, not here, since it needs a Qiskit BackendV2 object as input rather than a noise-specific dependency.

Everything in this package previously lived scattered across dense_evolution.circuits.registry (alongside unrelated hardware-detection code) and dense_evolution.mitigation.zne (alongside unrelated mitigation techniques, which only ever cancel noise, never generate it). Both modules re-export the relevant names from here for backward compatibility, but dense_evolution.noise is the canonical import path for new code.

NoiseModel

Stochastic single-qubit Kraus channels applied directly to a statevector. Each channel is a separate, importable module under dense_evolution. noise.kraus (dense_evolution.noise.kraus.depolarizing, etc.) -- this class is the shared dispatcher: RNG/key setup, the per-qubit loop, and final normalisation, common to every channel.

All channels are mathematically correct Kraus maps: - trace is preserved (normalisation enforced at the end) - phaseflip applies Z with probability p per qubit (non-deterministic) - amplitude_damping applies the correct K0/K1 Kraus operators - combined is a true worst-case NISQ mixture of all three Pauli errors plus amplitude damping

Supported models

'ideal' identity — no modification 'depolarizing' {√(1-p)I, √(p/3)X, √(p/3)Y, √(p/3)Z} 'bitflip' {√(1-p)I, √p·X} 'phaseflip' {√(1-p)I, √p·Z} ← was broken, now fixed 'amplitude_damping'{K0=diag(1,√(1-γ)), K1=[[0,√γ],[0,0]]} 'combined' depolarizing(p/2) + amplitude_damping(p/3), renormalised

Every channel draws one fire/no-fire decision per qubit per shot (plus one Pauli choice for depolarizing/combined's depolarizing sub-step), applied identically across the whole statevector -- the same single-Pauli-per-qubit-per-shot convention STIM's DEPOLARIZE1(p) uses. Prior to v8.1.57, every channel instead drew 2**(n-1) INDEPENDENT decisions per qubit per shot, one per amplitude pair (i.e. one per branch of the other n-1 qubits) -- inert on a product state, but on an entangled state it over-decohered any coherence-sensitive (off-diagonal) observable, up to hundreds of sigma vs the exact density-matrix Kraus-sum result on test cases (e.g. per-branch sampling dropped a measured value from 1.0 to 0.31 at p=0.15 on one such test -- see the v8.1.57 changelog entry for the full reproduction).

apply_to_sv staticmethod

apply_to_sv(
    sv: ndarray,
    n: int,
    model: str,
    p: float,
    rng: Optional[Generator] = None,
    qubits: Optional[List[int]] = None,
    jax_key: Optional[Any] = None,
) -> np.ndarray

Apply a stochastic Kraus channel to statevector sv in-place (numpy path) or via functional updates (JAX path).

Parameters:

Name	Type	Description	Default
sv	ndarray		required
n	int		required
model	str		required
p	float		required
rng	Optional[Generator]	*sv* is a NumPy array. When *sv* is a JAX array, `rng` used to be silently ignored in favor of `jax_key` (or a non-reproducible OS-entropy key if that was also None -- issue #7); it is now used to *derive* a reproducible jax_key (`rng.integers(...)` seeds `jax.random.PRNGKey`) whenever `jax_key` isn't given explicitly, so seeding `rng` has the effect a caller expects on both array types instead of only on one of them.	None
qubits	Optional[List[int]]		None
jax_key	optional JAX PRNGKey, only meaningful when *sv* is a JAX	array. Takes precedence over `rng` when both are given (explicit key beats a derived one). Created from OS entropy if neither `jax_key` nor `rng` is given.	None

Returns:

Type	Description
Normalised statevector (same array type as input).

Examples:

A real quantum computer is never perfect -- every gate has some chance of error. Once you have a statevector from running your own QASM circuit (the same circuit as the getting-started example), this function is how you find out what a noisy device would have actually given you instead.

Start from the circuit and statevector you already have:

>>> import numpy as np
>>> import dense_evolution as de
>>> qasm = 'OPENQASM 2.0; include "qelib1.inc"; qreg q[2]; creg c[2]; h q[0]; barrier q; cx q[0],q[1]; measure q -> c;'
>>> circuit = de.QASMParser().parse(qasm)
>>> sim = de.DenseSVSimulator(2)
>>> sim.run_circuit(circuit.to_tuples())
>>> sv = np.asarray(sim.get_statevector())

(the barrier is parsed and ignored -- it never becomes a gate tuple, so it has no effect on the statevector, only on how the circuit reads.)

Call NoiseModel.apply_to_sv on that same statevector, telling it the qubit count, which error model to simulate, and how strong it is:

>>> from dense_evolution.noise import NoiseModel
>>> rng = np.random.default_rng(0)
>>> sv_noisy = NoiseModel.apply_to_sv(sv.copy(), 2, 'depolarizing', 0.1, rng=rng)
>>> round(float(np.vdot(sv_noisy, sv_noisy).real), 4)  # still a valid, normalised state
1.0

'depolarizing' above is one of six models; pick any other one the same way, by name:

>>> NoiseModel.MODELS
['ideal', 'depolarizing', 'bitflip', 'phaseflip', 'amplitude_damping', 'combined']

p is that model's error probability (or damping rate for 'amplitude_damping') -- 0.1 above means each qubit has a 10% chance of a random Pauli error per call. Run it many times and average (see Density-matrix ZNE healing) to see what a real noisy device's typical output looks like, not just one random draw.

Source code in dense_evolution/noise/kraus_channels.py

@staticmethod
def apply_to_sv(
    sv:       np.ndarray,
    n:        int,
    model:    str,
    p:        float,
    rng:      Optional[np.random.Generator] = None,
    qubits:   Optional[List[int]] = None,
    jax_key:  Optional[Any] = None,
) -> np.ndarray:
    """
    Apply a stochastic Kraus channel to statevector *sv* in-place
    (numpy path) or via functional updates (JAX path).

    Parameters
    ----------
    sv      : complex statevector of length 2**n
    n       : number of qubits
    model   : one of NoiseModel.MODELS
    p       : error probability (or damping rate γ for amplitude_damping)
    rng     : optional pre-seeded numpy Generator. Used directly when
              *sv* is a NumPy array. When *sv* is a JAX array, `rng`
              used to be silently ignored in favor of `jax_key` (or a
              non-reproducible OS-entropy key if that was also None
              -- issue #7); it is now used to *derive* a reproducible
              jax_key (`rng.integers(...)` seeds `jax.random.PRNGKey`)
              whenever `jax_key` isn't given explicitly, so seeding
              `rng` has the effect a caller expects on both array
              types instead of only on one of them.
    qubits  : subset of qubits to apply the channel to; defaults to all
    jax_key : optional JAX PRNGKey, only meaningful when *sv* is a JAX
              array. Takes precedence over `rng` when both are given
              (explicit key beats a derived one). Created from OS
              entropy if neither `jax_key` nor `rng` is given.

    Returns
    -------
    Normalised statevector (same array type as input).

    Examples
    --------
    A real quantum computer is never perfect -- every gate has some chance of
    error. Once you have a statevector from running your own QASM circuit (the
    same circuit as the getting-started example), this function is how you find
    out what a noisy device would have actually given you instead.

    Start from the circuit and statevector you already have:

    >>> import numpy as np
    >>> import dense_evolution as de
    >>> qasm = 'OPENQASM 2.0; include "qelib1.inc"; qreg q[2]; creg c[2]; h q[0]; barrier q; cx q[0],q[1]; measure q -> c;'
    >>> circuit = de.QASMParser().parse(qasm)
    >>> sim = de.DenseSVSimulator(2)
    >>> sim.run_circuit(circuit.to_tuples())
    >>> sv = np.asarray(sim.get_statevector())

    (the `barrier` is parsed and ignored -- it never becomes a gate tuple, so it
    has no effect on the statevector, only on how the circuit reads.)

    Call `NoiseModel.apply_to_sv` on that same statevector, telling it the
    qubit count, which error model to simulate, and how strong it is:

    >>> from dense_evolution.noise import NoiseModel
    >>> rng = np.random.default_rng(0)
    >>> sv_noisy = NoiseModel.apply_to_sv(sv.copy(), 2, 'depolarizing', 0.1, rng=rng)
    >>> round(float(np.vdot(sv_noisy, sv_noisy).real), 4)  # still a valid, normalised state
    1.0

    `'depolarizing'` above is one of six models; pick any other one the same way,
    by name:

    >>> NoiseModel.MODELS
    ['ideal', 'depolarizing', 'bitflip', 'phaseflip', 'amplitude_damping', 'combined']

    `p` is that model's error probability (or damping rate for
    `'amplitude_damping'`) -- 0.1 above means each qubit has a 10% chance of a
    random Pauli error per call. Run it many times and average (see
    [Density-matrix ZNE healing](../examples.md#density-matrix-zne-healing))
    to see what a real noisy device's *typical* output looks like, not just one
    random draw.
    """
    if model == 'ideal':
        return sv
    try:
        if p <= 0.0:
            return sv
    except jax.errors.TracerBoolConversionError:
        # p is a traced value (e.g. flowing through jax.jit/vmap/grad
        # as a NoiseSpec pytree leaf) -- can't early-exit on a Python
        # bool of it. Falling through is still correct: every channel
        # below already reduces to a no-op at p=0 (`fire = r < p` is
        # always False), this skips only the eager-mode optimization,
        # not correctness.
        pass

    channel = _CHANNELS[model]
    is_jax = HAS_JAX and isinstance(sv, jnp.ndarray)
    dim    = len(sv)

    # ── RNG initialisation ────────────────────────────────────────
    if is_jax:
        if jax_key is not None:
            key = jax_key
        elif rng is not None:
            # Derive a reproducible JAX key from the caller's seeded
            # NumPy generator instead of silently ignoring it -- each
            # call advances `rng`'s state, so a fresh, identically-
            # seeded `rng` reproduces the exact same sequence of keys
            # across separate runs (same guarantee the NumPy path
            # already gives).
            key = jax.random.PRNGKey(int(rng.integers(0, 2**32 - 1)))
        else:
            seed_bytes = os.urandom(4)
            jax_seed   = int.from_bytes(seed_bytes, byteorder='big')
            jax_seed  ^= time.perf_counter_ns() & 0xFFFF_FFFF
            key = jax.random.PRNGKey(jax_seed)
    else:
        key = None
        if rng is None:
            rng = _fresh_rng()

    target_qubits = qubits if qubits is not None else list(range(n))
    sv_out = sv  # JAX: functional; NumPy: will be modified in-place copy

    if not is_jax:
        sv_out = sv.copy()  # never mutate the caller's array

    for q in target_qubits:
        idx_0, idx_1 = _qubit_index_pairs(dim, q)
        sv_out, key = channel.apply(sv_out, idx_0, idx_1, p, rng, key, is_jax)

    # ── normalise ─────────────────────────────────────────────────
    if is_jax:
        norm = jnp.linalg.norm(sv_out)
        return sv_out / (norm + 1e-15)
    else:
        norm = np.linalg.norm(sv_out)
        return sv_out / (norm + 1e-15)

kraus_description staticmethod

kraus_description(model: str) -> Dict

Human-readable Kraus-operator formula and physical meaning for one of NoiseModel.MODELS.

Examples:

>>> from dense_evolution.noise import NoiseModel
>>> NoiseModel.kraus_description('bitflip')['physical']
'Bit flip σ_x with probability p'

Source code in dense_evolution/noise/kraus_channels.py

@staticmethod
def kraus_description(model: str) -> Dict:
    """Human-readable Kraus-operator formula and physical meaning for
    one of `NoiseModel.MODELS`.

    Examples
    --------
    >>> from dense_evolution.noise import NoiseModel
    >>> NoiseModel.kraus_description('bitflip')['physical']
    'Bit flip σ_x with probability p'
    """
    desc = {
        'ideal': {
            'kraus': 1,
            'formula': 'K₀ = I',
            'physical': 'No noise',
        },
        'depolarizing': {
            'kraus': 4,
            'formula': 'K₀=√(1-p)I  K₁=√(p/3)X  K₂=√(p/3)Y  K₃=√(p/3)Z',
            'physical': 'Isotropic Pauli error — equiprobable X, Y, Z',
        },
        'bitflip': {
            'kraus': 2,
            'formula': 'K₀=√(1-p)I  K₁=√p·X',
            'physical': 'Bit flip σ_x with probability p',
        },
        'phaseflip': {
            'kraus': 2,
            'formula': 'K₀=√(1-p)I  K₁=√p·Z',
            'physical': 'Pure dephasing σ_z with probability p',
        },
        'amplitude_damping': {
            'kraus': 2,
            'formula': 'K₀=diag(1,√(1-γ))  K₁=[[0,√γ],[0,0]]',
            'physical': 'T₁ energy relaxation |1⟩→|0⟩ with rate γ',
        },
        'combined': {
            'kraus': 6,
            'formula': 'Depolarizing(p/2) ∘ AmplitudeDamping(p/3)',
            'physical': 'Worst-case NISQ: dephasing + relaxation',
        },
    }
    return desc.get(model, desc['ideal'])

NoiseSpec

NoiseSpec(
    model: str,
    p,
    jax_key,
    qubits: Optional[List[int]] = None,
)

Native JAX-differentiable representation of a noise configuration -- a real JAX PyTree, so noise parameters thread through jax.jit/jax.grad/ jax.vmap natively -- e.g. as the noise= argument to circuit_to_energy_fn's energy_fn -- instead of being applied as an external, Python-side step around the already-traced circuit (the old way: build sv, exit the trace, call apply_to_sv separately).

model/qubits are static (aux_data): they select which code path runs, not values to differentiate or batch over -- the same role static_argnames plays for a plain jax.jit function, but automatic here because it's part of the pytree structure. p/jax_key are pytree leaves (children): p can be a traced/differentiable value (e.g. optimizing noise strength itself), and jax_key flows through jit/vmap/scan the way any other JAX array does -- no external Python-level key management, no OS-entropy fallback (unlike apply_to_sv called standalone with jax_key=None), so a NoiseSpec's result is always reproducible from the key it was built with.

jax_key is required (not Optional) -- the whole point of wiring noise into the traced computation this way is to remove the need for an external, ad-hoc key-management workaround; a caller who wants a fresh key per call should split one themselves (jax.random.split) and build a fresh NoiseSpec, the same as any other JAX-idiomatic stateless-key pattern.

Examples:

>>> import jax
>>> from dense_evolution.noise import NoiseSpec
>>> key = jax.random.PRNGKey(0)
>>> spec = NoiseSpec(model="depolarizing", p=0.05, jax_key=key, qubits=[0, 1])
>>> spec
NoiseSpec(model='depolarizing', p=0.05, qubits=(0, 1))

Source code in dense_evolution/noise/differentiable.py

def __init__(self, model: str, p, jax_key, qubits: Optional[List[int]] = None):
    self.model = model
    self.p = p
    self.jax_key = jax_key
    self.qubits = tuple(qubits) if qubits is not None else None

apply_rz_all

apply_rz_all(sv0: ndarray, delta: ndarray) -> jnp.ndarray

Coherent per-qubit rz(delta_q) applied to every qubit of sv0 at once. rz gates are diagonal and all commute, so this is exact elementwise phase multiplication, not a per-gate circuit simulation -- and fully JAX-differentiable in delta.

Parameters:

Name	Type	Description	Default
sv0	statevector, length 2**n_qubits		required
delta	real array, length n_qubits -- rz angle for each qubit		required

Examples:

>>> import numpy as np
>>> import jax.numpy as jnp
>>> sv0 = jnp.array([1.0, 0.0], dtype=jnp.complex128)
>>> sv1 = apply_rz_all(sv0, jnp.array([np.pi]))
>>> round(float(jnp.abs(sv1[0]) ** 2), 6)
1.0

Source code in dense_evolution/noise/coherent_attack.py

def apply_rz_all(sv0: jnp.ndarray, delta: jnp.ndarray) -> jnp.ndarray:
    """Coherent per-qubit rz(delta_q) applied to every qubit of `sv0` at
    once. rz gates are diagonal and all commute, so this is exact
    elementwise phase multiplication, not a per-gate circuit simulation
    -- and fully JAX-differentiable in `delta`.

    Parameters
    ----------
    sv0 : statevector, length 2**n_qubits
    delta : real array, length n_qubits -- rz angle for each qubit

    Examples
    --------
    >>> import numpy as np
    >>> import jax.numpy as jnp
    >>> sv0 = jnp.array([1.0, 0.0], dtype=jnp.complex128)
    >>> sv1 = apply_rz_all(sv0, jnp.array([np.pi]))
    >>> round(float(jnp.abs(sv1[0]) ** 2), 6)
    1.0
    """
    n = delta.shape[0]
    dim = sv0.shape[0]
    idx = jnp.arange(dim, dtype=jnp.int32)
    bit_pos = jnp.arange(n - 1, -1, -1, dtype=jnp.int32)
    bits = (idx[:, None] >> bit_pos[None, :]) & 1
    s = (1 - 2 * bits).astype(jnp.float64)
    phase_arg = -0.5 * jnp.sum(delta[None, :] * s, axis=1)
    return sv0 * jnp.exp(1j * phase_arg)

x_stabilizer_leakage

x_stabilizer_leakage(
    delta: ndarray, sv0: ndarray, stabilizers
) -> jnp.ndarray

Total leakage of sv0, after a coherent apply_rz_all(sv0, delta) perturbation, out of the +1 joint eigenspace of stabilizers -- sum_i (1 - ) / 2, one term per generator. A smooth, bounded ([0, len(stabilizers)]) proxy for how much the coherent error disturbs the syndrome; NOT the actual decoder failure probability (see module docstring for the documented gap between the two). stabilizers should be X-type generators to pair with an rz (Z-type-diagonal) coherent error -- the same reasoning dense_evolution.qec.compute_syndrome uses generically for any Pauli-string stabilizer list.

Source code in dense_evolution/noise/coherent_attack.py

def x_stabilizer_leakage(delta: jnp.ndarray, sv0: jnp.ndarray, stabilizers) -> jnp.ndarray:
    """Total leakage of `sv0`, after a coherent `apply_rz_all(sv0, delta)`
    perturbation, out of the +1 joint eigenspace of `stabilizers` --
    sum_i (1 - <stabilizer_i>) / 2, one term per generator. A smooth,
    bounded ([0, len(stabilizers)]) proxy for how much the coherent
    error disturbs the syndrome; NOT the actual decoder failure
    probability (see module docstring for the documented gap between
    the two). `stabilizers` should be X-type generators to pair with an
    rz (Z-type-diagonal) coherent error -- the same reasoning
    `dense_evolution.qec.compute_syndrome` uses generically for any
    Pauli-string stabilizer list.
    """
    sv = apply_rz_all(sv0, delta)
    dim = sv0.shape[0]
    idx = jnp.arange(dim, dtype=jnp.int32)
    total = 0.0
    for g in stabilizers:
        mask = _flip_mask(g)
        src = idx ^ mask
        expectation = jnp.real(jnp.sum(jnp.conj(sv) * sv[src]))
        total = total + (1.0 - expectation) / 2.0
    return total

craft_adversarial_delta

craft_adversarial_delta(
    sv0: ndarray,
    stabilizers,
    epsilon: float,
    n_steps: int = 150,
    step_size: float = 0.05,
    seed: int = 0,
)

Gradient-ascent PGD on x_stabilizer_leakage, projected into the L2 epsilon-ball around delta=0 after every step. Returns (best_delta as numpy, best_leakage, leakage_history).

See the module docstring: on a distance-3 code this unconstrained search finds a direction with real decoder-failure rate 0 -- use craft_adversarial_delta_constrained for a genuine worst-case test.

Source code in dense_evolution/noise/coherent_attack.py

def craft_adversarial_delta(sv0: jnp.ndarray, stabilizers, epsilon: float,
                             n_steps: int = 150, step_size: float = 0.05, seed: int = 0):
    """Gradient-ascent PGD on `x_stabilizer_leakage`, projected into the L2
    epsilon-ball around delta=0 after every step. Returns (best_delta as
    numpy, best_leakage, leakage_history).

    See the module docstring: on a distance-3 code this unconstrained
    search finds a direction with real decoder-failure rate 0 -- use
    `craft_adversarial_delta_constrained` for a genuine worst-case test.
    """
    n_qubits = len(stabilizers[0])
    leakage_grad = jax.grad(lambda d, sv: x_stabilizer_leakage(d, sv, stabilizers), argnums=0)

    rng = np.random.default_rng(seed)
    init_dir = rng.normal(size=n_qubits)
    init_dir /= np.linalg.norm(init_dir)
    init_norm = min(epsilon, 1e-2)
    delta = jnp.array(init_dir * init_norm)

    best_delta = delta
    best_leakage = float(x_stabilizer_leakage(delta, sv0, stabilizers))
    history = [best_leakage]

    for _ in range(n_steps):
        grad = leakage_grad(delta, sv0)
        grad_norm = jnp.linalg.norm(grad)
        step = jnp.where(grad_norm > 1e-12, grad / grad_norm, jnp.zeros_like(grad))
        delta = delta + step_size * step
        delta_norm = jnp.linalg.norm(delta)
        delta = jnp.where(delta_norm > epsilon, delta / delta_norm * epsilon, delta)

        current = float(x_stabilizer_leakage(delta, sv0, stabilizers))
        history.append(current)
        if current > best_leakage:
            best_leakage = current
            best_delta = delta

    return np.asarray(best_delta), best_leakage, history

project_l2_linf

project_l2_linf(
    y: ndarray,
    epsilon: float,
    linf_cap: float,
    n_bisect: int = 60,
) -> np.ndarray

Exact projection of y onto the intersection of an L2 ball of radius epsilon and an L-infinity ball (box) of radius linf_cap -- not the same as clip-then-rescale, which can push coordinates back outside the box. Box-clips first; if that's already inside the L2 ball it's the exact answer, otherwise bisects the Lagrange multiplier on the L2 constraint until the clipped, rescaled point lands exactly on the L2 boundary.

Examples:

>>> import numpy as np
>>> project_l2_linf(np.array([10.0, 0.0]), epsilon=1.0, linf_cap=5.0).round(4)
array([1., 0.])
>>> project_l2_linf(np.array([0.1, 0.1]), epsilon=1.0, linf_cap=0.05).round(4)
array([0.05, 0.05])

Source code in dense_evolution/noise/coherent_attack.py

def project_l2_linf(y: np.ndarray, epsilon: float, linf_cap: float, n_bisect: int = 60) -> np.ndarray:
    """Exact projection of `y` onto the intersection of an L2 ball of
    radius `epsilon` and an L-infinity ball (box) of radius `linf_cap` --
    not the same as clip-then-rescale, which can push coordinates back
    outside the box. Box-clips first; if that's already inside the L2
    ball it's the exact answer, otherwise bisects the Lagrange
    multiplier on the L2 constraint until the clipped, rescaled point
    lands exactly on the L2 boundary.

    Examples
    --------
    >>> import numpy as np
    >>> project_l2_linf(np.array([10.0, 0.0]), epsilon=1.0, linf_cap=5.0).round(4)
    array([1., 0.])
    >>> project_l2_linf(np.array([0.1, 0.1]), epsilon=1.0, linf_cap=0.05).round(4)
    array([0.05, 0.05])
    """
    y = np.asarray(y, dtype=np.float64)
    box = np.clip(y, -linf_cap, linf_cap)
    if np.linalg.norm(box) <= epsilon + 1e-12:
        return box
    lo, hi = 0.0, 1e8
    for _ in range(n_bisect):
        mid = 0.5 * (lo + hi)
        z = np.clip(y / (1.0 + mid), -linf_cap, linf_cap)
        if np.linalg.norm(z) > epsilon:
            lo = mid
        else:
            hi = mid
    return np.clip(y / (1.0 + hi), -linf_cap, linf_cap)

craft_adversarial_delta_constrained

craft_adversarial_delta_constrained(
    sv0: ndarray,
    stabilizers,
    epsilon: float,
    linf_cap: float,
    n_steps: int = 150,
    step_size: float = 0.05,
    seed: int = 0,
)

Same PGD search as craft_adversarial_delta, but each step is projected into the L2-epsilon-ball INTERSECTED with an L-infinity box of radius linf_cap (via project_l2_linf) instead of the L2 ball alone. Capping the per-qubit angle forbids the degenerate one-qubit-takes-everything solution and forces the search to spread the budget across multiple qubits -- the regime that actually causes real decoder failures (see module docstring).

Source code in dense_evolution/noise/coherent_attack.py

def craft_adversarial_delta_constrained(sv0: jnp.ndarray, stabilizers, epsilon: float, linf_cap: float,
                                         n_steps: int = 150, step_size: float = 0.05, seed: int = 0):
    """Same PGD search as `craft_adversarial_delta`, but each step is
    projected into the L2-epsilon-ball INTERSECTED with an L-infinity
    box of radius `linf_cap` (via `project_l2_linf`) instead of the L2
    ball alone. Capping the per-qubit angle forbids the degenerate
    one-qubit-takes-everything solution and forces the search to spread
    the budget across multiple qubits -- the regime that actually causes
    real decoder failures (see module docstring)."""
    n_qubits = len(stabilizers[0])
    leakage = lambda d, sv: x_stabilizer_leakage(d, sv, stabilizers)
    leakage_grad = jax.grad(leakage, argnums=0)

    rng = np.random.default_rng(seed)
    init_dir = rng.normal(size=n_qubits)
    init_dir /= np.linalg.norm(init_dir)
    init_norm = min(epsilon, 1e-2)
    delta_np = project_l2_linf(init_dir * init_norm, epsilon, linf_cap)
    delta = jnp.array(delta_np)

    best_delta = delta
    best_leakage = float(leakage(delta, sv0))
    history = [best_leakage]

    for _ in range(n_steps):
        grad = leakage_grad(delta, sv0)
        grad_norm = jnp.linalg.norm(grad)
        step = jnp.where(grad_norm > 1e-12, grad / grad_norm, jnp.zeros_like(grad))
        delta_raw = np.asarray(delta) + step_size * np.asarray(step)
        delta_np = project_l2_linf(delta_raw, epsilon, linf_cap)
        delta = jnp.array(delta_np)

        current = float(leakage(delta, sv0))
        history.append(current)
        if current > best_leakage:
            best_leakage = current
            best_delta = delta

    return np.asarray(best_delta), best_leakage, history

decoder_failure_rate

decoder_failure_rate(
    delta_np: ndarray,
    sv0_np: ndarray,
    decode_fn,
    n_trials: int,
    rng: Generator,
) -> float

Real, discrete evaluation of how often a coherent error delta_np actually fools a decoder -- as opposed to x_stabilizer_leakage's smooth proxy. decode_fn(sv_noisy, rng) -> sv_corrected is any projective-measurement-based decoder (e.g. built around dense_evolution.qec.compute_syndrome); this function applies the coherent error once, then calls decode_fn n_trials times (the syndrome measurement that collapses the coherently-perturbed state is itself stochastic, so repeated trials are meaningful even though delta_np is fixed), and reports the fraction that fail to recover sv0_np exactly.

Source code in dense_evolution/noise/coherent_attack.py

def decoder_failure_rate(delta_np: np.ndarray, sv0_np: np.ndarray, decode_fn,
                          n_trials: int, rng: np.random.Generator) -> float:
    """Real, discrete evaluation of how often a coherent error `delta_np`
    actually fools a decoder -- as opposed to `x_stabilizer_leakage`'s
    smooth proxy. `decode_fn(sv_noisy, rng) -> sv_corrected` is any
    projective-measurement-based decoder (e.g. built around
    `dense_evolution.qec.compute_syndrome`); this function applies the
    coherent error once, then calls `decode_fn` `n_trials` times (the
    syndrome measurement that collapses the coherently-perturbed state
    is itself stochastic, so repeated trials are meaningful even though
    `delta_np` is fixed), and reports the fraction that fail to recover
    `sv0_np` exactly."""
    sv_delta = np.asarray(apply_rz_all(jnp.array(sv0_np), jnp.array(delta_np)))
    n_fail = 0
    for _ in range(n_trials):
        sv_corrected = decode_fn(sv_delta.copy(), rng)
        fidelity = np.abs(np.vdot(sv_corrected, sv0_np)) ** 2
        if fidelity < 1.0 - 1e-6:
            n_fail += 1
    return n_fail / n_trials

random_delta_failure_stats

random_delta_failure_stats(
    sv0_np: ndarray,
    decode_fn,
    epsilon: float,
    n_qubits: int,
    n_random: int,
    n_trials_each: int,
    rng: Generator,
) -> np.ndarray

Same evaluation as decoder_failure_rate, but over n_random random directions of L2 norm epsilon instead of one crafted delta -- the baseline craft_adversarial_delta's result should be compared against (see module docstring: the unconstrained crafted direction can score WORSE than this random baseline).

Source code in dense_evolution/noise/coherent_attack.py

def random_delta_failure_stats(sv0_np: np.ndarray, decode_fn, epsilon: float, n_qubits: int,
                                n_random: int, n_trials_each: int, rng: np.random.Generator) -> np.ndarray:
    """Same evaluation as `decoder_failure_rate`, but over `n_random`
    random directions of L2 norm `epsilon` instead of one crafted
    `delta` -- the baseline `craft_adversarial_delta`'s result should be
    compared against (see module docstring: the unconstrained crafted
    direction can score WORSE than this random baseline)."""
    rates = np.zeros(n_random)
    for i in range(n_random):
        d = rng.normal(size=n_qubits)
        d = d / np.linalg.norm(d) * epsilon
        rates[i] = decoder_failure_rate(d, sv0_np, decode_fn, n_trials_each, rng)
    return rates

global_depolarizing_channel

global_depolarizing_channel(
    rho: ndarray, p: float
) -> jnp.ndarray

Global n-qubit depolarizing channel, D_p(rho) = (1-p)rho + (p/dim)I.

Distinct from NoiseModel's 'depolarizing' model, which applies an independent PER-QUBIT local Kraus channel -- a different physical map from this GLOBAL channel, which mixes the whole dim-dimensional state toward the fully mixed state as one unit. Use this one when modeling e.g. state-prep/measurement (SPAM) error reported as a single joint depolarizing parameter over the whole register, not per-qubit gate noise (promoted from a real reproduction of arXiv:2608.16716's own SPAM model, Dense-Evolution-Discovery Experiment 33).

Source code in dense_evolution/noise/density_matrix_channels.py

def global_depolarizing_channel(rho: jnp.ndarray, p: float) -> jnp.ndarray:
    """Global n-qubit depolarizing channel, D_p(rho) = (1-p)*rho + (p/dim)*I.

    Distinct from `NoiseModel`'s `'depolarizing'` model, which applies an
    independent PER-QUBIT local Kraus channel -- a different physical map
    from this GLOBAL channel, which mixes the whole `dim`-dimensional state
    toward the fully mixed state as one unit. Use this one when modeling
    e.g. state-prep/measurement (SPAM) error reported as a single joint
    depolarizing parameter over the whole register, not per-qubit gate
    noise (promoted from a real reproduction of arXiv:2608.16716's own
    SPAM model, Dense-Evolution-Discovery Experiment 33).
    """
    rho = jnp.asarray(rho, dtype=jnp.complex128)
    dim = rho.shape[0]
    identity = jnp.eye(dim, dtype=jnp.complex128)
    return (1.0 - p) * rho + (p / dim) * identity

amplitude_damping_channel

amplitude_damping_channel(
    rho: ndarray, gamma: float
) -> jnp.ndarray

Single-qubit amplitude-damping channel: E0 @ rho @ E0.conj().T + E1 @ rho @ E1.conj().T, with E0=diag(1, sqrt(1-gamma)) and E1=[[0,sqrt(gamma)],[0,0]] -- population only ever moves |1>->|0>, never the reverse.

Distinct from global_depolarizing_channel (symmetric, mixes toward the fully-mixed state regardless of which state is |1> or |0>) -- this one is asymmetric by construction, the real signature of energy-relaxation (T1) processes and of quasiparticle poisoning (promoted from a real reproduction of arXiv:2104.05219's measured cosmic-ray-induced error bursts, Dense-Evolution-Discovery Experiment 34, where this asymmetry is exactly the mechanism's own reported signature: decay errors only, no excess excitation errors).

Single-qubit only (rho must be 2x2) -- unlike global_depolarizing_channel, this is not dimension-generic, since amplitude damping is inherently a per-qubit process, not a joint-register one.

Source code in dense_evolution/noise/density_matrix_channels.py

def amplitude_damping_channel(rho: jnp.ndarray, gamma: float) -> jnp.ndarray:
    """Single-qubit amplitude-damping channel: E0 @ rho @ E0.conj().T +
    E1 @ rho @ E1.conj().T, with E0=diag(1, sqrt(1-gamma)) and
    E1=[[0,sqrt(gamma)],[0,0]] -- population only ever moves |1>->|0>,
    never the reverse.

    Distinct from `global_depolarizing_channel` (symmetric, mixes toward
    the fully-mixed state regardless of which state is |1> or |0>) -- this
    one is asymmetric by construction, the real signature of energy-relaxation
    (T1) processes and of quasiparticle poisoning (promoted from a real
    reproduction of arXiv:2104.05219's measured cosmic-ray-induced error
    bursts, Dense-Evolution-Discovery Experiment 34, where this asymmetry is
    exactly the mechanism's own reported signature: decay errors only, no
    excess excitation errors).

    Single-qubit only (rho must be 2x2) -- unlike `global_depolarizing_channel`,
    this is not dimension-generic, since amplitude damping is inherently a
    per-qubit process, not a joint-register one.
    """
    rho = jnp.asarray(rho, dtype=jnp.complex128)
    e0 = jnp.array([[1.0, 0.0], [0.0, jnp.sqrt(1.0 - gamma)]], dtype=jnp.complex128)
    e1 = jnp.array([[0.0, jnp.sqrt(gamma)], [0.0, 0.0]], dtype=jnp.complex128)
    return e0 @ rho @ e0.conj().T + e1 @ rho @ e1.conj().T

cosmic_ray_burst_profile

cosmic_ray_burst_profile(
    time_us,
    baseline_gamma: float,
    ratio_intermediate: float = 2.5,
    ratio_peak: float = 3.75,
    tau1_us: float = 3.0,
    tau2_us: float = 300.0,
    tau_decay_ms: float = 25.0,
) -> jnp.ndarray

Time-dependent decay-probability profile for a cosmic-ray/gamma-ray- induced quasiparticle burst: a two-stage rise (fast to ratio_intermediatex baseline, slower to ratio_peakx baseline) times a single-exponential recovery, generalized out of a fixed, paper-number validation (Dense-Evolution-Discovery Experiment 34, reproducing arXiv:2104.05219's real measured event on a 26-qubit chip).

Feed the result to continuous_dissipative_evolve alongside amplitude_damping_channel (or any other single-time-varying-parameter channel) to inject a realistic burst into any circuit or QEC study, without re-deriving this shape by hand each time.

The default ratios/timescales are the paper's own real numbers -- see Experiment 34's docstring for exactly which are paper-fitted (the 25ms decay) versus chosen to match the paper's two described rise points (tau1/tau2). All are overridable for a different event severity or device generation; baseline_gamma is never derived here -- pass whatever per-slice decay probability corresponds to your own dt/T1 convention (see Experiment 34 for one worked example of that conversion).

Parameters:

Name	Type	Description	Default
time_us	array_like	Time since impact, in microseconds (t=0 is the impact instant).	required
baseline_gamma	float	Undisturbed per-slice decay probability; this profile scales it up, it does not derive it.	required
ratio_intermediate	float	Multiplier on baseline_gamma at the two described checkpoints (paper defaults 2.5=10/4, 3.75=15/4, from Fig. 3's ~10us/~1ms readings).	2.5
ratio_peak	float	Multiplier on baseline_gamma at the two described checkpoints (paper defaults 2.5=10/4, 3.75=15/4, from Fig. 3's ~10us/~1ms readings).	2.5
tau1_us	float	Rise timescales for the two saturating-exponential stages (paper defaults 3, 300 -- chosen to match its ~10us/~1ms descriptions, not fitted by the paper itself).	3.0
tau2_us	float	Rise timescales for the two saturating-exponential stages (paper defaults 3, 300 -- chosen to match its ~10us/~1ms descriptions, not fitted by the paper itself).	3.0
tau_decay_ms	float	Recovery time constant (paper default 25 -- its own fitted central value, real range 25-30ms across 415 events).	25.0

Returns:

Type	Description
ndarray	Per-slice decay probability at each entry of time_us.

Source code in dense_evolution/noise/cosmic_ray.py

def cosmic_ray_burst_profile(time_us, baseline_gamma: float, ratio_intermediate: float = 2.5,
                              ratio_peak: float = 3.75, tau1_us: float = 3.0,
                              tau2_us: float = 300.0, tau_decay_ms: float = 25.0) -> jnp.ndarray:
    """Time-dependent decay-probability profile for a cosmic-ray/gamma-ray-
    induced quasiparticle burst: a two-stage rise (fast to
    `ratio_intermediate`x baseline, slower to `ratio_peak`x baseline) times
    a single-exponential recovery, generalized out of a fixed, paper-number
    validation (Dense-Evolution-Discovery Experiment 34, reproducing
    arXiv:2104.05219's real measured event on a 26-qubit chip).

    Feed the result to `continuous_dissipative_evolve` alongside
    `amplitude_damping_channel` (or any other single-time-varying-parameter
    channel) to inject a realistic burst into any circuit or QEC study,
    without re-deriving this shape by hand each time.

    The default ratios/timescales are the paper's own real numbers -- see
    Experiment 34's docstring for exactly which are paper-fitted (the 25ms
    decay) versus chosen to match the paper's two described rise points
    (tau1/tau2). All are overridable for a different event severity or
    device generation; `baseline_gamma` is never derived here -- pass
    whatever per-slice decay probability corresponds to your own dt/T1
    convention (see Experiment 34 for one worked example of that
    conversion).

    Parameters
    ----------
    time_us : array_like
        Time since impact, in microseconds (t=0 is the impact instant).
    baseline_gamma : float
        Undisturbed per-slice decay probability; this profile scales it
        up, it does not derive it.
    ratio_intermediate, ratio_peak : float
        Multiplier on `baseline_gamma` at the two described checkpoints
        (paper defaults 2.5=10/4, 3.75=15/4, from Fig. 3's ~10us/~1ms
        readings).
    tau1_us, tau2_us : float
        Rise timescales for the two saturating-exponential stages (paper
        defaults 3, 300 -- chosen to match its ~10us/~1ms descriptions,
        not fitted by the paper itself).
    tau_decay_ms : float
        Recovery time constant (paper default 25 -- its own fitted central
        value, real range 25-30ms across 415 events).

    Returns
    -------
    jnp.ndarray
        Per-slice decay probability at each entry of `time_us`.
    """
    time_us = jnp.asarray(time_us)
    stage1 = (ratio_intermediate - 1.0) * (1.0 - jnp.exp(-time_us / tau1_us))
    stage2 = (ratio_peak - ratio_intermediate) * (1.0 - jnp.exp(-time_us / tau2_us))
    decay = jnp.exp(-time_us / (tau_decay_ms * 1000.0))
    scaling = 1.0 + (stage1 + stage2) * decay
    return baseline_gamma * scaling

oscillating_p_eff

oscillating_p_eff(
    base_p: float, factor: float, freq: float, amp: float
) -> jnp.ndarray

Effective noise probability that oscillates around base_p as a function of factor (e.g. a ZNE noise-scale factor), instead of scaling smoothly with it: base_p * (1 + amp * sin(factor * pi / freq)), clipped to [0.01, 0.5] so it always stays a valid probability.

Promoted from Dense-Evolution-Discovery's jsd_zne_oscillating_noise.py, where it was used to build a noise-vs-scale relationship deliberately NOT smooth/monotonic, to stress-test dense_evolution.mitigation.jsd_predictive_zne_density_matrix against noise models where plain Richardson extrapolation's smoothness assumption breaks down.

Examples:

>>> from dense_evolution.noise import oscillating_p_eff
>>> round(float(oscillating_p_eff(base_p=0.1, factor=0.0, freq=2.0, amp=0.5)), 4)
0.1
>>> round(float(oscillating_p_eff(base_p=0.1, factor=1.0, freq=2.0, amp=0.5)), 4)
0.15

Source code in dense_evolution/noise/oscillating.py

def oscillating_p_eff(base_p: float, factor: float, freq: float, amp: float) -> jnp.ndarray:
    """Effective noise probability that oscillates around `base_p` as a
    function of `factor` (e.g. a ZNE noise-scale factor), instead of
    scaling smoothly with it: `base_p * (1 + amp * sin(factor * pi /
    freq))`, clipped to `[0.01, 0.5]` so it always stays a valid
    probability.

    Promoted from Dense-Evolution-Discovery's jsd_zne_oscillating_noise.py,
    where it was used to build a noise-vs-scale relationship deliberately
    NOT smooth/monotonic, to stress-test
    `dense_evolution.mitigation.jsd_predictive_zne_density_matrix` against
    noise models where plain Richardson extrapolation's smoothness
    assumption breaks down.

    Examples
    --------
    >>> from dense_evolution.noise import oscillating_p_eff
    >>> round(float(oscillating_p_eff(base_p=0.1, factor=0.0, freq=2.0, amp=0.5)), 4)
    0.1
    >>> round(float(oscillating_p_eff(base_p=0.1, factor=1.0, freq=2.0, amp=0.5)), 4)
    0.15
    """
    p = base_p * (1.0 + amp * jnp.sin(factor * jnp.pi / freq))
    return jnp.clip(p, 0.01, 0.5)

See Also

• DenseSVSimulator, QASMParser — build the circuit and statevector every step above starts from.
• dense_evolution.mitigation — correcting noise after it's applied, instead of just simulating it.
• dense_evolution.qec — the decoder side of the code craft_adversarial_delta_constrained attacks.