diff --git a/pyproject.toml b/pyproject.toml index 5c90107b3..023af008e 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,6 +1,6 @@ [tool.poetry] name = "pyquil" -version = "4.20.0" +version = "4.21.0" description = "A Python library for creating Quantum Instruction Language (Quil) programs." authors = ["Rigetti Computing "] readme = "README.md" diff --git a/pyquil/noise/__init__.py b/pyquil/noise/__init__.py new file mode 100644 index 000000000..bace195bf --- /dev/null +++ b/pyquil/noise/__init__.py @@ -0,0 +1,76 @@ +"""pyquil.noise — Noise modeling for quantum simulators. + +This package provides: + +- **Noise model** (``_legacy_noise``): Kraus-map based noise construction + for the QVM, including ``KrausModel``, ``NoiseModel``, and decoherence + helpers. + +- **Channel classes** (``_channels``): quax-backed ``Channel``, + ``MeasurementChannel``, ``ResetChannel``, and ``CycleChannel`` dataclasses + for fine-grained noise modeling. These are private and not re-exported. + +- **New noise model** (``_noise_model``): The quax-based ``NoiseModel`` + container and program-level fidelity estimation utilities. These are + private and not re-exported; they will become the public API in the next + major version. +""" + +# ── Noise model (Kraus-map based) ─────────────────────────────────────── +from pyquil.noise._legacy_noise import ( + ANGLE_TOLERANCE, + INFINITY, + NO_NOISE, + KrausModel, + NoiseModel, + NoisyGateUndefined, + _bitstring_probs_by_qubit, # noqa: F401 + _check_kraus_ops, # noqa: F401 + _create_kraus_pragmas, # noqa: F401 + _decoherence_noise_model, # noqa: F401 + _get_program_gates, # noqa: F401 + _noise_model_program_header, # noqa: F401 + _run, # noqa: F401 + add_decoherence_noise, + append_kraus_to_gate, + apply_noise_model, + bitstring_probs_to_z_moments, + combine_kraus_maps, + correct_bitstring_probs, + corrupt_bitstring_probs, + damping_after_dephasing, + damping_kraus_map, + decoherence_noise_with_asymmetric_ro, + dephasing_kraus_map, + estimate_assignment_probs, + estimate_bitstring_probs, + get_noisy_gate, + pauli_kraus_map, + tensor_kraus_maps, +) + +__all__ = [ + # Noise model + "ANGLE_TOLERANCE", + "INFINITY", + "KrausModel", + "NO_NOISE", + "NoiseModel", + "NoisyGateUndefined", + "add_decoherence_noise", + "append_kraus_to_gate", + "apply_noise_model", + "bitstring_probs_to_z_moments", + "combine_kraus_maps", + "correct_bitstring_probs", + "corrupt_bitstring_probs", + "damping_after_dephasing", + "damping_kraus_map", + "decoherence_noise_with_asymmetric_ro", + "dephasing_kraus_map", + "estimate_assignment_probs", + "estimate_bitstring_probs", + "get_noisy_gate", + "pauli_kraus_map", + "tensor_kraus_maps", +] diff --git a/pyquil/noise/_channels.py b/pyquil/noise/_channels.py new file mode 100644 index 000000000..597057c67 --- /dev/null +++ b/pyquil/noise/_channels.py @@ -0,0 +1,1654 @@ +############################################################################## +# Copyright 2016-2026 Rigetti Computing +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +############################################################################## +"""Noise channel classes and gate-resolution utilities. + +This module defines ``Channel``, ``MeasurementChannel``, ``ResetChannel``, and +``CycleChannel`` dataclasses for representing noise in quantum circuits, along +with helper functions for resolving gate unitaries and extracting custom gate +definitions from Quil programs. +""" + +from __future__ import annotations + +import itertools +import json +import logging +from collections.abc import Callable +from dataclasses import dataclass, replace +from functools import cached_property, reduce +from itertools import product +from typing import TYPE_CHECKING, Any + +import jax.numpy as jnp +import numpy as np +import quax as qx +from jax import Array +from quil.expression import Expression as QuilExpression +from quil.program import Program as RSProgram +from scipy.linalg import expm as scipy_expm +from scipy.linalg import fractional_matrix_power +from scipy.linalg import logm as scipy_logm + +from pyquil.quilatom import Expression, FormalArgument, Parameter, substitute +from pyquil.quilbase import DefCircuit, DefGate, Gate, Measurement, Reset, ResetQubit + +if TYPE_CHECKING: + from plotly.graph_objs import Figure + + from pyquil import Program + +logger = logging.getLogger(__name__) + +# Type alias for the custom-gate lookup map used throughout the Channel constructors. +CustomGateMap = dict[str, qx.Unitary | Callable[..., qx.Unitary]] + + +def _parse_quil_instruction(quil_str: str) -> Gate | Measurement | Reset: + """Parse a single Quil instruction string into a pyquil instruction object. + + Uses the ``quil`` Rust parser directly, avoiding a dependency on ``pyquil.Program``. + """ + rs_inst = RSProgram.parse(quil_str).body_instructions[0] + if rs_inst.is_gate(): + return Gate._from_rs_gate(rs_inst.to_gate()) + elif rs_inst.is_measurement(): + return Measurement._from_rs_measurement(rs_inst.to_measurement()) + elif rs_inst.is_reset(): + reset = rs_inst.to_reset() + if reset.qubit is None: + return Reset._from_rs_reset(reset) + return ResetQubit._from_rs_reset(reset) + raise ValueError(f"Unsupported instruction type in: {quil_str}") + + +def _pack_complex_array(array: Array | np.ndarray) -> dict[str, Any]: + """Pack a complex array into JSON-compatible real/imaginary pairs.""" + np_array = np.asarray(array) + return { + "_complex_array": [[float(value.real), float(value.imag)] for value in np_array.flat], + "shape": list(np_array.shape), + } + + +def _unpack_complex_array(data: dict[str, Any]) -> Array: + """Unpack a complex array from :func:`_pack_complex_array` data.""" + shape = tuple(data["shape"]) + return jnp.array([complex(pair[0], pair[1]) for pair in data["_complex_array"]], dtype=complex).reshape(shape) + + +def _pack_dims(dims: tuple[tuple[int, ...], tuple[int, ...]]) -> list[list[int]]: + """Pack quax operator dims into JSON-compatible lists.""" + return [list(dims[0]), list(dims[1])] + + +def _unpack_dims(data: list[list[int]]) -> tuple[tuple[int, ...], tuple[int, ...]]: + """Unpack quax operator dims from JSON-compatible lists.""" + if len(data) != 2: + raise ValueError(f"Serialized operator dims must contain output and input dims, got {data}.") + return (tuple(int(dim) for dim in data[0]), tuple(int(dim) for dim in data[1])) + + +def _infer_legacy_qubit_dims(shape: tuple[int, ...], *, superoperator: bool) -> tuple[tuple[int, ...], tuple[int, ...]]: + """Infer qubit-only dims for data serialized before dims were stored.""" + hilbert_dim = int(round(np.sqrt(shape[0]))) if superoperator else int(shape[0]) + num_qubits = int(round(np.log2(hilbert_dim))) + if 2**num_qubits != hilbert_dim: + raise ValueError( + "Serialized operator data does not include dims and its shape is not compatible with qubit-only dims." + ) + return ((2,) * num_qubits, (2,) * num_qubits) + + +def _pack_operator(operator: qx.SuperOp | qx.Unitary | qx.Choi) -> dict[str, Any]: + """Pack a quax operator matrix with explicit dimension metadata.""" + data = _pack_complex_array(operator.matrix) + data["dims"] = _pack_dims(operator.dims) + return data + + +def _unpack_operator_dims( + data: dict[str, Any], shape: tuple[int, ...], *, superoperator: bool +) -> tuple[tuple[int, ...], tuple[int, ...]]: + """Read explicit dims, falling back to the legacy qubit-only inference.""" + if "dims" in data: + return _unpack_dims(data["dims"]) + return _infer_legacy_qubit_dims(shape, superoperator=superoperator) + + +def _resolve_params(params: list) -> list[float]: + """Resolve gate parameters to concrete float values. + + :param params: The gate parameters (may include symbolic Parameters or Expressions). + :return: A list of concrete float values. + :raises ValueError: If any parameter is symbolic and cannot be evaluated to a number. + """ + fixed_params = [] + for p in params: + if isinstance(p, (Parameter, Expression)): + evaluated = p._evaluate() + if isinstance(evaluated, (Parameter, Expression)): + raise ValueError( + f"Cannot resolve symbolic parameter {p}. Provide a gate with concrete numeric parameters." + ) + fixed_params.append(float(evaluated)) + elif isinstance(p, QuilExpression): + result = p.evaluate({}, {}) + fixed_params.append(float(result.to_number()) if hasattr(result, "to_number") else float(result)) # type: ignore[arg-type] + else: + fixed_params.append(float(p.real)) + return fixed_params + + +def get_custom_gates_from_program(program: Program) -> CustomGateMap: + """Extract custom gate definitions from a Quil program. + + Returns a dictionary mapping gate names to unitary matrices (for fixed gates) or callables + (for parametric gates). Does not include the standard gate set — use this to augment + the standard ``qx.gates.QUANTUM_GATES`` when resolving instructions with custom gates. + + :param program: A Quil program containing DefGate definitions. + :return: A dictionary of custom gate names to unitary matrices or callables. + """ + custom_gates: CustomGateMap = {} + for defgate in program.defined_gates: + if defgate.parameters: + + def parametric_gate(*args: float, defgate: DefGate = defgate) -> qx.Unitary: + parameter_map = {Parameter(p.name): arg for p, arg in zip(defgate.parameters, args, strict=False)} + matrix = jnp.asarray( + [[substitute(element, parameter_map) for element in row] for row in defgate.matrix], # type: ignore[arg-type] + dtype=complex, + ) + num_qubits = int(jnp.round(jnp.log2(matrix.shape[0]))) + return qx.Unitary.from_matrix(matrix, ((2,) * num_qubits, (2,) * num_qubits)) + + custom_gates[defgate.name] = parametric_gate + else: + matrix = jnp.asarray(defgate.matrix, dtype=complex) + num_qubits = int(jnp.round(jnp.log2(matrix.shape[0]))) + custom_gates[defgate.name] = qx.Unitary.from_matrix(matrix, ((2,) * num_qubits, (2,) * num_qubits)) + return custom_gates + + +def get_instruction_unitary( + inst: Gate, + custom_gates: CustomGateMap | None = None, +) -> qx.Unitary: + """Get the unitary matrix associated with a gate instruction. + + Looks up the gate by name — first in ``custom_gates`` (if provided), then in the + standard quax gate table ``qx.gates.QUANTUM_GATES``. Parametric gates are supported + provided all parameters are concrete numeric values. + + :param inst: The gate instruction. + :param custom_gates: Optional dictionary of additional gate definitions (e.g. from + :func:`get_custom_gates_from_program`). Takes precedence over the standard gate set. + :return: The unitary matrix. + :raises ValueError: If any gate parameter is symbolic. + :raises KeyError: If the gate name is not found in either the custom or standard gate set. + """ + name = inst.name + + # Look up gate definition: custom gates take precedence + if custom_gates is not None and name in custom_gates: + gate_def = custom_gates[name] + elif name in qx.gates.QUANTUM_GATES: + gate_def = qx.gates.QUANTUM_GATES[name] + else: + raise KeyError(f"Unknown gate '{name}'. Provide it via custom_gates (e.g. custom_gates={{'{name}': matrix}}).") + + if inst.params: + fixed_params = _resolve_params(list(inst.params)) + if not callable(gate_def): + raise ValueError(f"Gate '{name}' is not parametric but parameters were provided.") + result = gate_def(*fixed_params) + else: + if callable(gate_def): + result = gate_def() + else: + result = gate_def + + # quax parametric gates may return Operator instead of Unitary; wrap if needed + if not isinstance(result, qx.Unitary): + result = qx.Unitary.from_matrix(result.matrix, result.dims) + return result + + +@dataclass(frozen=True) +class Channel: + """A noise channel attaches a superoperator to a specific gate. + + The superoperator *includes* the gate unitary, so the channel replaces the gate + rather than being applied after it. + + The ``process`` field is a ``qx.SuperOp`` which can be converted to alternative + representations (Choi, Kraus, Pauli-Liouville) via ``quax``. + + Fidelity metrics are computed relative to the ideal gate unitary stored in + ``target_unitary``. For standard gates use the class methods (e.g. + :meth:`from_gate_fidelity`) which resolve the unitary automatically. + """ + + inst: Gate + """Quil gate to which the channel applies.""" + + process: qx.SuperOp + """The noisy process (superoperator) for the gate, including the gate unitary.""" + + target_unitary: qx.Unitary + """The noiseless unitary of the gate.""" + + @cached_property + def unitary(self) -> qx.Unitary: + """The noiseless unitary of the gate.""" + return self.target_unitary + + @cached_property + def qubits(self) -> list[int]: + """The qubits which the channel applies to.""" + return self.inst.get_qubit_indices() + + @cached_property + def num_qubits(self) -> int: + """The number of qubits the channel acts on.""" + return len(self.qubits) + + # ────────────────────────────────────────────── + # Constructors + # ────────────────────────────────────────────── + + @classmethod + def from_gate_fidelity( + cls: type[Channel], + inst: Gate, + fidelity: float, + custom_gates: CustomGateMap | None = None, + ) -> Channel: + r"""Create a depolarizing noise channel from an average gate fidelity. + + The resulting channel is the composition of the ideal gate unitary with a + depolarizing channel calibrated to the specified fidelity: + :math:`\\mathcal{E} = \\mathcal{D}_p \\circ \\mathcal{U}` + + :param inst: The gate to which the channel applies. + :param fidelity: The average gate fidelity, :math:`F_{\\mathrm{avg}} \\in [0, 1]`. + :param custom_gates: Optional dictionary of custom gate definitions. + :return: A Channel instance. + """ + unitary = get_instruction_unitary(inst, custom_gates) + p = qx.average_fidelity_to_depolarizing_constant(fidelity, unitary.dims[0]) + return cls.from_depolarizing_constant(inst, p, custom_gates) + + @classmethod + def from_pauli_fidelity( + cls: type[Channel], + inst: Gate, + pauli_fidelity: float, + custom_gates: CustomGateMap | None = None, + ) -> Channel: + r"""Create a depolarizing noise channel from a process (Pauli) fidelity. + + The process fidelity :math:`F_e` is related to the average gate fidelity by + :math:`F_{\\mathrm{avg}} = (d \\cdot F_e + 1) / (d + 1)`. + + :param inst: The gate to which the channel applies. + :param pauli_fidelity: The process fidelity (entanglement fidelity), :math:`F_e \\in [0, 1]`. + :param custom_gates: Optional dictionary of custom gate definitions. + :return: A Channel instance. + """ + unitary = get_instruction_unitary(inst, custom_gates) + p = qx.process_fidelity_to_depolarizing_constant(pauli_fidelity, unitary.dims[0]) + return cls.from_depolarizing_constant(inst, p, custom_gates) + + @classmethod + def from_depolarizing_constant( + cls: type[Channel], + inst: Gate, + depolarizing_constant: float, + custom_gates: CustomGateMap | None = None, + ) -> Channel: + r"""Create a depolarizing noise channel from a depolarization constant. + + The depolarizing constant :math:`p` parameterizes the channel as + :math:`\\mathcal{D}_p(\\rho) = p \\, \\rho + (1-p) \\, I/d`. + + :param inst: The gate to which the channel applies. + :param depolarizing_constant: The depolarization constant, e.g. 0.98 for 2% depolarization. + :param custom_gates: Optional dictionary of custom gate definitions. + :return: A Channel instance. + """ + unitary = get_instruction_unitary(inst, custom_gates) + depolarizing_superop = qx.depolarizing_channel_superoperator(1 - depolarizing_constant, unitary.dims[0]) + combined_superop = depolarizing_superop @ unitary + return cls(inst=inst, process=qx.to_superop(combined_superop), target_unitary=unitary) + + @classmethod + def from_pauli_noise( + cls: type[Channel], + inst: Gate, + pauli_noise: dict[str, float], + custom_gates: CustomGateMap | None = None, + ) -> Channel: + """Create a stochastic Pauli noise channel from Pauli error rates. + + The noise is specified as a dictionary mapping Pauli strings to error probabilities, + e.g. ``{"XX": 0.03, "ZI": 0.001}``. The probabilities must sum to at most 1.0; + any remainder is assigned to the identity (no-error) term. + + :param inst: The gate to which the channel applies. + :param pauli_noise: Pauli error rates, e.g. ``{"IX": 0.01, "ZZ": 0.02}``. + :param custom_gates: Optional dictionary of custom gate definitions. + :return: A Channel instance. + """ + unitary = get_instruction_unitary(inst, custom_gates) + num_qubits = len(unitary.dims[0]) + + total_error_rate = 0.0 + for pauli, error_rate in pauli_noise.items(): + if error_rate < 0.0: + raise ValueError(f"Pauli term '{pauli}' has negative error rate {error_rate}.") + total_error_rate += error_rate + if total_error_rate > 1.0: + raise ValueError(f"Pauli error rates must sum to at most 1.0, got {total_error_rate}.") + + for pauli in pauli_noise: + if len(pauli) != num_qubits: + raise ValueError(f"Pauli term '{pauli}' has length {len(pauli)}, expected {num_qubits}.") + + all_pauli_terms = tuple("".join(term) for term in product("IXYZ", repeat=num_qubits)) + + pauli_error_rates: list[float] = [] + for term in reversed(all_pauli_terms): + if term in pauli_noise: + error_rate = pauli_noise[term] + elif all(p == "I" for p in term): + error_rate = 1 - sum(pauli_error_rates) + else: + error_rate = 0 + pauli_error_rates.append(error_rate) + if not jnp.isclose(1.0, sum(pauli_error_rates)): + raise ValueError("Pauli error rates plus the implicit identity rate must sum to 1.0.") + pauli_error_rates = list(reversed(pauli_error_rates)) + + # Build the 4**num_qubits Pauli operators as tensor (Kronecker) products of the + # single-qubit Paulis, in lexicographic (I, X, Y, Z) order matching pauli_error_rates. + single_pauli_matrices = qx.ensembles.PAULIS.matrix # (4, 2, 2): I, X, Y, Z + pauli_op_matrices = jnp.stack( + [reduce(jnp.kron, paulis) for paulis in product(single_pauli_matrices, repeat=num_qubits)] + ) # (4**num_qubits, 2**num_qubits, 2**num_qubits) + + # Scale each Pauli by sqrt(probability) to form Kraus operators + coeffs = jnp.sqrt(jnp.array(pauli_error_rates, dtype=float)) + kraus_matrices = coeffs[:, None, None] * pauli_op_matrices + kraus_map = qx.KrausMap.from_matrix(kraus_matrices, unitary.dims) + + process_superop = qx.to_superop(kraus_map @ unitary) + return cls(inst=inst, process=process_superop, target_unitary=unitary) + + @classmethod + def from_random_coherent_error( + cls: type[Channel], + inst: Gate, + process_fidelity: float, + rng: np.random.Generator | None = None, + custom_gates: CustomGateMap | None = None, + ) -> Channel: + r"""Create a channel with a random coherent (unitary) error at the specified process fidelity. + + A random unitary close to identity is generated with the given process fidelity, + then composed with the ideal gate. + + :param inst: The gate to which the channel applies. + :param process_fidelity: The process fidelity of the coherent error, :math:`F_e \\in [0, 1]`. + :param rng: NumPy random number generator for reproducibility. + :param custom_gates: Optional dictionary of custom gate definitions. + :return: A Channel instance. + """ + if rng is None: + rng = np.random.default_rng() + + ideal = get_instruction_unitary(inst, custom_gates) + num_qubits = len(ideal.dims[0]) + d = 2**num_qubits + + # Generate a random unitary error with the specified process fidelity + # using Pauli generator decomposition + angle = jnp.arccos(2 * process_fidelity - 1) / (2 * jnp.pi) + id_coeff = 1 - float(angle) + coeffs = rng.random(4**num_qubits - 1) + coeffs = (1 - id_coeff) / np.sqrt(np.sum(np.square(coeffs))) * coeffs + + # Build Pauli generator sum using quax Pauli matrices + pauli_matrices = qx.ensembles.PAULIS.matrix # shape (4, 2, 2) + pauli_sum = jnp.eye(d, dtype=complex) * id_coeff + pauli_products = list(itertools.product(pauli_matrices, repeat=num_qubits))[1:] + for paulis, coefficient in zip(pauli_products, coeffs, strict=False): + pauli_sum = pauli_sum + reduce(jnp.kron, paulis) * coefficient + + from jax.scipy.linalg import expm as jax_expm + + error_unitary = jax_expm(-1j * jnp.pi * pauli_sum) + # Fix global phase + phase = jnp.exp(-1j * jnp.angle(error_unitary[0, 0])) + error_unitary = error_unitary * phase + + error_u = qx.Unitary.from_matrix(error_unitary, ideal.dims) + noisy_superop = qx.to_superop(error_u @ ideal) + return cls(inst=inst, process=noisy_superop, target_unitary=ideal) + + @classmethod + def from_mixture( + cls: type[Channel], + inst: Gate, + constituents: list[qx.Unitary], + probabilities: list[float], + custom_gates: CustomGateMap | None = None, + ) -> Channel: + r"""Create a mixture channel from a set of unitary errors with given probabilities. + + The channel is :math:`\\mathcal{E}(\\rho) = (1-\\sum p_i) U\\rho U^\\dagger + \\sum p_i V_i U \\rho U^\\dagger V_i^\\dagger` + where :math:`U` is the ideal gate and :math:`V_i` are the error unitaries. + + :param inst: The gate to which the channel applies. + :param constituents: Unitary error operators to mix. + :param probabilities: Probability of each unitary error. Must sum to at most 1.0. + :param custom_gates: Optional dictionary of custom gate definitions. + :return: A Channel instance. + """ + ideal = get_instruction_unitary(inst, custom_gates) + + if len(constituents) != len(probabilities): + raise ValueError("The number of constituents and probabilities must match.") + error_prob = sum(probabilities) + if error_prob > 1.0: + raise ValueError(f"The sum of probabilities ({error_prob}) must be at most 1.0.") + + # Build the mixture superop: (1-p_total) S(U) + sum p_i S(V_i @ U) + p0 = 1.0 - error_prob + noisy_superop_matrix = p0 * qx.to_superop(ideal).matrix + for p, v in zip(probabilities, constituents, strict=False): + composed = v @ ideal + noisy_superop_matrix = noisy_superop_matrix + p * qx.to_superop(composed).matrix + noisy_superop = qx.SuperOp.from_matrix(noisy_superop_matrix, ideal.dims) + return cls(inst=inst, process=noisy_superop, target_unitary=ideal) + + @classmethod + def from_coherence_times( + cls: type[Channel], + inst: Gate, + gate_duration: float, + t1s: list[float], + t2s: list[float] | None = None, + custom_gates: CustomGateMap | None = None, + ) -> Channel: + """Create a decoherence Channel based on the coherence times. + + In this construction, decoherence is applied _after_ the ideal gate unitary. + + :param inst: The target instruction. + :param gate_duration: The duration of the gate. + :param t1s: The t1 time(s) of the qubits + :param t2s: The t2 time(s) of the qubits. Default to 2*t1. + """ + unitary = get_instruction_unitary(inst, custom_gates) + qubits = inst.get_qubit_indices() + num_sys = len(qubits) + if num_sys != len(t1s): + raise ValueError(f"Expected {num_sys} T1 values for {inst.out()}, got {len(t1s)}.") + if t2s is None: + t2s = [2 * t1 for t1 in t1s] + else: + if num_sys != len(t2s): + raise ValueError(f"Expected {num_sys} T2 values for {inst.out()}, got {len(t2s)}.") + + t1_array = jnp.asarray(t1s) + tphi_array = 1 / (1 / jnp.asarray(t2s) - 1 / t1_array) + + choi = qx.thermal_relaxation_choi(t1s=t1_array, tphis=tphi_array, duration=gate_duration) + process = qx.to_superop(choi @ unitary) + return cls( + inst=inst, + process=process, + target_unitary=unitary, + ) + + @classmethod + def from_superoperator( + cls: type[Channel], + inst: Gate, + process: qx.SuperOp, + target_unitary: qx.Unitary | None = None, + custom_gates: CustomGateMap | None = None, + ) -> Channel: + """Create a Channel from a pre-built superoperator. + + If ``target_unitary`` is not provided it is inferred from the gate + instruction using the standard gate set (and ``custom_gates`` if given). + + :param inst: The gate to which the channel applies. + :param process: The noisy process superoperator (includes the gate unitary). + :param target_unitary: The ideal gate unitary. Resolved automatically + when omitted. + :param custom_gates: Optional dictionary of custom gate definitions, + used only when ``target_unitary`` is ``None``. + :return: A Channel instance. + """ + if target_unitary is None: + target_unitary = get_instruction_unitary(inst, custom_gates) + return cls(inst=inst, process=process, target_unitary=target_unitary) + + # ────────────────────────────────────────────── + # Cached representation conversions + # ────────────────────────────────────────────── + + @cached_property + def noise_process(self) -> qx.SuperOp: + r"""The noise-only channel with the ideal gate unitary factored out. + + If the full channel is :math:`\\mathcal{E} = \\Lambda \\circ \\mathcal{U}`, this + returns :math:`\\Lambda`. + """ + return qx.to_superop(self.process @ self.unitary.h) + + # ────────────────────────────────────────────── + # Fidelity properties + # ────────────────────────────────────────────── + + @cached_property + def fidelity(self) -> float: + r"""Average gate fidelity :math:`F_{\\mathrm{avg}}` of the channel relative to the ideal gate.""" + return float(qx.process_fidelity_to_average_fidelity(self.pauli_fidelity, dims=self.unitary.dims[0])) + + @cached_property + def infidelity(self) -> float: + r"""Average gate infidelity :math:`1 - F_{\\mathrm{avg}}`.""" + return 1.0 - self.fidelity + + @cached_property + def pauli_fidelity(self) -> float: + """Process fidelity (entanglement fidelity) :math:`F_e` relative to the ideal gate.""" + process, unitary = qx.promote_hilbert_space(self.process, qx.to_superop(self.unitary)) + return float(qx.process_fidelity(process, unitary)) + + @cached_property + def pauli_infidelity(self) -> float: + """Process infidelity :math:`1 - F_e`.""" + return 1.0 - self.pauli_fidelity + + @cached_property + def stochastic_infidelity(self) -> float: + """Stochastic (incoherent) component of the process infidelity.""" + return float(qx.stochastic_infidelity(self.noise_process)) + + @cached_property + def stochastic_fidelity(self) -> float: + """Stochastic fidelity :math:`1 - e_S`.""" + return 1.0 - self.stochastic_infidelity + + @cached_property + def coherent_infidelity(self) -> float: + """Coherent component of the process infidelity: :math:`e_C = e - e_S`.""" + return self.pauli_infidelity - self.stochastic_infidelity + + @cached_property + def coherent_fidelity(self) -> float: + """Coherent fidelity :math:`1 - e_C`.""" + return 1.0 - self.coherent_infidelity + + @cached_property + def unitarity(self) -> float: + """Unitarity of the channel.""" + return float(qx.unitarity(self.noise_process)) + + # ────────────────────────────────────────────── + # Channel analysis methods + # ────────────────────────────────────────────── + + def pauli_twirl(self) -> Channel: + """Return a Pauli-twirled version of this channel. + + Pauli twirling projects the channel onto the Pauli diagonal, eliminating + off-diagonal coherences in the Pauli-Liouville representation. The + resulting channel is a stochastic Pauli channel with the same diagonal + error rates. + """ + ptm = qx.to_pauli_liouville(self.process) + # Keep only the diagonal of the PTM + twirled_ptm_matrix = jnp.diag(jnp.diag(ptm.matrix)) + twirled_superop = qx.to_superop(qx.PauliLiouville.from_matrix(twirled_ptm_matrix, self.process.dims)) + return replace(self, process=twirled_superop) + + @cached_property + def _unitary_error_component(self) -> Array: + """Extract the dominant unitary from the noise-only channel. + + Uses eigendecomposition + SVD polar decomposition to find the closest + unitary to the noise channel. + """ + choi_matrix = qx.to_choi(self.noise_process).matrix + d = 2**self.num_qubits + + # Dominant eigenvector of the Choi matrix + eigenvalues, eigenvectors = jnp.linalg.eigh(choi_matrix) + dominant_eigenvector = eigenvectors[:, jnp.argmax(jnp.abs(eigenvalues))] + + # SVD polar decomposition to extract the closest unitary + u, _, vh = jnp.linalg.svd(dominant_eigenvector.reshape(d, d).T) + return u @ vh + + def to_coherent_channel(self) -> Channel: + """Isolate the coherent (unitary) component of the noise. + + Extracts the dominant unitary from the noise Choi matrix via polar + decomposition and returns a channel consisting of that unitary error + composed with the ideal gate. + """ + u_error = self._unitary_error_component + u_error_qx = qx.Unitary.from_matrix(u_error, self.process.dims) + coherent_superop = qx.to_superop(u_error_qx @ self.unitary) + return replace(self, process=coherent_superop) + + def to_stochastic_channel(self) -> Channel: + r"""Isolate the stochastic (incoherent) component of the noise. + + The full channel decomposes as + :math:`\\mathcal{E} = \\mathcal{S} \\circ \\mathcal{U}_{\\mathrm{err}} \\circ \\mathcal{U}_{\\mathrm{gate}}`. + This method factors out the coherent unitary error and returns + :math:`\\mathcal{S} \\circ \\mathcal{U}_{\\mathrm{gate}}`. + """ + u_error = self._unitary_error_component + # Get the noise-only superoperator and compose with U_err† + noise_superop = self.noise_process.matrix + u_err_inv_superop = jnp.kron(u_error.conj(), u_error.conj().T) + stochastic_noise_superop = noise_superop @ u_err_inv_superop + # Recompose with the ideal gate + ideal_superop = jnp.kron(self.unitary.matrix, self.unitary.matrix.conj()) + stochastic_superop = stochastic_noise_superop @ ideal_superop + return replace(self, process=qx.SuperOp.from_matrix(stochastic_superop, self.process.dims)) + + def is_pauli(self) -> bool: + """Check if the noise channel is a Pauli (stochastic Pauli) channel. + + A Pauli channel has a diagonal Pauli transfer matrix (noise-only part). + """ + ptm = qx.to_pauli_liouville(self.noise_process).matrix + mask = ~jnp.eye(ptm.shape[0], dtype=bool) + return bool(jnp.allclose(ptm[mask], 0)) + + def to_pauli_vector(self) -> Array: + """Convert the noise channel to a Pauli error probability vector. + + Returns the vector of probabilities for each Pauli error in lexicographic + order (II, IX, IY, IZ, XI, XX, ...). The vector sums to 1.0. + """ + noise_superop = self.noise_process.matrix + num_qubits = self.num_qubits + dim = noise_superop.shape[0] + + # Build all Pauli operators and their superoperators + pauli_matrices = qx.ensembles.PAULIS.matrix # (4, 2, 2): I, X, Y, Z + all_pauli_products = list(product(pauli_matrices, repeat=num_qubits)) + pauli_error_rates = [] + for pauli_tuple in all_pauli_products: + pauli_op = reduce(jnp.kron, pauli_tuple) + pauli_superop = jnp.kron(pauli_op, pauli_op.conj()) + rate = float(jnp.abs(jnp.trace(noise_superop @ pauli_superop) / dim)) + pauli_error_rates.append(rate) + + return jnp.array(pauli_error_rates, dtype=float) + + @cached_property + def pauli_vector(self) -> Array: + """The Pauli error probability vector of the noise channel.""" + return self.to_pauli_vector() + + # ────────────────────────────────────────────── + # Visualization + # ────────────────────────────────────────────── + + def plot(self, only_noise: bool = True, show_identity: bool = False) -> Figure: + """Plot the Pauli transfer matrix of the channel. + + :param only_noise: If True, plot the noise-only channel (gate unitary factored out). + If False, plot the full channel including the gate unitary. + :param show_identity: If True, include the identity component in the noise-only plot. + If False (default), visualize the generator of the noise channel via the matrix + logarithm of the PTM. For near-identity noise this approximates PTM - I, but + correctly captures the Lie-algebraic structure of the channel. + Only applies when ``only_noise=True``. + :return: A Plotly Figure. + """ + if only_noise: + channel = self.noise_process + if not show_identity: + ptm = qx.to_pauli_liouville(channel) + log_ptm = scipy_logm(np.asarray(ptm.matrix)) + channel = qx.PauliLiouville.from_matrix(jnp.array(log_ptm), channel.dims) + title_prefix = "Noise Channel" + else: + channel = self.process + title_prefix = "Full Channel" + + fig = qx.plot(channel) + fig.update_layout( + title=( + f"{title_prefix} for {self.inst.out()}
" + f"𝜀={self.pauli_infidelity * 100:.2f}%, " + f"𝜀u={self.coherent_infidelity * 100:.2f}%, " + f"𝜀s={self.stochastic_infidelity * 100:.2f}%" + ) + ) + return fig + + # ────────────────────────────────────────────── + # Serialization + # ────────────────────────────────────────────── + + def to_json(self) -> str: + """Serialize Channel to a JSON string. + + :return: JSON string representation. + """ + data = { + "schema_version": 1, + "inst": self.inst.out(), + "superop": _pack_operator(self.process), + } + data["target_unitary"] = _pack_operator(self.target_unitary) + + return json.dumps(data) + + @classmethod + def from_json(cls: type[Channel], json_str: str) -> Channel: + """Deserialize a Channel from a JSON string. + + :param json_str: JSON string as produced by :meth:`to_json`. + :return: Channel instance. + """ + data = json.loads(json_str) + inst = _parse_quil_instruction(data["inst"]) + if not isinstance(inst, Gate): + raise TypeError(f"Channel JSON must contain a gate instruction, got {type(inst).__name__}.") + + superop_data = data["superop"] + shape = tuple(superop_data["shape"]) + superop_array = _unpack_complex_array(superop_data) + dims = _unpack_operator_dims(superop_data, shape, superoperator=True) + superop = qx.SuperOp.from_matrix(superop_array, dims) + + if "target_unitary" in data: + u_data = data["target_unitary"] + u_shape = tuple(u_data["shape"]) + u_array = _unpack_complex_array(u_data) + u_dims = _unpack_operator_dims(u_data, u_shape, superoperator=False) + target_unitary = qx.Unitary.from_matrix(u_array, u_dims) + else: + target_unitary = get_instruction_unitary(inst) + + return cls(inst=inst, process=superop, target_unitary=target_unitary) + + # ────────────────────────────────────────────── + # Dunder methods + # ────────────────────────────────────────────── + + def __str__(self) -> str: + """Return a simplified string representation showing the gate and process fidelity.""" + return f"<{self.inst.out()} ~ ({100 * self.pauli_fidelity:.2f}%)>" + + def __eq__(self, other: object) -> bool: + """Check equality by instruction and exact process and ideal-gate matrices. + + Equality is exact (no fidelity tolerance): two channels are equal only if they + share the same instruction and bit-for-bit identical process and target-unitary + matrices. Making tolerance decisions on the user's behalf is deliberately avoided. + """ + if not isinstance(other, Channel): + return False + if self.inst != other.inst: + return False + return bool( + jnp.array_equal(self.process.matrix, other.process.matrix) + and jnp.array_equal(self.target_unitary.matrix, other.target_unitary.matrix) + ) + + __hash__ = None # type: ignore[assignment] + + def __matmul__(self, other: Channel) -> Channel: + r"""Compose two channels: ``channel_B @ channel_A``. + + Both channels share the same gate instruction. The composition factors + out one copy of the gate unitary so the result represents the sequential + application of the two noisy processes: + + :math:`\\mathcal{E}_B \\circ \\mathcal{U}^\\dagger \\circ \\mathcal{E}_A` + + This is the natural composition: if ``channel_A`` already includes the + gate, applying ``channel_B`` after it should not double-count the gate. + """ + if not isinstance(other, Channel): + return NotImplemented + if self.inst != other.inst: + raise ValueError(f"Cannot compose channels for different gates: {self.inst.out()} vs {other.inst.out()}") + # E_B @ U† @ E_A (factor out one gate unitary between the two channels) + u_dag_superop = qx.to_superop(self.unitary.h) + composed_superop = qx.to_superop(self.process @ u_dag_superop @ other.process) + return replace(self, process=composed_superop) + + def __pow__(self, power: float) -> Channel: + r"""Raise the channel's noise to a fractional ``power``, preserving the gate. + + With the channel written :math:`\\mathcal{E} = \\Lambda \\circ \\mathcal{U}`, this returns + :math:`\\Lambda^{power} \\circ \\mathcal{U}`: only the noise :math:`\\Lambda` is raised to the + fractional matrix power, while the ideal gate is kept. So ``power = 0`` yields the noiseless + gate, ``1`` leaves the channel unchanged, and ``> 1`` strengthens the noise. This is the knob + used to sweep noise strength. + """ + if not isinstance(power, (int, float)): + return NotImplemented + powered_noise_matrix = fractional_matrix_power(np.asarray(self.noise_process.matrix), power) + powered_noise = qx.SuperOp.from_matrix(jnp.asarray(powered_noise_matrix), self.noise_process.dims) + process = qx.to_superop(powered_noise @ qx.to_superop(self.unitary)) + return replace(self, process=process) + + def __or__(self, other: Channel | MeasurementChannel) -> CycleChannel: + """Tensor product of two channels on disjoint qubits, producing a CycleChannel. + + The result represents a cycle containing both operations acting in parallel + on disjoint qubits. The DefCircuit encodes the parallel operations as + formal instructions. + + :param other: Another Channel or MeasurementChannel on disjoint qubits. + :return: A CycleChannel representing the tensor product. + """ + if not isinstance(other, (Channel, MeasurementChannel)): + return NotImplemented + + # Validate disjoint qubits + self_qubits = set(self.qubits) + other_qubits = set(other.qubits) + if self_qubits & other_qubits: + raise ValueError(f"Cannot tensor channels with overlapping qubits: {self_qubits & other_qubits}") + + return _build_cycle_channel([self, other]) + + +@dataclass(frozen=True) +class MeasurementChannel: + """A measurement noise channel attaches a quantum instrument to a specific measurement operation. + + The ``process`` field is a ``qx.QuantumInstrument`` which models both classification + errors and post-measurement back-action. + """ + + inst: Measurement + """The measurement operation to which the channel applies.""" + + process: qx.QuantumInstrument + """A quantum instrument representation of the noisy measurement.""" + + @cached_property + def qubits(self) -> list[int]: + """The qubits which the measurement applies to.""" + qubit = self.inst.qubit + return [qubit.index if hasattr(qubit, "index") else int(qubit)] + + # ────────────────────────────────────────────── + # Constructors + # ────────────────────────────────────────────── + + @classmethod + def from_readout_fidelity( + cls: type[MeasurementChannel], + inst: Measurement, + fidelity: float, + asymmetry: float = 0.0, + dim: int = 2, + ) -> MeasurementChannel: + """Create a readout quantum instrument with optional asymmetry. + + Produces a perfectly QND measurement with the given classification fidelity. + Error is distributed only between adjacent levels: P(j+1|j) and P(j|j+1). + Non-adjacent confusion is zero. + + :param inst: The measurement instruction. + :param fidelity: The average readout fidelity. + :param asymmetry: Value between -1 and +1. Zero is symmetric. + Positive biases toward upward confusion P(j+1|j), negative toward downward P(j|j+1). + :param dim: The dimension of the measured system (2 for qubits, 3 for qutrits, etc.). + :return: A MeasurementChannel instance. + """ + # Compute per-pair error factor so that the average diagonal equals fidelity. + # Each adjacent pair (j, j+1) contributes error_factor*(1+a) + error_factor*(1-a) + # = 2*error_factor to total off-diagonal sum. With (dim-1) pairs, the average + # column error is 2*(dim-1)*error_factor/dim, which we set equal to (1-fidelity). + error_factor = dim * (1 - fidelity) / (2 * (dim - 1)) + + confusion = jnp.zeros((dim, dim)) + for j in range(dim - 1): + confusion = confusion.at[j + 1, j].set(error_factor * (1 + asymmetry)) + confusion = confusion.at[j, j + 1].set(error_factor * (1 - asymmetry)) + # Set diagonal so each column sums to 1 + col_sums = confusion.sum(axis=0) + confusion = confusion + jnp.diag(1 - col_sums) + + transition = jnp.eye(dim) + instrument = qx.instrument_from_confusion_and_transition( + confusion_matrix=confusion, + transition_matrix=transition, + dims=(dim,), + measured_qudits=(0,), + ) + return cls(inst=inst, process=instrument) + + @classmethod + def from_confusion_and_transition( + cls: type[MeasurementChannel], + inst: Measurement, + confusion_matrix: Array, + transition_matrix: Array, + ) -> MeasurementChannel: + """Create a MeasurementChannel from a confusion matrix and a transition matrix. + + Provides independent control over measurement classification accuracy + and post-measurement quantum state evolution. + + **Matrix Conventions (column-stochastic):** + + - ``confusion_matrix[i, j]``: P(outcome i | prepared j) + - ``transition_matrix[k, j]``: P(ending in k | input j) + - Columns sum to 1.0 + + :param inst: The measurement instruction. + :param confusion_matrix: A (d, d) classification matrix. + :param transition_matrix: A (d, d) post-measurement transition matrix. + :return: A MeasurementChannel instance. + """ + confusion = jnp.asarray(confusion_matrix) + dim = confusion.shape[0] + instrument = qx.instrument_from_confusion_and_transition( + confusion_matrix=confusion, + transition_matrix=jnp.asarray(transition_matrix), + dims=(dim,), + measured_qudits=(0,), + ) + return cls(inst=inst, process=instrument) + + @classmethod + def from_axis( + cls: type[MeasurementChannel], + inst: Measurement, + theta: float = 0.0, + phi: float = 0.0, + sharpness: float = 1.0, + ) -> MeasurementChannel: + """Create a MeasurementChannel from a Bloch sphere measurement axis. + + The angles refer to the standard Bloch sphere notation. + Theta=0, phi=0 is the Z axis (computational basis measurement). + + :param inst: The measurement instruction. + :param theta: The colatitude with respect to the z-axis. + :param phi: The longitude with respect to the x-axis. + :param sharpness: The sharpness of the measurement. 1.0 is projective, + 0.0 is no measurement. 0 < s < 1 is a weak measurement. + :return: A MeasurementChannel instance. + """ + instrument = qx.instrument_from_axis( + theta=theta, + phi=phi, + sharpness=sharpness, + ) + return cls(inst=inst, process=instrument) + + @classmethod + def from_binary_discriminator( + cls: type[MeasurementChannel], + inst: Measurement, + dim: int, + threshold: int, + fidelity: float = 1.0, + ) -> MeasurementChannel: + """Create a MeasurementChannel for a binary discriminator. + + Models a measurement that reports a single classical *bit* for a + ``dim``-level system by thresholding: levels ``[0, threshold)`` yield + outcome ``0`` and levels ``[threshold, dim)`` yield outcome ``1``. The + resulting instrument therefore always has exactly two outcomes, so leaked + levels are lumped in with whichever side of the threshold they fall on + (the usual case being a "dark" ground state vs. everything "bright"). + + For example, ``threshold=1, dim=2`` is an ordinary qubit readout + (``{0}`` -> 0, ``{1}`` -> 1); ``threshold=1, dim=3`` discriminates + ``{0}`` vs ``{1, 2}`` (ground vs. excited-or-leaked); ``threshold=2, + dim=3`` discriminates ``{0, 1}`` vs ``{2}`` (i.e. flags leakage only). + + An optional ``fidelity`` parameter degrades the ideal discriminator with + uniform classification noise. + + :param inst: The measurement instruction. + :param dim: The dimension of the measured system. + :param threshold: The split point: levels below it report 0, levels at or + above it report 1. Must satisfy ``1 <= threshold < dim``. + :param fidelity: Additional classification fidelity applied on top of the + discrimination (1.0 = perfect discriminator). + :return: A MeasurementChannel instance. + """ + if not (1 <= threshold < dim): + raise ValueError(f"threshold must satisfy 1 <= threshold < dim, got threshold={threshold}, dim={dim}") + + # Ideal two-outcome confusion matrix of shape (num_outcomes=2, dim): + # column j (prepared level j) puts all its weight on outcome 0 if + # j < threshold, else on outcome 1. Two rows so the instrument has + # exactly two outcomes (never a phantom, zero-probability outcome). + confusion = jnp.zeros((2, dim)) + for j in range(dim): + confusion = confusion.at[int(j >= threshold), j].set(1.0) + + # Optionally degrade with uniform noise across the two outcomes. + if fidelity < 1.0: + confusion = fidelity * confusion + (1 - fidelity) * jnp.ones((2, dim)) / 2 + + transition = jnp.eye(dim) + instrument = qx.instrument_from_confusion_and_transition( + confusion_matrix=confusion, + transition_matrix=transition, + dims=(dim,), + measured_qudits=(0,), + ) + return cls(inst=inst, process=instrument) + + # ────────────────────────────────────────────── + # Properties + # ────────────────────────────────────────────── + + @cached_property + def confusion_matrix(self) -> Array: + """The confusion matrix of the measurement. + + Shape ``(num_outcomes, d_measured)``. + Entry ``[i, j]`` is P(outcome i | prepared j). + """ + return self.process.confusion_matrix # type: ignore[no-any-return] + + @cached_property + def transition_matrix(self) -> Array: + """The post-measurement transition matrix. + + Shape ``(d, d)``. Entry ``[k, j]`` is P(ending in k | input j), + marginalized over all measurement outcomes. + """ + return self.process.transition_matrix # type: ignore[no-any-return] + + @cached_property + def non_demolition_fidelity(self) -> float: + """Quantum non-demolition (QND) fidelity. + + Measures how well the measurement preserves computational basis states, + averaged over outcomes and input states. + """ + return float(qx.non_demolition_fidelity(self.process)) + + @cached_property + def instrument_fidelity(self) -> float: + """Overall instrument fidelity w.r.t. ideal QND measurement. + + Accounts for both classification errors and post-measurement state disturbance. + """ + return float(qx.instrument_fidelity(self.process)) + + @cached_property + def classification_fidelity(self) -> float: + """Classification fidelity: average probability of correctly identifying the measurement outcome.""" + return float(qx.classification_fidelity(self.process)) + + # ────────────────────────────────────────────── + # Visualization + # ────────────────────────────────────────────── + + def plot(self) -> Figure: + """Plot the quantum instrument using the quax visualization. + + Shows per-outcome superoperator matrices and the total CPTP channel. + + :return: A Plotly Figure. + """ + fig = qx.plot(self.process) + fig.update_layout( + title=( + f"Quantum Instrument MEASURE {self.qubits[0]}
" + f"Cls: {100 * self.classification_fidelity:.2f}%, " + f"QND: {100 * self.non_demolition_fidelity:.2f}%, " + f"Instrument: {100 * self.instrument_fidelity:.2f}%" + ) + ) + return fig + + # ────────────────────────────────────────────── + # Serialization + # ────────────────────────────────────────────── + + def to_json(self) -> str: + """Serialize MeasurementChannel to a JSON string. + + :return: JSON string representation. + """ + # Store per-outcome superoperator matrices. + instrument_data = [] + for i in range(self.process.num_outcomes): + superop_i, _ = self.process.outcome_superop(i) + instrument_data.append(_pack_operator(superop_i)) + + data = { + "schema_version": 1, + "inst": self.inst.out(), + "instruments": instrument_data, + "measured_qudits": list(self.process.measured_qudits), + } + return json.dumps(data) + + @classmethod + def from_json(cls: type[MeasurementChannel], json_str: str) -> MeasurementChannel: + """Deserialize a MeasurementChannel from a JSON string. + + :param json_str: JSON string as produced by :meth:`to_json`. + :return: MeasurementChannel instance. + """ + data = json.loads(json_str) + inst = _parse_quil_instruction(data["inst"]) + if not isinstance(inst, Measurement): + raise TypeError( + f"MeasurementChannel JSON must contain a measurement instruction, got {type(inst).__name__}." + ) + measured_qudits = tuple(data["measured_qudits"]) + + superop_matrices = [] + instrument_dims = None + for inst_data in data["instruments"]: + shape = tuple(inst_data["shape"]) + arr = _unpack_complex_array(inst_data) + op_dims = _unpack_operator_dims(inst_data, shape, superoperator=True) + if instrument_dims is None: + instrument_dims = op_dims + elif instrument_dims != op_dims: + raise ValueError("All serialized measurement outcomes must have the same dims.") + superop_matrices.append(arr) + + if instrument_dims is None: + raise ValueError("MeasurementChannel JSON must contain at least one outcome superoperator.") + + instrument = qx.QuantumInstrument.from_matrix(jnp.stack(superop_matrices), instrument_dims, measured_qudits) + return cls(inst=inst, process=instrument) + + # ────────────────────────────────────────────── + # Dunder methods + # ────────────────────────────────────────────── + + def __str__(self) -> str: + """Return a simplified string representation.""" + return f"" + + def __eq__(self, other: object) -> bool: + """Check equality by instruction and exact instrument matrix (no tolerance).""" + if not isinstance(other, MeasurementChannel): + return False + if self.inst != other.inst: + return False + return bool(jnp.array_equal(self.process.matrix, other.process.matrix)) + + __hash__ = None # type: ignore[assignment] + + def __matmul__(self, other: MeasurementChannel) -> MeasurementChannel: + """Compose two measurement channels on the same qubit. + + Models sequential application: ``channel_B @ channel_A`` means + apply ``channel_A`` first, then ``channel_B``. + """ + if not isinstance(other, MeasurementChannel): + return NotImplemented + if self.inst != other.inst: + raise ValueError( + f"Cannot compose measurement channels for different qubits: {self.inst.out()} vs {other.inst.out()}" + ) + composed = self.process @ other.process + return replace(self, process=composed) + + def __pow__(self, power: float) -> MeasurementChannel: + r"""Raise the measurement's classification noise to a fractional ``power``. + + The confusion and transition matrices are column-stochastic, so a fractional power is + only well-defined through their *generator*: with :math:`M = e^{G}` (where :math:`G` has + zero column sums), the principled power is :math:`M^{p} = e^{p G}`, which is again + column-stochastic. This avoids the non-physical (e.g. negative) entries that a naive + ``fractional_matrix_power`` of a stochastic matrix can produce. + + So ``power = 0`` is an ideal (noiseless) measurement, ``1`` leaves it unchanged, and + ``> 1`` degrades it. Mirrors :meth:`Channel.__pow__` for sweeping readout-noise strength. + + :raises ValueError: If the matrix is not embeddable (no real generator) so that the + powered matrix is not a valid stochastic matrix. + """ + if not isinstance(power, (int, float)): + return NotImplemented + + def _powered_stochastic(matrix: Array) -> Array: + m = np.asarray(matrix, dtype=complex) + # Generator G = log(M); M**power = exp(power * G). Zero column sums of G are + # preserved under scaling, so exp(power * G) stays column-stochastic. + powered = scipy_expm(power * scipy_logm(m)) + if np.max(np.abs(powered.imag)) > 1e-9: + raise ValueError( + f"MeasurementChannel ** {power} has no real generator; the matrix is not " + "embeddable, so the powered measurement is not a valid stochastic matrix." + ) + powered = powered.real + if np.any(powered < -1e-9) or not np.allclose(powered.sum(axis=0), 1.0, atol=1e-6): + raise ValueError( + f"MeasurementChannel ** {power} is not a valid (non-negative, column-stochastic) " + "matrix; the underlying confusion/transition matrix is not embeddable for this power." + ) + # Clip away sub-tolerance numerical negatives validated above. + return jnp.asarray(np.clip(powered, 0.0, None)) + + confusion = _powered_stochastic(self.confusion_matrix) + transition = _powered_stochastic(self.transition_matrix) + return MeasurementChannel.from_confusion_and_transition(self.inst, confusion, transition) + + def __or__(self, other: Channel | MeasurementChannel) -> CycleChannel: + """Tensor product of two channels on disjoint qubits, producing a CycleChannel. + + :param other: Another Channel or MeasurementChannel on disjoint qubits. + :return: A CycleChannel representing the tensor product. + """ + if not isinstance(other, (Channel, MeasurementChannel)): + return NotImplemented + + self_qubits = set(self.qubits) + other_qubits = set(other.qubits) + if self_qubits & other_qubits: + raise ValueError(f"Cannot tensor channels with overlapping qubits: {self_qubits & other_qubits}") + + return _build_cycle_channel([self, other]) + + +@dataclass(frozen=True) +class ResetChannel: + """A reset noise channel attaches a superoperator to a specific reset operation. + + The ``process`` field is a ``qx.SuperOp`` which *includes* the ideal reset, so the channel + replaces the reset instruction rather than being applied after it. + """ + + inst: ResetQubit + """The reset operation to which the channel applies.""" + + process: qx.SuperOp + """A superoperator representation of the noisy reset (including ideal reset).""" + + def __post_init__(self) -> None: + """Validate that ResetChannel is attached to a targeted reset.""" + if not isinstance(self.inst, ResetQubit): + raise TypeError("ResetChannel only supports targeted ResetQubit instructions.") + + # ────────────────────────────────────────────── + # Constructors + # ────────────────────────────────────────────── + + @classmethod + def from_reset_fidelity( + cls: type[ResetChannel], + inst: ResetQubit, + fidelity: float, + dim: int = 2, + ) -> ResetChannel: + r"""Create a ResetChannel with depolarizing noise scaled to the given process fidelity. + + The ideal reset channel maps every state to :math:`|0\\rangle\\langle 0|`. Noise is + modelled as a depolarising channel applied after the ideal reset. + + :param inst: The reset instruction. + :param fidelity: Process fidelity of the reset channel, :math:`F \\in [0, 1]`. + 1.0 yields an ideal reset; values below 1 introduce depolarising noise. + :param dim: Hilbert-space dimension (2 for qubits). + :return: A ResetChannel instance. + """ + if not isinstance(inst, ResetQubit): + raise TypeError("ResetChannel only supports targeted ResetQubit instructions.") + + ideal_superop = qx.gates.RESET(dim=dim) + p = 1.0 - fidelity + d2 = dim * dim + # Depolarising channel in superop form: (1-p)*S_ideal + p*(I/d) for all inputs + # The completely depolarising superop maps everything to I/d: + # its rows are all zero except the diagonal entries corresponding to + # the trace extraction (maps vec(ρ) → vec(I/d) = vec(I)/d). + depol_superop_matrix = jnp.zeros((d2, d2), dtype=complex) + # vec(I/d) has value 1/d at positions 0, d+1, 2(d+1), ... i.e. diagonal entries + vec_identity_over_d = jnp.zeros(d2, dtype=complex) + for i in range(dim): + vec_identity_over_d = vec_identity_over_d.at[i * dim + i].set(1.0 / dim) + # The trace functional extracts sum of diagonal: positions 0, d+1, ... + trace_row = jnp.zeros(d2, dtype=complex) + for i in range(dim): + trace_row = trace_row.at[i * dim + i].set(1.0) + # Depolarising superop: each row of output is vec(I/d) * Tr(ρ) + depol_superop_matrix = jnp.outer(vec_identity_over_d, trace_row) + noisy_superop_matrix = (1.0 - p) * ideal_superop.matrix + p * depol_superop_matrix + noisy_superop = qx.SuperOp.from_matrix(noisy_superop_matrix, ideal_superop.dims) + return cls(inst=inst, process=noisy_superop) + + # ────────────────────────────────────────────── + # Properties + # ────────────────────────────────────────────── + + @cached_property + def qubits(self) -> list[int]: + """The qubit(s) that the reset applies to.""" + qubit = self.inst.qubit + if qubit is None: + return [] + return [qubit.index if hasattr(qubit, "index") else int(qubit)] + + @cached_property + def fidelity(self) -> float: + r"""Process fidelity of the reset channel relative to the ideal reset. + + Defined as :math:`F = \\mathrm{Tr}[\\Lambda_{\\mathrm{ideal}}^\\dagger \\Lambda] / d^2` + where :math:`\\Lambda` is the Choi matrix of the noisy channel and + :math:`\\Lambda_{\\mathrm{ideal}}` is the ideal-reset Choi. + """ + dim = self.process.dims[0][0] + ideal_choi = qx.to_choi(qx.gates.RESET(dim=dim)) + noisy_choi = qx.to_choi(self.process) + # Process fidelity = Tr[ideal_choi† @ noisy_choi] / d^2 + d2 = float(dim * dim) + return float(jnp.real(jnp.trace(ideal_choi.matrix.conj().T @ noisy_choi.matrix)) / d2) + + @cached_property + def noise_process(self) -> qx.SuperOp: + """The noise-only channel (ideal reset factored out). + + For a reset channel the noise framing is less natural than for unitary gates; + this property returns the full process superoperator. + """ + return self.process + + # ────────────────────────────────────────────── + # Visualization + # ────────────────────────────────────────────── + + def plot(self) -> Figure: + """Plot the Pauli transfer matrix of the reset channel. + + :return: A Plotly Figure. + """ + fig = qx.plot(self.process) + qubit_str = str(self.qubits[0]) if self.qubits else "?" + fig.update_layout(title=(f"Reset Channel RESET {qubit_str}
F_\u03c7={self.fidelity * 100:.2f}%")) + return fig + + # ────────────────────────────────────────────── + # Serialization + # ────────────────────────────────────────────── + + def to_json(self) -> str: + """Serialize ResetChannel to a JSON string. + + :return: JSON string representation. + """ + data = { + "schema_version": 1, + "inst": self.inst.out(), + "superop": _pack_operator(self.process), + } + return json.dumps(data) + + @classmethod + def from_json(cls: type[ResetChannel], json_str: str) -> ResetChannel: + """Deserialize a ResetChannel from a JSON string. + + :param json_str: JSON string as produced by :meth:`to_json`. + :return: ResetChannel instance. + """ + data = json.loads(json_str) + inst = _parse_quil_instruction(data["inst"]) + if not isinstance(inst, ResetQubit): + raise TypeError(f"ResetChannel JSON must contain a targeted reset instruction, got {type(inst).__name__}.") + superop_data = data["superop"] + shape = tuple(superop_data["shape"]) + arr = _unpack_complex_array(superop_data) + dims = _unpack_operator_dims(superop_data, shape, superoperator=True) + process = qx.SuperOp.from_matrix(arr, dims) + return cls(inst=inst, process=process) + + # ────────────────────────────────────────────── + # Dunder methods + # ────────────────────────────────────────────── + + def __str__(self) -> str: + """Return a simplified string representation.""" + qubit_str = str(self.qubits[0]) if self.qubits else "?" + return f"" + + def __eq__(self, other: object) -> bool: + """Check equality by instruction and exact process matrix (no tolerance).""" + if not isinstance(other, ResetChannel): + return False + if self.inst != other.inst: + return False + return bool(jnp.array_equal(self.process.matrix, other.process.matrix)) + + __hash__ = None # type: ignore[assignment] + + +@dataclass(frozen=True) +class CycleChannel: + """A cycle noise channel attaches superoperators to a specific cycle. + + Cycles can include gates and measurements. The constituent channels are stored + directly, allowing fidelity metrics and serialization to be derived from them. + """ + + inst: Gate + """The cycle to which the channel applies.""" + + defcircuit: DefCircuit + """The DefCircuit representing the logical cycle to which instruction represents.""" + + channels: tuple[Channel | MeasurementChannel, ...] + """Constituent channels (one per operation in the cycle) on disjoint qubits.""" + + def __post_init__(self) -> None: + """Validate that every instruction in the cycle body has a corresponding channel. + + Downstream consumers (the resolver, the stim converter) use only ``channels`` and + ignore ``defcircuit``; a missing channel would silently drop that operation's noise. + Operations are matched by identity (name, params, concrete qubits), independent of + the DefCircuit's formal-argument naming. + """ + if len(self.expanded_instructions) != len(self.channels): + raise ValueError( + "CycleChannel is incomplete: every instruction in the cycle's DefCircuit " + "body must have a corresponding channel. " + f"\nDefCircuit body: {self.expanded_instructions}" + f"\nChannels: {self.channels}" + ) + for instruction, channel in zip(self.expanded_instructions, self.channels, strict=True): + if str(instruction) != str(channel.inst): + raise ValueError( + "CycleChannel is incomplete: every instruction in the cycle's DefCircuit " + "body must have a corresponding channel. " + f"\nDefCircuit body: {instruction}" + f"\nChannels: {channel.inst}" + ) + + # ────────────────────────────────────────────── + # Derived properties + # ────────────────────────────────────────────── + + @cached_property + def expanded_instructions(self) -> list[Gate | Measurement | ResetQubit]: + """Return the expanded instructions of the defcircuit.""" + qarg_to_qubit = dict(zip(self.defcircuit.qubit_variables, self.inst.get_qubit_indices(), strict=False)) + instructions: list[Gate | Measurement | ResetQubit] = [] + for inst in self.defcircuit.instructions: + match inst: + case Measurement(): + instructions.append(Measurement(qubit=qarg_to_qubit[inst.qubit], classical_reg=inst.classical_reg)) # type: ignore[index] + case ResetQubit(): + instructions.append(ResetQubit(qarg_to_qubit[inst.qubit])) # type: ignore[index] + case Gate(): + instructions.append(Gate(inst.name, inst.params, [qarg_to_qubit[q] for q in inst.qubits])) # type: ignore[index] + case _: + raise TypeError(f"Unsupported instruction type in defcircuit: {type(inst).__name__}") + return instructions + + @cached_property + def operator(self) -> tuple[qx.SuperOp | qx.QuantumInstrument, ...]: + """Tuple of process superoperators, one per constituent channel.""" + return tuple(ch.process for ch in self.channels) + + @cached_property + def qubits(self) -> list[int]: + """All qubits in the cycle, derived from the instruction.""" + return self.inst.get_qubit_indices() + + @cached_property + def pauli_fidelity(self) -> float: + """Product of process (Pauli) fidelities over all gate channels in the cycle. + + Measurement channels do not contribute a gate fidelity and are skipped. + For near-ideal noise the product approximation is exact since constituent + channels act on disjoint subsystems. + """ + f = 1.0 + for ch in self.channels: + if isinstance(ch, Channel): + f *= ch.pauli_fidelity + return f + + @cached_property + def fidelity(self) -> float: + """Product of average gate fidelities over all gate channels in the cycle. + + Measurement channels do not contribute a gate fidelity and are skipped. + """ + f = 1.0 + for ch in self.channels: + if isinstance(ch, Channel): + f *= ch.fidelity + return f + + @cached_property + def infidelity(self) -> float: + """``1 - fidelity``.""" + return 1.0 - self.fidelity + + @cached_property + def pauli_infidelity(self) -> float: + """``1 - pauli_fidelity``.""" + return 1.0 - self.pauli_fidelity + + # ────────────────────────────────────────────── + # Serialization + # ────────────────────────────────────────────── + + def to_json(self) -> str: + """Serialize CycleChannel to a JSON string. + + :return: JSON string representation. + """ + ch_data = [] + for ch in self.channels: + ch_data.append({"type": type(ch).__name__, "data": ch.to_json()}) + data = { + "channels": ch_data, + } + return json.dumps(data) + + @classmethod + def from_json(cls: type[CycleChannel], json_str: str) -> CycleChannel: + """Deserialize a CycleChannel from a JSON string. + + The ``inst`` and ``defcircuit`` fields are reconstructed from the constituent + channels, consistent with how :func:`_build_cycle_channel` builds them. + + :param json_str: JSON string as produced by :meth:`to_json`. + :return: CycleChannel instance. + """ + data = json.loads(json_str) + _type_map: dict[str, type[Channel | MeasurementChannel]] = { + "Channel": Channel, + "MeasurementChannel": MeasurementChannel, + } + constituent_channels: list[Channel | MeasurementChannel] = [ + _type_map[ch_data["type"]].from_json(ch_data["data"]) for ch_data in data["channels"] + ] + return _build_cycle_channel(constituent_channels) + + # ────────────────────────────────────────────── + # Dunder methods + # ────────────────────────────────────────────── + + def __str__(self) -> str: + """Return a simplified string representation showing the gate and process fidelity.""" + return f"<{self.inst.out()} ~ ({100 * self.pauli_fidelity:.2f}%)>" + + def __eq__(self, other: object) -> bool: + """Check equality based on instruction and constituent channels.""" + if not isinstance(other, CycleChannel): + return False + if self.inst != other.inst: + return False + return self.channels == other.channels + + __hash__ = None # type: ignore[assignment] + + +def _channel_to_formal_inst(channel: Channel | MeasurementChannel) -> Gate | Measurement: + """Convert a channel's instruction to use formal arguments for DefCircuit.""" + if isinstance(channel, Channel): + inst = channel.inst + return Gate( + name=inst.name, + params=inst.params, + qubits=[FormalArgument(f"q{q}") for q in inst.get_qubit_indices()], + modifiers=inst.modifiers, # type: ignore[arg-type] + ) + elif isinstance(channel, MeasurementChannel): + qubit_idx = channel.qubits[0] + return Measurement( + qubit=FormalArgument(f"q{qubit_idx}"), + classical_reg=None, + ) + raise TypeError(f"Unsupported channel type: {type(channel)}") + + +def _build_cycle_channel( + channels: list[Channel | MeasurementChannel], +) -> CycleChannel: + """Build a CycleChannel from a list of Channel/MeasurementChannel on disjoint qubits.""" + all_qubits = sorted(q for ch in channels for q in ch.qubits) + cycle_name = "CYCLE" + formal_insts = [_channel_to_formal_inst(ch) for ch in channels] + + defcircuit = DefCircuit( + name=cycle_name, + parameters=[], + qubits=[FormalArgument(f"q{q}") for q in all_qubits], + instructions=list(formal_insts), + ) + inst = Gate(name=cycle_name, params=[], qubits=all_qubits) + return CycleChannel(inst=inst, defcircuit=defcircuit, channels=tuple(channels)) diff --git a/pyquil/noise.py b/pyquil/noise/_legacy_noise.py similarity index 99% rename from pyquil/noise.py rename to pyquil/noise/_legacy_noise.py index 9950e57a5..730a0b8a0 100644 --- a/pyquil/noise.py +++ b/pyquil/noise/_legacy_noise.py @@ -21,6 +21,7 @@ from typing import TYPE_CHECKING, Any, Optional, cast import numpy as np +from deprecated import deprecated from pyquil.external.rpcq import CompilerISA from pyquil.gates import MEASURE, RX, I @@ -41,6 +42,9 @@ class KrausModel(_KrausModel): """Encapsulate a single gate's noise model. + .. deprecated:: + Use :class:`pyquil.noise.Channel` for quax-based noise modeling. + :ivar str gate: The name of the gate. :ivar Sequence[float] params: Optional parameters for the gate. :ivar Sequence[int] targets: The target qubit ids. @@ -117,6 +121,11 @@ def __eq__(self, other: object) -> bool: _NoiseModel = namedtuple("_NoiseModel", ["gates", "assignment_probs"]) +@deprecated( + version="4.17.0", + reason="Use the quax-based noise model in pyquil.noise._noise_model instead. " + "This class will be removed in the next major version of pyquil.", +) class NoiseModel(_NoiseModel): """Encapsulate the QPU noise model containing information about the noisy gates. diff --git a/pyquil/noise/_noise_model.py b/pyquil/noise/_noise_model.py new file mode 100644 index 000000000..5bd1682ae --- /dev/null +++ b/pyquil/noise/_noise_model.py @@ -0,0 +1,489 @@ +############################################################################## +# Copyright 2016-2026 Rigetti Computing +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +############################################################################## +"""Noise model container and program-level fidelity estimation. + +This module defines: + +- ``NoiseModelLike``: A ``typing.Protocol`` defining the interface that all noise + models must satisfy (a single ``get_channel`` method). +- ``NoiseModel``: The primary concrete implementation — a frozen dataclass that + collects per-instruction noise channels. +- ``DepolarizingNoiseModel``: A convenience implementation that returns a + depolarizing channel for any gate. +- ``CompositeNoiseModel``: Chains multiple noise models, returning the first + non-None channel. +- Program-level fidelity estimation utilities. +""" + +from __future__ import annotations + +import json +import logging +from collections.abc import Iterable, Mapping, Sequence +from dataclasses import dataclass +from functools import reduce +from operator import mul +from types import MappingProxyType +from typing import ( + TYPE_CHECKING, + Protocol, + overload, + runtime_checkable, +) + +from pyquil.external.rpcq import CompilerISA +from pyquil.noise._channels import Channel, CycleChannel, MeasurementChannel, ResetChannel +from pyquil.quilbase import Gate, Measurement, ResetQubit + +if TYPE_CHECKING: + from pyquil import Program + from pyquil.paulis import PauliSum, PauliTerm + +logger = logging.getLogger(__name__) + +# ────────────────────────────────────────────────────────── +# Protocol +# ────────────────────────────────────────────────────────── + +# Channel union type returned by get_channel +ChannelType = Channel | MeasurementChannel | ResetChannel | CycleChannel +NoiseInstruction = Gate | Measurement | ResetQubit + + +@runtime_checkable +class NoiseModelLike(Protocol): + """Protocol defining the noise model interface. + + Any object that implements ``get_channel`` with the correct signature can be + used wherever a noise model is expected. This enables alternative noise model + strategies (depolarizing, dynamic, crosstalk-aware) without modifying + consumers. + + The standard concrete implementation is :class:`NoiseModel`. + """ + + @overload + def get_channel(self, inst: Gate) -> Channel | CycleChannel | None: ... + + @overload + def get_channel(self, inst: Measurement) -> MeasurementChannel | None: ... + + @overload + def get_channel(self, inst: ResetQubit) -> ResetChannel | None: ... + + def get_channel(self, inst: Gate | Measurement | ResetQubit) -> ChannelType | None: + """Retrieve the noise channel for a specific instruction. + + :param inst: A gate, measurement, or reset instruction. + :return: The associated noise channel, or ``None`` if the instruction + should be treated as ideal (noiseless). + """ + ... + + +@dataclass(frozen=True) +class NoiseModel: + """A noise model collects all the noise channels for a given quantum program. + + This includes gate channels, measurement channels, reset channels, and cycle channels. + + The constructor accepts a mapping from instruction to channel. Use + :meth:`from_channels` when constructing from a list, tuple, set, generator, + or other channel iterable. + """ + + channels: Mapping[NoiseInstruction, ChannelType] + """Immutable mapping from instruction to its noise channel.""" + + def __init__( + self, + channels: Mapping[NoiseInstruction, ChannelType] | None = None, + ) -> None: + if channels is None: + channels = {} + if not isinstance(channels, Mapping): + raise TypeError("NoiseModel channels must be a mapping. Use NoiseModel.from_channels(...) for iterables.") + + channel_map: dict[NoiseInstruction, ChannelType] = {} + for inst, channel in channels.items(): + if channel.inst != inst: + raise ValueError( + f"NoiseModel channel key {inst!r} does not match channel instruction {channel.inst!r}." + ) + channel_map[inst] = channel + object.__setattr__(self, "channels", MappingProxyType(channel_map)) + + def __getstate__(self) -> dict[str, dict[NoiseInstruction, ChannelType]]: + # ``channels`` is a MappingProxyType (unpicklable); serialize it as a plain dict. + return {"channels": dict(self.channels)} + + def __setstate__(self, state: Mapping[str, Mapping[NoiseInstruction, ChannelType]]) -> None: + object.__setattr__(self, "channels", MappingProxyType(dict(state["channels"]))) + + @overload + def get_channel(self, inst: Gate) -> Channel | CycleChannel | None: ... + + @overload + def get_channel(self, inst: Measurement) -> MeasurementChannel | None: ... + + @overload + def get_channel(self, inst: ResetQubit) -> ResetChannel | None: ... + + def get_channel( + self, inst: Gate | Measurement | ResetQubit + ) -> Channel | MeasurementChannel | ResetChannel | CycleChannel | None: + """Retrieve the noise channel associated with a specific instruction. + + :param inst: The instruction (gate, measurement, or reset) for which to retrieve the noise channel. + :return: The noise channel associated with the instruction, or None if no channel is found. + """ + return self.channels.get(inst) + + # ────────────────────────────────────────────── + # Constructors + # ────────────────────────────────────────────── + + @classmethod + def from_channels(cls: type[NoiseModel], channels: Iterable[ChannelType] = ()) -> NoiseModel: + """Create a noise model from an iterable of channels. + + :param channels: Noise channels to include in the model. + :return: A NoiseModel keyed by each channel's instruction. + :raises ValueError: If more than one channel targets the same instruction. + """ + channel_map: dict[NoiseInstruction, ChannelType] = {} + for channel in channels: + if channel.inst in channel_map: + raise ValueError(f"Duplicate noise channel for instruction {channel.inst!r}.") + channel_map[channel.inst] = channel + return cls(channels=channel_map) + + @classmethod + def from_isa(cls: type[NoiseModel], compiler_isa: CompilerISA) -> NoiseModel: + """Create a noise model from an instruction set architecture. + + Gate fidelities are converted to depolarizing channels and measurement + errors are symmetric. Only gates with concrete numeric parameters are + included. + + .. note:: + Two-qubit gate channels are keyed by the ISA edge's operand order + (e.g. ``CZ 0 1``). A program that issues the gate with the operands + reversed (``CZ 1 0``) will not match this channel, even for symmetric + gates. Issue gates in the ISA operand order, or build channels + explicitly for both orderings. + + :param compiler_isa: The compiler ISA. + :return: A NoiseModel with channels according to the provided fidelities. + """ + from pyquil.external.rpcq import GateInfo, MeasureInfo + from pyquil.quilatom import Qubit as QuilQubit + + channels: dict[NoiseInstruction, ChannelType] = {} + seen_measure_qubits: set[int] = set() + + for qubit_label, qubit in compiler_isa.qubits.items(): + for op_info in qubit.gates: + if isinstance(op_info, GateInfo): + qubits = [int(qubit_label)] + gate_name = op_info.operator + fidelity = op_info.fidelity + params = op_info.parameters + + if gate_name is None: + continue + # Skip gates with non-numeric parameters + if not all(isinstance(p, (float, int, complex)) for p in params): + continue + + numeric_params: list[float] = [float(p) for p in params if isinstance(p, (float, int, complex))] + inst = Gate(name=gate_name, params=numeric_params, qubits=qubits) + if fidelity is not None and fidelity < 1.0: + channels[inst] = Channel.from_gate_fidelity(inst=inst, fidelity=fidelity) + + elif isinstance(op_info, MeasureInfo): + if op_info.qubit is None: + continue + # Use qubit_label from the enclosing section when qubit is a wildcard + qubit_str = op_info.qubit if op_info.qubit != "_" else qubit_label + try: + qubit_idx = int(qubit_str) + except (ValueError, TypeError): + continue + fidelity = op_info.fidelity + if qubit_idx in seen_measure_qubits: + continue + if fidelity is None: + # Don't mark the qubit seen on a fidelity-less entry; a later + # MeasureInfo for the same qubit may carry a usable fidelity. + continue + seen_measure_qubits.add(qubit_idx) + m_inst = Measurement(qubit=QuilQubit(qubit_idx), classical_reg=None) + channels[m_inst] = MeasurementChannel.from_readout_fidelity(inst=m_inst, fidelity=fidelity) + + for edge_label, edge in compiler_isa.edges.items(): + for op_info in edge.gates: + if isinstance(op_info, GateInfo): + qubits = [int(q) for q in edge_label.split("-")] + gate_name = op_info.operator + fidelity = op_info.fidelity + params = op_info.parameters + + if gate_name is None: + continue + if not all(isinstance(p, (float, int, complex)) for p in params): + continue + + numeric_params = [float(p) for p in params if isinstance(p, (float, int, complex))] + inst = Gate(name=gate_name, params=numeric_params, qubits=qubits) + if fidelity is not None and fidelity < 1.0: + channels[inst] = Channel.from_gate_fidelity(inst=inst, fidelity=fidelity) + + return cls(channels=channels) + + # ────────────────────────────────────────────── + # Serialization + # ────────────────────────────────────────────── + + def to_json(self) -> str: + """Serialize NoiseModel to a JSON string. + + :return: JSON string representation. + """ + channel_data = [] + for ch in self.channels.values(): + if isinstance(ch, (Channel, MeasurementChannel, ResetChannel, CycleChannel)): + channel_data.append({"type": type(ch).__name__, "data": ch.to_json()}) + else: + logger.warning(f"Skipping serialization of {type(ch).__name__} (not yet supported).") + return json.dumps({"channels": channel_data}) + + @classmethod + def from_json(cls: type[NoiseModel], json_str: str) -> NoiseModel: + """Deserialize a NoiseModel from a JSON string. + + :param json_str: JSON string as produced by :meth:`to_json`. + :return: NoiseModel instance. + """ + data = json.loads(json_str) + _type_map = { + "Channel": Channel, + "MeasurementChannel": MeasurementChannel, + "ResetChannel": ResetChannel, + "CycleChannel": CycleChannel, + } + channels: list[Channel | MeasurementChannel | ResetChannel | CycleChannel] = [] + for ch_data in data["channels"]: + ch_cls = _type_map.get(ch_data["type"]) + if ch_cls is None: + raise ValueError(f"Unknown channel type: {ch_data['type']}") + channels.append(ch_cls.from_json(ch_data["data"])) # type: ignore[attr-defined] + return cls.from_channels(channels) + + # ────────────────────────────────────────────── + # Dunder methods + # ────────────────────────────────────────────── + + def __eq__(self, other: object) -> bool: + """Check if two NoiseModels contain equivalent channel maps.""" + if not isinstance(other, NoiseModel): + return False + return dict(self.channels) == dict(other.channels) + + # Unhashable: its channels hold jax arrays and are themselves unhashable, and an + # ``id``-based hash would be inconsistent with the value-based ``__eq__``. + __hash__ = None # type: ignore[assignment] + + def __add__(self, other: NoiseModel) -> NoiseModel: + """Combine two NoiseModels into their disjoint union. + + The two models must not both define a channel for the same instruction. + Addition is a union of channels, not a composition: overlapping channels are + a conflict, not an "addition". Compose channels explicitly + (``channel_a @ channel_b``) if that is what you intend. + + :raises ValueError: If both models define a channel for the same instruction. + """ + if not isinstance(other, NoiseModel): + return NotImplemented + + overlap = set(self.channels) & set(other.channels) + if overlap: + insts = ", ".join(repr(inst) for inst in overlap) + raise ValueError( + f"Cannot add NoiseModels: both define a channel for the same instruction(s): {insts}. " + "Addition is a disjoint union; compose the channels explicitly if that is intended." + ) + + return NoiseModel(channels={**self.channels, **other.channels}) + + def with_channels(self, channels: Iterable[ChannelType]) -> NoiseModel: + """Return a new model with additional channels. + + :param channels: New channels to add. + :return: A NoiseModel containing the existing and new channels. + :raises ValueError: If a new channel targets an instruction already in the model. + """ + combined = dict(self.channels) + for channel in channels: + if channel.inst in combined: + raise ValueError(f"Duplicate noise channel for instruction {channel.inst!r}.") + combined[channel.inst] = channel + return NoiseModel(channels=combined) + + +# ────────────────────────────────────────────────────────── +# Convenience NoiseModelLike implementations +# ────────────────────────────────────────────────────────── + + +@dataclass(frozen=True) +class DepolarizingNoiseModel: + r"""A noise model that applies uniform depolarizing noise to every gate. + + For any ``Gate`` instruction, returns a :class:`Channel` with the specified + depolarizing constant. Measurements and resets are treated as ideal. + + :param depolarizing_constant: The depolarization constant :math:`p` where + :math:`\\mathcal{D}_p(\\rho) = p \\, \\rho + (1-p) \\, I/d`. + A value of 1.0 means no noise; 0.0 means full depolarization. + """ + + depolarizing_constant: float + + @overload + def get_channel(self, inst: Gate) -> Channel | CycleChannel | None: ... + + @overload + def get_channel(self, inst: Measurement) -> MeasurementChannel | None: ... + + @overload + def get_channel(self, inst: ResetQubit) -> ResetChannel | None: ... + + def get_channel(self, inst: Gate | Measurement | ResetQubit) -> ChannelType | None: + """Return a depolarizing channel for gates; ``None`` for measurements/resets.""" + if isinstance(inst, Gate): + return Channel.from_depolarizing_constant(inst, self.depolarizing_constant) + return None + + +@dataclass(frozen=True) +class CompositeNoiseModel: + """A noise model that chains multiple models, returning the first non-None channel. + + Models are queried in order. The first model that returns a non-None channel + for a given instruction wins. + + :param models: Sequence of noise models to query in priority order. + """ + + models: tuple[NoiseModelLike, ...] + + def __init__(self, models: Sequence[NoiseModelLike]) -> None: + object.__setattr__(self, "models", tuple(models)) + + @overload + def get_channel(self, inst: Gate) -> Channel | CycleChannel | None: ... + + @overload + def get_channel(self, inst: Measurement) -> MeasurementChannel | None: ... + + @overload + def get_channel(self, inst: ResetQubit) -> ResetChannel | None: ... + + def get_channel(self, inst: Gate | Measurement | ResetQubit) -> ChannelType | None: + """Query each model in order, returning the first non-None result.""" + for model in self.models: + channel = model.get_channel(inst) + if channel is not None: + return channel + return None + + +# ────────────────────────────────────────────────────────── +# Program-level fidelity estimation +# ────────────────────────────────────────────────────────── + + +def estimate_program_fidelity(program: Program, noise_model: NoiseModelLike) -> float: + """Estimate the program fidelity for a given noise model. + + Works by multiplying the gate process fidelities together. Readout noise + is not considered. + + :param program: The program of interest. + :param noise_model: A noise model. + :return: The estimated process fidelity. + """ + gate_fidelities = [1.0] + for inst in program.instructions: + if isinstance(inst, Gate): + channel = noise_model.get_channel(inst) + if isinstance(channel, (Channel, CycleChannel)): + gate_fidelities.append(channel.pauli_fidelity) + + return reduce(mul, gate_fidelities) + + +def _light_cone_program(program: Program, qubits: list[int]) -> Program: + """Return a sub-program containing only gates in the backward light cone of *qubits*. + + Walks backward through the program's gate instructions. Any gate that + acts on a qubit currently in the light-cone set is included, and all of + its qubits are added to the set (because earlier gates on those qubits + are now causally relevant). + """ + from pyquil import Program as _Program + + gate_instructions = [inst for inst in program.instructions if isinstance(inst, Gate)] + relevant_qubits = set(qubits) + included: list[Gate] = [] + for inst in reversed(gate_instructions): + inst_qubits = {q.index for q in inst.qubits} # type: ignore[union-attr] + if inst_qubits & relevant_qubits: + included.append(inst) + relevant_qubits |= inst_qubits + reduced = _Program() + for inst in reversed(included): + reduced += inst + return reduced + + +def estimate_program_observable_fidelity( + program: Program, + noise_model: NoiseModelLike, + observable: PauliSum | PauliTerm, +) -> float: + """Estimate program fidelity restricted to the backward light cone of *observable*. + + Reduces the program to only the gates causally connected to the + observable qubits, then multiplies gate process fidelities together. + Readout noise is not considered. + + :param program: The program of interest. + :param noise_model: A noise model. + :param observable: A ``PauliTerm`` or ``PauliSum`` whose qubits define + the light cone. + :return: The estimated process fidelity for the light-cone-reduced program. + """ + from pyquil.paulis import PauliSum, PauliTerm + + if isinstance(observable, PauliTerm): + observable = PauliSum(terms=[observable]) + + qubits = [int(q) for term in observable.terms for q, _ in term.operations_as_set()] # type: ignore[arg-type] + reduced_program = _light_cone_program(program, qubits) + return estimate_program_fidelity(reduced_program, noise_model) diff --git a/pyquil/quilbase.py b/pyquil/quilbase.py index 0bbe39fb8..5f15345b6 100644 --- a/pyquil/quilbase.py +++ b/pyquil/quilbase.py @@ -56,10 +56,37 @@ if TYPE_CHECKING: # avoids circular import from pyquil.paulis import PauliSum +import math + import quil.expression as quil_rs_expr import quil.instructions as quil_rs +def _is_perfect_power(n: int) -> bool: + """Check if n is a prime power (p^k for prime p, k >= 1). + + This ensures the matrix dimension can be interpreted as k qudits of + dimension p. Composite non-prime-power dimensions like 6 = 2*3 are + ambiguous and rejected. + """ + if n < 2: + return False + # Find the smallest prime factor. + factor = 0 + for p in range(2, int(math.isqrt(n)) + 1): + if n % p == 0: + factor = p + break + if factor == 0: + # n is prime → n = n^1, valid single-qudit dimension. + return True + # Check that n is a power of this smallest prime factor. + val = factor + while val < n: + val *= factor + return val == n + + class _InstructionMeta(abc.ABCMeta): """A metaclass that allows us to group all instruction types from quil-rs and pyQuil as an `AbstractInstruction`. @@ -706,8 +733,10 @@ def _validate_matrix(matrix: list[list[Expression]] | np.ndarray | np.matrix, co else: raise TypeError("Matrix argument must be a list or NumPy array/matrix") - if 0 != rows & (rows - 1): - raise ValueError(f"Dimension of matrix must be a power of 2, got {rows}") + if not _is_perfect_power(rows): + raise ValueError( + f"Dimension of matrix must be a perfect power of an integer (e.g. 2, 3, 4, 8, 9, ...), got {rows}" + ) if not contains_parameters: np_matrix = np.asarray(matrix) @@ -732,9 +761,19 @@ def get_constructor(self) -> Callable[..., Gate] | Callable[..., Callable[..., G return lambda *qubits: Gate(name=self.name, params=[], qubits=list(map(unpack_qubit, qubits))) def num_args(self) -> int: - """Get the number of qubit arguments the gate takes.""" + """Get the number of qudit arguments the gate takes. + + For a matrix of dimension d^k, returns k where d is the smallest + integer base >= 2 such that rows = d^k. + """ rows = len(self.matrix) - return int(np.log2(rows)) + if rows < 2: + return 0 + for base in range(2, rows + 1): + k = int(round(math.log(rows, base))) + if base**k == rows: + return k + return 1 @property def matrix(self) -> np.ndarray: diff --git a/test/unit/conftest.py b/test/unit/conftest.py index fcb8d236a..445c6562c 100644 --- a/test/unit/conftest.py +++ b/test/unit/conftest.py @@ -1,6 +1,10 @@ import os from typing import Any, Dict +import jax + +jax.config.update("jax_enable_x64", True) + import numpy as np import pytest from qcs_sdk import QCSClient diff --git a/test/unit/test_noise.py b/test/unit/test_noise.py index 5a4861785..d3ab8255a 100644 --- a/test/unit/test_noise.py +++ b/test/unit/test_noise.py @@ -1,3 +1,7 @@ +# TODO(pyquil-5.0): Remove this file. These tests exercise the legacy KrausModel/NoiseModel +# API in pyquil.noise, which is superseded by pyquil.noise._channels / pyquil.noise._noise_model. +# Equivalent coverage lives in test_noise_model.py. + from collections import OrderedDict import numpy as np diff --git a/test/unit/test_noise_model.py b/test/unit/test_noise_model.py new file mode 100644 index 000000000..1da9848c6 --- /dev/null +++ b/test/unit/test_noise_model.py @@ -0,0 +1,644 @@ +# Copyright 2024-2026 Rigetti Computing +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +"""Unit tests for the quax-based noise model (Channel, MeasurementChannel, ResetChannel, NoiseModel).""" + +import jax.numpy as jnp +import numpy as np +import pytest +import quax as qx + +from pyquil.external.rpcq import CompilerISA +from pyquil.gates import CNOT, MEASURE, RESET, RX, RY, RZ, X +from pyquil.noise._channels import ( + Channel, + CycleChannel, + MeasurementChannel, + ResetChannel, + _build_cycle_channel, + get_instruction_unitary, +) +from pyquil.noise._noise_model import NoiseModel +from pyquil.quil import Program +from pyquil.quilatom import FormalArgument, Qubit +from pyquil.quilbase import DefCircuit, Gate, Measurement, ResetQubit + +# ────────────────────────────────────────────────────────── +# Channel tests +# ────────────────────────────────────────────────────────── + + +class TestChannel: + def test_from_depolarizing_constant(self): + """Channel.from_depolarizing_constant produces valid superoperator.""" + inst = RX(np.pi / 4, 0) + ch = Channel.from_depolarizing_constant(inst=inst, depolarizing_constant=0.98) + assert isinstance(ch.process, qx.SuperOp) + # Process fidelity should be close to the depolarizing constant + assert ch.fidelity < 1.0 + assert ch.fidelity > 0.95 + + def test_from_gate_fidelity(self): + """Channel.from_gate_fidelity produces correct fidelity.""" + inst = RX(np.pi / 2, 0) + ch = Channel.from_gate_fidelity(inst=inst, fidelity=0.99) + assert abs(ch.fidelity - 0.99) < 0.001 + + def test_from_pauli_fidelity(self): + """Channel.from_pauli_fidelity produces a valid channel.""" + inst = X(0) + ch = Channel.from_pauli_fidelity(inst=inst, pauli_fidelity=0.97) + assert isinstance(ch.process, qx.SuperOp) + assert ch.pauli_fidelity == pytest.approx(0.97, abs=0.001) + + def test_from_pauli_noise(self): + """Channel.from_pauli_noise produces a valid Pauli noise channel.""" + inst = RX(0.5, 0) + ch = Channel.from_pauli_noise(inst=inst, pauli_noise={"X": 0.01, "Z": 0.02}) + assert isinstance(ch.process, qx.SuperOp) + assert ch.fidelity < 1.0 + + def test_from_coherence_times(self): + """Channel.from_coherence_times produces a valid decoherence channel.""" + inst = RX(np.pi / 4, 0) + ch = Channel.from_coherence_times(inst=inst, gate_duration=40e-9, t1s=[30e-6], t2s=[20e-6]) + assert isinstance(ch.process, qx.SuperOp) + assert ch.fidelity < 1.0 + assert ch.fidelity > 0.99 # short gate relative to T1/T2 + + def test_qubits(self): + """Channel.qubits reflects the instruction's qubits.""" + ch = Channel.from_depolarizing_constant(inst=RX(0.1, 3), depolarizing_constant=0.99) + assert ch.qubits == [3] + + def test_num_qubits(self): + """Channel.num_qubits is correct for 2Q gates.""" + ch = Channel.from_depolarizing_constant(inst=CNOT(0, 1), depolarizing_constant=0.95) + assert ch.num_qubits == 2 + + def test_fidelity_properties(self): + """Fidelity, infidelity, pauli_fidelity, pauli_infidelity are consistent.""" + ch = Channel.from_gate_fidelity(inst=RX(0.3, 0), fidelity=0.98) + assert ch.fidelity + ch.infidelity == pytest.approx(1.0) + assert ch.pauli_fidelity + ch.pauli_infidelity == pytest.approx(1.0) + assert ch.stochastic_fidelity + ch.stochastic_infidelity == pytest.approx(1.0) + assert ch.coherent_fidelity + ch.coherent_infidelity == pytest.approx(1.0) + + def test_noise_process(self): + """noise_process factors out the ideal gate unitary.""" + ch = Channel.from_depolarizing_constant(inst=RX(np.pi / 4, 0), depolarizing_constant=0.99) + noise = ch.noise_process + assert isinstance(noise, qx.SuperOp) + + def test_is_pauli(self): + """A depolarizing channel on a Clifford gate should be a Pauli channel.""" + ch = Channel.from_depolarizing_constant(inst=X(0), depolarizing_constant=0.98) + assert ch.is_pauli() + + def test_pauli_twirl(self): + """Pauli twirl of a channel on a Clifford gate should be a Pauli channel.""" + ch = Channel.from_random_coherent_error(inst=X(0), process_fidelity=0.97, rng=np.random.default_rng(42)) + twirled = ch.pauli_twirl() + assert twirled.is_pauli() + + def test_unitarity(self): + """A depolarizing channel should have unitarity < 1.""" + ch = Channel.from_depolarizing_constant(inst=RX(0.5, 0), depolarizing_constant=0.95) + assert 0 < ch.unitarity < 1.0 + + def test_pauli_vector_sums_to_one(self): + """Pauli error probability vector should sum to 1.""" + ch = Channel.from_depolarizing_constant(inst=X(0), depolarizing_constant=0.97) + pv = ch.pauli_vector + assert float(jnp.sum(pv)) == pytest.approx(1.0, abs=1e-8) + + def test_perfect_channel(self): + """A depolarizing constant of 1.0 should give fidelity 1.0.""" + ch = Channel.from_depolarizing_constant(inst=RX(np.pi, 0), depolarizing_constant=1.0) + assert ch.fidelity == pytest.approx(1.0, abs=1e-10) + + def test_from_pauli_noise_rejects_invalid_probabilities(self): + """Pauli error rates must be probabilities with total error no greater than 1.""" + with pytest.raises(ValueError, match="negative"): + Channel.from_pauli_noise(inst=RX(0.5, 0), pauli_noise={"X": -0.1}) + + with pytest.raises(ValueError, match="at most 1.0"): + Channel.from_pauli_noise(inst=RX(0.5, 0), pauli_noise={"X": 0.6, "Z": 0.5}) + + def test_from_pauli_noise_two_qubit(self): + """from_pauli_noise builds the correct 16-term Pauli channel for a 2Q gate (regression).""" + pauli_noise = {"IX": 0.01, "XI": 0.005, "ZZ": 0.02} + ch = Channel.from_pauli_noise(inst=CNOT(0, 1), pauli_noise=pauli_noise) + pv = np.asarray(ch.pauli_vector) + assert pv.size == 16 + assert float(jnp.sum(ch.pauli_vector)) == pytest.approx(1.0, abs=1e-3) + terms = [a + b for a in "IXYZ" for b in "IXYZ"] + rates = dict(zip(terms, pv, strict=True)) + for term, rate in pauli_noise.items(): + assert rates[term] == pytest.approx(rate, abs=1e-3) + assert rates["II"] == pytest.approx(1.0 - sum(pauli_noise.values()), abs=1e-3) + + def test_pow_scales_noise(self): + """Channel ** power scales the noise while preserving the gate.""" + ch = Channel.from_depolarizing_constant(inst=RX(np.pi / 2, 0), depolarizing_constant=0.99) + assert (ch**0.0).pauli_infidelity == pytest.approx(0.0, abs=1e-3) + assert (ch**1.0).pauli_infidelity == pytest.approx(ch.pauli_infidelity, abs=1e-3) + assert (ch**2.0).pauli_infidelity > ch.pauli_infidelity + # The ideal gate is preserved. + assert (ch**2.0).qubits == ch.qubits + assert jnp.allclose((ch**2.0).target_unitary.matrix, ch.target_unitary.matrix) + + def test_json_roundtrip_preserves_qutrit_dims(self): + """Channel JSON includes explicit dims for non-qubit operators.""" + qutrit_x = jnp.array( + [ + [0.0, 0.0, 1.0], + [1.0, 0.0, 0.0], + [0.0, 1.0, 0.0], + ], + dtype=complex, + ) + target_unitary = qx.Unitary.from_matrix(qutrit_x, ((3,), (3,))) + channel = Channel( + inst=Gate("TX", [], [0]), process=qx.to_superop(target_unitary), target_unitary=target_unitary + ) + + restored = Channel.from_json(channel.to_json()) + + assert restored.inst == channel.inst + assert restored.process.dims == ((3,), (3,)) + assert restored.target_unitary.dims == ((3,), (3,)) + assert jnp.allclose(restored.process.matrix, channel.process.matrix) + + +# ────────────────────────────────────────────────────────── +# MeasurementChannel tests +# ────────────────────────────────────────────────────────── + + +class TestMeasurementChannel: + def test_from_readout_fidelity(self): + """MeasurementChannel.from_readout_fidelity produces a valid quantum instrument.""" + # Use the pyquil MEASURE gate to get qubit + prog = Program(MEASURE(0, None)) + meas_inst = [i for i in prog if isinstance(i, Measurement)][0] + ch = MeasurementChannel.from_readout_fidelity(inst=meas_inst, fidelity=0.95) + assert isinstance(ch.process, qx.QuantumInstrument) + + def test_from_readout_fidelity_with_asymmetry(self): + """MeasurementChannel with asymmetry produces asymmetric confusion.""" + prog = Program(MEASURE(0, None)) + meas_inst = [i for i in prog if isinstance(i, Measurement)][0] + ch = MeasurementChannel.from_readout_fidelity(inst=meas_inst, fidelity=0.95, asymmetry=0.5) + assert isinstance(ch.process, qx.QuantumInstrument) + + def test_qubits(self): + """MeasurementChannel.qubits returns the correct qubit.""" + prog = Program(MEASURE(5, None)) + meas_inst = [i for i in prog if isinstance(i, Measurement)][0] + ch = MeasurementChannel.from_readout_fidelity(inst=meas_inst, fidelity=0.99) + assert ch.qubits == [5] + + @pytest.mark.parametrize("asymmetry", [0.0, 0.5]) + def test_pow_scales_readout_noise(self, asymmetry): + """MeasurementChannel ** power scales readout noise via the stochastic generator.""" + prog = Program(MEASURE(0, None)) + meas_inst = [i for i in prog if isinstance(i, Measurement)][0] + ch = MeasurementChannel.from_readout_fidelity(inst=meas_inst, fidelity=0.95, asymmetry=asymmetry) + + def bitflip(channel): + cm = np.asarray(channel.confusion_matrix) + return 1.0 - 0.5 * (float(cm[0, 0]) + float(cm[1, 1])) + + assert bitflip(ch**0.0) == pytest.approx(0.0, abs=1e-3) + assert bitflip(ch**1.0) == pytest.approx(bitflip(ch), abs=1e-3) + assert bitflip(ch**2.0) > bitflip(ch) + # The generator construction keeps the result exactly column-stochastic and non-negative. + powered = np.asarray((ch**1.5).confusion_matrix) + assert np.all(powered >= -1e-9) + assert np.allclose(powered.sum(axis=0), 1.0, atol=1e-6) + + def test_pow_rejects_non_embeddable_measurement(self): + """A confusion matrix with no real generator cannot be fractionally powered.""" + prog = Program(MEASURE(0, None)) + meas_inst = [i for i in prog if isinstance(i, Measurement)][0] + # Near-complete bit flip: eigenvalue ~ -0.8, so the matrix is not embeddable. + confusion = jnp.array([[0.1, 0.9], [0.9, 0.1]]) + ch = MeasurementChannel.from_confusion_and_transition(meas_inst, confusion, jnp.eye(2)) + with pytest.raises(ValueError, match="not embeddable|not a valid"): + _ = ch**0.5 + + def test_from_binary_discriminator_qubit_is_faithful_readout(self): + """Regression: dim=2/threshold=1 is a real qubit readout, not an always-0 collapse. + + The previous implementation mapped both |0> and |1> to outcome 0 for a qubit, + silently erasing all measurement information. + """ + prog = Program(MEASURE(0, None)) + meas_inst = [i for i in prog if isinstance(i, Measurement)][0] + ch = MeasurementChannel.from_binary_discriminator(inst=meas_inst, dim=2, threshold=1) + cm = np.asarray(ch.confusion_matrix) + assert cm.shape == (2, 2) # two reachable outcomes + assert np.allclose(cm, np.eye(2)) # |0> -> 0, |1> -> 1 + + def test_from_binary_discriminator_qutrit_split(self): + """dim=3 splits into exactly two outcomes at the threshold.""" + prog = Program(MEASURE(0, None)) + meas_inst = [i for i in prog if isinstance(i, Measurement)][0] + # threshold=2: {0,1} -> 0, {2} -> 1 (flag leakage only) + cm2 = np.asarray(MeasurementChannel.from_binary_discriminator(inst=meas_inst, dim=3, threshold=2).confusion_matrix) + assert np.allclose(cm2, [[1, 1, 0], [0, 0, 1]]) + # threshold=1: {0} -> 0, {1,2} -> 1 (ground vs excited-or-leaked) + cm1 = np.asarray(MeasurementChannel.from_binary_discriminator(inst=meas_inst, dim=3, threshold=1).confusion_matrix) + assert np.allclose(cm1, [[1, 0, 0], [0, 1, 1]]) + + def test_from_binary_discriminator_fidelity_degrades(self): + """Sub-unit fidelity stays column-stochastic and keeps both outcomes reachable.""" + prog = Program(MEASURE(0, None)) + meas_inst = [i for i in prog if isinstance(i, Measurement)][0] + cm = np.asarray( + MeasurementChannel.from_binary_discriminator(inst=meas_inst, dim=2, threshold=1, fidelity=0.9).confusion_matrix + ) + assert cm.shape == (2, 2) + assert np.allclose(cm.sum(axis=0), 1.0) + assert cm[1, 1] > cm[0, 1] # |1> still most likely reads as outcome 1 + + def test_from_binary_discriminator_rejects_bad_threshold(self): + """threshold must satisfy 1 <= threshold < dim.""" + prog = Program(MEASURE(0, None)) + meas_inst = [i for i in prog if isinstance(i, Measurement)][0] + with pytest.raises(ValueError, match="threshold"): + MeasurementChannel.from_binary_discriminator(inst=meas_inst, dim=2, threshold=2) + + def test_json_roundtrip_preserves_qutrit_dims(self): + """MeasurementChannel JSON includes explicit dims for non-qubit instruments.""" + prog = Program(MEASURE(0, None)) + meas_inst = [i for i in prog if isinstance(i, Measurement)][0] + channel = MeasurementChannel.from_readout_fidelity(inst=meas_inst, fidelity=0.95, dim=3) + + restored = MeasurementChannel.from_json(channel.to_json()) + + assert restored.inst == channel.inst + assert restored.process.dims == channel.process.dims + assert restored.process.measured_qudits == channel.process.measured_qudits + assert jnp.allclose(restored.process.matrix, channel.process.matrix) + + +# ────────────────────────────────────────────────────────── +# NoiseModel tests +# ────────────────────────────────────────────────────────── + + +class TestNoiseModel: + def test_empty_model(self): + """An empty NoiseModel has no channels.""" + nm = NoiseModel() + assert nm.get_channel(RX(0.5, 0)) is None + + def test_constructor_accepts_instruction_mapping(self): + """NoiseModel stores channels keyed by instruction.""" + inst = RX(np.pi / 4, 0) + ch = Channel.from_depolarizing_constant(inst=inst, depolarizing_constant=0.98) + nm = NoiseModel(channels={inst: ch}) + assert nm.channels[inst] is ch + assert nm.get_channel(inst) is ch + + def test_constructor_rejects_channel_iterable(self): + """Sequence construction should go through from_channels.""" + inst = RX(np.pi / 4, 0) + ch = Channel.from_depolarizing_constant(inst=inst, depolarizing_constant=0.98) + with pytest.raises(TypeError, match="from_channels"): + NoiseModel(channels=[ch]) # type: ignore[arg-type] + + def test_pickle_roundtrip(self): + """NoiseModel survives pickling (its MappingProxyType channels would otherwise block it). + + This is what lets a model be shipped to multiprocessing workers. + """ + import pickle + + gate = RX(np.pi / 4, 0) + prog = Program(MEASURE(0, None)) + meas_inst = [i for i in prog if isinstance(i, Measurement)][0] + nm = NoiseModel.from_channels( + [ + Channel.from_depolarizing_constant(inst=gate, depolarizing_constant=0.98), + MeasurementChannel.from_readout_fidelity(inst=meas_inst, fidelity=0.95), + ] + ) + restored = pickle.loads(pickle.dumps(nm)) + assert set(restored.channels) == set(nm.channels) + assert isinstance(restored.get_channel(gate), Channel) + # Channels survive the round-trip intact (exercises value-based __eq__). + assert restored == nm + + def test_constructor_rejects_mismatched_mapping_key(self): + """Mapping keys must match the instruction stored on each channel.""" + inst = RX(np.pi / 4, 0) + ch = Channel.from_depolarizing_constant(inst=inst, depolarizing_constant=0.98) + with pytest.raises(ValueError, match="does not match"): + NoiseModel(channels={RY(np.pi / 2, 0): ch}) + + def test_from_channels_rejects_duplicates(self): + """Duplicate instruction channels are ambiguous and rejected.""" + inst = RX(np.pi / 4, 0) + ch1 = Channel.from_depolarizing_constant(inst=inst, depolarizing_constant=0.98) + ch2 = Channel.from_depolarizing_constant(inst=inst, depolarizing_constant=0.97) + with pytest.raises(ValueError, match="Duplicate noise channel"): + NoiseModel.from_channels([ch1, ch2]) + + def test_get_channel_gate(self): + """NoiseModel.get_channel returns the correct Channel for a gate.""" + inst = RX(np.pi / 4, 0) + ch = Channel.from_depolarizing_constant(inst=inst, depolarizing_constant=0.98) + nm = NoiseModel.from_channels([ch]) + retrieved = nm.get_channel(inst) + assert retrieved is ch + + def test_get_channel_returns_none_for_missing(self): + """get_channel returns None for instructions not in the model.""" + inst = RX(np.pi / 4, 0) + ch = Channel.from_depolarizing_constant(inst=inst, depolarizing_constant=0.98) + nm = NoiseModel.from_channels([ch]) + other_inst = RY(np.pi / 2, 1) + assert nm.get_channel(other_inst) is None + + def test_multiple_channels(self): + """NoiseModel with multiple channels retrieves each correctly.""" + inst1 = RX(0.5, 0) + inst2 = RY(0.3, 1) + inst3 = CNOT(0, 1) + ch1 = Channel.from_depolarizing_constant(inst=inst1, depolarizing_constant=0.99) + ch2 = Channel.from_depolarizing_constant(inst=inst2, depolarizing_constant=0.97) + ch3 = Channel.from_depolarizing_constant(inst=inst3, depolarizing_constant=0.95) + nm = NoiseModel.from_channels([ch1, ch2, ch3]) + assert nm.get_channel(inst1) is ch1 + assert nm.get_channel(inst2) is ch2 + assert nm.get_channel(inst3) is ch3 + + def test_json_roundtrip(self): + """NoiseModel JSON keeps the existing channel-list wire format.""" + ch = Channel.from_depolarizing_constant(inst=RX(np.pi / 4, 0), depolarizing_constant=0.98) + meas_ch = MeasurementChannel.from_readout_fidelity(inst=MEASURE(1, None), fidelity=0.95) + nm = NoiseModel.from_channels([ch, meas_ch]) + + restored = NoiseModel.from_json(nm.to_json()) + + assert restored == nm + assert set(restored.channels) == {ch.inst, meas_ch.inst} + + def test_add_combines_disjoint_channels(self): + """NoiseModel addition preserves disjoint channels from both operands.""" + ch1 = Channel.from_depolarizing_constant(inst=RX(np.pi / 4, 0), depolarizing_constant=0.98) + ch2 = Channel.from_depolarizing_constant(inst=RY(np.pi / 4, 1), depolarizing_constant=0.97) + + combined = NoiseModel.from_channels([ch1]) + NoiseModel.from_channels([ch2]) + + assert combined.get_channel(ch1.inst) == ch1 + assert combined.get_channel(ch2.inst) == ch2 + + def test_add_rejects_overlapping_channels(self): + """Addition is a disjoint union; a shared instruction is a conflict, not a composition.""" + inst = RX(np.pi / 4, 0) + ch1 = Channel.from_depolarizing_constant(inst=inst, depolarizing_constant=0.98) + ch2 = Channel.from_depolarizing_constant(inst=inst, depolarizing_constant=0.97) + + with pytest.raises(ValueError, match="same instruction"): + _ = NoiseModel.from_channels([ch1]) + NoiseModel.from_channels([ch2]) + + def test_with_channels_returns_extended_model(self): + """with_channels returns a new model and rejects duplicate instructions.""" + ch1 = Channel.from_depolarizing_constant(inst=RX(np.pi / 4, 0), depolarizing_constant=0.98) + ch2 = Channel.from_depolarizing_constant(inst=RY(np.pi / 4, 1), depolarizing_constant=0.97) + nm = NoiseModel.from_channels([ch1]) + + extended = nm.with_channels([ch2]) + + assert nm.get_channel(ch2.inst) is None + assert extended.get_channel(ch1.inst) is ch1 + assert extended.get_channel(ch2.inst) is ch2 + with pytest.raises(ValueError, match="Duplicate noise channel"): + nm.with_channels([ch1]) + + def test_noise_model_is_unhashable(self): + """NoiseModel is unhashable (consistent with value-based equality).""" + nm = NoiseModel.from_channels([Channel.from_depolarizing_constant(RX(0.5, 0), 0.98)]) + with pytest.raises(TypeError): + hash(nm) + + +# ────────────────────────────────────────────────────────── +# get_instruction_unitary tests +# ────────────────────────────────────────────────────────── + + +class TestGetInstructionUnitary: + def test_standard_gate(self): + """get_instruction_unitary resolves standard gates.""" + u = get_instruction_unitary(RX(np.pi / 2, 0)) + assert isinstance(u, qx.Unitary) + assert u.matrix.shape == (2, 2) + + def test_two_qubit_gate(self): + """get_instruction_unitary resolves 2Q gates.""" + u = get_instruction_unitary(CNOT(0, 1)) + assert isinstance(u, qx.Unitary) + assert u.matrix.shape == (4, 4) + + def test_custom_gate(self): + """get_instruction_unitary resolves custom gates from custom_gates dict.""" + custom_matrix = np.array([[0, 1], [1, 0]], dtype=complex) + inst = Gate("MY_GATE", [], [0]) + u = get_instruction_unitary(inst, custom_gates={"MY_GATE": qx.Unitary.from_matrix(custom_matrix, ((2,), (2,)))}) + assert isinstance(u, qx.Unitary) + assert np.allclose(np.asarray(u.matrix), custom_matrix) + + +# ────────────────────────────────────────────────────────── +# ResetChannel tests +# ────────────────────────────────────────────────────────── + + +class TestResetChannel: + def test_from_reset_fidelity_perfect(self): + """A perfect reset (fidelity=1.0) should produce a valid superoperator.""" + inst = RESET(0) + ch = ResetChannel.from_reset_fidelity(inst=inst, fidelity=1.0) + assert isinstance(ch.process, qx.SuperOp) + + def test_from_reset_fidelity_noisy(self): + """A noisy reset should have lower fidelity than a perfect one.""" + inst = RESET(0) + ch_perfect = ResetChannel.from_reset_fidelity(inst=inst, fidelity=1.0) + ch_noisy = ResetChannel.from_reset_fidelity(inst=inst, fidelity=0.95) + assert isinstance(ch_noisy.process, qx.SuperOp) + assert ch_noisy.fidelity < ch_perfect.fidelity + + def test_qubits(self): + """ResetChannel.qubits returns the correct qubit.""" + inst = RESET(3) + ch = ResetChannel.from_reset_fidelity(inst=inst, fidelity=0.99) + assert ch.qubits == [3] + + def test_global_reset_channel_rejected(self): + """ResetChannel is intentionally scoped to targeted resets.""" + with pytest.raises(TypeError, match="targeted"): + ResetChannel.from_reset_fidelity(inst=RESET(), fidelity=1.0) + + def test_json_roundtrip_preserves_qutrit_dims(self): + """ResetChannel JSON includes explicit dims for non-qubit resets.""" + channel = ResetChannel.from_reset_fidelity(inst=ResetQubit(0), fidelity=0.9, dim=3) + + restored = ResetChannel.from_json(channel.to_json()) + + assert restored.inst == channel.inst + assert restored.process.dims == ((3,), (3,)) + assert jnp.allclose(restored.process.matrix, channel.process.matrix) + + +# ────────────────────────────────────────────────────────── +# Channel equality / hashing semantics +# ────────────────────────────────────────────────────────── + + +class TestChannelEqualityAndHashing: + def test_channels_are_unhashable(self): + """Channels hold jax arrays and are intentionally unhashable.""" + ch = Channel.from_depolarizing_constant(RX(0.5, 0), 0.98) + with pytest.raises(TypeError): + hash(ch) + + def test_channel_equality_is_exact(self): + """Channel equality is exact: identical builds are equal, different ones are not.""" + ch1 = Channel.from_depolarizing_constant(RX(0.5, 0), 0.98) + ch2 = Channel.from_depolarizing_constant(RX(0.5, 0), 0.98) + ch3 = Channel.from_depolarizing_constant(RX(0.5, 0), 0.97) + assert ch1 == ch2 + assert ch1 != ch3 + assert ch1 != "not a channel" + + def test_channel_inequality_on_different_instruction(self): + """Channels on different instructions are never equal.""" + ch_a = Channel.from_depolarizing_constant(RX(0.5, 0), 0.98) + ch_b = Channel.from_depolarizing_constant(RX(0.5, 1), 0.98) + assert ch_a != ch_b + + +# ────────────────────────────────────────────────────────── +# Channel construction / analysis coverage +# ────────────────────────────────────────────────────────── + + +class TestChannelAnalysis: + def test_from_mixture(self): + """from_mixture builds a noisy channel from unitary errors with probabilities.""" + z = qx.Unitary.from_matrix(jnp.array([[1, 0], [0, -1]], dtype=complex), ((2,), (2,))) + ch = Channel.from_mixture(X(0), constituents=[z], probabilities=[0.1]) + assert isinstance(ch.process, qx.SuperOp) + assert ch.pauli_infidelity > 0.0 + + def test_coherent_and_stochastic_decomposition(self): + """to_coherent_channel / to_stochastic_channel split the noise into components.""" + ch = Channel.from_random_coherent_error(X(0), process_fidelity=0.95, rng=np.random.default_rng(0)) + coherent = ch.to_coherent_channel() + stochastic = ch.to_stochastic_channel() + assert isinstance(coherent, Channel) + assert isinstance(stochastic, Channel) + # Coherent + stochastic infidelity decomposition is non-negative and finite. + assert np.isfinite(ch.coherent_infidelity) + assert np.isfinite(ch.stochastic_infidelity) + + def test_pauli_twirl_is_pauli(self): + """Twirling a coherent-error channel yields a stochastic Pauli channel.""" + ch = Channel.from_random_coherent_error(X(0), process_fidelity=0.95, rng=np.random.default_rng(1)) + twirled = ch.pauli_twirl() + assert twirled.is_pauli() + + +# ────────────────────────────────────────────────────────── +# NoiseModel.from_isa +# ────────────────────────────────────────────────────────── + + +class TestFromIsa: + @staticmethod + def _isa() -> CompilerISA: + return CompilerISA.parse_obj( + { + "1Q": { + "0": { + "id": 0, + "gates": [ + { + "operator_type": "gate", + "operator": "RX", + "parameters": [1.5707963267948966], + "arguments": ["_"], + "fidelity": 0.99, + }, + # A fidelity-less measurement entry must not mask the real one below. + {"operator_type": "measure", "qubit": "0", "fidelity": None}, + {"operator_type": "measure", "qubit": "0", "fidelity": 0.95}, + ], + }, + "1": {"id": 1, "gates": []}, + }, + "2Q": { + "0-1": { + "ids": [0, 1], + "gates": [ + { + "operator_type": "gate", + "operator": "CZ", + "parameters": [], + "arguments": ["_", "_"], + "fidelity": 0.9, + } + ], + } + }, + } + ) + + def test_builds_gate_and_edge_channels(self): + nm = NoiseModel.from_isa(self._isa()) + assert isinstance(nm.get_channel(Gate("RX", [1.5707963267948966], [0])), Channel) + assert isinstance(nm.get_channel(Gate("CZ", [], [0, 1])), Channel) + + def test_measurement_dedup_prefers_real_fidelity(self): + """A None-fidelity measure entry must not block a later usable one (dedup ordering).""" + nm = NoiseModel.from_isa(self._isa()) + channel = nm.get_channel(Measurement(qubit=Qubit(0), classical_reg=None)) + assert isinstance(channel, MeasurementChannel) + + +class TestCycleChannel: + def test_complete_cycle_constructs(self): + """A CycleChannel whose channels cover every DefCircuit body instruction is valid.""" + channels = tuple( + Channel.from_depolarizing_constant(inst, depolarizing_constant=0.99) for inst in (RX(0.1, 0), RZ(0.2, 1)) + ) + cycle = _build_cycle_channel(list(channels)) + assert cycle.channels == channels + + def test_incomplete_cycle_rejected(self): + """A body instruction with no corresponding channel is a footgun and must raise.""" + q0, q1 = FormalArgument("q0"), FormalArgument("q1") + # DefCircuit body has two gates, but only one channel is supplied for q0. + defcircuit = DefCircuit("CYCLE", [], [q0, q1], [RX(0.1, q0), RZ(0.2, q1)]) + cycle_inst = Gate("CYCLE", [], [Qubit(0), Qubit(1)]) + channels = (Channel.from_depolarizing_constant(RX(0.1, 0), depolarizing_constant=0.99),) + + with pytest.raises(ValueError, match="incomplete"): + _ = CycleChannel(inst=cycle_inst, defcircuit=defcircuit, channels=channels)