Move MIS sampler into backend implementation

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,

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