Bosonic engine results (#546)

* Bosonic backend fixes

- Allows for multiple shots for measurements (new test)
- Raises an error if symbolic non-Gaussian state prep is attempted

* Bosonic engine

Defines the BosonicEngine class

Co-Authored-By: ilan-tz <57886357+ilan-tz@users.noreply.github.com>
Co-Authored-By: Josh Izaac <josh@iza.ac>
Co-Authored-By: Nicolas Quesada <991946+nquesada@users.noreply.github.com>

* Adds ancilla samples to Results class

Ancilla samples are produced by ancillary modes used for measurement-based gates

* Integration tests

Updates/adds new integration tests to apply to the bosonic backend

* Codefactor changes

* Codefactor changes

* Test for parameter_checker function

boosts code coverage

* Update engine

Deletes the run method from the BosonicEngine, since it can be inherited from LocalEngine with a minimal change.

* Fix empty circuit bug

* disable too-many-branches for run_prog

* Ran black on test_bosonic_backend

* Apply suggestions from code review

Co-authored-by: antalszava <antalszava@gmail.com>

* Implemented changes from code review

- discovered and fixed bug in array dimension of returned samples
- makes tests more readable
- adds ancilla information to string of Results __repr__
- adds BosonicEngine to docs

Co-Authored-By: antalszava <24476053+antalszava@users.noreply.github.com>

* Code factor + test tidying

* Use isinstance and Iterable to check type

Co-Authored-By: Theodor <6934626+thisac@users.noreply.github.com>

* Apply suggestions from code review

Co-authored-by: antalszava <antalszava@gmail.com>

* Applied suggestions from code review

Co-Authored-By: antalszava <24476053+antalszava@users.noreply.github.com>

* Apply suggestions from code review

Co-authored-by: ilan-tz <57886357+ilan-tz@users.noreply.github.com>
Co-authored-by: Josh Izaac <josh@iza.ac>
Co-authored-by: Nicolas Quesada <991946+nquesada@users.noreply.github.com>
Co-authored-by: antalszava <antalszava@gmail.com>
Co-authored-by: antalszava <24476053+antalszava@users.noreply.github.com>
Co-authored-by: Theodor <6934626+thisac@users.noreply.github.com>
Co-authored-by: Theodor <theodor@xanadu.ai>

strawberryfields/api/result.py

@@ -32,6 +32,9 @@ class Result:
 
     * ``results.samples``: Measurement samples from any measurements performed.
 
+    * ``results.ancilla_samples``: Measurement samples from any ancillary states
+      used for measurement-based gates.
+
     **Example:**
 
     The following example runs an existing Strawberry Fields
@@ -58,11 +61,12 @@ class Result:
         but the return value of ``Result.state`` will be ``None``.
     """
 
-    def __init__(self, samples, all_samples=None, is_stateful=True):
+    def __init__(self, samples, all_samples=None, is_stateful=True, ancilla_samples=None):
         self._state = None
         self._is_stateful = is_stateful
         self._samples = samples
         self._all_samples = all_samples
+        self._ancilla_samples = ancilla_samples
 
     @property
     def samples(self):
@@ -90,6 +94,21 @@ def all_samples(self):
         """
         return self._all_samples
 
+    @property
+    def ancilla_samples(self):
+        """All measurement samples from ancillary modes used for measurement-based
+        gates.
+
+        Returns a dictionary which associates each mode (keys) with the
+        list of measurements outcomes (values) from all the ancilla-assisted
+        gates applied to that mode.
+
+        Returns:
+            dict[int, list]: mode index associated with the list of ancilla
+            measurement outcomes
+        """
+        return self._ancilla_samples
+
     @property
     def state(self):
         """The quantum state object.
@@ -119,6 +138,16 @@ def __repr__(self):
         if self.samples.ndim == 2:
             # if samples has dim 2, assume they're from a standard Program
             shots, modes = self.samples.shape
+
+            if self.ancilla_samples is not None:
+                ancilla_modes = 0
+                for i in self.ancilla_samples.keys():
+                    ancilla_modes += len(self.ancilla_samples[i])
+                return (
+                    f"<Result: shots={shots}, num_modes={modes}, num_ancillae={ancilla_modes}, "
+                    f"contains state={self._is_stateful}>"
+                )
+
             return "<Result: shots={}, num_modes={}, contains state={}>".format(
                 shots, modes, self._is_stateful
             )

strawberryfields/backends/bosonicbackend/backend.py

@@ -16,6 +16,7 @@
 """Bosonic backend"""
 import itertools as it
 from functools import reduce
+from collections.abc import Iterable
 
 import numpy as np
 
@@ -28,13 +29,28 @@
 
 from strawberryfields.backends.bosonicbackend.bosoniccircuit import BosonicModes
 from strawberryfields.backends.base import NotApplicableError
+from strawberryfields.program_utils import CircuitError
+import sympy
 
 
 def kron_list(l):
     """Take Kronecker products of a list of lists."""
     return reduce(np.kron, l)
 
 
+def parameter_checker(parameters):
+    """Checks if any items in an iterable are sympy objects."""
+    for item in parameters:
+        if isinstance(item, sympy.Expr):
+            return True
+
+        # This checks all the nested items if item is an iterable
+        if isinstance(item, Iterable) and not isinstance(item, str):
+            if parameter_checker(item):
+                return True
+    return False
+
+
 class BosonicBackend(BaseBosonic):
     r"""The BosonicBackend implements a simulation of quantum optical circuits
     in NumPy by representing states as linear combinations of Gaussian functions
@@ -69,6 +85,7 @@ def __init__(self):
         self.circuit = None
         self.ancillae_samples_dict = {}
 
+    # pylint: disable=too-many-branches
     # pylint: disable=import-outside-toplevel
     def run_prog(self, prog, **kwargs):
         """Runs a strawberryfields program using the bosonic backend.
@@ -95,8 +112,9 @@ def run_prog(self, prog, **kwargs):
             _New_modes,
         )
 
-        # Initialize the circuit. This applies all non-Gaussian state-prep
-        self.init_circuit(prog)
+        # If a circuit exists, initialize the circuit. This applies all non-Gaussian state-prep
+        if prog.circuit:
+            self.init_circuit(prog)
 
         # Apply operations to circuit. For now, copied from LocalEngine;
         # only change is to ignore preparation classes and ancilla-assisted gates
@@ -124,9 +142,9 @@ def run_prog(self, prog, **kwargs):
                 if val is not None:
                     for i, r in enumerate(cmd.reg):
                         if r.ind not in self.ancillae_samples_dict.keys():
-                            self.ancillae_samples_dict[r.ind] = [val[:, i]]
+                            self.ancillae_samples_dict[r.ind] = [val]
                         else:
-                            self.ancillae_samples_dict[r.ind].append(val[:, i])
+                            self.ancillae_samples_dict[r.ind].append(val)
 
                 applied.append(cmd)
 
@@ -162,7 +180,6 @@ def run_prog(self, prog, **kwargs):
 
         return applied, samples_dict, all_samples
 
-    # pylint: disable=too-many-branches
     # pylint: disable=import-outside-toplevel
     def init_circuit(self, prog):
         """Instantiate the circuit and initialize weights, means, and covs
@@ -173,6 +190,8 @@ def init_circuit(self, prog):
 
         Raises:
             NotImplementedError: if ``Ket`` or ``DensityMatrix`` preparation used
+            CircuitError: if any of the parameters for non-Gaussian state preparation
+                are symbolic
         """
         from strawberryfields.ops import (
             Bosonic,
@@ -184,7 +203,10 @@ def init_circuit(self, prog):
             _New_modes,
         )
 
-        non_gauss_preps = (Bosonic, Catstate, DensityMatrix, Fock, GKP, Ket)
+        # _New_modes is what gets checked when New() is called in a program circuit.
+        # It is included here since it could be used to instantiate a mode for non-Gaussian
+        # state preparation, and it's best to initialize any new modes from the outset.
+        non_gauss_preps = (Bosonic, Catstate, DensityMatrix, Fock, GKP, Ket, _New_modes)
         nmodes = prog.init_num_subsystems
         self.begin_circuit(nmodes)
         # Dummy initial weights, means and covs
@@ -204,6 +226,12 @@ def init_circuit(self, prog):
             if np.any(isitnew):
                 # Operation parameters
                 pars = cmd.op.p
+                # Check if any of the preparations rely on symbolic quantities
+                if isinstance(cmd.op, non_gauss_preps) and parameter_checker(pars):
+                    raise CircuitError(
+                        "Symbolic non-Gaussian preparations have not been implemented "
+                        "in the bosonic backend."
+                    )
                 for reg in labels:
                     # All the possible preparations should go in this loop
                     if isinstance(cmd.op, Bosonic):
@@ -361,7 +389,7 @@ def prepare_cat(self, alpha, phi, representation, cutoff, D):
         :math:`\ket{\text{cat}(\alpha)} = \frac{1}{N} (\ket{\alpha} +e^{i\phi\pi} \ket{-\alpha})`.
 
         Args:
-            alpha (float): alpha value of cat state
+            alpha (complex): alpha value of cat state
             phi (float): phi value of cat state
             representation (str): whether to use the ``'real'`` or ``'complex'`` representation
             cutoff (float): if using the ``'real'`` representation, this determines
@@ -422,7 +450,7 @@ def prepare_cat_real_rep(self, alpha, phi, cutoff, D):
         For this representation, weights, means and covariances are real-valued.
 
         Args:
-            alpha (float): alpha value of cat state
+            alpha (complex): alpha value of cat state
             phi (float): phi value of cat state
             cutoff (float): this determines how many terms to keep
             D (float): quality parameter of approximation
@@ -739,17 +767,18 @@ def measure_homodyne(self, phi, mode, shots=1, select=None, **kwargs):
         self.circuit.phase_shift(-phi, mode)
 
         if select is None:
-            val = self.circuit.homodyne(mode, **kwargs)[0, 0]
+            val = self.circuit.homodyne(mode, shots=shots, **kwargs)[:, 0]
         else:
             val = select * 2 / np.sqrt(2 * self.circuit.hbar)
-            self.circuit.post_select_homodyne(mode, val, **kwargs)
+            self.circuit.post_select_homodyne(mode, val)
+            val = np.array([val])
 
-        return np.array([[val * np.sqrt(2 * self.circuit.hbar) / 2]])
+        return np.array([val]).T * np.sqrt(2 * self.circuit.hbar) / 2
 
     def measure_heterodyne(self, mode, shots=1, select=None):
         if select is None:
             res = 0.5 * self.circuit.heterodyne(mode, shots=shots)
-            return np.array([[res[:, 0] + 1j * res[:, 1]]])
+            return np.array([res[:, 0] + 1j * res[:, 1]]).T
 
         res = select
         self.circuit.post_select_heterodyne(mode, select)

strawberryfields/backends/bosonicbackend/bosoniccircuit.py

@@ -353,7 +353,7 @@ def mb_squeeze_single_shot(self, k, r, phi, r_anc, eta_anc):
         self.beamsplitter(theta, 0, k, new_mode)
         self.loss(eta_anc, new_mode)
         self.phase_shift(np.pi / 2, new_mode)
-        val = self.homodyne(new_mode)
+        val = self.homodyne(new_mode)[0][0]
 
         # Delete all record of ancilla mode
         self.del_mode(new_mode)
@@ -367,7 +367,7 @@ def mb_squeeze_single_shot(self, k, r, phi, r_anc, eta_anc):
 
         # Feedforward displacement
         prefac = -np.tan(theta) / np.sqrt(2 * self.hbar * eta_anc)
-        self.displace(prefac * val[0][0], np.pi / 2, k)
+        self.displace(prefac * val, np.pi / 2, k)
 
         self.phase_shift(phi / 2, k)
 

strawberryfields/engine.py

@@ -33,7 +33,7 @@
 from .backends.base import BaseBackend, NotApplicableError
 
 # for automodapi, do not include the classes that should appear under the top-level strawberryfields namespace
-__all__ = ["BaseEngine", "LocalEngine"]
+__all__ = ["BaseEngine", "LocalEngine", "BosonicEngine"]
 
 
 class BaseEngine(abc.ABC):
@@ -317,6 +317,14 @@ class LocalEngine(BaseEngine):
         backend_options (None, Dict[str, Any]): keyword arguments to be passed to the backend
     """
 
+    def __new__(cls, backend, *, backend_options=None):
+        if backend == "bosonic":
+            bos_eng = super().__new__(BosonicEngine)
+            bos_eng.__init__(backend, backend_options=backend_options)
+            return bos_eng
+
+        return super().__new__(cls)
+
     def __str__(self):
         return self.__class__.__name__ + "({})".format(self.backend_name)
 
@@ -501,7 +509,8 @@ def run(self, program, *, args=None, compile_options=None, **kwargs):
             # session and feed_dict are needed by TF backend both during simulation (if program
             # contains measurements) and state object construction.
             result._state = self.backend.state(**temp_run_options)
-
+            if self.backend_name == "bosonic":
+                result._ancilla_samples = self.backend.ancillae_samples_dict.copy()
         return result
 
 
@@ -723,3 +732,17 @@ class Engine(LocalEngine):
 
     # alias for backwards compatibility
     __doc__ = LocalEngine.__doc__
+
+
+class BosonicEngine(LocalEngine):
+    """Local quantum program executor engine for programs executed on the bosonic backend.
+
+    The BosonicEngine is used to execute :class:`.Program` instances on the bosonic backend,
+    and makes the results available via :class:`.Result`.
+    """
+
+    def _run_program(self, prog, **kwargs):
+        # Custom Bosonic run code
+        applied, samples_dict, all_samples = self.backend.run_prog(prog, **kwargs)
+        samples = self._combine_and_sort_samples(samples_dict)
+        return applied, samples, all_samples

tests/backend/test_bosonic_backend.py

@@ -20,6 +20,7 @@
 import strawberryfields as sf
 import strawberryfields.backends.bosonicbackend.backend as bosonic
 import pytest
+import sympy
 
 pytestmark = pytest.mark.bosonic
 
@@ -29,6 +30,7 @@
 r_fock = 0.05
 EPS_VALS = np.array([0.01, 0.05, 0.1, 0.5])
 R_VALS = np.linspace(-1, 1, 5)
+SHOTS_VALS = np.array([1, 19, 100])
 
 
 class TestKronList:
@@ -41,6 +43,24 @@ def test_kron_list(self):
         assert np.allclose(list_compare, bosonic.kron_list([l1, l2]))
 
 
+class TestParameterChecker:
+    """Test parameter_checker function from the bosonic backend."""
+
+    def test_parameter_checker(self):
+        symbolic_param = sympy.Expr()
+        params = []
+        assert not bosonic.parameter_checker(params)
+
+        params = [1, "real", 3.0, [4 + 1j, 5], [[3.5, 4.8, "complex"]], np.array([7, 9]), range(3)]
+        assert not bosonic.parameter_checker(params)
+
+        params = [1, "real", symbolic_param]
+        assert bosonic.parameter_checker(params)
+
+        params = [1, [3, symbolic_param]]
+        assert bosonic.parameter_checker(params)
+
+
 class TestBosonicCatStates:
     r"""Tests cat state method of the BosonicBackend class."""
 
@@ -493,7 +513,7 @@ def test_measurement(self, alpha, r):
             sf.ops.Catstate(alpha) | q[0]
             sf.ops.Squeezed(r) | q[1]
             sf.ops.MeasureX | q[0]
-            sf.ops.MeasureX | q[1]
+            sf.ops.MeasureHD | q[1]
         backend = bosonic.BosonicBackend()
         applied, samples, all_samples = backend.run_prog(prog)
         state = backend.state()
@@ -506,6 +526,26 @@ def test_measurement(self, alpha, r):
             assert i in samples.keys()
             assert samples[i].shape == (1,)
 
+    @pytest.mark.parametrize("alpha", ALPHA_VALS)
+    @pytest.mark.parametrize("shots", SHOTS_VALS)
+    def test_measurement_many_shots(self, alpha, shots):
+        """Runs a program with measurements."""
+        prog = sf.Program(1)
+        with prog.context as q:
+            sf.ops.Catstate(alpha) | q[0]
+            sf.ops.MeasureHomodyne(0) | q[0]
+
+        backend = bosonic.BosonicBackend()
+        applied, samples, all_samples = backend.run_prog(prog, shots=shots)
+        state = backend.state()
+
+        # Check output is vacuum since everything was measured
+        assert np.allclose(state.fidelity_vacuum(), 1)
+
+        # Check samples
+        assert 0 in samples.keys()
+        assert samples[0].shape == (int(shots),)
+
     @pytest.mark.parametrize("alpha", ALPHA_VALS)
     @pytest.mark.parametrize("r", R_VALS)
     def test_mb_gates(self, alpha, r):