diff --git a/thewalrus/quantum.py b/thewalrus/quantum.py index 8be31c20e..545e57e53 100644 --- a/thewalrus/quantum.py +++ b/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: diff --git a/thewalrus/tests/test_integration.py b/thewalrus/tests/test_integration.py new file mode 100644 index 000000000..db575dc5c --- /dev/null +++ b/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) diff --git a/thewalrus/tests/test_quantum.py b/thewalrus/tests/test_quantum.py index b7d1036bf..1c9b6d2e7 100644 --- a/thewalrus/tests/test_quantum.py +++ b/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) ] )