Generaldyne sampling in gaussian (#614)

Fixes the formula used for sampling generaldyne outcomes in the gaussian backend.

The covariance matrix of the sampling distribution needed also to depend on the covariance of the state onto which the target state was begin measured.

Updated the test in test_heterodyne to account for this.

strawberryfields/backends/gaussianbackend/gaussiancircuit.py

@@ -478,7 +478,7 @@ def measure_dyne(self, covmat, indices, shots=1):
 
         r = self.smean()
         (va, vc) = ops.chop_in_blocks_vector(r, expind)
-        vm = np.random.multivariate_normal(vc, C, size=shots)
+        vm = np.random.multivariate_normal(vc, C + covmat, size=shots)
         # The next line is a hack in that it only updates conditioned on the first samples value
         # should still work if shots = 1
         va = va + np.dot(np.dot(B, np.linalg.inv(C + covmat)), vm[0] - vc)

tests/backend/test_heterodyne.py

@@ -20,10 +20,11 @@
 phase_alphas = np.linspace(0, 2 * np.pi, 7, endpoint=False)
 squeeze_val = np.arcsinh(1.0)
 
-n_meas = 300
+n_meas = 500
 disp_val = 1.0 + 1j * 1.0
 num_stds = 10.0
 std_10 = num_stds / np.sqrt(n_meas)
+R_VALS = [-1.2,0,1,1.4]
 
 
 @pytest.mark.backends("gaussian","bosonic")
@@ -65,18 +66,21 @@ def test_mean_coherent(self, setup_backend, pure, tol):
 
         assert np.allclose(x.mean(), disp_val, atol=std_10 + tol, rtol=0)
 
-    def test_std_vacuum(self, setup_backend, pure, tol):
-        """Test heterodyne provides the correct standard deviation of the vacuum"""
+    @pytest.mark.parametrize("r", R_VALS)
+    def test_std(self, r, setup_backend, pure, tol):
+        """Test heterodyne provides the correct standard deviation for heterodyne outcomes on 
+        squeezed states"""
         backend = setup_backend(1)
         x = np.empty(0)
 
         for i in range(n_meas):
             backend.reset(pure=pure)
+            backend.prepare_squeezed_state(r, 0, 0)
             meas_result = backend.measure_heterodyne(0)
             x = np.append(x, meas_result)
 
         xr = x.real
         xi = x.imag
-        xvar = xi.std() ** 2 + xr.std() ** 2
 
-        assert np.allclose(np.sqrt(xvar), np.sqrt(0.5), atol=std_10 + tol, rtol=0)
+        assert np.allclose(xr.std(), np.sqrt(1 + np.exp(-2 * r)) / 2, atol=std_10 + tol, rtol=0)
+        assert np.allclose(xi.std(), np.sqrt(1 + np.exp(2 * r)) / 2, atol=std_10 + tol, rtol=0)