Update decomposition function

strawberryfields/decompositions.py

@@ -612,7 +612,7 @@ def triangular(V, tol=1e-11):
             :math:`|VV^\dagger-I| \leq` tol
 
     Returns:
-        tuple[array]: returns a tuple of the form ``(tlist,np.diag(localV), None)``
+        tuple[array]: returns a tuple of the form ``(None,np.diag(localV), tlist)``
             where:
 
             * ``tlist``: list containing ``[n,m,theta,phi,n_size]`` of the T unitaries needed
@@ -630,7 +630,7 @@ def triangular(V, tol=1e-11):
             tlist.append(nullT(nsize - j - 1, nsize - i - 2, localV))
             localV = T(*tlist[-1]) @ localV
 
-    return list(reversed(tlist)), np.diag(localV), None
+    return None, np.diag(localV), tlist
 
 
 def M(n, sigma, delta, m):

strawberryfields/ops.py

@@ -2664,26 +2664,29 @@ def _decompose(self, reg, **kwargs):
             decomp_fn = getattr(dec, mesh)
             BS1, R, BS2 = decomp_fn(self.p[0], tol=tol)
 
-            for n, m, theta, phi, _ in BS1:
-                theta = theta if np.abs(theta) >= _decomposition_tol else 0
-                phi = phi if np.abs(phi) >= _decomposition_tol else 0
 
-                if "symmetric" in mesh:
-                    # Mach-Zehnder interferometers
-                    cmds.append(
-                        Command(
-                            MZgate(np.mod(theta, 2 * np.pi), np.mod(phi, 2 * np.pi)),
-                            (reg[n], reg[m]),
+            if BS1 is not None:
+
+                for n, m, theta, phi, _ in BS1:
+                    theta = theta if np.abs(theta) >= _decomposition_tol else 0
+                    phi = phi if np.abs(phi) >= _decomposition_tol else 0
+
+                    if "symmetric" in mesh:
+                        # Mach-Zehnder interferometers
+                        cmds.append(
+                            Command(
+                                MZgate(np.mod(theta, 2 * np.pi), np.mod(phi, 2 * np.pi)),
+                                (reg[n], reg[m]),
+                            )
                         )
-                    )
 
-                else:
-                    # Clements style beamsplitters
-                    if not (drop_identity and phi == 0):
-                        cmds.append(Command(Rgate(phi), reg[n]))
+                    else:
+                        # Clements style beamsplitters
+                        if not (drop_identity and phi == 0):
+                            cmds.append(Command(Rgate(phi), reg[n]))
 
-                    if not (drop_identity and theta == 0):
-                        cmds.append(Command(BSgate(theta, 0), (reg[n], reg[m])))
+                        if not (drop_identity and theta == 0):
+                            cmds.append(Command(BSgate(theta, 0), (reg[n], reg[m])))
 
             for n, expphi in enumerate(R):
                 # local phase shifts

tests/frontend/test_decompositions.py

@@ -351,11 +351,11 @@ def test_identity(self, tol):
         n = 20
         U = np.identity(n)
 
-        tlist, diags, _ = dec.triangular(U)
+        _, diags, tlist = dec.triangular(U)
 
         qrec = np.diag(diags)
 
-        for i in tlist:
+        for i in reversed(tlist):
             qrec = dec.Ti(*i) @ qrec
 
         assert np.allclose(U, qrec, atol=tol, rtol=0)
@@ -373,7 +373,7 @@ def test_random_unitary(self, tol):
         n = 20
         U = haar_measure(n)
 
-        tlist, diags, _ = dec.triangular(U)
+        _, diags, tlist = dec.triangular(U)
 
         qrec = np.diag(diags)
 

tests/frontend/test_interferometer_triangular.py

@@ -0,0 +1,52 @@
+import strawberryfields as sf
+from strawberryfields.ops import *
+import numpy as np
+
+"""
+Test to check that the triangular and rectangular meshes give the same results in the
+interferometer opreation.
+The tests creates a random M by M unitary matrix and applies it to a M-mode interferometer
+with the rectangular and triangular meshes. The test checks that the output probabilities.
+"""
+
+M = 5
+# Create a random complex matrix
+random_matrix = np.random.rand(M, M) + 1j * np.random.rand(M, M)
+# Perform QR decomposition to get a unitary matrix
+q, r = np.linalg.qr(random_matrix)
+# Normalize to ensure unitary property (q should already be unitary, but we'll double-check)
+unitary_matrix = q / np.linalg.norm(q, axis=0)
+
+boson_sampling_triangular = sf.Program(M)
+with boson_sampling_triangular.context as q:
+    # prepare the input fock states
+    for i in range(M):
+        if i % 2 == 0:
+            Fock(1) | q[i]
+        else:
+            Vac | q[i]
+    Interferometer(unitary_matrix, mesh = 'triangular', tol = 1e-10) | q
+
+boson_sampling_triangular.compile(compiler="fock")
+eng = sf.Engine(backend="fock", backend_options={"cutoff_dim": 5})
+results = eng.run(boson_sampling_triangular)
+# extract the joint Fock probabilities
+probs_triangular = results.state.all_fock_probs()
+
+#now try with mesh = 'rectangular'
+boson_sampling_rectangular = sf.Program(M)
+
+with boson_sampling_rectangular.context as q:
+    for i in range(M):
+        if i % 2 == 0:
+            Fock(1) | q[i]
+        else:
+            Vac | q[i]
+    Interferometer(unitary_matrix, mesh = 'rectangular', tol = 1e-10) | q
+
+boson_sampling_rectangular.compile(compiler="fock")
+eng = sf.Engine(backend="fock", backend_options={"cutoff_dim": 5})
+results = eng.run(boson_sampling_rectangular)
+# extract the joint Fock probabilities
+probs_rect = results.state.all_fock_probs()
+assert(np.allclose(probs_triangular, probs_rect))