Integration tests

thewalrus/quantum.py

@@ -331,7 +331,6 @@ def density_matrix_element(mu, cov, i, j, include_prefactor=True, tol=1e-10, hba
             haf = hafnian_repeated(A, rpt)
     else:
         # replace the diagonal of A with gamma
-        # gamma = X @ np.linalg.inv(Q).conj() @ beta
         gamma = beta.conj() - A @ beta
         if np.prod([k + 1 for k in rpt]) ** (1 / len(rpt)) < 3:
             A_rpt = reduction(A, rpt)
@@ -383,11 +382,12 @@ def density_matrix(mu, cov, post_select=None, normalize=False, cutoff=5, hbar=2)
     pref = prefactor(mu, cov, hbar=hbar)
 
     if post_select is None:
-        A = Amat(cov, hbar=hbar)
+        A = Amat(cov, hbar=hbar).conj()
         if np.allclose(mu, np.zeros_like(mu)):
             return pref * hermite_multidimensional(-A, cutoff, renorm=True)
         try:
-            y = np.linalg.inv(A) @ np.linalg.inv(Qmat(cov)) @ Beta(mu, hbar=hbar)
+            beta = Beta(mu)
+            y = np.linalg.inv(A) @ (beta - A @ beta.conj())
             return pref * hermite_multidimensional(-A, cutoff, y=-y, renorm=True)
         except np.linalg.LinAlgError:
             pass
@@ -454,7 +454,7 @@ def pure_state_amplitude(mu, cov, i, include_prefactor=True, tol=1e-10, hbar=2,
     A = Amat(cov, hbar=hbar)
     (n, _) = cov.shape
     N = n // 2
-    B = A[0:N, 0:N]
+    B = A[0:N, 0:N].conj()
     alpha = beta[0:N]
 
     if np.linalg.norm(alpha) < tol:
@@ -465,19 +465,16 @@ def pure_state_amplitude(mu, cov, i, include_prefactor=True, tol=1e-10, hbar=2,
         else:
             haf = hafnian_repeated(B, rpt)
     else:
-        # replace the diagonal of A with gamma
-        # gamma = X @ np.linalg.inv(Q).conj() @ beta
-        zeta = alpha - B @ (alpha.conj())
-
+        gamma = alpha - B @ np.conj(alpha)
         if np.prod([k + 1 for k in rpt]) ** (1 / len(rpt)) < 3:
             B_rpt = reduction(B, rpt)
-            np.fill_diagonal(B_rpt, reduction(zeta, rpt))
+            np.fill_diagonal(B_rpt, reduction(gamma, rpt))
             haf = hafnian(B_rpt, loop=True)
         else:
-            haf = hafnian_repeated(B, rpt, mu=zeta, loop=True)
+            haf = hafnian_repeated(B, rpt, mu=gamma, loop=True)
 
     if include_prefactor:
-        pref = np.exp(-0.5 * (np.linalg.norm(alpha) ** 2 - alpha.conj() @ B @ alpha.conj()))
+        pref = np.exp(-0.5 * (np.linalg.norm(alpha) ** 2 - alpha @ B @ alpha))
         haf *= pref
 
     return haf / np.sqrt(np.prod(fac(rpt)) * np.sqrt(np.linalg.det(Q)))
@@ -527,12 +524,13 @@ def state_vector(mu, cov, post_select=None, normalize=False, cutoff=5, hbar=2, c
 
     B = A[0:N, 0:N]
     alpha = beta[0:N]
-    pref = np.exp(-0.5 * (np.linalg.norm(alpha) ** 2 - alpha.conj() @ B @ alpha.conj()))
-
+    gamma = np.conj(alpha) - B @ alpha
+    pref = np.exp(-0.5 * (np.linalg.norm(alpha) ** 2 - alpha @ B @ alpha))
     if post_select is None:
+
         psi = (
             pref
-            * hafnian_batched(B, cutoff, mu=alpha, renorm=True)
+            * hafnian_batched(B.conj(), cutoff, mu=gamma.conj(), renorm=True)
             / np.sqrt(np.sqrt(np.linalg.det(Q).real))
         )
     else:

thewalrus/tests/test_integration.py

@@ -0,0 +1,51 @@
+# 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.
+"""Tests for The Walrus quantum functions"""
+# pylint: disable=no-self-use,redefined-outer-name
+
+import numpy as np
+from thewalrus.quantum import (
+    density_matrix,
+    state_vector,
+)
+
+
+def test_cubic_phase():
+    """Test that all the possible ways of obtaining a cubic phase state using the different methods agree"""
+    mu = np.array([-0.50047867, 0.37373598, 0.01421683, 0.26999427, 0.04450994, 0.01903583])
+
+    cov = np.array(
+        [
+            [1.57884241, 0.81035494, 1.03468307, 1.14908791, 0.09179507, -0.11893174],
+            [0.81035494, 1.06942863, 0.89359234, 0.20145142, 0.16202296, 0.4578259],
+            [1.03468307, 0.89359234, 1.87560498, 0.16915661, 1.0836528, -0.09405278],
+            [1.14908791, 0.20145142, 0.16915661, 2.37765137, -0.93543385, -0.6544286],
+            [0.09179507, 0.16202296, 1.0836528, -0.93543385, 2.78903152, -0.76519088],
+            [-0.11893174, 0.4578259, -0.09405278, -0.6544286, -0.76519088, 1.51724222],
+        ]
+    )
+
+    cutoff = 7
+    # the Fock state measurement of mode 0 to be post-selected
+    m1 = 1
+    # the Fock state measurement of mode 1 to be post-selected
+    m2 = 2
+
+    psi = state_vector(mu, cov, post_select={0: m1, 1: m2}, cutoff=cutoff, hbar=2)
+    psi_c = state_vector(mu, cov, cutoff=cutoff, hbar=2)[m1, m2, :]
+    rho = density_matrix(mu, cov, post_select={0: m1, 1: m2}, cutoff=cutoff, hbar=2)
+    rho_c = density_matrix(mu, cov, cutoff=cutoff, hbar=2)[m1, m2, :, m1, m2, :]
+    assert np.allclose(np.outer(psi, psi.conj()), rho)
+    assert np.allclose(np.outer(psi_c, psi_c.conj()), rho)
+    assert np.allclose(rho_c, rho)

thewalrus/tests/test_quantum.py

@@ -372,7 +372,6 @@ def test_coherent_squeezed():
          [0.07862323 + 0.00868528j, 0.01274241 + 0.00614023j, -0.02127257 - 0.00122123j, -0.00624626 - 0.00288134j, 0.00702606]]
     )
     # fmt:on
-
     assert np.allclose(res, expected)
 
 
@@ -529,7 +528,7 @@ def test_is_pure_cov_thermal(nbar):
 
 @pytest.mark.parametrize("i", [0, 1, 2, 3, 4])
 @pytest.mark.parametrize("j", [0, 1, 2, 3, 4])
-def test_pure_state_amplitude_two_mode_squezed(i, j):
+def test_pure_state_amplitude_two_mode_squeezed(i, j):
     """ Tests pure state amplitude for a two mode squeezed vacuum state """
     nbar = 1.0
     phase = np.pi / 8
@@ -539,8 +538,9 @@ def test_pure_state_amplitude_two_mode_squezed(i, j):
     if i != j:
         exact = 0.0
     else:
-        exact = np.exp(-1j * i * phase) * (nbar / (1.0 + nbar)) ** (i / 2) / np.sqrt(1.0 + nbar)
+        exact = np.exp(1j * i * phase) * (nbar / (1.0 + nbar)) ** (i / 2) / np.sqrt(1.0 + nbar)
     num = pure_state_amplitude(mu, cov, [i, j])
+
     assert np.allclose(exact, num)
 
 
@@ -587,7 +587,7 @@ def test_state_vector_two_mode_squeezed():
     mu = np.zeros([4], dtype=np.complex)
     exact = np.array(
         [
-            (np.exp(-1j * i * phase) * (nbar / (1.0 + nbar)) ** (i / 2) / np.sqrt(1.0 + nbar))
+            (np.exp(1j * i * phase) * (nbar / (1.0 + nbar)) ** (i / 2) / np.sqrt(1.0 + nbar))
             for i in range(cutoff)
         ]
     )
@@ -607,7 +607,7 @@ def test_state_vector_two_mode_squeezed_post():
     exact = np.diag(
         np.array(
             [
-                (np.exp(-1j * i * phase) * (nbar / (1.0 + nbar)) ** (i / 2) / np.sqrt(1.0 + nbar))
+                (np.exp(1j * i * phase) * (nbar / (1.0 + nbar)) ** (i / 2) / np.sqrt(1.0 + nbar))
                 for i in range(cutoff)
             ]
         )
@@ -647,7 +647,7 @@ def test_state_vector_two_mode_squeezed_post_normalize():
     exact = np.diag(
         np.array(
             [
-                (np.exp(-1j * i * phase) * (nbar / (1.0 + nbar)) ** (i / 2) / np.sqrt(1.0 + nbar))
+                (np.exp(1j * i * phase) * (nbar / (1.0 + nbar)) ** (i / 2) / np.sqrt(1.0 + nbar))
                 for i in range(cutoff)
             ]
         )