diff --git a/perceval/backends/__init__.py b/perceval/backends/__init__.py index 84e2736e..2920be3f 100644 --- a/perceval/backends/__init__.py +++ b/perceval/backends/__init__.py @@ -29,6 +29,7 @@ from ._abstract_backends import ABackend, ASamplingBackend, AStrongSimulationBackend, IFFBackend from ._clifford2017 import Clifford2017Backend +from ._mis import MISBackend from ._mps import MPSBackend from ._naive import NaiveBackend from ._naive_approx import NaiveApproxBackend @@ -38,6 +39,7 @@ BACKEND_LIST = { "CliffordClifford2017": Clifford2017Backend, + "MIS": MISBackend, "MPS": MPSBackend, "Naive": NaiveBackend, "NaiveApprox": NaiveApproxBackend, diff --git a/perceval/backends/_mis.py b/perceval/backends/_mis.py new file mode 100644 index 00000000..40c9cdc0 --- /dev/null +++ b/perceval/backends/_mis.py @@ -0,0 +1,140 @@ +# MIT License +# +# Copyright (c) 2022 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. +# +# As a special exception, the copyright holders of exqalibur library give you +# permission to combine exqalibur with code included in the standard release of +# Perceval under the MIT license (or modified versions of such code). You may +# copy and distribute such a combined system following the terms of the MIT +# license for both exqalibur and Perceval. This exception for the usage of +# exqalibur is limited to the python bindings used by Perceval. +# +# 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. + +from __future__ import annotations + +import numpy as np + +from perceval.utils import FockState +from perceval.components import ACircuit + +from ._abstract_backends import ASamplingBackend +from ._slos import SLOSBackend + + +class MISBackend(ASamplingBackend): + """ + Metropolised Independence Sampling (MIS) backend for boson sampling. + + The proposal distribution is generated by sampling each photon independently + from |U|^2 (distinguishable photons), while acceptance uses the exact + probability from a strong simulator. + + :param burn_in: number of initial steps discarded when initializing the chain + :param seed: optional RNG seed for reproducibility + """ + + def __init__(self, burn_in: int = 0, seed: int | None = None): + super().__init__() + if burn_in < 0: + raise ValueError("burn_in must be a non-negative integer") + self._burn_in = int(burn_in) + self._rng = np.random.default_rng(seed) + self._transition_probs = None + self._input_modes = None + self._current_state = None + self._current_g = None + self._current_p = None + self._burned_in = False + self._slos = SLOSBackend() + + @property + def name(self) -> str: + return "MIS" + + def set_circuit(self, circuit: ACircuit): + super().set_circuit(circuit) + self._transition_probs = np.abs(np.asarray(self._umat)) ** 2 + self._current_state = None + self._current_g = None + self._current_p = None + self._burned_in = False + self._slos.set_circuit(circuit) + + def set_input_state(self, input_state: FockState): + super().set_input_state(input_state) + self._input_modes = [mode for mode, count in enumerate(input_state) for _ in range(count)] + self._current_state = None + self._current_g = None + self._current_p = None + self._burned_in = False + self._slos.set_input_state(input_state) + + def _distinguishable_proposal(self) -> tuple[FockState, float]: + counts = [0] * self._circuit.m + prob = 1.0 + for mode in self._input_modes: + probs = self._transition_probs[:, mode] + chosen = self._rng.choice(self._circuit.m, p=probs) + counts[chosen] += 1 + prob *= probs[chosen] + return FockState(counts), prob + + def _target_probability(self, output_state: FockState) -> float: + return self._slos.probability(output_state) + + def _ensure_chain(self): + if self._current_state is None: + self._current_state, self._current_g = self._distinguishable_proposal() + self._current_p = self._target_probability(self._current_state) + + if not self._burned_in and self._burn_in > 0: + for _ in range(self._burn_in): + self._step_chain() + self._burned_in = True + + def _step_chain(self): + proposal_state, proposal_g = self._distinguishable_proposal() + proposal_p = self._target_probability(proposal_state) + + if self._current_p == 0 and proposal_p > 0: + accept = True + elif proposal_p == 0: + accept = False + else: + ratio = (proposal_p * self._current_g) / (self._current_p * proposal_g) + accept = self._rng.random() < min(1.0, ratio) + + if accept: + self._current_state = proposal_state + self._current_g = proposal_g + self._current_p = proposal_p + + def sample(self): + self._ensure_chain() + self._step_chain() + return self._current_state + + def samples(self, count: int): + self._ensure_chain() + results = [] + for _ in range(count): + self._step_chain() + results.append(self._current_state) + return results