diff --git a/merlin/builder/circuit_builder.py b/merlin/builder/circuit_builder.py index 94c6cd115..cc5a3aa37 100644 --- a/merlin/builder/circuit_builder.py +++ b/merlin/builder/circuit_builder.py @@ -575,6 +575,7 @@ def add_entangling_layer( name: str | None = None, trainable_inner: bool | None = None, trainable_outer: bool | None = None, + seed: int | None = None, ) -> "CircuitBuilder": """Add an entangling layer spanning a range of modes. @@ -595,6 +596,9 @@ def add_entangling_layer( Override for the internal phase shifters. trainable_outer : bool | None Override for the output phase shifters. + seed : int | None + Optional seed for the random phases assigned to non-trainable + phase shifters. ``None`` draws fresh entropy each call. Returns ------- @@ -648,6 +652,7 @@ def add_entangling_layer( model=normalized_model, trainable_inner=trainable_inner, trainable_outer=trainable_outer, + seed=seed, ) self.circuit.add(component) @@ -939,16 +944,19 @@ def _mzi_factory( component, "trainable_outer", component.trainable ), base: str = prefix, + block: GenericInterferometer = component, ): """Build a Mach-Zehnder interferometer optionally parameterised per index.""" if inner_trainable: phi_inner = pcvl_module.P(f"{base}_li{i}") else: - phi_inner = 0.0 + fixed = block.fixed_inner_values + phi_inner = fixed[i] if i < len(fixed) else 0.0 if outer_trainable: phi_outer = pcvl_module.P(f"{base}_lo{i}") else: - phi_outer = 0.0 + fixed = block.fixed_outer_values + phi_outer = fixed[i] if i < len(fixed) else 0.0 return ( pcvl_module.BS() // pcvl_module.PS(phi_inner) @@ -975,16 +983,19 @@ def _bell_factory( component, "trainable_outer", component.trainable ), base: str = prefix, + block: GenericInterferometer = component, ): """Build a Mach-Zehnder interferometer optionally parameterised per index.""" if inner_trainable: phi_inner = pcvl_module.P(f"{base}_li{i}") else: - phi_inner = 0.0 + fixed = block.fixed_inner_values + phi_inner = fixed[i] if i < len(fixed) else 0.0 if outer_trainable: phi_outer = pcvl_module.P(f"{base}_lo{i}") else: - phi_outer = 0.0 + fixed = block.fixed_outer_values + phi_outer = fixed[i] if i < len(fixed) else 0.0 circuit = pcvl_module.Circuit(2) circuit.add(0, pcvl_module.BS()) diff --git a/merlin/core/components.py b/merlin/core/components.py index 855ac10d2..909073504 100644 --- a/merlin/core/components.py +++ b/merlin/core/components.py @@ -25,7 +25,9 @@ Components are platform-agnostic and focus on intent rather than implementation. """ -from dataclasses import dataclass +import math +import random +from dataclasses import dataclass, field from enum import Enum from typing import Any @@ -189,6 +191,27 @@ class GenericInterferometer: Whether inner phase shifters are trainable. trainable_outer : bool | None Whether outer phase shifters are trainable. + seed : int | None + Optional seed controlling the random phases drawn for non-trainable + inner/outer phase shifters. ``None`` draws fresh entropy without + touching global random state. + fixed_inner_values : list[float] + Fixed (non-trainable) values for the inner phase shifters, drawn + randomly unless explicitly provided. When provided, must contain + exactly ``span * (span - 1) // 2`` entries. + fixed_outer_values : list[float] + Fixed (non-trainable) values for the outer phase shifters, drawn + randomly unless explicitly provided. When provided, must contain + exactly ``span * (span - 1) // 2`` entries. + + Raises + ------ + TypeError + If ``model`` is not a string. + ValueError + If ``model`` is not ``"mzi"`` or ``"bell"``, or if a provided + ``fixed_inner_values``/``fixed_outer_values`` list does not contain + exactly one entry per phase shifter. """ start_mode: int @@ -198,6 +221,9 @@ class GenericInterferometer: model: str = "mzi" trainable_inner: bool | None = None trainable_outer: bool | None = None + seed: int | None = None + fixed_inner_values: list[float] = field(default_factory=list) + fixed_outer_values: list[float] = field(default_factory=list) def __post_init__(self): """Validate and normalize the interferometer configuration.""" @@ -219,6 +245,35 @@ def __post_init__(self): # Normalise the aggregate flag so downstream logic can rely on it self.trainable = bool(self.trainable_inner or self.trainable_outer) + # Non-trainable phase shifters must still be initialised to random + # phases (not a constant 0.0), otherwise a "fixed" interferometer + # degenerates into a deterministic swap/identity instead of the + # random fixed unitary that untrained entangling layers are meant to + # provide (e.g. for reservoir-computing style architectures). + count = self.span * (self.span - 1) // 2 if self.span > 1 else 0 + # A partially filled list would silently fall back to 0.0 phases for + # the missing indices downstream, so explicit values must cover every + # phase shifter. + for label, values in ( + ("fixed_inner_values", self.fixed_inner_values), + ("fixed_outer_values", self.fixed_outer_values), + ): + if values and len(values) != count: + raise ValueError( + f"GenericInterferometer {label} must contain exactly " + f"{count} entries for span {self.span}, got {len(values)}" + ) + if count > 0: + rng = random.Random(self.seed) + if not self.trainable_inner and not self.fixed_inner_values: + self.fixed_inner_values = [ + rng.uniform(0, 2 * math.pi) for _ in range(count) + ] + if not self.trainable_outer and not self.fixed_outer_values: + self.fixed_outer_values = [ + rng.uniform(0, 2 * math.pi) for _ in range(count) + ] + def get_params(self) -> dict[str, Any]: """Return placeholder names for every internal interferometer parameter. diff --git a/merlin/pcvl_pytorch/noisy_slos.py b/merlin/pcvl_pytorch/noisy_slos.py index bc5e7f183..eaab5a0ee 100644 --- a/merlin/pcvl_pytorch/noisy_slos.py +++ b/merlin/pcvl_pytorch/noisy_slos.py @@ -109,6 +109,11 @@ def __init__( for n_i in range(1, (2 * self.n_photons) + 1) ], ) + # Kept for the extra-photon sectors: photon-number-indexed SLOS + # graphs spanning up to 2*n_photons, reused by + # _augmented_obb_probs to simulate both the (grown) coherent cell + # and the one-hot distinguishable cells of each OBB partition. + self._regular_slos_graphs = regular_slos_graphs self._slos_graphs = [ NoisySLOSComputeGraph( noise_groups=noise_groups, @@ -193,6 +198,98 @@ def _get_extra_photon_combinations( return output + def _augmented_obb_probs( + self, + unitary: torch.Tensor, + obb: _InputStateNoisySLOSComputeGraph, + extra_vec: torch.Tensor, + n_augmented: int, + ) -> torch.Tensor: + """Compute one g2 extra-photon sector for a fixed emission combination. + + Perceval's source model draws one independent coherent-vs-distinguishable + outcome per *original* input photon only -- the same Orthogonal Bad Bits + partitions ``obb`` already built for the noiseless ``n_photons`` state -- + never an extra independent draw per g2-emitted sibling. Whatever the fate + of the original photons in a given partition cell, every g2 sibling + always joins that cell's coherent (non-excluded) group. This reuses + ``obb``'s cached partitions/weights, growing only the coherent cell's + Fock state by ``extra_vec`` before convolving with the (unchanged) + one-hot distinguishable cells. + + Parameters + ---------- + unitary : torch.Tensor + Circuit unitary, batched ``[batch, m, m]``. + obb : _InputStateNoisySLOSComputeGraph + Cached OBB graph for the base (non-augmented) input state. + extra_vec : torch.Tensor + Per-mode count of g2-emitted photons for this combination. + n_augmented : int + Total photon number of this sector (``n_photons`` + emitted). + + Returns + ------- + torch.Tensor + Probabilities over the ``n_augmented``-photon Fock basis, shape + ``[batch, n_states]``. + """ + batch_size = unitary.size(0) + fock_keys = [ + tuple(row) for row in self._fock_states_per_n[n_augmented].tolist() + ] + key_to_idx = {key: idx for idx, key in enumerate(fock_keys)} + output_probs = torch.zeros( + batch_size, len(fock_keys), device=self.device, dtype=self.dtype + ) + + for order, (cells, counts) in enumerate(obb._partitions): + bit_weight = obb._weights[order] + for cell, count in zip(cells, counts, strict=True): + if order == obb.n_photons: + base_state = extra_vec + bad_states = cell + else: + base_state = cell[0] + extra_vec + bad_states = cell[1:] + + base_n = int(base_state.sum().item()) + _, base_probs = self._regular_slos_graphs[base_n - 1].compute_probs( + unitary, base_state + ) + if base_probs.ndim == 1: + base_probs = base_probs.unsqueeze(0) + cell_keys = [self._fock_states_per_n[base_n]] + cell_probs = [base_probs] + for bad_state in bad_states: + _, bad_probs = self._regular_slos_graphs[0].compute_probs( + unitary, bad_state + ) + if bad_probs.ndim == 1: + bad_probs = bad_probs.unsqueeze(0) + cell_keys.append(self._fock_states_per_n[1]) + cell_probs.append(bad_probs) + + conv_keys, conv_probs = convolve_distributions(cell_keys, *cell_probs) + if conv_probs.ndim == 1: + conv_probs = conv_probs.unsqueeze(0) + conv_keys_list = ( + [tuple(k.tolist()) for k in conv_keys] + if isinstance(conv_keys, torch.Tensor) + else [tuple(k) for k in conv_keys] + ) + + for local_idx, key in enumerate(conv_keys_list): + idx = key_to_idx.get(key) + if idx is not None: + output_probs[:, idx] = ( + output_probs[:, idx] + + bit_weight * count.item() * conv_probs[:, local_idx] + ) + + output_probs = output_probs / output_probs.sum(dim=1, keepdim=True) + return output_probs + def compute_probs( self, unitary: torch.Tensor, @@ -314,14 +411,19 @@ def compute_probs( probs = reordered_probs else: - input_state_to_run = list(input_state) + extra_vec = torch.zeros( + self.m, dtype=torch.int32, device=self.device + ) for photon in combination: - input_state_to_run[photon] += 1 - probs = cast( - torch.Tensor, - slos_graphs_list[num_photons_added].compute_probs( - unitary, input_state_to_run - ), + extra_vec[photon] += 1 + obb = slos_graphs_list[0]._slos_graph_per_input[ + tuple(input_state) + ] + probs = self._augmented_obb_probs( + unitary, + obb, + extra_vec, + self.n_photons + num_photons_added, ) sector.tensor = sector.tensor + weight_k * probs @@ -367,6 +469,9 @@ def to(self, device: str | torch.device) -> NoisyG2SLOSComputeGraph: else: for graph in self._slos_graphs: graph.to(self.device) + if hasattr(self, "_regular_slos_graphs"): + for regular_graph in self._regular_slos_graphs: + regular_graph.to(self.device) return self diff --git a/tests/builder/test_circuit_builder.py b/tests/builder/test_circuit_builder.py index 32ebf96c6..42fc83a0b 100644 --- a/tests/builder/test_circuit_builder.py +++ b/tests/builder/test_circuit_builder.py @@ -1,8 +1,10 @@ from __future__ import annotations +import math import os from pathlib import Path +import numpy as np import perceval as pcvl import pytest import torch @@ -378,6 +380,96 @@ def test_entangling_layer_mode_range_and_non_trainable(): assert "block" in builder.trainable_parameter_prefixes +def test_non_trainable_entangling_layer_uses_random_phases_not_zero(): + """A non-trainable entangling layer must still carry random fixed phases. + + Regression test: previously every phase shifter of a ``trainable=False`` + entangling layer was silently pinned to 0.0, collapsing the block into a + deterministic swap/identity instead of a random fixed unitary (breaking, + e.g., reservoir-computing style architectures that rely on + ``add_entangling_layer(trainable=False)`` to inject a fixed random mix). + """ + builder = CircuitBuilder(n_modes=4) + builder.add_entangling_layer(trainable=False, name="pre_mix") + + component = builder.circuit.components[-1] + assert isinstance(component, GenericInterferometer) + assert component.trainable is False + + # Phases must be populated and not trivially all-zero. + assert component.fixed_inner_values + assert component.fixed_outer_values + # random.uniform() draws from a continuous distribution; all values being + # exactly 0.0 has measure-zero probability and is effectively impossible. + # These assertions guard against the old bug where phases were hardcoded to 0.0. + assert any(v != 0.0 for v in component.fixed_inner_values) + assert any(v != 0.0 for v in component.fixed_outer_values) + for v in component.fixed_inner_values + component.fixed_outer_values: + assert 0.0 <= v < 2 * math.pi + + pcvl_circuit = builder.to_pcvl_circuit(pcvl) + + # No trainable parameters should have been registered. + assert pcvl_circuit.get_parameters() == [] + + unitary = np.array(pcvl_circuit.compute_unitary()) + identity = np.eye(4) + swap = np.array( + [[0, 0, 0, 1], [0, 0, 1, 0], [0, 1, 0, 0], [1, 0, 0, 0]], dtype=complex + ) + assert not np.allclose(unitary, identity, atol=1e-6) + assert not np.allclose(np.abs(unitary), np.abs(swap), atol=1e-6) + + +def test_non_trainable_entangling_layer_is_reproducible_and_varies_by_instance(): + """Repeated conversion is stable; separate builders draw fresh phases.""" + builder = CircuitBuilder(n_modes=4) + builder.add_entangling_layer(trainable=False, name="pre_mix") + + unitary_first = np.array(builder.to_pcvl_circuit(pcvl).compute_unitary()) + unitary_second = np.array(builder.to_pcvl_circuit(pcvl).compute_unitary()) + assert np.allclose(unitary_first, unitary_second) + + other_builder = CircuitBuilder(n_modes=4) + other_builder.add_entangling_layer(trainable=False, name="pre_mix") + unitary_other = np.array(other_builder.to_pcvl_circuit(pcvl).compute_unitary()) + assert not np.allclose(unitary_first, unitary_other) + + +def test_non_trainable_entangling_layer_seed_is_reproducible(): + builder_a = CircuitBuilder(n_modes=4) + builder_a.add_entangling_layer(trainable=False, name="pre_mix", seed=1234) + + builder_b = CircuitBuilder(n_modes=4) + builder_b.add_entangling_layer(trainable=False, name="pre_mix", seed=1234) + + unitary_a = np.array(builder_a.to_pcvl_circuit(pcvl).compute_unitary()) + unitary_b = np.array(builder_b.to_pcvl_circuit(pcvl).compute_unitary()) + assert np.allclose(unitary_a, unitary_b) + + +def test_non_trainable_entangling_layer_partial_trainable_mix(): + """Only the trainable side should expose parameters; the other side stays random.""" + builder = CircuitBuilder(n_modes=4) + builder.add_entangling_layer( + name="mix", trainable_inner=True, trainable_outer=False + ) + + component = builder.circuit.components[-1] + assert component.trainable_inner is True + assert component.trainable_outer is False + assert component.fixed_inner_values == [] + assert component.fixed_outer_values + # random.uniform() draws from a continuous distribution; all values being + # exactly 0.0 has measure-zero probability and is effectively impossible. + assert any(v != 0.0 for v in component.fixed_outer_values) + + pcvl_circuit = builder.to_pcvl_circuit(pcvl) + params = pcvl_circuit.get_parameters() + assert any(p.name.startswith("mix_li") for p in params) + assert not any(p.name.startswith("mix_lo") for p in params) + + def test_entangling_layer_invalid_modes(): builder = CircuitBuilder(n_modes=4) diff --git a/tests/core/test_components.py b/tests/core/test_components.py index 96f07c396..421a3d2fc 100644 --- a/tests/core/test_components.py +++ b/tests/core/test_components.py @@ -1,10 +1,15 @@ from __future__ import annotations +import math + +import pytest + import merlin.core.components as components_mod Rotation = components_mod.Rotation BeamSplitter = components_mod.BeamSplitter EntanglingBlock = components_mod.EntanglingBlock +GenericInterferometer = components_mod.GenericInterferometer ParameterRole = components_mod.ParameterRole @@ -34,3 +39,83 @@ def test_beam_splitter_get_params_exposes_non_fixed_names(): def test_entangling_block_exposes_no_parameters(): block = EntanglingBlock(targets=[0, 1], depth=2, trainable=True) assert block.get_params() == {} + + +def test_generic_interferometer_non_trainable_gets_random_phases(): + """Non-trainable phase shifters must draw random phases, not default to 0.0.""" + gi = GenericInterferometer(start_mode=0, span=4, trainable=False) + + assert gi.trainable is False + count = gi.span * (gi.span - 1) // 2 + assert len(gi.fixed_inner_values) == count + assert len(gi.fixed_outer_values) == count + # random.uniform() draws from a continuous distribution; all values being + # exactly 0.0 has measure-zero probability and is effectively impossible. + # These assertions guard against the old bug where phases were hardcoded to 0.0. + assert any(v != 0.0 for v in gi.fixed_inner_values) + assert any(v != 0.0 for v in gi.fixed_outer_values) + for value in gi.fixed_inner_values + gi.fixed_outer_values: + assert 0.0 <= value < 2 * math.pi + + +def test_generic_interferometer_trainable_gets_no_fixed_phases(): + gi = GenericInterferometer(start_mode=0, span=4, trainable=True) + assert gi.fixed_inner_values == [] + assert gi.fixed_outer_values == [] + + +def test_generic_interferometer_seed_is_reproducible(): + gi_a = GenericInterferometer(start_mode=0, span=3, trainable=False, seed=99) + gi_b = GenericInterferometer(start_mode=0, span=3, trainable=False, seed=99) + assert gi_a.fixed_inner_values == gi_b.fixed_inner_values + assert gi_a.fixed_outer_values == gi_b.fixed_outer_values + + +def test_generic_interferometer_no_seed_varies_between_instances(): + gi_a = GenericInterferometer(start_mode=0, span=3, trainable=False) + gi_b = GenericInterferometer(start_mode=0, span=3, trainable=False) + assert gi_a.fixed_inner_values != gi_b.fixed_inner_values + + +@pytest.mark.parametrize("field_name", ["fixed_inner_values", "fixed_outer_values"]) +def test_generic_interferometer_rejects_short_fixed_values(field_name): + """A provided fixed-values list shorter than the shifter count must raise. + + Without this check, the circuit builder silently substitutes 0.0 for the + missing indices, reintroducing the zero-phase bug that random + initialization of non-trainable shifters is meant to prevent. + """ + # span=4 -> 6 phase shifters; provide only 3 values + with pytest.raises(ValueError, match=f"{field_name} must contain exactly 6"): + GenericInterferometer( + start_mode=0, + span=4, + trainable=False, + **{field_name: [0.1, 0.2, 0.3]}, + ) + + +@pytest.mark.parametrize("field_name", ["fixed_inner_values", "fixed_outer_values"]) +def test_generic_interferometer_rejects_overlong_fixed_values(field_name): + with pytest.raises(ValueError, match=f"{field_name} must contain exactly 1"): + GenericInterferometer( + start_mode=0, + span=2, + trainable=False, + **{field_name: [0.1, 0.2]}, + ) + + +def test_generic_interferometer_accepts_complete_fixed_values(): + """Explicit lists covering every shifter are preserved as-is.""" + inner = [0.1, 0.2, 0.3] + outer = [0.4, 0.5, 0.6] + gi = GenericInterferometer( + start_mode=0, + span=3, + trainable=False, + fixed_inner_values=inner, + fixed_outer_values=outer, + ) + assert gi.fixed_inner_values == inner + assert gi.fixed_outer_values == outer \ No newline at end of file diff --git a/tests/core/test_merlin_processor_unit.py b/tests/core/test_merlin_processor_unit.py index b640fb11f..886423144 100644 --- a/tests/core/test_merlin_processor_unit.py +++ b/tests/core/test_merlin_processor_unit.py @@ -25,6 +25,7 @@ from merlin.algorithms.module import MerlinModule from merlin.builder.circuit_builder import CircuitBuilder from merlin.core.circuit import Circuit +from merlin.core.components import GenericInterferometer from merlin.core.computation_space import ComputationSpace from merlin.core.merlin_processor import ( BackendCapabilities, @@ -1527,7 +1528,9 @@ class ReuploadInput(torch.nn.Module): def forward(self, input_tensor: torch.Tensor) -> torch.Tensor: return torch.cat((input_tensor, input_tensor), dim=-1) - def make_builder_layer(prefixes: Sequence[str]) -> QuantumLayer: + def make_builder_layer( + prefixes: Sequence[str], + ) -> tuple[QuantumLayer, list[GenericInterferometer]]: builder = CircuitBuilder(n_modes=n_modes) builder.add_entangling_layer(trainable=False, name=f"{prefixes[0]}_pre") for index, prefix in enumerate(prefixes): @@ -1539,7 +1542,12 @@ def make_builder_layer(prefixes: Sequence[str]) -> QuantumLayer: trainable=False, name=entangler_name, ) - return QuantumLayer( + entanglers = [ + component + for component in builder.circuit.components + if isinstance(component, GenericInterferometer) + ] + layer = QuantumLayer( input_size=n_modes * len(prefixes), builder=builder, input_state=[1, 0, 0], @@ -1548,14 +1556,27 @@ def make_builder_layer(prefixes: Sequence[str]) -> QuantumLayer: ), dtype=torch.float64, ).eval() + return layer, entanglers def make_builder_equivalent_perceval_circuit( - prefixes: Sequence[str], + prefixes: Sequence[str], entanglers: Sequence[GenericInterferometer] ) -> pcvl.Circuit: - def fixed_mzi(_index: int) -> pcvl.Circuit: - return pcvl.BS() // pcvl.PS(0.0) // pcvl.BS() // pcvl.PS(0.0) + # Non-trainable entangling layers get random (not zero) fixed phases; + # reuse each layer's actual drawn values so this reference circuit + # matches the builder-produced one exactly instead of assuming 0.0. + entangler_iter = iter(entanglers) def add_fixed_entangler(circuit: pcvl.Circuit) -> None: + entangler = next(entangler_iter) + + def fixed_mzi(index: int) -> pcvl.Circuit: + return ( + pcvl.BS() + // pcvl.PS(entangler.fixed_inner_values[index]) + // pcvl.BS() + // pcvl.PS(entangler.fixed_outer_values[index]) + ) + circuit.add( 0, pcvl.GenericInterferometer( @@ -1603,10 +1624,12 @@ def run_perceval_probabilities( output[row_index, state_to_index[str(state)]] = float(probability) return output - first_circuit = make_builder_equivalent_perceval_circuit(["a"]) - second_circuit = make_builder_equivalent_perceval_circuit(["b", "c"]) - first_layer = make_builder_layer(["a"]) - second_layer = make_builder_layer(["b", "c"]) + first_layer, first_entanglers = make_builder_layer(["a"]) + second_layer, second_entanglers = make_builder_layer(["b", "c"]) + first_circuit = make_builder_equivalent_perceval_circuit(["a"], first_entanglers) + second_circuit = make_builder_equivalent_perceval_circuit( + ["b", "c"], second_entanglers + ) model = torch.nn.Sequential(first_layer, ReuploadInput(), second_layer).eval() input_tensor = torch.tensor( [[0.1, 0.7, 1.4], [1.2, 0.3, 0.6], [2.4, 1.1, 0.2]], diff --git a/tests/pcvl_pytorch/test_noisy_g2_slos.py b/tests/pcvl_pytorch/test_noisy_g2_slos.py index 5b6081516..e496bec79 100644 --- a/tests/pcvl_pytorch/test_noisy_g2_slos.py +++ b/tests/pcvl_pytorch/test_noisy_g2_slos.py @@ -1109,3 +1109,94 @@ def test_g2_layer_to_moves_per_sector_transforms() -> None: layer = _g2_layer_with_loss() layer.to("cpu") + +class TestG2AugmentedObbRegression: + """Pin the indistinguishable extra-photon sector probabilities. + + Regression guard for ``NoisyG2SLOSComputeGraph._augmented_obb_probs``: + with partial indistinguishability, each extra-photon sector must reuse the + OBB partitions already built for the base ``n_photons`` input state and + fold every g2-emitted sibling into its partition cell's coherent group. + An earlier implementation instead grew a whole new noisy SLOS graph on the + augmented input state, re-drawing an independent + coherent-vs-distinguishable outcome for the g2 sibling; on this fixture + that skews sector 1 by up to 6.8e-3 per entry (e.g. the (1, 1, 1) entry + comes out ~0.031014 instead of the correct 0.024246), far beyond the pin + tolerance below. + + Golden values were produced by the fixed implementation on a + deterministic 3-mode Fourier interferometer in float64 and match Perceval + 1.2.4 ``Simulator.probs_svd`` with ``NoiseModel(g2=0.3, + g2_distinguishable=False, indistinguishability=0.6)`` to below 1e-15. + """ + + # Fock basis order is Combinadics("fock", n, m).enumerate_states(). + EXPECTED_SECTOR_1 = [ + 0.030686209636, # (3, 0, 0) + 0.038768111696, # (2, 1, 0) + 0.038768111696, # (2, 0, 1) + 0.038768111696, # (1, 2, 0) + 0.024245706182, # (1, 1, 1) + 0.038768111696, # (1, 0, 2) + 0.030686209636, # (0, 3, 0) + 0.038768111696, # (0, 2, 1) + 0.038768111696, # (0, 1, 2) + 0.030686209636, # (0, 0, 3) + ] + EXPECTED_SECTOR_2 = [ + 0.002972157995, # (4, 0, 0) + 0.003723146497, # (3, 1, 0) + 0.003723146497, # (3, 0, 1) + 0.004881425266, # (2, 2, 0) + 0.001597365052, # (2, 1, 1) + 0.004881425266, # (2, 0, 2) + 0.003723146497, # (1, 3, 0) + 0.001597365052, # (1, 2, 1) + 0.001597365052, # (1, 1, 2) + 0.003723146497, # (1, 0, 3) + 0.002972157995, # (0, 4, 0) + 0.003723146497, # (0, 3, 1) + 0.004881425266, # (0, 2, 2) + 0.003723146497, # (0, 1, 3) + 0.002972157995, # (0, 0, 4) + ] + + def test_extra_photon_sectors_match_golden_values(self): + """Extra-photon sectors equal the Perceval-verified golden values.""" + m = 3 + n_photons = 2 + modes = np.arange(m) + fourier = np.exp(2j * np.pi * np.outer(modes, modes) / m) / np.sqrt(m) + unitary = torch.tensor(fourier, dtype=torch.complex128) + + groups = NoiseGroups( + source={ + "g2": 0.3, + "g2_distinguishable": False, + "indistinguishability": 0.6, + }, + circuit=None, + post_measurement=None, + ) + noisy_slos = NoisyG2SLOSComputeGraph( + groups, + m=m, + n_photons=n_photons, + computation_space=ComputationSpace.FOCK, + dtype=torch.float64, + ) + result = noisy_slos.compute_probs(unitary, [1, 1, 0]) + + assert isinstance(result, SectoredDistribution) + assert len(result.sectors) == n_photons + 1 + + for sector_idx, expected in ( + (1, self.EXPECTED_SECTOR_1), + (2, self.EXPECTED_SECTOR_2), + ): + actual = result.sectors[sector_idx].tensor.squeeze() + expected_tensor = torch.tensor(expected, dtype=torch.float64) + assert torch.allclose(actual, expected_tensor, atol=1e-8), ( + f"Sector {sector_idx} deviates from golden values by " + f"{(actual - expected_tensor).abs().max().item():.3e}" + )