Adds Torontonian sampling support to the Gaussian backend

.github/CHANGELOG.md

@@ -17,6 +17,10 @@
   allowing a bipartite graph to be embedded on a device that allows for
   initial two-mode squeezed states, and block diagonal unitaries.
 
+* Added support for threshold measurement in the Gaussian backend, via the new backend API
+  method `measure_threshold`.
+  [#152](https://github.com/XanaduAI/strawberryfields/pull/152)
+
 ### API Changes
 
 * The `strawberryfields.ops.Measure` shorthand has been deprecated in favour

doc/conf.py

@@ -49,10 +49,10 @@ class TypeMock(type):
     'numbers',
     'blackbird',
     'blackbird.utils',
-    'hafnian',
-    'hafnian.lib',
-    'hafnian.samples',
-    'hafnian.quantum',
+    'thewalrus',
+    'thewalrus.libwalrus',
+    'thewalrus.samples',
+    'thewalrus.quantum',
     ]
 
 np_math_fns = ['abs',

doc/gallery/scattershot-boson-sampling/scattershot-bs.ipynb

@@ -292,7 +292,7 @@
     "\n",
     "Using that, we can now perform the computation.\n",
     "\n",
-    "First, the permanent of the matrix can be calculated via the [Hafnian](https://hafnian.readthedocs.io) library:"
+    "First, the permanent of the matrix can be calculated via [The Walrus](https://the-walrus.readthedocs.io) library:"
    ]
   },
   {
@@ -301,7 +301,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "from hafnian import perm"
+    "from thewalrus import perm"
    ]
   },
   {

doc/installing.rst

@@ -21,7 +21,7 @@ as well as the following Python packages:
 * `SciPy <http://scipy.org/>`_  >=1.0.0
 * `NetworkX <http://networkx.github.io/>`_ >=2.0
 * `Blackbird <https://quantum-blackbird.readthedocs.io>`_ >=0.2.0
-* `Hafnian <https://hafnian.readthedocs.io>`_ >=0.6.0
+* `The Walrus <https://the-walrus.readthedocs.io>`_ >=0.7.0
 * `toml <https://pypi.org/project/toml/>`_
 * `appdirs <https://pypi.org/project/appdirs/>`_
 

doc/tutorials/tutorial_boson_sampling.rst

@@ -180,9 +180,9 @@ For this example, we'll consider the output state :math:`\ket{2,0,0,1}`. Extract
 >>> probs[2,0,0,1]
 0.10644192724642336
 
-Before we can calculate the right hand side of equation, we need a method of calculating the permanent. Since the permanent is classically hard to compute, it is not provided in either NumPy *or* SciPy, so we will use `hafnian <https://hafnian.readthedocs.io>`_ package, installed alongside Strawberry Fields:
+Before we can calculate the right-hand-side of the equation, we need a method of calculating the permanent. Since the permanent is classically hard to compute, it is not provided in either NumPy *or* SciPy, so we will use `The Walrus <https://the-walrus.readthedocs.io>`_ library, installed alongside Strawberry Fields:
 
->>> from hafnian import perm
+>>> from thewalrus import perm
 
 Finally, how do we determine the submatrix :math:`U_{st}`? This isn't too hard :cite:`tillmann2013`; first, we calculate the submatrix :math:`U_s` by taking :math:`m_k` copies of the :math:`k\text{th}` **columns** of :math:`U`, where :math:`m_k` is the photon number of the :math:`k\text{th}` input state.
 

doc/tutorials/tutorial_gaussian_boson_sampling.rst

@@ -122,9 +122,9 @@ Now that we have the interferometer unitary transformation :math:`U`, as well as
 Calculating the hafnian
 ------------------------
 
-Before we can calculate the right hand side of the Gaussian boson sampling equation, we need a method of calculating the hafnian. Since the hafnian is classically hard to compute, it is not provided in either NumPy *or* SciPy, so we will use `hafnian <https://hafnian.readthedocs.io>`_ package, installed alongside Strawberry Fields:
+Before we can calculate the right hand side of the Gaussian boson sampling equation, we need a method of calculating the hafnian. Since the hafnian is classically hard to compute, it is not provided in either NumPy *or* SciPy, so we will use `The Walrus <https://the-walrus.readthedocs.io>`_ library, installed alongside Strawberry Fields:
 
->>> from hafnian import haf
+>>> from thewalrus import haf
 
 Now, for the right hand side numerator, we first calculate the submatrix :math:`[(UU^T\tanh(r))]_{st}`:
 

requirements.txt

@@ -4,6 +4,6 @@ scipy>=1.0.0
 tensorflow==1.3
 tensorflow-tensorboard>=0.1.8
 quantum-blackbird
-hafnian>=0.6
+thewalrus>=0.7
 toml
 appdirs

setup.py

@@ -28,7 +28,7 @@
     "scipy>=1.0.0",
     "networkx>=2.0",
     "quantum-blackbird>=0.2.0",
-    "hafnian>=0.6",
+    "thewalrus>=0.7",
     "toml",
     "appdirs"
 ]

strawberryfields/__init__.py

@@ -98,7 +98,7 @@ def about():
     import os
     import numpy
     import scipy
-    import hafnian
+    import thewalrus
     import blackbird
 
     # a QuTiP-style infobox
@@ -111,7 +111,7 @@ def about():
     print('Strawberry Fields version: {}'.format(__version__))
     print('Numpy version:             {}'.format(numpy.__version__))
     print('Scipy version:             {}'.format(scipy.__version__))
-    print('Hafnian version:           {}'.format(hafnian.__version__))
+    print('The Walrus version:           {}'.format(thewalrus.__version__))
     print('Blackbird version:         {}'.format(blackbird.__version__))
 
     try:

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 thewalrus.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/backends/gaussianbackend/ops.py

@@ -17,7 +17,7 @@
 
 from scipy.linalg import sqrtm
 
-from hafnian.quantum import density_matrix_element, is_pure_cov, pure_state_amplitude, state_vector, density_matrix
+from thewalrus.quantum import density_matrix_element, is_pure_cov, pure_state_amplitude, state_vector, density_matrix
 
 
 

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/decompositions.py

@@ -47,7 +47,7 @@
 
 import numpy as np
 from scipy.linalg import block_diag, sqrtm, polar, schur
-from hafnian.quantum import find_scaling_adjacency_matrix
+from thewalrus.quantum import find_scaling_adjacency_matrix
 
 from .backends.shared_ops import sympmat, changebasis
 

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):
+    """Measures a set of modes with thresholded Fock-state measurements, i.e.,
+    measuring whether a mode contain 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,85 @@
+# 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)
+
+
+    def test_binary_outcome(self, setup_backend, pure):
+        """Test that the outcomes of a threshold measurement is zero or one."""
+        num_modes = 2
+        for _ in range(NUM_REPEATS):
+            backend = setup_backend(num_modes)
+            backend.reset(pure=pure)
+
+            r = 0.5
+            backend.squeeze(r, 0)
+            backend.beamsplitter(np.sqrt(0.5), -np.sqrt(0.5), 0, 1)
+            meas_modes = [0, 1]
+            meas_results = backend.measure_threshold(meas_modes)
+
+            for i in range(num_modes):
+                assert meas_results[i] == 0 or meas_results[i] == 1
+
+
+

tests/frontend/test_about.py

@@ -36,7 +36,7 @@ def test_about(capfd):
 
     assert "Numpy version" in out
     assert "Scipy version" in out
-    assert "Hafnian version" in out
+    assert "The Walrus version" in out
     assert "Blackbird version" in out
 
 

tests/frontend/test_ops_decompositions.py

@@ -17,7 +17,7 @@
 pytestmark = pytest.mark.frontend
 
 import numpy as np
-from hafnian.quantum import Amat
+from thewalrus.quantum import Amat
 
 import strawberryfields as sf
 from strawberryfields import decompositions as dec