diff --git a/docs/source/quantum_expert_area/noisy_simulations.rst b/docs/source/quantum_expert_area/noisy_simulations.rst index bf646d37..f4b0a793 100644 --- a/docs/source/quantum_expert_area/noisy_simulations.rst +++ b/docs/source/quantum_expert_area/noisy_simulations.rst @@ -103,6 +103,28 @@ phase shifter, builds one unitary, computes probabilities, then averages the probability distributions. Amplitudes are not averaged. Use ``torch.manual_seed(...)`` to make the sampled perturbations reproducible. +Black-box unitaries under circuit phase noise +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +``pcvl.Unitary`` components contain no explicit phase shifters, but on hardware +they are implemented by a programmable interferometer. When ``phase_imprecision`` +or ``phase_error`` is active, Merlin decomposes every ``pcvl.Unitary`` component +into a rectangular (Clements) mesh of Mach-Zehnder interferometers. The mesh +phases are represented by ``pcvl.PS`` components and receive phase noise. The +fixed 50:50 ``pcvl.BS`` (beam splitter) parameters carry no phase noise. +``pcvl.PERM`` components are waveguide crossings without programmable phases +and are never decomposed. + +The decomposition runs once, when the layer is built (about 0.2 s for a +10-mode unitary), and only when circuit phase noise is configured; without +noise, each unitary stays precomputed as a constant tensor, exactly as before. +The numerical fit reproduces the target unitary with a fidelity error of about +``1e-6``, which corresponds to matrix entries agreeing to about ``1e-3``. The +fit's arbitrary global phase is compensated on the mesh's output phase layer, +so unitaries acting on a subset of modes keep their relative phase with the +rest of the circuit. If the fit does not converge, construction raises an +error instead of silently changing the mesh topology. + Coherent and incoherent error ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ diff --git a/merlin/pcvl_pytorch/locirc_to_tensor.py b/merlin/pcvl_pytorch/locirc_to_tensor.py index ab3d7a9f..4252db21 100644 --- a/merlin/pcvl_pytorch/locirc_to_tensor.py +++ b/merlin/pcvl_pytorch/locirc_to_tensor.py @@ -22,8 +22,10 @@ from __future__ import annotations +import copy import warnings +import numpy as np import torch from multipledispatch import dispatch from perceval.components import ( @@ -36,6 +38,8 @@ Circuit, Unitary, ) +from perceval.utils import Matrix +from perceval.utils.algorithms.circuit_optimizer import CircuitOptimizer from ..utils.dtypes import resolve_float_complex @@ -50,6 +54,45 @@ Barrier: Synchronization barrier (removed during compilation) """ +_DECOMPOSITION_CACHE: dict[bytes, Circuit] = {} +"""Cache mapping target-matrix bytes to its alpha-corrected Clements mesh. + +Keys are ``ndarray.tobytes()`` of the complex128 target matrix (so identical +unitaries — same bytes — always hit the same entry). Values are the completed +(fixed-parameter, alpha-corrected) ``Circuit`` objects from the first +successful ``CircuitOptimizer().optimize_rectangle`` call for that matrix. + +Consequences: + +* Each target matrix is decomposed at most once regardless of how many + ``CircuitConverter`` instances are built for circuits that contain it. +* The global NumPy RNG state is perturbed at most once per unique matrix, + so repeated or rebuild operations are deterministic from the caller's + perspective and do not accumulate RNG drift. +* ``phase_imprecision`` quantization, being a discontinuous function of the + mesh phases, now also produces consistent results across rebuilds because + the mesh phases are identical every time. + +Thread safety: the cache is populated lazily without a lock. Concurrent +writes for the same key both produce valid (equivalent) results, so the +only observable consequence is a redundant optimizer call during a race — +not incorrect output. +""" + +_DEFAULT_OPTIMIZE_RECTANGLE = CircuitOptimizer.optimize_rectangle +"""Optimizer implementation used to create entries in the decomposition cache.""" + +_FIT_TOLERANCE = 1e-3 +"""Tolerance for unitary decomposition fit quality. + +After fitting a unitary to an MZI mesh, the overlap magnitude + abs(trace(target^H @ fitted) / dimension) +should be close to 1.0. If the overlap falls below (1 - tolerance), +the fit is considered to have failed and a RuntimeError is raised. +With default CircuitOptimizer settings (fidelity error ~1e-6), +the overlap should be > 0.999 in all normal cases. +""" + def _phase_shifter_max_error(component: AComponent) -> float: """Return the local Perceval phase-error half-width for a phase shifter. @@ -86,6 +129,246 @@ def _circuit_has_local_phase_error(circuit: Circuit) -> bool: return any(_phase_shifter_max_error(component) > 0.0 for _, component in circuit) +def _decompose_unitaries( + circuit: Circuit, + phase_imprecision: float = 0.0, + phase_error: float = 0.0, +) -> Circuit: + """Replace black-box ``Unitary`` blocks with Clements MZI meshes. + + Called when circuit phase noise is configured (``phase_imprecision > 0`` + or ``phase_error > 0``). Each non-``PERM`` ``Unitary`` component is fitted + to a rectangular (Clements) mesh of MZIs with all parameters fixed as + constants. This makes the unitary's structure explicit. The mesh phases + are represented by ``PS`` components and receive phase noise; the fixed + ``BS`` parameters do not. + ``PERM`` components (waveguide crossings with no programmable phases) are + left untouched. + A 1-mode ``Unitary`` (a bare global phase) is replaced by a single ``PS`` + without running the optimizer. + + The fit reproduces each target unitary within the optimizer threshold + (fidelity error ~1e-6, matrix entries within ~1e-3). Its arbitrary global + phase is folded into the mesh's output phase layer so that blocks acting + on a subset of modes keep their relative phase with the rest of the + circuit. + + Decomposition results are memoized in ``_DECOMPOSITION_CACHE`` keyed on + the matrix bytes. After the first construction for a given matrix, every + subsequent ``CircuitConverter`` build reuses the cached mesh (deep-copied + for independence), so: (a) the global NumPy RNG is perturbed at most once + per unique matrix; (b) rebuilding or reloading a model always produces the + same mesh phases; and (c) ``phase_imprecision`` quantization, which is a + discontinuous function of the individual mesh phases, gives consistent + results across builds. + + .. note:: + Decomposition complexity is ``O(m^3)`` in the optimizer iterations and + typically takes ~0.2 s for 10 modes. For ``m >= 16`` the first build + may take several seconds; consider building the converter once and + reusing it. + + When phase noise is configured, a warning is emitted if the decomposition + fit error is comparable to or exceeds the configured noise scale. In this + case, the fit error becomes an uncontrolled artifact that can dominate + over the modeled physical noise. + + Parameters + ---------- + circuit : Circuit + Perceval circuit to transform. + phase_imprecision : float + Deterministic phase quantization step in radians. Used to estimate + the physical noise scale for warning purposes. + phase_error : float + Stochastic phase perturbation half-width in radians. Used to estimate + the physical noise scale for warning purposes. + + Returns + ------- + Circuit + The original circuit object when it contains no non-``PERM`` + ``Unitary`` component, otherwise a new circuit with each such + component replaced by an equivalent MZI mesh. + + Raises + ------ + ValueError + If a ``Unitary`` component uses polarization. + RuntimeError + If the Clements fit does not converge for a component. + """ + if not any( + isinstance(component, Unitary) and not isinstance(component, PERM) + for _, component in circuit + ): + return circuit + + new_circuit = Circuit(circuit.m) + for r, component in circuit: + if not isinstance(component, Unitary) or isinstance(component, PERM): + new_circuit.add(r, component) + continue + if component.requires_polarization: + raise ValueError( + f"Circuit phase noise cannot be applied to polarized unitary " + f"'{component.name}' on modes {tuple(r)}: decomposition into " + f"an MZI mesh is not supported for polarized components." + ) + + target = np.asarray(component.compute_unitary(), dtype=complex) + + # A 1-mode Unitary is a bare global phase exp(i*phi) with no beam + # splitters to place. The MZI optimizer has no defined behaviour for + # a 1x1 matrix, so handle it directly by emitting a single PS. + if component.m == 1: + ps_phase = float(np.angle(target[0, 0])) + new_circuit.add(r[0], PS(ps_phase)) + continue + + # Warn for large mode counts where the optimizer cost is superlinear + # and can take several seconds on the first build. + if component.m >= 16: + warnings.warn( + f"Decomposing a {component.m}-mode unitary '{component.name}' " + f"into a Clements mesh. The optimizer cost is O(m^3) in " + f"iteration count; for m >= 16 this may take several seconds. " + f"The result is cached after the first build for each unique matrix.", + UserWarning, + stacklevel=3, + ) + + # Check cache before running the optimizer. The cache key is the + # bytes representation of the complex128 target matrix, which + # uniquely identifies the unitary up to floating-point equality. + cache_key = target.tobytes() + # Do not use entries created by the default optimizer when the method + # has been replaced. Besides keeping cache contents tied to the + # optimizer that created them, this ensures validation is performed + # against the active implementation. + optimizer_is_default = ( + CircuitOptimizer.optimize_rectangle is _DEFAULT_OPTIMIZE_RECTANGLE + ) + if cache_key in _DECOMPOSITION_CACHE and optimizer_is_default: + # Reuse a deep copy of the cached mesh so each CircuitConverter + # gets an independent set of components. + cached_mesh = copy.deepcopy(_DECOMPOSITION_CACHE[cache_key]) + cached_fitted = np.asarray(cached_mesh.compute_unitary(), dtype=complex) + if not np.allclose( + cached_fitted.conj().T @ cached_fitted, + np.eye(component.m), + atol=1e-10, + ): + raise RuntimeError( + f"Cached Clements decomposition is not unitary for " + f"unitary '{component.name}' on modes {tuple(r)}." + ) + cached_overlap = np.trace(target.conj().T @ cached_fitted) / component.m + if abs(cached_overlap) < 1.0 - _FIT_TOLERANCE: + raise RuntimeError( + f"Cached Clements decomposition failed the fit quality check " + f"for unitary '{component.name}' on modes {tuple(r)}: " + f"overlap magnitude = {abs(cached_overlap):.6e}." + ) + cached_parameters = cached_mesh.get_parameters(all_params=True) + cached_output_phases = [ + parameter + for parameter in cached_parameters + if parameter.name.startswith("phL") + ] + if len(cached_output_phases) == component.m: + new_circuit.add(r[0], cached_mesh, merge=True) + continue + + # Entries created by an older optimizer/template are not valid + # for the current output-phase contract. Treat them as misses. + del _DECOMPOSITION_CACHE[cache_key] + + try: + mesh = CircuitOptimizer().optimize_rectangle(Matrix(target)) + except RuntimeError as exc: + raise RuntimeError( + f"Clements decomposition did not converge for unitary " + f"'{component.name}' on modes {tuple(r)}: {exc}" + ) from exc + + # The fit reproduces the target only up to a global phase e^{i*alpha}. + # On a block covering a subset of modes that phase is physical (it is + # relative to the untouched modes), so fold -alpha into the mesh's + # output phase layer, which has exactly one PS per mode. + fitted = np.asarray(mesh.compute_unitary(), dtype=complex) + overlap = np.trace(target.conj().T @ fitted) / component.m + + # Verify fit quality: if overlap magnitude is too small, the fit failed + # and alpha becomes meaningless (trace near zero makes phase arbitrary). + # So failures are raised explicitly and not hidden. + if abs(overlap) < 1.0 - _FIT_TOLERANCE: + raise RuntimeError( + f"Clements decomposition fit quality check failed for unitary " + f"'{component.name}' on modes {tuple(r)}: overlap magnitude " + f"|trace(target^H @ fitted) / m| = {abs(overlap):.6e} < " + f"{1.0 - _FIT_TOLERANCE:.6e}. The fitted mesh does not adequately " + f"reproduce the target unitary. This may indicate a bug in " + f"CircuitOptimizer or insufficient optimization iterations." + ) + + alpha = float(np.angle(overlap)) + + # Warn if decomposition fit error is comparable to configured phase noise. + # For unitary matrices, the fit quality is measured by the overlap magnitude: + # overlap = trace(target^H @ fitted) / m, which is bounded in [0, 1]. + # The fit error is 1 - |overlap|. With default CircuitOptimizer settings + # (fidelity error ~1e-6), we expect |overlap| > 0.999, i.e., fit_error < 0.001. + # + # For comparison, a phase error of delta_phi rad induces a unitary entry + # change ~delta_phi, so the noise scale is approximately + # max(phase_imprecision, phase_error). If fit_error > 0.1 * noise_scale, + # the fit error likely dominates the physical noise being modeled. + if phase_imprecision > 0.0 or phase_error > 0.0: + noise_scale = max(phase_imprecision, phase_error) + # Remove the arbitrary global phase before measuring the matrix + # residual. This detects entry-wise fit errors even when the + # trace overlap happens to cancel them for a particular target. + phase_aligned_fitted = fitted * np.exp(-1j * np.angle(overlap)) + fit_error = float( + np.linalg.norm(target - phase_aligned_fitted, ord="fro") + / np.sqrt(component.m) + ) + if fit_error > 0.1 * noise_scale: + warnings.warn( + f"Unitary decomposition phase-aligned fit residual " + f"({fit_error:.3e}) is " + f"comparable to or exceeds the configured phase noise scale " + f"(max(phase_imprecision, phase_error)={noise_scale:.3e}). " + f"For component '{component.name}' on modes {tuple(r)}, the " + f"decomposition artifact may dominate over the modeled physical noise. " + f"Consider using smaller phase noise values or an exact decomposition method.", + UserWarning, + stacklevel=3, + ) + + mesh_parameters = mesh.get_parameters(all_params=False) + output_phases = [p for p in mesh_parameters if p.name.startswith("phL")] + if len(output_phases) != component.m: + raise RuntimeError( + f"Unexpected mesh structure from CircuitOptimizer for unitary " + f"'{component.name}': expected {component.m} output phases " + f"('phL*'), found {len(output_phases)}." + ) + for parameter in mesh_parameters: + value = float(parameter) + if parameter.name.startswith("phL"): + value -= alpha + parameter.fix_value(value) + + # Store the completed mesh in the cache before adding to the circuit, + # so that future builds for the same matrix skip the optimizer entirely. + new_circuit.add(r[0], mesh, merge=True) + if optimizer_is_default: + _DECOMPOSITION_CACHE[cache_key] = copy.deepcopy(mesh) + return new_circuit + + class CircuitConverter: """Convert a parameterized Perceval circuit into a differentiable PyTorch unitary matrix. @@ -175,6 +458,17 @@ class CircuitConverter: The quantization uses nearest-grid rounding through :func:`torch.round`; it is not floor or truncation. + Black-Box Unitaries Under Circuit Noise: + When ``phase_imprecision > 0`` or ``phase_error > 0``, every + non-``PERM`` ``Unitary`` component is automatically decomposed at + construction into a rectangular (Clements) mesh of MZIs with fixed + numeric phases plus an output PS layer. Phase noise is then applied + to the PS components of the mesh (fidelity error ~1e-6, global phase + compensated), matching Perceval's convention that ``phase_imprecision`` + targets phase shifters only. Without circuit noise the fast path is + kept: ``Unitary`` components stay precomputed constant tensors. + ``PERM`` components (waveguide crossings) are never decomposed. + Example: Basic usage with a single phase shifter: @@ -311,6 +605,16 @@ def __init__( assert isinstance(circuit, Circuit), ( f"Expected a Perceval LO circuit, but got {type(circuit).__name__}" ) + # Circuit phase noise must reach the phases inside black-box Unitary + # blocks, so decompose them into MZI meshes before building the + # parameter mapping. Without noise the fast path (precomputed constant + # tensor per Unitary) is kept. + if self._phase_imprecision > 0.0 or self._phase_error > 0.0: + circuit = _decompose_unitaries( + circuit, + phase_imprecision=self._phase_imprecision, + phase_error=self._phase_error, + ) self.circuit = circuit if self._phase_error > 0.0 and _circuit_has_local_phase_error(self.circuit): warnings.warn( @@ -447,11 +751,15 @@ def _compile_circuit(self): ) if isinstance(c, Barrier): continue - # Check if this PS component requires dynamic handling due to phase_error. - # Phase shifters with active phase_error must remain dynamic (not precomputed) - # because fresh perturbations must be drawn on each call to to_tensor(). - # In contrast, phase_imprecision-only PS can still be precomputed since - # quantization is deterministic. + # A component must remain dynamic (not precomputed as a tensor) when + # it carries phases that receive stochastic perturbations, because + # fresh samples must be drawn on every call to to_tensor(). + # Deterministic phase_imprecision-only components can still be + # precomputed: _compute_tensor applies quantization at that point, + # and the result is the same on every call. + # + # PS: also sensitive when the component carries a per-component + # max_error (local phase error override). is_phase_error_sensitive = isinstance(c, PS) and ( self._phase_error > 0.0 or _phase_shifter_max_error(c) > 0.0 ) @@ -654,6 +962,55 @@ def to_tensor( return converted_tensor + def _apply_phase_noise( + self, + phase: torch.Tensor, + phase_error_half_width: float, + ) -> torch.Tensor: + """Apply configured phase noise to a phase tensor in-place (returns new tensor). + + Applies, in order: + + 1. Deterministic nearest-grid quantization (STE) when + ``self._phase_imprecision > 0``. + 2. Stochastic uniform perturbation when ``self._apply_phase_error`` is + ``True`` and ``phase_error_half_width > 0``. + + This helper is used by the :meth:`_compute_tensor` overload for ``PS``. + Phase noise is not applied to ``BS`` parameters, matching Perceval's + ``phase_imprecision`` noise model which targets only phase shifters. + + Parameters + ---------- + phase : torch.Tensor + Phase value(s) as a real-valued tensor. The tensor may be scalar, + 1-D ``(batch_size,)``, or any shape that is compatible with the + calling code. + phase_error_half_width : float + Half-width of the ``Uniform(-w, w)`` perturbation to apply when + ``self._apply_phase_error`` is ``True``. The caller selects the + per-component ``max_error`` or falls back to ``self._phase_error``. + + Returns + ------- + torch.Tensor + Phase tensor with quantization and/or perturbation applied. The + returned tensor shares the device and dtype of the input. + """ + if self._phase_imprecision > 0.0: + step = phase.new_tensor(self._phase_imprecision) + phase_quantized = torch.round(phase / step) * step + # Straight-through estimator: autograd sees d phase / d commanded = 1. + phase = phase + (phase_quantized - phase).detach() + + if self._apply_phase_error and phase_error_half_width > 0.0: + noise = torch.empty_like(phase).uniform_( + -phase_error_half_width, phase_error_half_width + ) + phase = phase + noise + + return phase + @dispatch((Unitary, PERM)) def _compute_tensor(self, comp: AComponent) -> torch.Tensor: """Compute tensor for Unitary and Permutation components. @@ -678,29 +1035,38 @@ def _compute_tensor(self, comp: AComponent) -> torch.Tensor: # type: ignore[no- """Compute tensor for Beam Splitter component. Handles different BS conventions (Rx, Ry, H) and processes 5 parameters: - theta, phi_tl, phi_bl, phi_tr, phi_br + theta, phi_tl, phi_bl, phi_tr, phi_br. - Args: - comp: BS component with parameters + Phase noise (``phase_imprecision`` and ``phase_error``) is not applied to + BS parameters. Noise is applied exclusively to PS (phase shifter) + components, matching Perceval's ``phase_imprecision`` noise model. - Returns: - Batched 2x2 unitary tensor of shape (batch_size, 2, 2) + Parameters + ---------- + comp : AComponent + ``BS`` component with parameters. + + Returns + ------- + torch.Tensor + Batched 2×2 unitary tensor of shape ``(batch_size, 2, 2)``. - Raises: - NotImplementedError: If BS convention is not supported + Raises + ------ + NotImplementedError + If the BS convention is not ``Rx``, ``Ry``, or ``H``. """ param_values = [] for _index, param in enumerate(comp.get_parameters(all_params=True)): if param.is_variable: tensor_id, idx_in_tensor = self.param_mapping[param.name] - param_values.append(self.torch_params[tensor_id][..., idx_in_tensor]) + raw = self.torch_params[tensor_id][..., idx_in_tensor] else: - param_values.append( - torch.tensor( - float(param), dtype=self.tensor_fdtype, device=self.device - ) + raw = torch.tensor( + float(param), dtype=self.tensor_fdtype, device=self.device ) + param_values.append(raw) cos_theta = torch.cos(param_values[0] / 2) sin_theta = torch.sin(param_values[0] / 2) @@ -822,35 +1188,15 @@ def _compute_tensor(self, comp: AComponent) -> torch.Tensor: # type: ignore[no- if phase.ndim == 0 and self.batch_size > 1: phase = phase.expand(self.batch_size) - # Apply finite phase resolution with nearest-grid quantization. This is - # not truncation: a phase is mapped to - # round(phase / phase_imprecision) * phase_imprecision. torch.round - # decides exact half-step ties, so pi/8 with a pi/4 step maps to 0. - # The STE keeps gradients attached to the commanded phase. - if self._phase_imprecision > 0.0: - phase_imprecision = phase.new_tensor(self._phase_imprecision) - phase_quantized = torch.round(phase / phase_imprecision) * phase_imprecision - # Straight-through estimator: adding a detached delta makes the - # forward value equal to phase_quantized, while autograd sees - # d phase / d commanded_phase = 1 because the delta is constant. - phase = phase + (phase_quantized - phase).detach() - - # Apply stochastic phase perturbation after quantization. Each call - # draws fresh Uniform(-phase_error, phase_error) samples. The samples - # are noise, not trainable parameters; gradients flow only through the - # commanded phase value. + # Compute phase error half-width: prefer component-specific max error, + # then fall back to the global phase_error component_phase_error = _phase_shifter_max_error(comp) phase_error_half_width = ( component_phase_error if component_phase_error > 0.0 else self._phase_error ) - if self._apply_phase_error and phase_error_half_width > 0.0: - # Use torch.empty_like() to ensure perturbations follow the same - # device, dtype, and batch structure as the phase tensor. - phase_error = torch.empty_like(phase).uniform_( - -phase_error_half_width, - phase_error_half_width, - ) - phase = phase + phase_error + + # Apply quantization and perturbation via the common helper + phase = self._apply_phase_noise(phase, phase_error_half_width) unitary_tensor = torch.exp(1j * phase.to(self.tensor_cdtype)).reshape( -1, 1 diff --git a/tests/algorithms/test_noise_contract.py b/tests/algorithms/test_noise_contract.py index 0fe29a2a..ed665696 100644 --- a/tests/algorithms/test_noise_contract.py +++ b/tests/algorithms/test_noise_contract.py @@ -1574,3 +1574,81 @@ def test_g2_output_keys_match_tensor_order_with_forward(): assert perceval_states_pl.issubset(set(keys_g2_pl)), ( f"Some Perceval states not in Merlin keys. Perceval: {perceval_states_pl}, Merlin: {set(keys_g2_pl)}" ) + + +# Black-box Unitary blocks under circuit phase noise (reservoir-style circuits) +def _reservoir_unitary_circuit(m: int = 4, n_features: int = 2) -> pcvl.Circuit: + """Create a reservoir circuit: random Unitary / px encoder / random Unitary.""" + circuit = pcvl.Circuit(m) + circuit.add(0, pcvl.Unitary(pcvl.Matrix.random_unitary(m))) + for i in range(n_features): + circuit.add(i, pcvl.PS(pcvl.P(f"px{i + 1}"))) + circuit.add(0, pcvl.Unitary(pcvl.Matrix.random_unitary(m))) + return circuit + + +def _reservoir_layer( + circuit: pcvl.Circuit, noise: pcvl.NoiseModel | None = None, **kwargs +) -> ml.QuantumLayer: + return ml.QuantumLayer( + input_size=2, + circuit=circuit, + input_parameters=["px"], + input_state=[1, 0, 1, 0], + n_photons=2, + noise=noise, + measurement_strategy=ml.MeasurementStrategy.probs( + computation_space=ml.ComputationSpace.FOCK + ), + dtype=torch.float64, + **kwargs, + ) + + +def test_reservoir_unitary_layer_with_phase_error_forwards_and_backwards(): + np.random.seed(42) + circuit = _reservoir_unitary_circuit() + noisy_layer = _reservoir_layer( + deepcopy(circuit), + noise=pcvl.NoiseModel(phase_error=0.1), + n_phase_error_samples=3, + ) + noiseless_layer = _reservoir_layer(deepcopy(circuit)) + + x = torch.tensor([[0.3, 0.7]], dtype=torch.float64, requires_grad=True) + noisy_output = noisy_layer(x) + + # 10 = C(4 + 2 - 1, 2) Fock states for 2 photons in 4 modes + _assert_normalized_distribution(noisy_output, 10) + noisy_output.pow(2).sum().backward() + assert x.grad is not None + assert torch.isfinite(x.grad).all() + assert not torch.allclose(noisy_output, noiseless_layer(x.detach())) + + +def test_reservoir_unitary_layer_with_imprecision_only_shifts_output(): + np.random.seed(43) + circuit = _reservoir_unitary_circuit() + quantized_layer = _reservoir_layer( + deepcopy(circuit), noise=pcvl.NoiseModel(phase_imprecision=0.8) + ) + noiseless_layer = _reservoir_layer(deepcopy(circuit)) + + x = torch.tensor([[0.3, 0.7]], dtype=torch.float64) + quantized_output = quantized_layer(x) + + _assert_normalized_distribution(quantized_output, 10) + assert not torch.allclose(quantized_output, noiseless_layer(x)) + + +def test_reservoir_unitary_layer_neutral_noise_matches_no_noise(): + # Fast-path regression guard: a neutral NoiseModel must not trigger the + # Clements decomposition of Unitary blocks. + np.random.seed(44) + circuit = _reservoir_unitary_circuit() + neutral_layer = _reservoir_layer(deepcopy(circuit), noise=pcvl.NoiseModel()) + plain_layer = _reservoir_layer(deepcopy(circuit)) + + x = torch.tensor([[0.3, 0.7]], dtype=torch.float64) + + assert torch.allclose(neutral_layer(x), plain_layer(x)) diff --git a/tests/models/test_reservoir_classifier.py b/tests/models/test_reservoir_classifier.py index b8ff4f91..4357c7ba 100644 --- a/tests/models/test_reservoir_classifier.py +++ b/tests/models/test_reservoir_classifier.py @@ -905,6 +905,277 @@ def test_layer_noise_model_rebuilds_layer_and_invalidates_fit(): assert model._fit_quantum_cache is None +def test_reservoir_classifier_with_phase_imprecision_and_phase_error(): + """Test ReservoirClassifier with circuit phase noise (phase_imprecision + phase_error). + + Verifies that: + - NoiseModel with phase_imprecision and phase_error can be set + - Model fits and predicts successfully with phase noise + - Phase noise invalidates the fitted state (triggering refitting) + - Multiple predictions differ due to phase_error stochastic sampling (cache disabled) + """ + X, y = _toy_data() + model = ReservoirClassifier( + in_features=4, + out_features=2, + n_photons=1, + reduction=PCA(n_components=2), + cache=False, # Disable cache to test noise variability + ) + + # Fit without noise + model.fit_reservoir(X) + dataset_before = model.make_dataset(X, y) + features_before, _ = dataset_before.tensors + + # Create NoiseModel with both phase_imprecision and phase_error + noise_model = pcvl.NoiseModel( + phase_imprecision=0.1, # ~5.7° quantization step + phase_error=0.05 # ±2.9° stochastic noise + ) + model.layer.noise = noise_model + + # Verify noise is set + assert model.layer.noise is noise_model + assert model._quantum_layer.noise is noise_model + + # Verify fit is invalidated + assert model._is_fitted is False + assert model._fit_quantum_cache is None + + # Refit with noise and verify it works + model.fit_reservoir(X) + dataset_with_noise = model.make_dataset(X, y) + features_with_noise, _ = dataset_with_noise.tensors + + assert features_with_noise.shape == features_before.shape + assert features_with_noise.dtype == features_before.dtype + + # Verify predictions work with noise + logits = model.predict(X) + assert logits.shape == (len(X), 2) + + # With cache=False, phase_error causes stochastic sampling, so predictions differ + logits1 = model.predict(X) + logits2 = model.predict(X) + assert not torch.allclose(logits1, logits2, rtol=0.1), ( + "With cache=False and phase_error, predictions should differ due to stochastic sampling" + ) + + # Note: Features will differ due to phase noise, but shouldn't be drastically different + # (they should remain within reasonable bounds for classification to still work) + feature_diff = torch.abs(features_with_noise - features_before).mean() + # Allow for variation but ensure noise didn't completely break the features + assert feature_diff < 1.0 # Reasonable tolerance for noise impact + + +def test_reservoir_classifier_with_phase_imprecision_only(): + """Test ReservoirClassifier with deterministic phase imprecision only. + + Verifies that: + - Phase imprecision alone (without phase_error) produces deterministic results + - Model fits and predicts successfully + """ + X, y = _toy_data() + model = ReservoirClassifier( + in_features=4, + out_features=2, + n_photons=1, + reduction=PCA(n_components=2), + ) + model.fit_reservoir(X) + + # Set phase imprecision only (no stochastic phase_error) + noise_model = pcvl.NoiseModel(phase_imprecision=0.1) + model.layer.noise = noise_model + model.fit_reservoir(X) + + # Multiple calls should produce identical results (deterministic) + features1 = model.transform_reservoir(X) + features2 = model.transform_reservoir(X) + + assert torch.allclose(features1, features2, rtol=1e-5) + + # Predictions should still work + logits1 = model.predict(X) + logits2 = model.predict(X) + assert torch.allclose(logits1, logits2, rtol=1e-5) + + +def test_reservoir_classifier_with_phase_error_only(): + """Test ReservoirClassifier with stochastic phase_error only. + + Verifies that: + - Phase error (without quantization) can be configured + - With cache=False, phase_error triggers stochastic sampling on every forward pass + - Each prediction produces different results due to stochastic sampling + - Model fits and predicts successfully + + Note: phase_error automatically triggers stochastic sampling via + ComputationProcess._compute_phase_error_probabilities() on every forward pass. + With caching enabled (default), identical input X returns the cached result. + With cache=False, each call re-runs the quantum circuit with new noise samples. + """ + X, y = _toy_data() + model = ReservoirClassifier( + in_features=4, + out_features=2, + n_photons=1, + reduction=PCA(n_components=2), + cache=False, # Disable cache to observe stochastic phase_error sampling + ) + model.fit_reservoir(X) + + # Set phase_error only (no deterministic quantization) + noise_model = pcvl.NoiseModel(phase_error=0.1) + model.layer.noise = noise_model + model.fit_reservoir(X) + + # With cache=False, multiple calls trigger stochastic sampling and produce DIFFERENT results + features1 = model.transform_reservoir(X) + features2 = model.transform_reservoir(X) + + # Should differ due to stochastic phase_error sampling + assert not torch.allclose(features1, features2, rtol=0.1), ( + "With cache=False and phase_error, features should differ due to stochastic sampling" + ) + + # Predictions should also differ + logits1 = model.predict(X) + logits2 = model.predict(X) + assert not torch.allclose(logits1, logits2, rtol=0.1), ( + "With cache=False and phase_error, predictions should differ due to stochastic sampling" + ) + + # Verify model can be moved and still works + assert model.layer.noise is noise_model + assert model._is_fitted is True + + +def test_reservoir_classifier_with_phase_imprecision_and_error_combined(): + """Test ReservoirClassifier with both phase_imprecision and phase_error. + + Verifies that: + - Combined noise (quantization + stochastic error) works correctly + - With cache=True (default), identical input returns cached result (deterministic) + - With cache=False, phase_error causes stochastic variation on each forward pass + - Model fits and predicts successfully + """ + X, y = _toy_data() + model = ReservoirClassifier( + in_features=4, + out_features=2, + n_photons=1, + reduction=PCA(n_components=2), + cache=False, # Disable cache to observe combined noise effects + ) + model.fit_reservoir(X) + + # Set both phase_imprecision AND phase_error + noise_model = pcvl.NoiseModel( + phase_imprecision=0.05, # Deterministic quantization step + phase_error=0.02, # Stochastic noise sampled on every forward pass + ) + model.layer.noise = noise_model + model.fit_reservoir(X) + + # With cache=False, phase_error causes stochastic sampling, so multiple calls differ + features1 = model.transform_reservoir(X) + features2 = model.transform_reservoir(X) + + # Should differ due to phase_error stochastic sampling + assert not torch.allclose(features1, features2, rtol=0.1), ( + "With cache=False and phase_error, features should differ due to stochastic sampling" + ) + + # Predictions should also differ + logits1 = model.predict(X) + logits2 = model.predict(X) + assert not torch.allclose(logits1, logits2, rtol=0.1), ( + "With cache=False and phase_error, predictions should differ due to stochastic sampling" + ) + + # Verify the noise affected output compared to no-noise case + model_clean = ReservoirClassifier( + in_features=4, + out_features=2, + n_photons=1, + reduction=PCA(n_components=2), + ) + model_clean.fit_reservoir(X) + features_clean = model_clean.transform_reservoir(X) + + # Noisy and clean should differ (due to quantization) + feature_diff = torch.abs(features1 - features_clean).mean() + # Should have some difference due to phase_imprecision + assert feature_diff > 0.0 + + +@pytest.mark.skipif(not torch.cuda.is_available(), reason="CUDA not available") +def test_reservoir_classifier_with_phase_noise_on_gpu(): + """Test ReservoirClassifier with phase noise on GPU. + + Verifies that: + - Phase noise configuration works on GPU device + - Model can be moved to CUDA and still apply noise correctly + - Forward pass (fit + predict) works end-to-end on GPU + - Multiple forward passes with noise produce different results (cache disabled) + - Tensor outputs remain on GPU + """ + X, y = _toy_data() + model = ReservoirClassifier( + in_features=4, + out_features=2, + n_photons=1, + reduction=PCA(n_components=2), + cache=False, # Disable cache to allow re-running quantum layer with noise + ) + + # Move model to GPU + model.to("cuda") + assert model.device.type == "cuda" + + # Fit on GPU without noise + model.fit_reservoir(X) + + # Apply phase noise on GPU + noise_model = pcvl.NoiseModel( + phase_imprecision=0.1, + phase_error=0.05, + ) + model.layer.noise = noise_model + + # Verify noise is set + assert model.layer.noise is noise_model + assert model._quantum_layer.noise is noise_model + + # Refit with noise on GPU + model.fit_reservoir(X) + + # Make dataset on GPU + dataset = model.make_dataset(X, y) + features, targets = dataset.tensors + + # Verify tensors are on CPU (dataset tensors are returned on CPU by default) + assert targets.device.type == "cpu" + + # Get predictions on GPU + logits = model.predict(X) + assert logits.shape == (len(X), 2) + assert logits.device.type == "cpu" # Predictions are moved back to CPU + + # Verify readout layer weights are on GPU + assert model.readout.weight.device.type == "cuda" + + # Test multiple forward passes with different noise samples work + logits1 = model.predict(X) + logits2 = model.predict(X) + + # With phase_error active, results should differ slightly + assert not torch.allclose(logits1, logits2, rtol=0.1) + assert logits1.shape == logits2.shape + + def test_reservoir_memory_cache_stays_bounded(monkeypatch): # Use a larger synthetic dataset than the unit tests above so the cache path # is exercised with a non-trivial number of samples. diff --git a/tests/pcvl_pytorch/test_unitary_decomposition_noise.py b/tests/pcvl_pytorch/test_unitary_decomposition_noise.py new file mode 100644 index 00000000..33b125ea --- /dev/null +++ b/tests/pcvl_pytorch/test_unitary_decomposition_noise.py @@ -0,0 +1,512 @@ +# MIT License +# +# Copyright (c) 2026 Quandela +# +# Permission is hereby granted, free of charge, to any person obtaining a copy +# of this software and associated documentation files (the "Software"), to deal +# in the Software without restriction, including without limitation the rights +# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +# copies of the Software, and to permit persons to whom the Software is +# furnished to do so, subject to the following conditions: +# +# The above copyright notice and this permission notice shall be included in all +# copies or substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +# SOFTWARE. + +"""Tests for automatic Clements decomposition of Unitary blocks under phase noise. + +The Clements fit reproduces each target unitary to ~1e-6 fidelity error, which +translates to matrix entries agreeing to ~1e-3, so unitaries are compared with +``_ELEMENTWISE_ATOL`` (the fit's global phase is compensated in production +code, making elementwise comparison valid) — never through raw phase values, +which differ between fits. +""" + +from __future__ import annotations + +import numpy as np +import perceval as pcvl +import pytest +import torch +from perceval.components import PS, Unitary + +import merlin.pcvl_pytorch.locirc_to_tensor as locirc_to_tensor +from merlin.pcvl_pytorch.locirc_to_tensor import ( + CircuitConverter, + _DECOMPOSITION_CACHE, + _decompose_unitaries, +) + +_ELEMENTWISE_ATOL = 1e-2 +_FIDELITY_ATOL = 1e-5 + + +@pytest.fixture(autouse=True) +def clear_decomposition_cache(): + """Keep process-global decomposition state isolated between tests.""" + locirc_to_tensor._DECOMPOSITION_CACHE.clear() + yield + locirc_to_tensor._DECOMPOSITION_CACHE.clear() + + +def _circuit_unitary(circuit: pcvl.Circuit) -> np.ndarray: + return np.asarray(circuit.compute_unitary(), dtype=complex) + + +def _fidelity(u: np.ndarray, v: np.ndarray) -> float: + return float(np.abs(np.trace(u.conj().T @ v)) / u.shape[-1]) + + +def _has_black_box_unitary(circuit: pcvl.Circuit) -> bool: + return any( + isinstance(c, Unitary) and not isinstance(c, pcvl.PERM) for _, c in circuit + ) + + +def _single_unitary_circuit(m: int = 4) -> pcvl.Circuit: + return pcvl.Circuit(m) // pcvl.Unitary(pcvl.Matrix.random_unitary(m)) + + +def _reservoir_circuit(m: int = 4, n_features: int = 2) -> pcvl.Circuit: + circuit = pcvl.Circuit(m) + circuit.add(0, pcvl.Unitary(pcvl.Matrix.random_unitary(m))) + for i in range(n_features): + circuit.add(i, pcvl.PS(pcvl.P(f"px{i + 1}"))) + circuit.add(0, pcvl.Unitary(pcvl.Matrix.random_unitary(m))) + return circuit + + +def test_decompose_replaces_unitary_with_equivalent_mesh(): + np.random.seed(0) + circuit = _single_unitary_circuit(4) + + decomposed = _decompose_unitaries(circuit) + + assert decomposed is not circuit + assert not _has_black_box_unitary(decomposed) + assert decomposed.get_parameters() == [] + original = _circuit_unitary(circuit) + rebuilt = _circuit_unitary(decomposed) + assert _fidelity(original, rebuilt) == pytest.approx(1.0, abs=_FIDELITY_ATOL) + assert np.allclose(rebuilt, original, atol=_ELEMENTWISE_ATOL) + + +def test_decompose_returns_same_object_without_black_box_unitary(): + circuit = pcvl.Circuit(2) + circuit.add((0, 1), pcvl.BS.H()) + circuit.add(0, pcvl.PS(pcvl.P("phi"))) + circuit.add((0, 1), pcvl.PERM([1, 0])) + + assert _decompose_unitaries(circuit) is circuit + + +def test_decompose_compensates_global_phase_of_sub_mode_blocks(): + # A block on a subset of modes must not pick up the fit's arbitrary + # global phase: relative to the untouched modes it would be physical. + # Elementwise comparison of the full-circuit unitary catches it. + np.random.seed(1) + circuit = pcvl.Circuit(4) + circuit.add(1, pcvl.Unitary(pcvl.Matrix.random_unitary(2))) + + decomposed = _decompose_unitaries(circuit) + + assert np.allclose( + _circuit_unitary(decomposed), _circuit_unitary(circuit), atol=_ELEMENTWISE_ATOL + ) + + +def test_decomposition_cache_reuses_fixed_mesh(monkeypatch): + """Reuse a cached mesh after all of its parameters are fixed.""" + np.random.seed(101) + circuit = _single_unitary_circuit(4) + _DECOMPOSITION_CACHE.clear() + + first = _decompose_unitaries(circuit) + first_mesh = next(component for _, component in first) + first_parameters = [ + float(parameter) + for parameter in first_mesh.get_parameters(all_params=True) + if parameter.name.startswith("phL") + ] + + optimizer_call_count = 0 + original_optimizer = locirc_to_tensor.CircuitOptimizer.optimize_rectangle + + def count_optimizer_calls(self, target): + nonlocal optimizer_call_count + optimizer_call_count += 1 + return original_optimizer(self, target) + + monkeypatch.setattr( + locirc_to_tensor.CircuitOptimizer, + "optimize_rectangle", + count_optimizer_calls, + ) + monkeypatch.setattr(locirc_to_tensor, "_DEFAULT_OPTIMIZE_RECTANGLE", count_optimizer_calls) + second = _decompose_unitaries(circuit) + second_mesh = next(component for _, component in second) + second_parameters = [ + float(parameter) + for parameter in second_mesh.get_parameters(all_params=True) + if parameter.name.startswith("phL") + ] + + assert optimizer_call_count == 0 + assert second_parameters == first_parameters + + +def test_no_noise_converter_keeps_unitary_fast_path(): + np.random.seed(2) + circuit = _single_unitary_circuit(4) + + converter = CircuitConverter(circuit, dtype=torch.float64) + + assert converter.circuit is circuit + assert len(converter.list_rct) == 1 + assert isinstance(converter.list_rct[0][1], torch.Tensor) + + +def test_phase_error_triggers_decomposition_and_noiseless_eval_matches(): + np.random.seed(3) + circuit = _single_unitary_circuit(4) + + converter = CircuitConverter(circuit, dtype=torch.float64, phase_error=0.1) + + assert converter.circuit is not circuit + assert any(isinstance(c, PS) for _, c in converter.list_rct) + unitary = converter.to_tensor().numpy() + original = _circuit_unitary(circuit) + assert _fidelity(original, unitary) == pytest.approx(1.0, abs=_FIDELITY_ATOL) + assert np.allclose(unitary, original, atol=_ELEMENTWISE_ATOL) + + +def test_decomposed_mesh_keeps_all_phase_shifters_noise_sensitive(): + """Test that every phase shifter in a Clements mesh receives phase noise.""" + np.random.seed(31) + mode_count = 4 + circuit = _single_unitary_circuit(mode_count) + + converter = CircuitConverter( + circuit, + dtype=torch.float64, + phase_error=0.1, + ) + + expected_phase_shifter_count = ( + 2 * (mode_count * (mode_count - 1) // 2) + mode_count + ) + noisy_phase_shifters = [ + component for _, component in converter.list_rct if isinstance(component, PS) + ] + + assert len(noisy_phase_shifters) == expected_phase_shifter_count + + +def test_phase_error_sampling_varies_and_is_reproducible(): + np.random.seed(4) + circuit = _single_unitary_circuit(4) + converter = CircuitConverter(circuit, dtype=torch.float64, phase_error=0.1) + + torch.manual_seed(123) + first = converter.to_tensor(apply_phase_error=True) + second = converter.to_tensor(apply_phase_error=True) + torch.manual_seed(123) + replayed = converter.to_tensor(apply_phase_error=True) + + assert not torch.allclose(first, second) + assert torch.allclose(first, replayed) + + +def test_imprecision_only_triggers_decomposition_and_quantizes(): + np.random.seed(5) + circuit = _single_unitary_circuit(4) + + converter = CircuitConverter(circuit, dtype=torch.float64, phase_imprecision=0.5) + + assert converter.circuit is not circuit + # Mesh phases are quantized with a coarse step, so the effective unitary + # must deviate from the target well beyond the decomposition error. + unitary = converter.to_tensor().numpy() + assert not np.allclose(unitary, _circuit_unitary(circuit), atol=_ELEMENTWISE_ATOL) + + +def test_reservoir_circuit_keeps_input_parameters_and_gradients(): + np.random.seed(6) + circuit = _reservoir_circuit(4, 2) + + converter = CircuitConverter( + circuit, input_specs=["px"], dtype=torch.float64, phase_error=0.05 + ) + + assert converter.spec_mappings["px"] == ["px1", "px2"] + params = torch.tensor( + [[0.1, 0.2], [0.3, 0.4]], dtype=torch.float64, requires_grad=True + ) + unitary = converter.to_tensor(params) + assert unitary.shape == (2, 4, 4) + identity = torch.eye(4, dtype=unitary.dtype).expand(2, 4, 4) + assert torch.allclose(unitary @ unitary.mH, identity, atol=1e-6) + unitary.real.sum().backward() + assert params.grad is not None + assert torch.isfinite(params.grad).all() + + +def test_perm_is_not_decomposed(): + np.random.seed(7) + circuit = pcvl.Circuit(4) + circuit.add(0, pcvl.PERM([1, 0])) + circuit.add(2, pcvl.Unitary(pcvl.Matrix.random_unitary(2))) + perm = next(c for _, c in circuit if isinstance(c, pcvl.PERM)) + + decomposed = _decompose_unitaries(circuit) + + assert any(c is perm for _, c in decomposed) + assert not _has_black_box_unitary(decomposed) + + +def test_polarized_unitary_raises_value_error(): + circuit = pcvl.Circuit(1) + circuit.add(0, pcvl.Unitary(pcvl.Matrix.random_unitary(2), use_polarization=True)) + + with pytest.raises(ValueError, match="polarized unitary"): + _decompose_unitaries(circuit) + + +def test_non_convergence_error_names_component(monkeypatch): + def _fail(self, target): + raise RuntimeError("optimization above threshold") + + monkeypatch.setattr(locirc_to_tensor.CircuitOptimizer, "optimize_rectangle", _fail) + circuit = _single_unitary_circuit(4) + + with pytest.raises(RuntimeError, match=r"did not converge .* modes \(0, 1, 2, 3\)"): + _decompose_unitaries(circuit) + + +def test_one_mode_unitary_creates_single_ps(): + """Test that a 1-mode unitary is replaced by a single PS without optimizer.""" + np.random.seed(2) + # Create a 1-mode circuit with a bare global phase + circuit = pcvl.Circuit(1) + phase_value = np.pi / 4 + u = pcvl.Matrix.random_unitary(1) + circuit.add(0, pcvl.Unitary(u)) + + # Decompose with phase noise enabled + decomposed = _decompose_unitaries(circuit, phase_imprecision=0.1, phase_error=0.05) + + # Should replace Unitary with PS, not an MZI mesh + assert decomposed is not circuit + components = [c for _, c in decomposed] + assert len(components) == 1 + assert isinstance(components[0], PS) + + +def test_slow_decomposition_warning_for_large_unitary(): + """Test that a warning is emitted for m >= 16 mode unitaries.""" + np.random.seed(3) + circuit = pcvl.Circuit(16) + circuit.add((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15), + pcvl.Unitary(pcvl.Matrix.random_unitary(16))) + + # Should emit warning about slow decomposition + with pytest.warns(UserWarning, match=r"16-mode.*Clements mesh.*O\(m\^3\)"): + _decompose_unitaries(circuit, phase_imprecision=0.1) + + +def test_fit_error_vs_noise_scale_warning(monkeypatch): + """Test that a warning is emitted when fit error exceeds 0.1 * noise_scale.""" + np.random.seed(4) + circuit = _single_unitary_circuit(4) + + # Mock CircuitOptimizer.optimize_rectangle to return a mesh with controlled fit error + class MockMesh(pcvl.Circuit): + def __init__(self, target): + super().__init__(4) + self.target = target + for index in range(4): + self.add(index, PS(pcvl.P(f"phL{index}"))) + + def get_parameters(self, all_params=False): + return [ + type( + "MockParameter", + (), + { + "name": f"phL{index}", + "__float__": lambda self: 0.1, + "fix_value": lambda self, value: None, + }, + )() + for index in range(4) + ] + + def compute_unitary(self): + # Apply a small unitary perturbation that is larger than the + # configured noise threshold after global-phase alignment. + phase_rotation = np.diag(np.exp(1j * np.array([0.002, 0, 0, 0]))) + return self.target @ phase_rotation + + def mock_optimize(self, target): + return MockMesh(target) + + monkeypatch.setattr(locirc_to_tensor.CircuitOptimizer, "optimize_rectangle", mock_optimize) + + # Use phase_error small enough that 0.1 * noise_scale < fit_error (≈0.001) + # With phase_error = 0.005, noise_scale = 0.005, threshold = 0.0005 < 0.001 ✓ + with pytest.warns(UserWarning, match=r"fit residual.*comparable.*phase noise scale"): + _decompose_unitaries(circuit, phase_error=0.005) + + +def test_fit_quality_error_if_overlap_too_small(monkeypatch): + """Test that RuntimeError is raised if fit overlap magnitude is too small.""" + np.random.seed(5) + circuit = _single_unitary_circuit(4) + + # Mock CircuitOptimizer.optimize_rectangle to return a mesh with terrible fit + class BadMesh: + def __init__(self, target): + self.target = target + + def get_parameters(self): + return [] + + def compute_unitary(self): + # Return an orthogonal matrix unrelated to target + # This gives overlap ≈ 0, triggering the fit quality error + return np.array([[0, 1, 0, 0], + [-1, 0, 0, 0], + [0, 0, 0, 1], + [0, 0, -1, 0]], dtype=complex) + + def mock_optimize(self, target): + return BadMesh(target) + + monkeypatch.setattr(locirc_to_tensor.CircuitOptimizer, "optimize_rectangle", mock_optimize) + + with pytest.raises(RuntimeError, match=r"fit quality check failed.*overlap magnitude"): + _decompose_unitaries(circuit) + + +def test_phL_count_mismatch_error(monkeypatch): + """Test that RuntimeError is raised if mesh has wrong number of output phases.""" + np.random.seed(6) + circuit = _single_unitary_circuit(4) + + # Mock CircuitOptimizer.optimize_rectangle to return a mesh with wrong structure + class BadStructureMesh: + def __init__(self, target): + self.target = target + + def get_parameters(self, all_params=False): + # Return parameters with no phL entries (missing output layer) + # This simulates a malformed mesh structure + return [type('param', (), {'name': f'other{i}', '__float__': lambda self: 0.1})() + for i in range(5)] + + def compute_unitary(self): + # Return the target itself so overlap check passes (overlap = 1) + # Then the phL count check fails + return self.target + + def mock_optimize(self, target): + return BadStructureMesh(target) + + monkeypatch.setattr(locirc_to_tensor.CircuitOptimizer, "optimize_rectangle", mock_optimize) + + with pytest.raises(RuntimeError, match=r"Unexpected mesh structure.*expected 4 output phases.*found 0"): + _decompose_unitaries(circuit) + + +@pytest.mark.skipif(not torch.cuda.is_available(), reason="CUDA not available") +def test_circuit_converter_with_phase_noise_on_gpu(): + """Test CircuitConverter unitary decomposition and phase noise on GPU. + + Verifies that: + - Unitary decomposition works on GPU device + - Phase imprecision and phase_error produce correct outputs on GPU + - Tensors are created on the specified GPU device + - Stochastic sampling works correctly with fresh samples per call + - Results match CPU baseline (fidelity check) + """ + np.random.seed(42) + circuit = _single_unitary_circuit(4) + + # Create converter with phase noise on GPU + converter_gpu = CircuitConverter( + circuit, + dtype=torch.float64, + device="cuda", + phase_imprecision=0.1, + phase_error=0.05, + ) + + assert str(converter_gpu.device) == "cuda:0" or str(converter_gpu.device) == "cuda" + + # Get unitary on GPU + unitary_gpu = converter_gpu.to_tensor() + assert unitary_gpu.device.type == "cuda" + assert unitary_gpu.shape == (4, 4) + + # Verify decomposition occurred + assert converter_gpu.circuit is not circuit + assert any(isinstance(c, PS) for _, c in converter_gpu.list_rct) + + # Decomposed mesh PS phases receive noise; BS parameters do not. + # The converter's output should still be a valid unitary (U†U ≈ I). + # Verify this rather than comparing against the unquantized original. + unitary_gpu_numpy = unitary_gpu.cpu().numpy() + identity = np.eye(unitary_gpu_numpy.shape[0], dtype=complex) + assert np.allclose( + unitary_gpu_numpy.conj().T @ unitary_gpu_numpy, identity, atol=1e-6 + ) + + # Phase error sampling should vary on GPU + torch.manual_seed(123) + first_sample_gpu = converter_gpu.to_tensor(apply_phase_error=True) + second_sample_gpu = converter_gpu.to_tensor(apply_phase_error=True) + + assert first_sample_gpu.device.type == "cuda" + assert second_sample_gpu.device.type == "cuda" + assert not torch.allclose(first_sample_gpu, second_sample_gpu) + + # Verify reproducibility with same seed + torch.manual_seed(456) + third_sample_gpu = converter_gpu.to_tensor(apply_phase_error=True) + torch.manual_seed(456) + replayed_sample_gpu = converter_gpu.to_tensor(apply_phase_error=True) + + assert torch.allclose(third_sample_gpu, replayed_sample_gpu) + + # Compare CPU and GPU decomposition without quantization. Quantization is + # discontinuous at grid boundaries, and CPU/CUDA can round a boundary + # value differently even when the source mesh parameters are identical. + torch.manual_seed(789) + converter_cpu = CircuitConverter( + circuit, + dtype=torch.float64, + device="cpu", + phase_imprecision=0.0, + phase_error=0.05, + ) + unitary_cpu = converter_cpu.to_tensor() + + # The GPU result above includes quantization; compare fresh unquantized + # GPU output with the CPU baseline. + converter_gpu_no_quantization = CircuitConverter( + circuit, + dtype=torch.float64, + device="cuda", + phase_error=0.05, + ) + unitary_gpu_no_quantization = converter_gpu_no_quantization.to_tensor() + assert torch.allclose( + unitary_gpu_no_quantization.cpu(), unitary_cpu, atol=1e-3 + )