Adds Torontonian sampling support to the Gaussian backend

strawberryfields/backends/base.py

@@ -506,6 +506,25 @@ def measure_fock(self, modes, shots=1, select=None, **kwargs):
         """
         raise NotImplementedError
 
+    def measure_threshold(self, modes, shots=1, select=None, **kwargs):
+        """Measure the given modes in the thresholded Fock basis, i.e., zero or nonzero photons).
+
+        .. note::
+
+            When :code:``shots == 1``, updates the current system state to the conditional state of that
+            measurement result. When :code:``shots > 1``, the system state is not updated.
+
+        Args:
+            modes (Sequence[int]): which modes to measure
+            shots (int): number of measurement samples to obtain
+            select (None or Sequence[int]): If not None: desired values of the measurement results.
+                Enables post-selection on specific measurement results instead of random sampling.
+                ``len(select) == len(modes)`` is required.
+        Returns:
+            tuple[int]: measurement results
+        """
+        raise NotImplementedError
+
     def is_vacuum(self, tol=0.0, **kwargs):
         r"""Test whether the current circuit state is vacuum (up to given tolerance).
 

strawberryfields/backends/gaussianbackend/backend.py

@@ -30,7 +30,7 @@
     ix_,
 )
 from numpy.linalg import inv
-from hafnian.samples import hafnian_sample_state
+from hafnian.samples import hafnian_sample_state, torontonian_sample_state
 
 from strawberryfields.backends import BaseGaussian
 from strawberryfields.backends.shared_ops import changebasis
@@ -196,7 +196,7 @@ def measure_fock(self, modes, shots=1, select=None):
         if shots != 1:
             if select is not None:
                 raise NotImplementedError("Gaussian backend currently does not support "
-                                          "postselection if shots != 1 for Fock measurement")
+                                          "postselection")
             warnings.warn("Cannot simulate non-Gaussian states. "
                           "Conditional state after Fock measurement has not been updated.")
 
@@ -217,6 +217,33 @@ def measure_fock(self, modes, shots=1, select=None):
             samples = samples.reshape((len(modes),))
         return samples
 
+
+    def measure_threshold(self, modes, shots=1, select=None):
+        if shots != 1:
+            if select is not None:
+                raise NotImplementedError("Gaussian backend currently does not support "
+                                          "postselection")
+            warnings.warn("Cannot simulate non-Gaussian states. "
+                          "Conditional state after Threshold measurement has not been updated.")
+
+        mu = self.circuit.mean
+        cov = self.circuit.scovmatxp()
+        # check we are sampling from a gaussian state with zero mean
+        if not allclose(mu, zeros_like(mu)):
+            raise NotImplementedError("Threshold measurement is only supported for "
+                                      "Gaussian states with zero mean")
+        x_idxs = array(modes)
+        p_idxs = x_idxs + len(mu)
+        modes_idxs = concatenate([x_idxs, p_idxs])
+        reduced_cov = cov[ix_(modes_idxs, modes_idxs)]
+        samples = torontonian_sample_state(reduced_cov, shots)
+        # for backward compatibility with previous measurement behaviour,
+        # if only one shot, then we drop the shots axis
+        if shots == 1:
+            samples = samples.reshape((len(modes),))
+        return samples
+
+
     def state(self, modes=None, **kwargs):
         """Returns the state of the quantum simulation.
 

strawberryfields/circuitspecs/gaussian.py

@@ -40,6 +40,7 @@ class GaussianSpecs(CircuitSpecs):
         "MeasureHomodyne",
         "MeasureHeterodyne",
         "MeasureFock",
+        "MeasureThreshold",
         # channels
         "LossChannel",
         "ThermalLossChannel",

strawberryfields/circuitspecs/gbs.py

@@ -39,6 +39,7 @@ class GBSSpecs(GaussianSpecs):
         "MeasureHomodyne",
         "MeasureHeterodyne",
         "MeasureFock",
+        "MeasureThreshold",
         # channels
         "LossChannel",
         "ThermalLossChannel",

strawberryfields/ops.py

@@ -190,6 +190,7 @@ def square(q):
 
 .. autosummary::
     MeasureFock
+    MeasureThreshold
     MeasureHomodyne
     MeasureHeterodyne
 
@@ -1052,6 +1053,23 @@ def _apply(self, reg, backend, shots=1, **kwargs):
         return backend.measure_fock(reg, shots=shots, select=self.select, **kwargs)
 
 
+class MeasureThreshold(Measurement):
+    """:ref:`photo_detection`: measures a set of modes in the thresholded Fock basis, i.e.
+    zero or nonzero photons.
+
+    After measurement, the modes are reset to the vacuum state.
+    """
+    ns = None
+
+    def __init__(self, select=None):
+        if select is not None and not isinstance(select, Sequence):
+            select = [select]
+        super().__init__([], select)
+
+    def _apply(self, reg, backend, shots=1, **kwargs):
+        return backend.measure_threshold(reg, shots=shots, select=self.select, **kwargs)
+
+
 class MeasureHomodyne(Measurement):
     r"""Performs a :ref:`homodyne measurement <homodyne>`, measures one quadrature of a mode.
 
@@ -2131,7 +2149,7 @@ def _decompose(self, reg, **kwargs):
 simple_state_preparations = (Vacuum, Coherent, Squeezed, DisplacedSqueezed, Fock, Catstate, Thermal)  # have __init__ methods with default arguments
 state_preparations = simple_state_preparations + (Ket, DensityMatrix)
 
-measurements = (MeasureFock, MeasureHomodyne, MeasureHeterodyne)
+measurements = (MeasureFock, MeasureHomodyne, MeasureHeterodyne, MeasureThreshold)
 
 decompositions = (Interferometer, GraphEmbed, GaussianTransform, Gaussian)
 

tests/backend/test_threshold_measurement.py

@@ -0,0 +1,65 @@
+# Copyright 2019 Xanadu Quantum Technologies Inc.
+
+# 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.
+
+r"""Unit tests for measurements in the Fock basis"""
+import pytest
+
+import numpy as np
+
+
+NUM_REPEATS = 50
+
+
+@pytest.mark.backends("gaussian")
+class TestGaussianRepresentation:
+    """Tests that make use of the Fock basis representation."""
+
+    def measure_threshold_gaussian_warning(self, setup_backend):
+        """Tests the warning message when MeasureThreshold is called."""
+
+        backend = setup_backend(3)
+
+        with pytest.warns(Warning, match="Cannot simulate non-Gaussian states. Conditional state after "
+                                         "Threshold measurement has not been updated."):
+            backend.measure_threshold([0, 1], shots=5)
+
+
+@pytest.mark.backends("gaussian")
+class TestRepresentationIndependent:
+    """Basic implementation-independent tests."""
+
+    def test_two_mode_squeezed_measurements(self, setup_backend, pure):
+        """Tests Threshold measurement on the two mode squeezed vacuum state."""
+        for _ in range(NUM_REPEATS):
+            backend = setup_backend(2)
+            backend.reset(pure=pure)
+
+            r = 0.25
+            # Circuit to prepare two mode squeezed vacuum
+            backend.squeeze(-r, 0)
+            backend.squeeze(r, 1)
+            backend.beamsplitter(np.sqrt(0.5), -np.sqrt(0.5), 0, 1)
+            meas_modes = [0, 1]
+            meas_results = backend.measure_threshold(meas_modes)
+            assert np.all(meas_results[0] == meas_results[1])
+
+    def test_vacuum_measurements(self, setup_backend, pure):
+        """Tests Threshold measurement on the vacuum state."""
+        backend = setup_backend(3)
+
+        for _ in range(NUM_REPEATS):
+            backend.reset(pure=pure)
+
+            meas = backend.measure_threshold([0, 1, 2])[0]
+            assert np.all(np.array(meas) == 0)