new passive circuit compiler (#600)

* new passive circuit compiler

* add PassiveChannel as valid gate for passive compiler

* generalise channel merging to multimode channels (ie PassiveChannel)

* PassiveChannel tests

* tests for passive compiler

* tests for passive compiler auxillary functions

* code factor stuff

* flip the conj in the apply_u method

* now featuring relevant docstrings

* dev version of thewalrus

* run black

* remove dependency on dev version of thewalrus

* fix docstring

* remove expand_passive from backend

* trailing whitespace

* remove numba

* remove numba imports

* minor comments

* minor comments

* minor comments

* black

* black

* declare hbar

* more functional gate application

* hbar in test args

* local tests are more happy

* run black

* better test coverage

* run black

Co-authored-by: Nicolas Quesada <zeitus@gmail.com>

strawberryfields/backends/gaussianbackend/backend.py

@@ -15,17 +15,7 @@
 """Gaussian backend"""
 import warnings
 
-from numpy import (
-    empty,
-    concatenate,
-    array,
-    identity,
-    sqrt,
-    vstack,
-    zeros_like,
-    allclose,
-    ix_,
-)
+from numpy import empty, concatenate, array, identity, sqrt, vstack, zeros_like, allclose, ix_
 from thewalrus.samples import hafnian_sample_state, torontonian_sample_state
 from thewalrus.symplectic import xxpp_to_xpxp
 
@@ -208,6 +198,14 @@ def is_vacuum(self, tol=0.0, **kwargs):
     def loss(self, T, mode):
         self.circuit.loss(T, mode)
 
+    def passive(self, T, *modes):
+        """
+        linear optical passive transformations
+        """
+        T_expand = identity(self.circuit.nlen, dtype=T.dtype)
+        T_expand[ix_(modes, modes)] = T
+        self.circuit.passive(T_expand)
+
     def thermal_loss(self, T, nbar, mode):
         self.circuit.thermal_loss(T, nbar, mode)
 

strawberryfields/backends/gaussianbackend/gaussiancircuit.py

@@ -540,6 +540,10 @@ def post_select_heterodyne(self, n, alpha_val):
 
     def apply_u(self, U):
         """Transforms the state according to the linear optical unitary that maps a[i] \to U[i, j]^*a[j]"""
-        self.mean = np.dot(np.conj(U), self.mean)
-        self.nmat = np.dot(np.dot(U, self.nmat), np.conj(np.transpose(U)))
-        self.mmat = np.dot(np.dot(np.conj(U), self.mmat), np.conj(np.transpose(U)))
+        self.mean = U @ self.mean
+        self.nmat = U.conj() @ self.nmat @ U.T
+        self.mmat = U @ self.mmat @ U.T
+
+    def passive(self, T):
+        """Transforms the state according to the arbitrary linear transformation that maps a[i] \to T[i, j]^*a[j]"""
+        self.apply_u(T)

strawberryfields/compilers/__init__.py

@@ -46,8 +46,20 @@
 from .gbs import GBS
 from .gaussian_unitary import GaussianUnitary
 from .gaussian_merge import GaussianMerge
+from .passive import Passive
 
-compilers = (Fock, Gaussian, Bosonic, GBS, GaussianUnitary, Xcov, Xstrict, Xunitary, GaussianMerge)
+compilers = (
+    Fock,
+    Gaussian,
+    Bosonic,
+    GBS,
+    GaussianUnitary,
+    Xcov,
+    Xstrict,
+    Xunitary,
+    GaussianMerge,
+    Passive,
+)
 
 compiler_db = {c.short_name: c for c in compilers}
 """dict[str, ~strawberryfields.compilers.Compiler]: Map from compiler name to the corresponding

strawberryfields/compilers/gaussian.py

@@ -41,6 +41,7 @@ class Gaussian(Compiler):
         # channels
         "LossChannel",
         "ThermalLossChannel",
+        "PassiveChannel",
         # single mode gates
         "Dgate",
         "Sgate",

strawberryfields/compilers/passive.py

@@ -0,0 +1,209 @@
+# 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.
+"""This module contains a compiler to reduce a sequence of passive gates into a single multimode linear passive operation"""
+
+import numpy as np
+from strawberryfields.program_utils import Command
+from strawberryfields import ops
+from strawberryfields.parameters import par_evaluate
+
+from .compiler import Compiler
+
+
+def _apply_one_mode_gate(G, T, i):
+    """In-place applies a one mode gate G into the process matrix T in mode i
+
+    Args:
+        G (complex/float): one mode gate
+        T (array): passive transformation
+        i (int): index of one mode gate
+    """
+
+    T[i] *= G
+    return T
+
+
+def _apply_two_mode_gate(G, T, i, j):
+    """In-place applies a two mode gate G into the process matrix T in modes i and j
+
+    Args:
+        G (array): 2x2 matrix for two mode gate
+        T (array): passive transformation
+        i (int): index of first mode of gate
+        j (int): index of second mode of gate
+    """
+    (T[i], T[j]) = (G[0, 0] * T[i] + G[0, 1] * T[j], G[1, 0] * T[i] + G[1, 1] * T[j])
+    return T
+
+
+def _beam_splitter_passive(theta, phi):
+    """Beam-splitter.
+
+    Args:
+        theta (float): transmissivity parameter
+        phi (float): phase parameter
+
+    Returns:
+        array: unitary 2x2 transformation matrix of an interferometer with angles theta and phi
+    """
+    ct = np.cos(theta)
+    st = np.sin(theta)
+    eip = np.cos(phi) + 1j * np.sin(phi)
+    U = np.array(
+        [
+            [ct, -np.conj(eip) * st],
+            [eip * st, ct],
+        ]
+    )
+    return U
+
+
+class Passive(Compiler):
+    """Compiler to write a sequence of passive operations as a single passive operation
+
+    This compiler checks whether the circuit can be implemented as a sequence of
+    passive operations. If so, it arranges them in a single matrix, T. It then returns an PassiveChannel
+    operation which can act this transformation.
+
+    This compiler can be accessed by calling :meth:`.Program.compile` with `'passive'` specified.
+
+    **Example:**
+
+    Consider the following Strawberry Fields program, compiled using the `'passive'` compiler:
+
+    .. code-block:: python3
+
+        from strawberryfields.ops import BSgate, LossChannel, Rgate
+        import strawberryfields as sf
+
+        circuit = sf.Program(2)
+        with circuit.context as q:
+            Rgate(np.pi) | q[0]
+            BSgate(0.25 * np.pi, 0) | (q[0], q[1])
+            LossChannel(0.9) | q[1]
+
+        compiled_circuit = circuit.compile(compiler="passive")
+
+    We can now print the compiled circuit, consisting a single
+    :class:`~.PassiveChannel`:
+
+    >>> compiled_circuit.print()
+    PassiveChannel([[-0.7071+8.6596e-17j -0.7071+0.0000e+00j]
+     [-0.6708+8.2152e-17j  0.6708+0.0000e+00j]]) | (q[0], q[1])
+    """
+
+    short_name = "passive"
+    interactive = True
+
+    primitives = {
+        # meta operations
+        "All",
+        "_New_modes",
+        "_Delete",
+        # single mode gates
+        "Rgate",
+        "LossChannel",
+        # multi mode gates
+        "MZgate",
+        "sMZgate",
+        "BSgate",
+        "Interferometer",  # Note that interferometer is accepted as a primitive
+        "PassiveChannel",  # and PassiveChannels!
+    }
+
+    decompositions = {}
+    # pylint: disable=too-many-branches, too-many-statements
+    def compile(self, seq, registers):
+        """Try to arrange a passive circuit into a single multimode passive operation
+
+        This method checks whether the circuit can be implemented as a sequence of passive gates.
+        If the answer is yes it arranges them into a single operation.
+
+        Args:
+            seq (Sequence[Command]): passive quantum circuit to modify
+            registers (Sequence[RegRefs]): quantum registers
+        Returns:
+            List[Command]: compiled circuit
+        Raises:
+            CircuitError: the circuit does not correspond to a passive unitary
+        """
+
+        # Check which modes are actually being used
+        used_modes = []
+        for operations in seq:
+            modes = [modes_label.ind for modes_label in operations.reg]
+            used_modes.append(modes)
+
+        used_modes = list(set(item for sublist in used_modes for item in sublist))
+
+        # dictionary mapping the used modes to consecutive non-negative integers
+        dict_indices = {used_modes[i]: i for i in range(len(used_modes))}
+        nmodes = len(used_modes)
+
+        # We start with an identity then sequentially update with the gate transformations
+        T = np.identity(nmodes, dtype=np.complex128)
+
+        # Now we will go through each operation in the sequence `seq` and apply it to T
+        for operations in seq:
+            name = operations.op.__class__.__name__
+            params = par_evaluate(operations.op.p)
+            modes = [modes_label.ind for modes_label in operations.reg]
+            if name == "Rgate":
+                G = np.exp(1j * params[0])
+                T = _apply_one_mode_gate(G, T, dict_indices[modes[0]])
+            elif name == "LossChannel":
+                G = np.sqrt(params[0])
+                T = _apply_one_mode_gate(G, T, dict_indices[modes[0]])
+            elif name == "Interferometer":
+                U = params[0]
+                if U.shape == (1, 1):
+                    T = _apply_one_mode_gate(U[0, 0], T, dict_indices[modes[0]])
+                elif U.shape == (2, 2):
+                    T = _apply_two_mode_gate(U, T, dict_indices[modes[0]], dict_indices[modes[1]])
+                else:
+                    modes = [dict_indices[mode] for mode in modes]
+                    U_expand = np.eye(nmodes, dtype=np.complex128)
+                    U_expand[np.ix_(modes, modes)] = U
+                    T = U_expand @ T
+            elif name == "PassiveChannel":
+                T0 = params[0]
+                if T0.shape == (1, 1):
+                    T = _apply_one_mode_gate(T0[0, 0], T, dict_indices[modes[0]])
+                elif T0.shape == (2, 2):
+                    T = _apply_two_mode_gate(T0, T, dict_indices[modes[0]], dict_indices[modes[1]])
+                else:
+                    modes = [dict_indices[mode] for mode in modes]
+                    T0_expand = np.eye(nmodes, dtype=np.complex128)
+                    T0_expand[np.ix_(modes, modes)] = T0
+                    T = T0_expand @ T
+            elif name == "BSgate":
+                G = _beam_splitter_passive(params[0], params[1])
+                T = _apply_two_mode_gate(G, T, dict_indices[modes[0]], dict_indices[modes[1]])
+            elif name == "MZgate":
+                v = np.exp(1j * params[0])
+                u = np.exp(1j * params[1])
+                U = 0.5 * np.array([[u * (v - 1), 1j * (1 + v)], [1j * u * (1 + v), 1 - v]])
+                T = _apply_two_mode_gate(U, T, dict_indices[modes[0]], dict_indices[modes[1]])
+            elif name == "sMZgate":
+                exp_sigma = np.exp(1j * (params[0] + params[1]) / 2)
+                delta = (params[0] - params[1]) / 2
+                U = exp_sigma * np.array(
+                    [[np.sin(delta), np.cos(delta)], [np.cos(delta), -np.sin(delta)]]
+                )
+                T = _apply_two_mode_gate(U, T, dict_indices[modes[0]], dict_indices[modes[1]])
+
+        ord_reg = [r for r in list(registers) if r.ind in used_modes]
+        ord_reg = sorted(list(ord_reg), key=lambda x: x.ind)
+
+        return [Command(ops.PassiveChannel(T), ord_reg)]

strawberryfields/ops.py

@@ -391,9 +391,11 @@ def merge(self, other):
         # channels can be merged if they are the same class and share all the other parameters
         if self.p[1:] == other.p[1:]:
             # determine the combined first parameter
-            T = self.p[0] * other.p[0]
+
+            T = np.dot(other.p[0], self.p[0])
             # if one, replace with the identity
-            if T == 1:
+            T_arr = np.atleast_2d(T)
+            if np.allclose(T_arr, np.eye(T_arr.shape[0])):
                 return None
 
             # return a copy
@@ -1411,6 +1413,29 @@ def _apply(self, reg, backend, **kwargs):
         return ancilla_val / s
 
 
+class PassiveChannel(Channel):
+    r"""Perform an arbitrary multimode passive operation
+
+    Args:
+        T (array): an NxN matrix acting on a N mode state
+
+    .. details::
+
+        Acts the following transformation on the state:
+
+        .. math::
+            a^{\dagger}_i \to \sum_j T_{ij} a^{\dagger}j
+    """
+
+    def __init__(self, T):
+        super().__init__([T])
+        self.ns = T.shape[0]
+
+    def _apply(self, reg, backend, **kwargs):
+        p = par_evaluate(self.p)
+        backend.passive(p[0], *reg)
+
+
 # ====================================================================
 # Unitary gates
 # ====================================================================