From bec434106b259af98be17a6f2dfdec1df2cae5a2 Mon Sep 17 00:00:00 2001 From: Theodor Isacsson Date: Fri, 4 Sep 2020 17:18:33 -0400 Subject: [PATCH 01/17] Refactor thewalrus.quantum --- thewalrus/quantum/__init__.py | 189 ++++++ thewalrus/quantum/adjacency_matrices.py | 185 ++++++ thewalrus/quantum/covariance_matrices.py | 288 ++++++++ thewalrus/quantum/fock_states_and_tensors.py | 623 ++++++++++++++++++ thewalrus/quantum/means_and_variances.py | 236 +++++++ .../quantum/photon_number_distributions.py | 94 +++ 6 files changed, 1615 insertions(+) create mode 100644 thewalrus/quantum/__init__.py create mode 100644 thewalrus/quantum/adjacency_matrices.py create mode 100644 thewalrus/quantum/covariance_matrices.py create mode 100644 thewalrus/quantum/fock_states_and_tensors.py create mode 100644 thewalrus/quantum/means_and_variances.py create mode 100644 thewalrus/quantum/photon_number_distributions.py diff --git a/thewalrus/quantum/__init__.py b/thewalrus/quantum/__init__.py new file mode 100644 index 000000000..915c7cb8e --- /dev/null +++ b/thewalrus/quantum/__init__.py @@ -0,0 +1,189 @@ +# Copyright 2019-2020 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. +""" +Quantum algorithms +================== + +.. currentmodule:: thewalrus.quantum + +This submodule provides access to various utility functions that act on Gaussian +quantum states. + +For more details on how the hafnian relates to various properties of Gaussian quantum +states, see: + +* Kruse, R., Hamilton, C. S., Sansoni, L., Barkhofen, S., Silberhorn, C., & Jex, I. + "Detailed study of Gaussian boson sampling." `Physical Review A 100, 032326 (2019) + `_ + +* Hamilton, C. S., Kruse, R., Sansoni, L., Barkhofen, S., Silberhorn, C., & Jex, I. + "Gaussian boson sampling." `Physical Review Letters, 119(17), 170501 (2017) + `_ + +* Quesada, N. + "Franck-Condon factors by counting perfect matchings of graphs with loops." + `Journal of Chemical Physics 150, 164113 (2019) `_ + +* Quesada, N., Helt, L. G., Izaac, J., Arrazola, J. M., Shahrokhshahi, R., Myers, C. R., & Sabapathy, K. K. + "Simulating realistic non-Gaussian state preparation." `Physical Review A 100, 022341 (2019) + `_ + + +Fock states and tensors +----------------------- + +.. autosummary:: + + pure_state_amplitude + state_vector + density_matrix_element + density_matrix + fock_tensor + probabilities + update_probabilities_with_loss + update_probabilities_with_noise + normal_ordered_expectation + fidelity + +Details +^^^^^^^ + +.. autofunction:: + pure_state_amplitude + +.. autofunction:: + state_vector + +.. autofunction:: + density_matrix_element + +.. autofunction:: + density_matrix + +.. autofunction:: + fock_tensor + +.. autofunction:: + probabilities + +.. autofunction:: + update_probabilities_with_loss + +.. autofunction:: + update_probabilities_with_noise + +.. autofunction:: + normal_ordered_expectation + +.. autofunction:: + fidelity + + +Utility functions +----------------- + +.. autosummary:: + + reduced_gaussian + Xmat + Qmat + Covmat + Amat + Beta + Means + prefactor + find_scaling_adjacency_matrix + mean_number_of_clicks + find_scaling_adjacency_matrix_torontonian + gen_Qmat_from_graph + photon_number_mean + photon_number_mean_vector + photon_number_covar + photon_number_covmat + is_valid_cov + is_pure_cov + is_classical_cov + total_photon_num_dist_pure_state + gen_single_mode_dist + gen_multi_mode_dist + normal_ordered_complex_cov + + +Details +^^^^^^^ +""" + +from .fock_states_and_tensors import ( + pure_state_amplitude, + state_vector, + density_matrix_element, + density_matrix, + fock_tensor, + probabilities, + update_probabilities_with_loss, + update_probabilities_with_noise, + fidelity, + prefactor, +) + +from .adjacency_matrices import ( + mean_number_of_clicks, + find_scaling_adjacency_matrix, + find_scaling_adjacency_matrix_torontonian, + gen_Qmat_from_graph, +) + +from .covariance_matrices import ( + Xmat, + Qmat, + Covmat, + Amat, + normal_ordered_expectation, + normal_ordered_complex_cov, + Beta, + Means, + + is_valid_cov, + is_pure_cov, + is_classical_cov +) + +from .means_and_variances import ( + reduced_gaussian, + photon_number_mean, + photon_number_mean_vector, + photon_number_covar, + photon_number_covmat, + photon_number_expectation, + photon_number_squared_expectation, +) + +from .photon_number_distributions import ( + total_photon_num_dist_pure_state, + gen_single_mode_dist, + gen_multi_mode_dist, +) + +__all__ = [ + "pure_state_amplitude", + "state_vector", + "density_matrix_element", + "density_matrix", + "fock_tensor", + "probabilities", + "update_probabilities_with_loss", + "update_probabilities_with_noise", + "normal_ordered_expectation", + "fidelity", +] \ No newline at end of file diff --git a/thewalrus/quantum/adjacency_matrices.py b/thewalrus/quantum/adjacency_matrices.py new file mode 100644 index 000000000..f661e76da --- /dev/null +++ b/thewalrus/quantum/adjacency_matrices.py @@ -0,0 +1,185 @@ +# Copyright 2019-2020 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. +""" +Functions for rescaling adjacency matrices as well as for calculating the mean +number of clicks for an adjacency matrix. +""" + +import numpy as np +from scipy.optimize import root_scalar + +from .covariance_matrices import Xmat + +def mean_number_of_clicks(A): + r""" Given an adjacency matrix this function calculates the mean number of clicks. + For this to make sense the user must provide a matrix with singular values + less than or equal to one. See Appendix A.3 of `_ + by Banchi et al. + + Args: + A (array): rescaled adjacency matrix + + Returns: + float: mean number of clicks + """ + n, _ = A.shape + idn = np.identity(n) + X = np.block([[0 * idn, idn], [idn, 0 * idn]]) + B = np.block([[A, 0 * A], [0 * A, np.conj(A)]]) + Q = np.linalg.inv(np.identity(2 * n) - X @ B) + meanc = 1.0 * n + + for i in range(n): + det_val = np.real(Q[i, i]*Q[i+n, i+n] - Q[i+n, i]*Q[i, i+n]) + meanc -= 1.0 / np.sqrt(det_val) + return meanc + + +################################################################################ +# Rescale adjacency matrices +################################################################################ + + +def find_scaling_adjacency_matrix(A, n_mean): + r""" Returns the scaling parameter by which the adjacency matrix A + should be rescaled so that the Gaussian state that endodes it has + a total mean photon number n_mean. + + Args: + A (array): Adjacency matrix + n_mean (float): Mean photon number of the Gaussian state + + Returns: + float: Scaling parameter + """ + eps = 1e-10 + ls = np.linalg.svd(A, compute_uv=False) + max_sv = ls[0] + a_lim = 0.0 + b_lim = 1.0 / (eps + max_sv) + x_init = 0.5 * b_lim + + if 1000 * eps >= max_sv: + raise ValueError("The singular values of the matrix A are too small.") + + def mean_photon_number(x, vals): + r""" Returns the mean number of photons in the Gaussian state that + encodes the adjacency matrix x*A where vals are the singular values of A + + Args: + x (float): Scaling parameter + vals (array): Singular values of the matrix A + + Returns: + n_mean: Mean photon number in the Gaussian state + """ + vals2 = (x * vals) ** 2 + n = np.sum(vals2 / (1.0 - vals2)) + return n + + # The following function is implicitly tested in test_find_scaling_adjacency_matrix + def grad_mean_photon_number(x, vals): # pragma: no cover + r""" Returns the gradient od the mean number of photons in the Gaussian state that + encodes the adjacency matrix x*A with respect to x. + vals are the singular values of A + + Args: + x (float): Scaling parameter + vals (array): Singular values of the matrix A + + Returns: + d_n_mean: Derivative of the mean photon number in the Gaussian state + with respect to x + """ + vals1 = vals * x + dn = (2.0 / x) * np.sum((vals1 / (1 - vals1 ** 2)) ** 2) + return dn + + f = lambda x: mean_photon_number(x, ls) - n_mean + df = lambda x: grad_mean_photon_number(x, ls) + res = root_scalar(f, fprime=df, x0=x_init, bracket=(a_lim, b_lim)) + + if not res.converged: + raise ValueError("The search for a scaling value failed") + + return res.root + + +def find_scaling_adjacency_matrix_torontonian(A, c_mean): + r""" Returns the scaling parameter by which the adjacency matrix A + should be rescaled so that the Gaussian state that encodes it has + give a mean number of clicks equal to ``c_mean`` when measured with + threshold detectors. + + Args: + A (array): adjacency matrix + c_mean (float): mean photon number of the Gaussian state + + Returns: + float: scaling parameter + """ + n, _ = A.shape + if c_mean < 0 or c_mean > n: + raise ValueError("The mean number of clicks should be smaller than the number of modes") + + vals = np.linalg.svd(A, compute_uv=False) + localA = A / vals[0] # rescale the matrix so that the singular values are between 0 and 1. + + def cost(x): + r""" Cost function giving the difference between the wanted number of clicks and the number + of clicks at a given scaling value. It assumes that the adjacency matrix has been rescaled + so that it has singular values between 0 and 1. + + Args: + x (float): scaling value + + Return: + float: difference between desired and obtained mean number of clicks + """ + if x >= 1.0: + return c_mean - n + if x <= 0: + return c_mean + return c_mean - mean_number_of_clicks(x * localA) + + res = root_scalar(cost, x0=0.5, bracket=(0.0, 1.0)) # Do the optimization + + if not res.converged: + raise ValueError("The search for a scaling value failed") + return res.root / vals[0] + + +def gen_Qmat_from_graph(A, n_mean): + r""" Returns the Qmat xp-covariance matrix associated to a graph with + adjacency matrix :math:`A` and with mean photon number :math:`n_{mean}`. + + Args: + A (array): a :math:`N\times N` ``np.float64`` (symmetric) adjacency matrix + n_mean (float): mean photon number of the Gaussian state + + Returns: + array: the :math:`2N\times 2N` Q matrix. + """ + n, m = A.shape + + if n != m: + raise ValueError("Matrix must be square.") + + sc = find_scaling_adjacency_matrix(A, n_mean) + Asc = sc * A + A = np.block([[Asc, 0 * Asc], [0 * Asc, Asc.conj()]]) + I = np.identity(2 * n) + X = Xmat(n) + Q = np.linalg.inv(I - X @ A) + return Q diff --git a/thewalrus/quantum/covariance_matrices.py b/thewalrus/quantum/covariance_matrices.py new file mode 100644 index 000000000..d5426ef00 --- /dev/null +++ b/thewalrus/quantum/covariance_matrices.py @@ -0,0 +1,288 @@ +# Copyright 2019-2020 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. +""" +Functions for transforming one type of covariance-matrix-like object into +another as well as various property tests for covariance matrices. +""" + +import numpy as np + +from ..symplectic import sympmat +from .._hafnian import hafnian, reduction + +################################################################################ +# Transform one type of covariance-matrix-like object into another +################################################################################ + + +def Xmat(N): + r"""Returns the matrix :math:`X_n = \begin{bmatrix}0 & I_n\\ I_n & 0\end{bmatrix}` + + Args: + N (int): positive integer + + Returns: + array: :math:`2N\times 2N` array + """ + I = np.identity(N) + O = np.zeros_like(I) + X = np.block([[O, I], [I, O]]) + return X + + +def Qmat(cov, hbar=2): + r"""Returns the :math:`Q` Husimi matrix of the Gaussian state. + + Args: + cov (array): :math:`2N\times 2N xp-` Wigner covariance matrix + hbar (float): the value of :math:`\hbar` in the commutation + relation :math:`[\x,\p]=i\hbar`. + + Returns: + array: the :math:`Q` matrix. + """ + # number of modes + N = len(cov) // 2 + I = np.identity(N) + + x = cov[:N, :N] * 2 / hbar + xp = cov[:N, N:] * 2 / hbar + p = cov[N:, N:] * 2 / hbar + # the (Hermitian) matrix elements + aidaj = (x + p + 1j * (xp - xp.T) - 2 * I) / 4 + # the (symmetric) matrix elements + aiaj = (x - p + 1j * (xp + xp.T)) / 4 + + # calculate the covariance matrix sigma_Q appearing in the Q function: + # Q(alpha) = exp[-(alpha-beta).sigma_Q^{-1}.(alpha-beta)/2]/|sigma_Q| + Q = np.block([[aidaj, aiaj.conj()], [aiaj, aidaj.conj()]]) + np.identity(2 * N) + return Q + + +def Covmat(Q, hbar=2): + r"""Returns the Wigner covariance matrix in the :math:`xp`-ordering of the Gaussian state. + This is the inverse function of Qmat. + + Args: + Q (array): :math:`2N\times 2N` Husimi Q matrix + hbar (float): the value of :math:`\hbar` in the commutation + relation :math:`[\x,\p]=i\hbar`. + + Returns: + array: the :math:`xp`-ordered covariance matrix in the xp-ordering. + """ + # number of modes + n = len(Q) // 2 + I = np.identity(n) + N = Q[0:n, 0:n] - I + M = Q[n : 2 * n, 0:n] + mm11a = 2 * (N.real + M.real) + np.identity(n) + mm22a = 2 * (N.real - M.real) + np.identity(n) + mm12a = 2 * (M.imag + N.imag) + cov = np.block([[mm11a, mm12a], [mm12a.T, mm22a]]) + + return (hbar / 2) * cov + + +def Amat(cov, hbar=2, cov_is_qmat=False): + r"""Returns the :math:`A` matrix of the Gaussian state whose hafnian gives the photon number probabilities. + + Args: + cov (array): :math:`2N\times 2N` covariance matrix + hbar (float): the value of :math:`\hbar` in the commutation + relation :math:`[\x,\p]=i\hbar`. + cov_is_qmat (bool): if ``True``, it is assumed that ``cov`` is in fact the Q matrix. + + Returns: + array: the :math:`A` matrix. + """ + # number of modes + N = len(cov) // 2 + X = Xmat(N) + + # inverse Q matrix + if cov_is_qmat: + Q = cov + else: + Q = Qmat(cov, hbar=hbar) + + Qinv = np.linalg.inv(Q) + + # calculate Hamilton's A matrix: A = X.(I-Q^{-1})* + A = X @ (np.identity(2 * N) - Qinv).conj() + return A + + +def normal_ordered_expectation(mu, cov, rpt, hbar=2): + r"""Calculates the expectation value of the normal ordered product + :math:`\prod_{i=0}^{N-1} a_i^{\dagger n_i} \prod_{j=0}^{N-1} a_j^{m_j}` with respect to an N-mode Gaussian state, + where :math:`\text{rpt}=(n_0, n_1, \ldots, n_{N-1}, m_0, m_1, \ldots, m_{N-1})`. + + Args: + mu (array): length-:math:`2N` means vector in xp-ordering. + cov (array): :math:`2N\times 2N` covariance matrix in xp-ordering. + rpt (list): integers specifying the terms to calculate. + hbar (float): value of hbar in the uncertainty relation. + + Returns: + (float): expectation value of the normal ordered product of operators + """ + alpha = Beta(mu, hbar=hbar) + V = normal_ordered_complex_cov(cov, hbar=hbar) + A = reduction(V, rpt) + if np.allclose(mu, 0): + res = np.conj(hafnian(A)) + else: + np.fill_diagonal(A, reduction(np.conj(alpha), rpt)) + res = np.conj(hafnian(A, loop=True)) + return np.conj(res) + + +def normal_ordered_complex_cov(cov, hbar=2): + r"""Calculates the normal ordered covariance matrix in the complex basis. + + Args: + cov (array): xp-covariance matrix. + hbar (float): value of hbar in the uncertainty relation. + + Returns: + (array): covariance matrix in the creation/annihilation operator basis. + """ + + n, _ = cov.shape + n_modes = n // 2 + cov = cov / (hbar / 2) + A = cov[:n_modes, :n_modes] + B = cov[:n_modes, n_modes:] + C = cov[n_modes:, n_modes:] + N = 0.25 * (A + C + 1j * (B - B.T) - 2 * np.identity(n_modes)) + M = 0.25 * (A - C + 1j * (B + B.T)) + mat = np.block([[M.conj(), N], [N.T, M]]) + return mat + + +def Beta(mu, hbar=2): + r"""Returns the vector of complex displacements and conjugate displacements. + + Args: + mu (array): length-:math:`2N` means vector + hbar (float): the value of :math:`\hbar` in the commutation + relation :math:`[\x,\p]=i\hbar`. + + Returns: + array: the expectation values + :math:`[\langle a_1\rangle, \langle a_2\rangle,\dots,\langle a_N\rangle, \langle a^\dagger_1\rangle, \dots, \langle a^\dagger_N\rangle]` + """ + N = len(mu) // 2 + # mean displacement of each mode + alpha = (mu[:N] + 1j * mu[N:]) / np.sqrt(2 * hbar) + # the expectation values (, ,...,, , ..., ) + return np.concatenate([alpha, alpha.conj()]) + + +def Means(beta, hbar=2): + r"""Returns the vector of real quadrature displacements. + + Args: + beta (array): length-:math:`2N` means bivector + hbar (float): the value of :math:`\hbar` in the commutation + relation :math:`[\x,\p]=i\hbar`. + + Returns: + array: the quadrature expectation values + :math:`[\langle q_1\rangle, \langle q_2\rangle,\dots,\langle q_N\rangle, \langle p_1\rangle, \dots, \langle p_N\rangle]` + """ + + N = len(beta) // 2 + alpha = beta[0:N] + return np.sqrt(2 * hbar) * np.concatenate([alpha.real, alpha.imag]) + +################################################################################ +# Test properties of covariance matrices +################################################################################ + + +def is_valid_cov(cov, hbar=2, rtol=1e-05, atol=1e-08): + r""" Checks if the covariance matrix is a valid quantum covariance matrix. + + Args: + cov (array): a covariance matrix + hbar (float): value of hbar in the uncertainty relation + rtol (float): the relative tolerance parameter used in `np.allclose` + atol (float): the absolute tolerance parameter used in `np.allclose` + + Returns: + (bool): whether the given covariance matrix is a valid covariance matrix + """ + (n, m) = cov.shape + if n != m: + # raise ValueError("The input matrix must be square") + return False + if not np.allclose(cov, np.transpose(cov), rtol=rtol, atol=atol): + # raise ValueError("The input matrix is not symmetric") + return False + if n % 2 != 0: + # raise ValueError("The input matrix is of even dimension") + return False + + nmodes = n // 2 + vals = np.linalg.eigvalsh(cov + 0.5j * hbar * sympmat(nmodes)) + vals[np.abs(vals) < atol] = 0.0 + if np.all(vals >= 0): + # raise ValueError("The input matrix violates the uncertainty relation") + return True + + return False + + +def is_pure_cov(cov, hbar=2, rtol=1e-05, atol=1e-08): + r""" Checks if the covariance matrix is a valid quantum covariance matrix + that corresponds to a quantum pure state + + Args: + cov (array): a covariance matrix + hbar (float): value of hbar in the uncertainty relation + rtol (float): the relative tolerance parameter used in `np.allclose` + atol (float): the absolute tolerance parameter used in `np.allclose` + + Returns: + (bool): whether the given covariance matrix corresponds to a pure state + """ + if is_valid_cov(cov, hbar=hbar, rtol=rtol, atol=atol): + purity = 1 / np.sqrt(np.linalg.det(2 * cov / hbar)) + if np.allclose(purity, 1.0, rtol=rtol, atol=atol): + return True + + return False + + +def is_classical_cov(cov, hbar=2, atol=1e-08): + r""" Checks if the covariance matrix can be efficiently sampled. + + Args: + cov (array): a covariance matrix + hbar (float): value of hbar in the uncertainty relation + + Returns: + (bool): whether the given covariance matrix corresponds to a classical state + """ + + if is_valid_cov(cov, hbar=hbar, atol=atol): + (n, _) = cov.shape + vals = np.linalg.eigvalsh(cov - 0.5 * hbar * np.identity(n)) + vals[np.abs(vals) < atol] = 0.0 + + if np.all(vals >= 0): + return True + return False diff --git a/thewalrus/quantum/fock_states_and_tensors.py b/thewalrus/quantum/fock_states_and_tensors.py new file mode 100644 index 000000000..904d88a3a --- /dev/null +++ b/thewalrus/quantum/fock_states_and_tensors.py @@ -0,0 +1,623 @@ +# Copyright 2019-2020 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. +""" +Set of functions for calculating various state representations, probabilites and +fidelities of Gaussian states. +""" + +from itertools import count, product, chain + +import numpy as np +import dask + +from scipy.special import factorial as fac +from scipy.linalg import sqrtm +from numba import jit + +from ..symplectic import expand, sympmat, is_symplectic +from ..libwalrus import interferometer, interferometer_real + +from .._hafnian import hafnian, hafnian_repeated, reduction +from .._hermite_multidimensional import hermite_multidimensional, hafnian_batched + +from .covariance_matrices import ( + is_pure_cov, + Amat, + Beta, + Qmat, +) + + +# pylint: disable=too-many-arguments +def pure_state_amplitude(mu, cov, i, include_prefactor=True, tol=1e-10, hbar=2, check_purity=True): + r"""Returns the :math:`\langle i | \psi\rangle` element of the state ket + of a Gaussian state defined by covariance matrix cov. + + + Args: + mu (array): length-:math:`2N` quadrature displacement vector + cov (array): length-:math:`2N` covariance matrix + i (list): list of amplitude elements + include_prefactor (bool): if ``True``, the prefactor is automatically calculated + used to scale the result. + tol (float): tolerance for determining if displacement is negligible + hbar (float): the value of :math:`\hbar` in the commutation + relation :math:`[\x,\p]=i\hbar`. + check_purity (bool): if ``True``, the purity of the Gaussian state is checked + before calculating the state vector. + + Returns: + complex: the pure state amplitude + """ + if check_purity: + if not is_pure_cov(cov, hbar=hbar, rtol=1e-05, atol=1e-08): + raise ValueError("The covariance matrix does not correspond to a pure state") + + rpt = i + beta = Beta(mu, hbar=hbar) + Q = Qmat(cov, hbar=hbar) + A = Amat(cov, hbar=hbar) + (n, _) = cov.shape + N = n // 2 + B = A[0:N, 0:N].conj() + alpha = beta[0:N] + + if np.linalg.norm(alpha) < tol: + # no displacement + if np.prod([k + 1 for k in rpt]) ** (1 / len(rpt)) < 3: + B_rpt = reduction(B, rpt) + haf = hafnian(B_rpt) + else: + haf = hafnian_repeated(B, rpt) + else: + 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(gamma, rpt)) + haf = hafnian(B_rpt, loop=True) + else: + haf = hafnian_repeated(B, rpt, mu=gamma, loop=True) + + if include_prefactor: + 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))) + + +def state_vector( + mu, cov, post_select=None, normalize=False, cutoff=5, hbar=2, check_purity=True, **kwargs +): + r"""Returns the state vector of a (PNR post-selected) Gaussian state. + + The resulting density matrix will have shape + + .. math:: \underbrace{D\times D \times \cdots \times D}_M + + where :math:`D` is the Fock space cutoff, and :math:`M` is the + number of *non* post-selected modes, i.e. ``M = len(mu)//2 - len(post_select)``. + + If post_select is None then the density matrix elements are calculated using + the multidimensional Hermite polynomials which provide a significantly faster + evaluation. + + + Args: + mu (array): length-:math:`2N` means vector in xp-ordering + cov (array): :math:`2N\times 2N` covariance matrix in xp-ordering + post_select (dict): dictionary containing the post-selected modes, of + the form ``{mode: value}``. + normalize (bool): If ``True``, a post-selected density matrix is re-normalized. + cutoff (dim): the final length (i.e., Hilbert space dimension) of each + mode in the density matrix. + hbar (float): the value of :math:`\hbar` in the commutation + relation :math:`[\x,\p]=i\hbar`. + check_purity (bool): if ``True``, the purity of the Gaussian state is checked + before calculating the state vector. + + Keyword Args: + choi_r (float or None): Value of the two-mode squeezing parameter used in Choi-Jamiolkoski + trick in :func:`~.fock_tensor`. This keyword argument should only be used when ``state_vector`` + is called by :func:`~.fock_tensor`. + + Returns: + np.array[complex]: the state vector of the Gaussian state + """ + if check_purity: + if not is_pure_cov(cov, hbar=hbar, rtol=1e-05, atol=1e-08): + raise ValueError("The covariance matrix does not correspond to a pure state") + + beta = Beta(mu, hbar=hbar) + A = Amat(cov, hbar=hbar) + Q = Qmat(cov, hbar=hbar) + + (n, _) = cov.shape + N = n // 2 + + B = A[0:N, 0:N] + alpha = beta[0:N] + gamma = np.conj(alpha) - B @ alpha + prefexp = -0.5 * (np.linalg.norm(alpha) ** 2 - alpha @ B @ alpha) + pref = np.exp(prefexp.conj()) + if post_select is None: + choi_r = kwargs.get("choi_r", None) + if choi_r is None: + denom = np.sqrt(np.sqrt(np.linalg.det(Q).real)) + else: + rescaling = np.concatenate( + [np.ones([N // 2]), (1.0 / np.tanh(choi_r)) * np.ones([N // 2])] + ) + B = np.diag(rescaling) @ B @ np.diag(rescaling) + gamma = rescaling * gamma + denom = np.sqrt(np.sqrt(np.linalg.det(Q / np.cosh(choi_r)).real)) + + psi = ( + np.real_if_close(pref) + * hafnian_batched(B.conj(), cutoff, mu=gamma.conj(), renorm=True) + / denom + ) + else: + M = N - len(post_select) + psi = np.zeros([cutoff] * (M), dtype=np.complex128) + + for idx in product(range(cutoff), repeat=M): + el = [] + + counter = count(0) + modes = (np.arange(N)).tolist() + el = [post_select[i] if i in post_select else idx[next(counter)] for i in modes] + psi[idx] = pure_state_amplitude( + mu, cov, el, check_purity=False, include_prefactor=False, hbar=hbar + ) + + psi = psi * pref + + if normalize: + norm = np.sqrt(np.sum(np.abs(psi) ** 2)) + psi = psi / norm + + return psi + + +def density_matrix_element(mu, cov, i, j, include_prefactor=True, tol=1e-10, hbar=2): + r"""Returns the :math:`\langle i | \rho | j \rangle` element of the density matrix + of a Gaussian state defined by covariance matrix cov. + + Args: + mu (array): length-:math:`2N` quadrature displacement vector + cov (array): length-:math:`2N` covariance matrix + i (list): list of density matrix rows + j (list): list of density matrix columns + include_prefactor (bool): if ``True``, the prefactor is automatically calculated + used to scale the result. + tol (float): tolerance for determining if displacement is negligible + hbar (float): the value of :math:`\hbar` in the commutation + relation :math:`[\x,\p]=i\hbar`. + + Returns: + complex: the density matrix element + """ + rpt = i + j + beta = Beta(mu, hbar=hbar) + A = Amat(cov, hbar=hbar) + if np.linalg.norm(beta) < tol: + # no displacement + if np.prod([k + 1 for k in rpt]) ** (1 / len(rpt)) < 3: + A_rpt = reduction(A, rpt) + haf = hafnian(A_rpt) + else: + haf = hafnian_repeated(A, rpt) + else: + # replace the diagonal of A with gamma + gamma = beta.conj() - A @ beta + if np.prod([k + 1 for k in rpt]) ** (1 / len(rpt)) < 3: + A_rpt = reduction(A, rpt) + np.fill_diagonal(A_rpt, reduction(gamma, rpt)) + haf = hafnian(A_rpt, loop=True) + else: + haf = hafnian_repeated(A, rpt, mu=gamma, loop=True) + + if include_prefactor: + haf *= prefactor(mu, cov, hbar=hbar) + + return haf / np.sqrt(np.prod(fac(rpt))) + + +def density_matrix(mu, cov, post_select=None, normalize=False, cutoff=5, hbar=2): + r"""Returns the density matrix of a (PNR post-selected) Gaussian state. + + The resulting density matrix will have shape + + .. math:: \underbrace{D\times D \times \cdots \times D}_{2M} + + where :math:`D` is the Fock space cutoff, and :math:`M` is the + number of *non* post-selected modes, i.e. ``M = len(mu)//2 - len(post_select)``. + + Note that we use the Strawberry Fields convention for indexing the density + matrix; the first two dimensions correspond to subsystem 1, the second two + dimensions correspond to subsystem 2, etc. + If post_select is None then the density matrix elements are calculated using + the multidimensional Hermite polynomials which provide a significantly faster + evaluation. + + Args: + mu (array): length-:math:`2N` means vector in xp-ordering + cov (array): :math:`2N\times 2N` covariance matrix in xp-ordering + post_select (dict): dictionary containing the post-selected modes, of + the form ``{mode: value}``. If post_select is None the whole non post-selected density matrix + is calculated directly using (multidimensional) Hermite polynomials, which is significantly faster + than calculating one hafnian at a time. + normalize (bool): If ``True``, a post-selected density matrix is re-normalized. + cutoff (dim): the final length (i.e., Hilbert space dimension) of each + mode in the density matrix. + hbar (float): the value of :math:`\hbar` in the commutation + relation :math:`[\x,\p]=i\hbar`. + + Returns: + np.array[complex]: the density matrix of the Gaussian state + """ + N = len(mu) // 2 + pref = prefactor(mu, cov, hbar=hbar) + + if post_select is None: + A = Amat(cov, hbar=hbar).conj() + sf_order = tuple(chain.from_iterable([[i, i + N] for i in range(N)])) + + if np.allclose(mu, np.zeros_like(mu)): + tensor = np.real_if_close(pref) * hermite_multidimensional( + -A, cutoff, renorm=True, modified=True + ) + return tensor.transpose(sf_order) + beta = Beta(mu, hbar=hbar) + y = beta - A @ beta.conj() + tensor = np.real_if_close(pref) * hermite_multidimensional( + -A, cutoff, y=y, renorm=True, modified=True + ) + return tensor.transpose(sf_order) + + M = N - len(post_select) + rho = np.zeros([cutoff] * (2 * M), dtype=np.complex128) + + for idx in product(range(cutoff), repeat=2 * M): + el = [] + + counter = count(0) + modes = (np.arange(2 * N) % N).tolist() + el = [post_select[i] if i in post_select else idx[next(counter)] for i in modes] + + el = np.array(el).reshape(2, -1) + el0 = el[0].tolist() + el1 = el[1].tolist() + + sf_idx = np.array(idx).reshape(2, -1) + sf_el = tuple(sf_idx[::-1].T.flatten()) + + rho[sf_el] = density_matrix_element(mu, cov, el0, el1, include_prefactor=False, hbar=hbar) + + rho *= pref + + if normalize: + # construct the standard 2D density matrix, and take the trace + new_ax = np.arange(2 * M).reshape([M, 2]).T.flatten() + tr = np.trace(rho.transpose(new_ax).reshape([cutoff ** M, cutoff ** M])).real + # renormalize + rho /= tr + + return rho + + +def fock_tensor( + S, + alpha, + cutoff, + choi_r=np.arcsinh(1.0), + check_symplectic=True, + sf_order=False, + rtol=1e-05, + atol=1e-08, +): + r""" + Calculates the Fock representation of a Gaussian unitary parametrized by + the symplectic matrix S and the displacements alpha up to cutoff in Fock space. + + Args: + S (array): symplectic matrix + alpha (array): complex vector of displacements + cutoff (int): cutoff in Fock space + choi_r (float): squeezing parameter used for the Choi expansion + check_symplectic (boolean): checks whether the input matrix is symplectic + sf_order (boolean): reshapes the tensor so that it follows the sf ordering of indices + rtol (float): the relative tolerance parameter used in `np.allclose` + atol (float): the absolute tolerance parameter used in `np.allclose` + + Return: + (array): Tensor containing the Fock representation of the Gaussian unitary + """ + # Check the matrix is symplectic + if check_symplectic: + if not is_symplectic(S, rtol=rtol, atol=atol): + raise ValueError("The matrix S is not symplectic") + + # And that S and alpha have compatible dimensions + m, _ = S.shape + l = m // 2 + if l != len(alpha): + raise ValueError( + "The matrix S and the vector alpha do not have compatible dimensions" + ) + # Check if S corresponds to an interferometer, if so use optimized routines + if np.allclose(S @ S.T, np.identity(m), rtol=rtol, atol=atol) and np.allclose( + alpha, 0, rtol=rtol, atol=atol + ): + reU = S[:l, :l] + imU = S[:l, l:] + if np.allclose(imU, 0, rtol=rtol, atol=atol): + Ub = np.block([[0 * reU, -reU], [-reU.T, 0 * reU]]) + tensor = interferometer_real(Ub, cutoff) + else: + U = reU - 1j * imU + Ub = np.block([[0 * U, -U], [-U.T, 0 * U]]) + tensor = interferometer(Ub, cutoff) + else: + # Construct the covariance matrix of l two-mode squeezed vacua pairing modes i and i+l + ch = np.cosh(choi_r) * np.identity(l) + sh = np.sinh(choi_r) * np.identity(l) + zh = np.zeros([l, l]) + Schoi = np.block( + [[ch, sh, zh, zh], [sh, ch, zh, zh], [zh, zh, ch, -sh], [zh, zh, -sh, ch]] + ) + # And then its Choi expanded symplectic + S_exp = expand(S, list(range(l)), 2 * l) @ Schoi + # And this is the corresponding covariance matrix + cov = S_exp @ S_exp.T + alphat = np.array(list(alpha) + ([0] * l)) + x = 2 * alphat.real + p = 2 * alphat.imag + mu = np.concatenate([x, p]) + + tensor = state_vector( + mu, + cov, + normalize=False, + cutoff=cutoff, + hbar=2, + check_purity=False, + choi_r=choi_r, + ) + + if sf_order: + sf_indexing = tuple(chain.from_iterable([[i, i + l] for i in range(l)])) + return tensor.transpose(sf_indexing) + + return tensor + + +def probabilities(mu, cov, cutoff, parallel=False, hbar=2.0, rtol=1e-05, atol=1e-08): + r"""Generate the Fock space probabilities of a Gaussian state up to a Fock space cutoff. + + .. note:: + + Individual density matrix elements are computed using multithreading by OpenMP. + Setting ``parallel=True`` will further result in *multiple* density matrix elements + being computed in parallel. + + When setting ``parallel=True``, OpenMP will need to be turned off by setting the + environment variable ``OMP_NUM_THREADS=1`` (forcing single threaded use for individual + matrix elements). Remove the environment variable or set it to ``OMP_NUM_THREADS=''`` + to again use multithreading with OpenMP. + + Args: + mu (array): vector of means of length ``2*n_modes`` + cov (array): covariance matrix of shape ``[2*n_modes, 2*n_modes]`` + cutoff (int): cutoff in Fock space + parallel (bool): if ``True``, uses ``dask`` for parallelization instead of OpenMP + hbar (float): value of :math:`\hbar` in the commutation relation :math;`[\hat{x}, \hat{p}]=i\hbar` + rtol (float): the relative tolerance parameter used in ``np.allclose`` + atol (float): the absolute tolerance parameter used in ``np.allclose`` + + Returns: + (array): Fock space probabilities up to cutoff. The shape of this tensor is ``[cutoff]*num_modes``. + """ + if is_pure_cov(cov, hbar=hbar, rtol=rtol, atol=atol): # Check if the covariance matrix cov is pure + return np.abs(state_vector(mu, cov, cutoff=cutoff, hbar=hbar, check_purity=False)) ** 2 + num_modes = len(mu) // 2 + + if parallel: + compute_list = [] + # create a list of parallelizable computations + for i in product(range(cutoff), repeat=num_modes): + compute_list.append(dask.delayed(density_matrix_element)(mu, cov, i, i, hbar=hbar)) + + probs = np.maximum( + 0.0, np.real_if_close(dask.compute(*compute_list, scheduler="processes")) + ).reshape([cutoff] * num_modes) + # maximum is needed because sometimes a probability is very close to zero from below + else: + probs = np.zeros([cutoff] * num_modes) + for i in product(range(cutoff), repeat=num_modes): + probs[i] = np.maximum( + 0.0, np.real_if_close(density_matrix_element(mu, cov, i, i, hbar=hbar)) + ) + # maximum is needed because sometimes a probability is very close to zero from below + return probs + +@jit(nopython=True) +def _loss_mat(eta, cutoff): # pragma: no cover + r"""Constructs a binomial loss matrix with transmission eta up to n photons. + + Args: + eta (float): Transmission coefficient. ``eta=0.0`` corresponds to complete loss and ``eta=1.0`` corresponds to no loss. + cutoff (int): cutoff in Fock space. + + Returns: + array: :math:`n\times n` matrix representing the loss. + """ + # If full transmission return the identity + + if eta < 0.0 or eta > 1.0: + raise ValueError("The transmission parameter eta should be a number between 0 and 1.") + + if eta == 1.0: + return np.identity(cutoff) + + # Otherwise construct the matrix elements recursively + lm = np.zeros((cutoff, cutoff)) + mu = 1.0 - eta + lm[:, 0] = mu ** (np.arange(cutoff)) + for i in range(cutoff): + for j in range(1, i + 1): + lm[i, j] = lm[i, j - 1] * (eta / mu) * (i - j + 1) / (j) + return lm + +def update_probabilities_with_loss(etas, probs): + """Given a list of transmissivities a tensor of probabilitites, calculate + an updated tensor of probabilities after loss is applied. + + Args: + etas (list): List of transmissitivities describing the loss in each of the modes + probs (array): Array of probabilitites in the different modes + + Returns: + array: List of loss-updated probabilities with the same shape as probs. + """ + + probs_shape = probs.shape + if len(probs_shape) != len(etas): + raise ValueError("The list of transmission etas and the tensor of probabilities probs have incompatible dimensions.") + + alphabet = "abcdefghijklmnopqrstuvwxyz" + cutoff = probs_shape[0] + for i, eta in enumerate(etas): + einstrings = "ij,{}i...->{}j...".format(alphabet[:i], alphabet[:i]) + + qein = np.zeros_like(probs) + qein = np.einsum(einstrings, _loss_mat(eta, cutoff), probs) + probs = np.copy(qein) + return qein + +@jit(nopython=True) +def _update_1d(probs, one_d, cutoff): # pragma: no cover + """ Performs a convolution of the two arrays. The first one does not need to be one dimensional, which is why we do not use ``np.convolve``. + + Args: + probs (array): (multidimensional) array + one_d (array): one dimensional array + cutoff (int): cutoff in Fock space for the first array + + Returns: + (array): the convolution of the two arrays, with the same shape as ``probs``. + """ + new_d = np.zeros_like(probs) + for i in range(cutoff): + for j in range(min(i + 1, len(one_d))): + new_d[i] += probs[i - j] * one_d[j] + return new_d + +def update_probabilities_with_noise(probs_noise, probs): + """Given a list of noise probability distributions for each of the modes and a tensor of + probabilitites, calculate an updated tensor of probabilities after noise is applied. + + Args: + probs_noise (list): List of probability distributions describing the noise in each of the modes + probs (array): Array of probabilitites in the different modes + + Returns: + array: List of noise-updated probabilities with the same shape as probs. + """ + probs_shape = probs.shape + num_modes = len(probs_shape) + cutoff = probs_shape[0] + if num_modes != len(probs_noise): + raise ValueError( + "The list of probability distributions probs_noise and the tensor of probabilities probs have incompatible dimensions." + ) + + for k in range(num_modes): #update one mode at a time + perm = np.arange(num_modes) + perm[0] = k + perm[k] = 0 + one_d = probs_noise[k] + probs_masked = np.transpose(probs, axes=perm) + probs_masked = _update_1d(probs_masked, one_d, cutoff) + probs = np.transpose(probs_masked, axes=perm) + return probs + + +def fidelity(mu1, cov1, mu2, cov2, hbar=2, rtol=1e-05, atol=1e-08): + """Calculates the fidelity between two Gaussian quantum states. + + Note that if the covariance matrices correspond to pure states this + function reduces to the modulus square of the overlap of their state vectors. + For the derivation see `'Quantum Fidelity for Arbitrary Gaussian States', Banchi et al. <10.1103/PhysRevLett.115.260501>`_. + + Args: + mu1 (array): vector of means of the first state + cov1 (array): covariance matrix of the first state + mu2 (array): vector of means of the second state + cov2 (array): covariance matrix of the second state + hbar (float): value of hbar in the uncertainty relation + rtol (float): the relative tolerance parameter used in `np.allclose` + atol (float): the absolute tolerance parameter used in `np.allclose` + + Returns: + (float): value of the fidelity between the two states + """ + + n0, n1 = cov1.shape + m0, m1 = cov2.shape + (l0,) = mu1.shape + (l1,) = mu1.shape + if not n0 == n1 == m0 == m1 == l0 == l1: + raise ValueError("The inputs have incompatible shapes") + + v1 = cov1 / hbar + v2 = cov2 / hbar + deltar = (mu1 - mu2) / np.sqrt(hbar / 2) + n0, n1 = cov1.shape + n = n0 // 2 + W = sympmat(n) + + si12 = np.linalg.inv(v1 + v2) + vaux = W.T @ si12 @ (0.25 * W + v2 @ W @ v1) + p1 = vaux @ W + p1 = p1 @ p1 + p1 = np.identity(2 * n) + 0.25 * np.linalg.inv(p1) + if np.allclose(p1, 0, rtol=rtol, atol=atol): + p1 = np.zeros_like(p1) + else: + p1 = sqrtm(p1) + p1 = 2 * (p1 + np.identity(2 * n)) + p1 = p1 @ vaux + f = np.sqrt(np.linalg.det(si12) * np.linalg.det(p1)) * np.exp( + -0.25 * deltar @ si12 @ deltar + ) + return f + + +def prefactor(mu, cov, hbar=2): + r"""Returns the prefactor. + + .. math:: prefactor = \frac{e^{-\beta Q^{-1}\beta^*/2}}{n_1!\cdots n_m! \sqrt{|Q|}} + + Args: + mu (array): length-:math:`2N` vector of mean values :math:`[\alpha,\alpha^*]` + cov (array): length-:math:`2N` `xp`-covariance matrix + + Returns: + float: the prefactor + """ + Q = Qmat(cov, hbar=hbar) + beta = Beta(mu, hbar=hbar) + Qinv = np.linalg.inv(Q) + return np.exp(-0.5 * beta @ Qinv @ beta.conj()) / np.sqrt(np.linalg.det(Q)) diff --git a/thewalrus/quantum/means_and_variances.py b/thewalrus/quantum/means_and_variances.py new file mode 100644 index 000000000..b53a74900 --- /dev/null +++ b/thewalrus/quantum/means_and_variances.py @@ -0,0 +1,236 @@ +# Copyright 2019-2020 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. +""" +Functions for constructing/calculating the means, variances and covariances of +Gaussian states and photon number distributions of Gaussian states. +""" + +from itertools import product + +import numpy as np + +from .covariance_matrices import normal_ordered_expectation + +################################################################################ +# Construct the reduced means and cov of a Gaussian state +################################################################################ + + +def reduced_gaussian(mu, cov, modes): + r""" Returns the vector of means and the covariance matrix of the specified modes. + + Args: + mu (array): a length-:math:`2N` ``np.float64`` vector of means. + cov (array): a :math:`2N\times 2N` ``np.float64`` covariance matrix + representing an :math:`N` mode quantum state. + modes (int of Sequence[int]): indices of the requested modes + + Returns: + tuple (means, cov): where means is an array containing the vector of means, + and cov is a square array containing the covariance matrix. + """ + N = len(mu) // 2 + + # reduce rho down to specified subsystems + if isinstance(modes, int): + modes = [modes] + + if np.any(np.array(modes) > N): + raise ValueError("Provided mode is larger than the number of subsystems.") + + if len(modes) == N: + # reduced state is full state + return mu, cov + + ind = np.concatenate([np.array(modes), np.array(modes) + N]) + rows = ind.reshape(-1, 1) + cols = ind.reshape(1, -1) + + return mu[ind], cov[rows, cols] + + +################################################################################ +# Calculate means or variances of photon number +################################################################################ + +def photon_number_mean(mu, cov, j, hbar=2): + r""" Calculate the mean photon number of mode j of a Gaussian state. + + Args: + mu (array): vector of means of the Gaussian state using the ordering + :math:`[q_1, q_2, \dots, q_n, p_1, p_2, \dots, p_n]` + cov (array): the covariance matrix of the Gaussian state + j (int): the j :sup:`th` mode + hbar (float): the ``hbar`` convention used in the commutation + relation :math:`[q, p]=i\hbar` + + Returns: + float: the mean photon number in mode :math:`j`. + """ + num_modes = len(mu) // 2 + return ( + mu[j] ** 2 + + mu[j + num_modes] ** 2 + + cov[j, j] + + cov[j + num_modes, j + num_modes] + - hbar + ) / (2 * hbar) + + +def photon_number_mean_vector(mu, cov, hbar=2): + r""" Calculate the mean photon number of each of the modes in a Gaussian state + + Args: + mu (array): vector of means of the Gaussian state using the ordering + :math:`[q_1, q_2, \dots, q_n, p_1, p_2, \dots, p_n]` + cov (array): the covariance matrix of the Gaussian state + hbar (float): the ``hbar`` convention used in the commutation + relation :math:`[q, p]=i\hbar` + + Returns: + array: the vector of means of the photon number distribution + """ + + N = len(mu) // 2 + return np.array([photon_number_mean(mu, cov, j, hbar=hbar) for j in range(N)]) + + +def photon_number_covar(mu, cov, j, k, hbar=2): + r""" Calculate the variance/covariance of the photon number distribution + of a Gaussian state. + + Implements the covariance matrix of the photon number distribution of a + Gaussian state according to the Last two eq. of Part II. in + `'Multidimensional Hermite polynomials and photon distribution for polymode + mixed light', Dodonov et al. `_ + + .. math:: + \sigma_{n_j n_j} &= \frac{1}{2}\left(T_j^2 - 2d_j - \frac{1}{2}\right) + + \left<\mathbf{Q}_j\right>\mathcal{M}_j\left<\mathbf{Q}_j\right>, \\ + \sigma_{n_j n_k} &= \frac{1}{2}\mathrm{Tr}\left(\Lambda_j \mathbf{M} \Lambda_k \mathbf{M}\right) + + \frac{1}{2}\left<\mathbf{Q}\right>\Lambda_j \mathbf{M} \Lambda_k\left<\mathbf{Q}\right>, + + where :math:`T_j` and :math:`d_j` are the trace and the determinant of + :math:`2 \times 2` matrix :math:`\mathcal{M}_j` whose elements coincide + with the nonzero elements of matrix :math:`\mathbf{M}_j = \Lambda_j \mathbf{M} \Lambda_k` + while the two-vector :math:`\mathbf{Q}_j` has the components :math:`(q_j, p_j)`. + :math:`2N \times 2N` projector matrix :math:`\Lambda_j` has only two nonzero + elements: :math:`\left(\Lambda_j\right)_{jj} = \left(\Lambda_j\right)_{j+N,j+N} = 1`. + Note that the convention for ``mu`` used here differs from the one used in Dodonov et al., + They both provide the same results in this particular case. + + Args: + mu (array): vector of means of the Gaussian state using the ordering + :math:`[q_1, q_2, \dots, q_n, p_1, p_2, \dots, p_n]` + cov (array): the covariance matrix of the Gaussian state + j (int): the j :sup:`th` mode + k (int): the k :sup:`th` mode + hbar (float): the ``hbar`` convention used in the commutation + relation :math:`[q, p]=i\hbar` + + Returns: + float: the covariance for the photon numbers at modes :math:`j` and :math:`k`. + """ + # renormalise the covariance matrix + cov = cov / hbar + + N = len(mu) // 2 + mu = np.array(mu) / np.sqrt(hbar) + + lambda_1 = np.zeros((2 * N, 2 * N)) + lambda_1[j, j] = lambda_1[j + N, j + N] = 1 + + lambda_2 = np.zeros((2 * N, 2 * N)) + lambda_2[k, k] = lambda_2[k + N, k + N] = 1 + + if j == k: + idxs = ((j, j, j + N, j + N), (j, j + N, j, j + N)) + M = (lambda_1 @ cov @ lambda_2)[idxs].reshape(2, 2) + + term_1 = (np.trace(M) ** 2 - 2 * np.linalg.det(M) - 0.5) / 2 + term_2 = mu[[j, j + N]] @ M @ mu[[j, j + N]] + else: + term_1 = np.trace(lambda_1 @ cov @ lambda_2 @ cov) / 2 + term_2 = (mu @ lambda_1 @ cov @ lambda_2 @ mu) / 2 + + return term_1 + term_2 + + +def photon_number_covmat(mu, cov, hbar=2): + r""" Calculate the covariance matrix of the photon number distribution of a + Gaussian state. + + Args: + mu (array): vector of means of the Gaussian state using the ordering + :math:`[q_1, q_2, \dots, q_n, p_1, p_2, \dots, p_n]` + cov (array): the covariance matrix of the Gaussian state + hbar (float): the ``hbar`` convention used in the commutation + relation :math:`[q, p]=i\hbar` + + Returns: + array: the covariance matrix of the photon number distribution + """ + N = len(mu) // 2 + pnd_cov = np.zeros((N, N)) + for i in range(N): + for j in range(i+1): + pnd_cov[i][j] = photon_number_covar(mu, cov, i, j, hbar=hbar) + pnd_cov[j][i] = pnd_cov[i][j] + return pnd_cov + + +def photon_number_expectation(mu, cov, modes, hbar=2): + r"""Calculates the expectation value of the product of the number operator of the modes in a Gaussian state. + + Args: + mu (array): length-:math:`2N` means vector in xp-ordering. + cov (array): :math:`2N\times 2N` covariance matrix in xp-ordering. + modes (list): list of modes + hbar (float): value of hbar in the uncertainty relation. + + Returns: + (float): expectation value of the product of the number operators of the modes. + """ + n, _ = cov.shape + n_modes = n // 2 + rpt = np.zeros([n], dtype=int) + for i in modes: + rpt[i] = 1 + rpt[i + n_modes] = 1 + + return normal_ordered_expectation(mu, cov, rpt, hbar=hbar) + + +def photon_number_squared_expectation(mu, cov, modes, hbar=2): + r"""Calculates the expectation value of the square of the product of the number operator of the modes in + a Gaussian state. + + Args: + mu (array): length-:math:`2N` means vector in xp-ordering. + cov (array): :math:`2N\times 2N` covariance matrix in xp-ordering. + modes (list): list of modes + hbar (float): value of hbar in the uncertainty relation. + + Returns: + (float): expectation value of the square of the product of the number operator of the modes. + """ + n_modes = len(modes) + + mu_red, cov_red = reduced_gaussian(mu, cov, modes) + result = 0 + for item in product([1, 2], repeat=n_modes): + rpt = item + item + term = normal_ordered_expectation(mu_red, cov_red, rpt, hbar=hbar) + result += term + return result diff --git a/thewalrus/quantum/photon_number_distributions.py b/thewalrus/quantum/photon_number_distributions.py new file mode 100644 index 000000000..526bd48bf --- /dev/null +++ b/thewalrus/quantum/photon_number_distributions.py @@ -0,0 +1,94 @@ +# Copyright 2019-2020 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. +""" +Functions for calculating the photon number distributions of various states. +""" + +import numpy as np +from scipy.stats import nbinom + +from .covariance_matrices import Amat, is_pure_cov + +################################################################################ +# Calculate the total photon number distribution +################################################################################ + + +def total_photon_num_dist_pure_state(cov, cutoff=50, hbar=2, padding_factor=2): + r""" Calculates the total photon number distribution of a pure state + with zero mean. + + Args: + cov (array): :math:`2N\times 2N` covariance matrix in xp-ordering + cutoff (int): Fock cutoff + tol (float): tolerance for determining if displacement is negligible + hbar (float): the value of :math:`\hbar` in the commutation + padding_factor (int): expanded size of the photon distribution to avoid accumulation of errors + + Returns: + (array): Total photon number distribution + """ + if is_pure_cov(cov): + A = Amat(cov, hbar=hbar) + (n, _) = A.shape + N = n // 2 + B = A[0:N, 0:N] + rs = np.arctanh(np.linalg.svd(B, compute_uv=False)) + return gen_multi_mode_dist(rs, cutoff=cutoff, padding_factor=padding_factor)[0:cutoff] + raise ValueError("The Gaussian state is not pure") + + +def gen_single_mode_dist(s, cutoff=50, N=1): + """Generate the photon number distribution of :math:`N` identical single mode squeezed states. + + Args: + s (float): squeezing parameter + cutoff (int): Fock cutoff + N (float): number of squeezed states + + Returns: + (array): Photon number distribution + """ + r = 0.5 * N + q = 1.0 - np.tanh(s) ** 2 + N = cutoff // 2 + ps_tot = np.zeros(cutoff) + if cutoff % 2 == 0: + ps = nbinom.pmf(np.arange(N), p=q, n=r) + ps_tot[0::2] = ps + else: + ps = nbinom.pmf(np.arange(N + 1), p=q, n=r) + ps_tot[0:-1][0::2] = ps[0:-1] + ps_tot[-1] = ps[-1] + + return ps_tot + + +def gen_multi_mode_dist(s, cutoff=50, padding_factor=2): + """Generates the total photon number distribution of single mode squeezed states with different squeezing values. + + Args: + s (array): array of squeezing parameters + cutoff (int): Fock cutoff + + Returns: + (array[int]): total photon number distribution + """ + scale = padding_factor + cutoff_sc = scale * cutoff + ps = np.zeros(cutoff_sc) + ps[0] = 1.0 + for s_val in s: + ps = np.convolve(ps, gen_single_mode_dist(s_val, cutoff_sc))[0:cutoff_sc] + return ps From 0ca5a00103937b3a21b6e590b97a2a07e04d32dd Mon Sep 17 00:00:00 2001 From: Theodor Isacsson Date: Fri, 4 Sep 2020 18:05:25 -0400 Subject: [PATCH 02/17] Fix imports --- thewalrus/quantum/__init__.py | 11 +- thewalrus/quantum/adjacency_matrices.py | 156 ++++++++++++------- thewalrus/quantum/covariance_matrices.py | 66 +++++--- thewalrus/quantum/fock_states_and_tensors.py | 33 +++- thewalrus/quantum/means_and_variances.py | 37 ----- 5 files changed, 185 insertions(+), 118 deletions(-) diff --git a/thewalrus/quantum/__init__.py b/thewalrus/quantum/__init__.py index 915c7cb8e..d3c9fa1e4 100644 --- a/thewalrus/quantum/__init__.py +++ b/thewalrus/quantum/__init__.py @@ -105,6 +105,7 @@ prefactor find_scaling_adjacency_matrix mean_number_of_clicks + mean_number_of_clicks_graph find_scaling_adjacency_matrix_torontonian gen_Qmat_from_graph photon_number_mean @@ -114,6 +115,7 @@ is_valid_cov is_pure_cov is_classical_cov + find_classical_subsystem total_photon_num_dist_pure_state gen_single_mode_dist gen_multi_mode_dist @@ -134,23 +136,29 @@ update_probabilities_with_loss, update_probabilities_with_noise, fidelity, + + find_classical_subsystem, prefactor, + loss_mat, ) from .adjacency_matrices import ( mean_number_of_clicks, + mean_number_of_clicks_graph, + variance_number_of_clicks, find_scaling_adjacency_matrix, find_scaling_adjacency_matrix_torontonian, gen_Qmat_from_graph, ) from .covariance_matrices import ( + reduced_gaussian, + Xmat, Qmat, Covmat, Amat, normal_ordered_expectation, - normal_ordered_complex_cov, Beta, Means, @@ -160,7 +168,6 @@ ) from .means_and_variances import ( - reduced_gaussian, photon_number_mean, photon_number_mean_vector, photon_number_covar, diff --git a/thewalrus/quantum/adjacency_matrices.py b/thewalrus/quantum/adjacency_matrices.py index f661e76da..882895415 100644 --- a/thewalrus/quantum/adjacency_matrices.py +++ b/thewalrus/quantum/adjacency_matrices.py @@ -19,9 +19,66 @@ import numpy as np from scipy.optimize import root_scalar -from .covariance_matrices import Xmat +from .covariance_matrices import Xmat, Covmat, Qmat, reduced_gaussian -def mean_number_of_clicks(A): +def mean_number_of_clicks(cov, hbar=2): + r""" Calculates the total mean number of clicks when a zero-mean gaussian state + is measured using threshold detectors. + + Args + cov (array): :math:`2N\times 2N` covariance matrix in xp-ordering + hbar (float): the value of :math:`\hbar` in the commutation relation :math:`[\x,\p]=i\hbar` + + Returns + float: mean number of clicks + """ + n, _ = cov.shape + nmodes = n // 2 + Q = Qmat(cov, hbar=hbar) + meanc = 1.0 * nmodes + + for i in range(nmodes): + det_val = np.real(Q[i, i] * Q[i + nmodes, i + nmodes] - Q[i + nmodes, i] * Q[i, i + nmodes]) + meanc -= 1.0 / np.sqrt(det_val) + return meanc + + +def variance_number_of_clicks(cov, hbar=2): + r""" Calculates the variance of the total number of clicks when a zero-mean gaussian state + is measured using threshold detectors. + + Args + cov (array): :math:`2N\times 2N` covariance matrix in xp-ordering + hbar (float): the value of :math:`\hbar` in the commutation relation :math:`[\x,\p]=i\hbar` + + Returns + float: variance in the total number of clicks + """ + n, _ = cov.shape + means = np.zeros([n]) + nmodes = n // 2 + Q = Qmat(cov, hbar=hbar) + vac_probs = np.array( + [ + np.real(Q[i, i] * Q[i + nmodes, i + nmodes] - Q[i + nmodes, i] * Q[i, i + nmodes]) + for i in range(nmodes) + ] + ) + vac_probs = np.sqrt(vac_probs) + vac_probs = 1 / vac_probs + term1 = np.sum(vac_probs * (1 - vac_probs)) + term2 = 0 + for i in range(nmodes): + for j in range(i): + _, Qij = reduced_gaussian(means, Q, [i, j]) + prob_vac_ij = np.linalg.det(Qij).real + prob_vac_ij = 1.0 / np.sqrt(prob_vac_ij) + term2 += prob_vac_ij - vac_probs[i] * vac_probs[j] + + return term1 + 2 * term2 + + +def mean_number_of_clicks_graph(A): r""" Given an adjacency matrix this function calculates the mean number of clicks. For this to make sense the user must provide a matrix with singular values less than or equal to one. See Appendix A.3 of `_ @@ -38,12 +95,7 @@ def mean_number_of_clicks(A): X = np.block([[0 * idn, idn], [idn, 0 * idn]]) B = np.block([[A, 0 * A], [0 * A, np.conj(A)]]) Q = np.linalg.inv(np.identity(2 * n) - X @ B) - meanc = 1.0 * n - - for i in range(n): - det_val = np.real(Q[i, i]*Q[i+n, i+n] - Q[i+n, i]*Q[i, i+n]) - meanc -= 1.0 / np.sqrt(det_val) - return meanc + return mean_number_of_clicks(Covmat(Q)) ################################################################################ @@ -51,6 +103,50 @@ def mean_number_of_clicks(A): ################################################################################ +def find_scaling_adjacency_matrix_torontonian(A, c_mean): + r""" Returns the scaling parameter by which the adjacency matrix A + should be rescaled so that the Gaussian state that encodes it has + give a mean number of clicks equal to ``c_mean`` when measured with + threshold detectors. + + Args: + A (array): adjacency matrix + c_mean (float): mean photon number of the Gaussian state + + Returns: + float: scaling parameter + """ + n, _ = A.shape + if c_mean < 0 or c_mean > n: + raise ValueError("The mean number of clicks should be smaller than the number of modes") + + vals = np.linalg.svd(A, compute_uv=False) + localA = A / vals[0] # rescale the matrix so that the singular values are between 0 and 1. + + def cost(x): + r""" Cost function giving the difference between the wanted number of clicks and the number + of clicks at a given scaling value. It assumes that the adjacency matrix has been rescaled + so that it has singular values between 0 and 1. + + Args: + x (float): scaling value + + Return: + float: difference between desired and obtained mean number of clicks + """ + if x >= 1.0: + return c_mean - n + if x <= 0: + return c_mean + return c_mean - mean_number_of_clicks_graph(x * localA) + + res = root_scalar(cost, x0=0.5, bracket=(0.0, 1.0)) # Do the optimization + + if not res.converged: + raise ValueError("The search for a scaling value failed") + return res.root / vals[0] + + def find_scaling_adjacency_matrix(A, n_mean): r""" Returns the scaling parameter by which the adjacency matrix A should be rescaled so that the Gaussian state that endodes it has @@ -116,50 +212,6 @@ def grad_mean_photon_number(x, vals): # pragma: no cover return res.root -def find_scaling_adjacency_matrix_torontonian(A, c_mean): - r""" Returns the scaling parameter by which the adjacency matrix A - should be rescaled so that the Gaussian state that encodes it has - give a mean number of clicks equal to ``c_mean`` when measured with - threshold detectors. - - Args: - A (array): adjacency matrix - c_mean (float): mean photon number of the Gaussian state - - Returns: - float: scaling parameter - """ - n, _ = A.shape - if c_mean < 0 or c_mean > n: - raise ValueError("The mean number of clicks should be smaller than the number of modes") - - vals = np.linalg.svd(A, compute_uv=False) - localA = A / vals[0] # rescale the matrix so that the singular values are between 0 and 1. - - def cost(x): - r""" Cost function giving the difference between the wanted number of clicks and the number - of clicks at a given scaling value. It assumes that the adjacency matrix has been rescaled - so that it has singular values between 0 and 1. - - Args: - x (float): scaling value - - Return: - float: difference between desired and obtained mean number of clicks - """ - if x >= 1.0: - return c_mean - n - if x <= 0: - return c_mean - return c_mean - mean_number_of_clicks(x * localA) - - res = root_scalar(cost, x0=0.5, bracket=(0.0, 1.0)) # Do the optimization - - if not res.converged: - raise ValueError("The search for a scaling value failed") - return res.root / vals[0] - - def gen_Qmat_from_graph(A, n_mean): r""" Returns the Qmat xp-covariance matrix associated to a graph with adjacency matrix :math:`A` and with mean photon number :math:`n_{mean}`. diff --git a/thewalrus/quantum/covariance_matrices.py b/thewalrus/quantum/covariance_matrices.py index d5426ef00..917b5d637 100644 --- a/thewalrus/quantum/covariance_matrices.py +++ b/thewalrus/quantum/covariance_matrices.py @@ -21,6 +21,45 @@ from ..symplectic import sympmat from .._hafnian import hafnian, reduction + +################################################################################ +# Construct the reduced means and cov of a Gaussian state +################################################################################ + + +def reduced_gaussian(mu, cov, modes): + r""" Returns the vector of means and the covariance matrix of the specified modes. + + Args: + mu (array): a length-:math:`2N` ``np.float64`` vector of means. + cov (array): a :math:`2N\times 2N` ``np.float64`` covariance matrix + representing an :math:`N` mode quantum state. + modes (int of Sequence[int]): indices of the requested modes + + Returns: + tuple (means, cov): where means is an array containing the vector of means, + and cov is a square array containing the covariance matrix. + """ + N = len(mu) // 2 + + # reduce rho down to specified subsystems + if isinstance(modes, int): + modes = [modes] + + if np.any(np.array(modes) > N): + raise ValueError("Provided mode is larger than the number of subsystems.") + + if len(modes) == N: + # reduced state is full state + return mu, cov + + ind = np.concatenate([np.array(modes), np.array(modes) + N]) + rows = ind.reshape(-1, 1) + cols = ind.reshape(1, -1) + + return mu[ind], cov[rows, cols] + + ################################################################################ # Transform one type of covariance-matrix-like object into another ################################################################################ @@ -139,7 +178,8 @@ def normal_ordered_expectation(mu, cov, rpt, hbar=2): (float): expectation value of the normal ordered product of operators """ alpha = Beta(mu, hbar=hbar) - V = normal_ordered_complex_cov(cov, hbar=hbar) + n = len(cov) + V = (Qmat(cov, hbar=hbar) - np.identity(n)) @ Xmat(n // 2) A = reduction(V, rpt) if np.allclose(mu, 0): res = np.conj(hafnian(A)) @@ -149,29 +189,6 @@ def normal_ordered_expectation(mu, cov, rpt, hbar=2): return np.conj(res) -def normal_ordered_complex_cov(cov, hbar=2): - r"""Calculates the normal ordered covariance matrix in the complex basis. - - Args: - cov (array): xp-covariance matrix. - hbar (float): value of hbar in the uncertainty relation. - - Returns: - (array): covariance matrix in the creation/annihilation operator basis. - """ - - n, _ = cov.shape - n_modes = n // 2 - cov = cov / (hbar / 2) - A = cov[:n_modes, :n_modes] - B = cov[:n_modes, n_modes:] - C = cov[n_modes:, n_modes:] - N = 0.25 * (A + C + 1j * (B - B.T) - 2 * np.identity(n_modes)) - M = 0.25 * (A - C + 1j * (B + B.T)) - mat = np.block([[M.conj(), N], [N.T, M]]) - return mat - - def Beta(mu, hbar=2): r"""Returns the vector of complex displacements and conjugate displacements. @@ -273,6 +290,7 @@ def is_classical_cov(cov, hbar=2, atol=1e-08): Args: cov (array): a covariance matrix hbar (float): value of hbar in the uncertainty relation + atol (float): the absolute tolerance parameter used in `np.allclose` Returns: (bool): whether the given covariance matrix corresponds to a classical state diff --git a/thewalrus/quantum/fock_states_and_tensors.py b/thewalrus/quantum/fock_states_and_tensors.py index 904d88a3a..cbaabf23b 100644 --- a/thewalrus/quantum/fock_states_and_tensors.py +++ b/thewalrus/quantum/fock_states_and_tensors.py @@ -32,14 +32,15 @@ from .._hermite_multidimensional import hermite_multidimensional, hafnian_batched from .covariance_matrices import ( + is_classical_cov, is_pure_cov, Amat, Beta, Qmat, + reduced_gaussian, ) -# pylint: disable=too-many-arguments def pure_state_amplitude(mu, cov, i, include_prefactor=True, tol=1e-10, hbar=2, check_purity=True): r"""Returns the :math:`\langle i | \psi\rangle` element of the state ket of a Gaussian state defined by covariance matrix cov. @@ -453,7 +454,7 @@ def probabilities(mu, cov, cutoff, parallel=False, hbar=2.0, rtol=1e-05, atol=1e return probs @jit(nopython=True) -def _loss_mat(eta, cutoff): # pragma: no cover +def loss_mat(eta, cutoff): # pragma: no cover r"""Constructs a binomial loss matrix with transmission eta up to n photons. Args: @@ -502,7 +503,7 @@ def update_probabilities_with_loss(etas, probs): einstrings = "ij,{}i...->{}j...".format(alphabet[:i], alphabet[:i]) qein = np.zeros_like(probs) - qein = np.einsum(einstrings, _loss_mat(eta, cutoff), probs) + qein = np.einsum(einstrings, loss_mat(eta, cutoff), probs) probs = np.copy(qein) return qein @@ -605,6 +606,32 @@ def fidelity(mu1, cov1, mu2, cov2, hbar=2, rtol=1e-05, atol=1e-08): return f +def find_classical_subsystem(cov, hbar=2, atol=1e-08): + """Find the largest integer ``k`` so that subsystem in modes ``[0,1,...,k-1]`` is a classical state. + + + Args: + cov (array): a covariance matrix + hbar (float): value of hbar in the uncertainty relation + atol (float): the absolute tolerance parameter used when determining if the state is classical + + Returns: + int: the largest k so that modes ``[0,1,...,k-1]`` are in a classical state. + """ + n, _ = cov.shape + nmodes = n // 2 + if is_classical_cov(cov, hbar=hbar, atol=atol): + return nmodes + k = 0 + mu = np.zeros(n) + is_classical = True + while is_classical: + _, Vk = reduced_gaussian(mu, cov, list(range(k + 1))) + is_classical = is_classical_cov(Vk, hbar=hbar, atol=atol) + k += 1 + return k - 1 + + def prefactor(mu, cov, hbar=2): r"""Returns the prefactor. diff --git a/thewalrus/quantum/means_and_variances.py b/thewalrus/quantum/means_and_variances.py index b53a74900..4ed183b3a 100644 --- a/thewalrus/quantum/means_and_variances.py +++ b/thewalrus/quantum/means_and_variances.py @@ -22,43 +22,6 @@ from .covariance_matrices import normal_ordered_expectation -################################################################################ -# Construct the reduced means and cov of a Gaussian state -################################################################################ - - -def reduced_gaussian(mu, cov, modes): - r""" Returns the vector of means and the covariance matrix of the specified modes. - - Args: - mu (array): a length-:math:`2N` ``np.float64`` vector of means. - cov (array): a :math:`2N\times 2N` ``np.float64`` covariance matrix - representing an :math:`N` mode quantum state. - modes (int of Sequence[int]): indices of the requested modes - - Returns: - tuple (means, cov): where means is an array containing the vector of means, - and cov is a square array containing the covariance matrix. - """ - N = len(mu) // 2 - - # reduce rho down to specified subsystems - if isinstance(modes, int): - modes = [modes] - - if np.any(np.array(modes) > N): - raise ValueError("Provided mode is larger than the number of subsystems.") - - if len(modes) == N: - # reduced state is full state - return mu, cov - - ind = np.concatenate([np.array(modes), np.array(modes) + N]) - rows = ind.reshape(-1, 1) - cols = ind.reshape(1, -1) - - return mu[ind], cov[rows, cols] - ################################################################################ # Calculate means or variances of photon number From 4f448a4c1c37ac0d01376441db33c7a8e0a97f82 Mon Sep 17 00:00:00 2001 From: Theodor Isacsson Date: Fri, 4 Sep 2020 18:05:37 -0400 Subject: [PATCH 03/17] Remove old quantum.py --- thewalrus/quantum.py | 1471 ------------------------------------------ 1 file changed, 1471 deletions(-) delete mode 100644 thewalrus/quantum.py diff --git a/thewalrus/quantum.py b/thewalrus/quantum.py deleted file mode 100644 index 16850fb67..000000000 --- a/thewalrus/quantum.py +++ /dev/null @@ -1,1471 +0,0 @@ -# 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. -""" -Quantum algorithms -================== - -.. currentmodule:: thewalrus.quantum - -This submodule provides access to various utility functions that act on Gaussian -quantum states. - -For more details on how the hafnian relates to various properties of Gaussian quantum -states, see: - -* Kruse, R., Hamilton, C. S., Sansoni, L., Barkhofen, S., Silberhorn, C., & Jex, I. - "Detailed study of Gaussian boson sampling." `Physical Review A 100, 032326 (2019) - `_ - -* Hamilton, C. S., Kruse, R., Sansoni, L., Barkhofen, S., Silberhorn, C., & Jex, I. - "Gaussian boson sampling." `Physical Review Letters, 119(17), 170501 (2017) - `_ - -* Quesada, N. - "Franck-Condon factors by counting perfect matchings of graphs with loops." - `Journal of Chemical Physics 150, 164113 (2019) `_ - -* Quesada, N., Helt, L. G., Izaac, J., Arrazola, J. M., Shahrokhshahi, R., Myers, C. R., & Sabapathy, K. K. - "Simulating realistic non-Gaussian state preparation." `Physical Review A 100, 022341 (2019) - `_ - - -Fock states and tensors ------------------------ - -.. autosummary:: - - pure_state_amplitude - state_vector - density_matrix_element - density_matrix - fock_tensor - probabilities - update_probabilities_with_loss - update_probabilities_with_noise - normal_ordered_expectation - fidelity - -Details -^^^^^^^ - -.. autofunction:: - pure_state_amplitude - -.. autofunction:: - state_vector - -.. autofunction:: - density_matrix_element - -.. autofunction:: - density_matrix - -.. autofunction:: - fock_tensor - -.. autofunction:: - probabilities - -.. autofunction:: - update_probabilities_with_loss - -.. autofunction:: - update_probabilities_with_noise - -.. autofunction:: - normal_ordered_expectation - -.. autofunction:: - fidelity - - -Utility functions ------------------ - -.. autosummary:: - - reduced_gaussian - Xmat - Qmat - Covmat - Amat - Beta - Means - prefactor - find_scaling_adjacency_matrix - mean_number_of_clicks - find_scaling_adjacency_matrix_torontonian - gen_Qmat_from_graph - photon_number_mean - photon_number_mean_vector - photon_number_covar - photon_number_covmat - is_valid_cov - is_pure_cov - is_classical_cov - find_classical_subsystem - total_photon_num_dist_pure_state - gen_single_mode_dist - gen_multi_mode_dist - variance_number_of_clicks - mean_number_of_clicks_graph - -Details -^^^^^^^ -""" -# pylint: disable=too-many-arguments -from itertools import count, product, chain - -import numpy as np -import dask - -from scipy.optimize import root_scalar -from scipy.special import factorial as fac -from scipy.stats import nbinom -from scipy.linalg import sqrtm -from numba import jit - -from thewalrus.symplectic import expand, sympmat, is_symplectic -from thewalrus.libwalrus import interferometer, interferometer_real - -from ._hafnian import hafnian, hafnian_repeated, reduction -from ._hermite_multidimensional import hermite_multidimensional, hafnian_batched - - -def reduced_gaussian(mu, cov, modes): - r""" Returns the vector of means and the covariance matrix of the specified modes. - - Args: - mu (array): a length-:math:`2N` ``np.float64`` vector of means. - cov (array): a :math:`2N\times 2N` ``np.float64`` covariance matrix - representing an :math:`N` mode quantum state. - modes (int of Sequence[int]): indices of the requested modes - - Returns: - tuple (means, cov): where means is an array containing the vector of means, - and cov is a square array containing the covariance matrix. - """ - N = len(mu) // 2 - - # reduce rho down to specified subsystems - if isinstance(modes, int): - modes = [modes] - - if np.any(np.array(modes) > N): - raise ValueError("Provided mode is larger than the number of subsystems.") - - if len(modes) == N: - # reduced state is full state - return mu, cov - - ind = np.concatenate([np.array(modes), np.array(modes) + N]) - rows = ind.reshape(-1, 1) - cols = ind.reshape(1, -1) - - return mu[ind], cov[rows, cols] - - -def Xmat(N): - r"""Returns the matrix :math:`X_n = \begin{bmatrix}0 & I_n\\ I_n & 0\end{bmatrix}` - - Args: - N (int): positive integer - - Returns: - array: :math:`2N\times 2N` array - """ - I = np.identity(N) - O = np.zeros_like(I) - X = np.block([[O, I], [I, O]]) - return X - - -def Qmat(cov, hbar=2): - r"""Returns the :math:`Q` Husimi matrix of the Gaussian state. - - Args: - cov (array): :math:`2N\times 2N xp-` Wigner covariance matrix - hbar (float): the value of :math:`\hbar` in the commutation - relation :math:`[\x,\p]=i\hbar`. - - Returns: - array: the :math:`Q` matrix. - """ - # number of modes - N = len(cov) // 2 - I = np.identity(N) - - x = cov[:N, :N] * 2 / hbar - xp = cov[:N, N:] * 2 / hbar - p = cov[N:, N:] * 2 / hbar - # the (Hermitian) matrix elements - aidaj = (x + p + 1j * (xp - xp.T) - 2 * I) / 4 - # the (symmetric) matrix elements - aiaj = (x - p + 1j * (xp + xp.T)) / 4 - - # calculate the covariance matrix sigma_Q appearing in the Q function: - # Q(alpha) = exp[-(alpha-beta).sigma_Q^{-1}.(alpha-beta)/2]/|sigma_Q| - Q = np.block([[aidaj, aiaj.conj()], [aiaj, aidaj.conj()]]) + np.identity(2 * N) - return Q - - -def Covmat(Q, hbar=2): - r"""Returns the Wigner covariance matrix in the :math:`xp`-ordering of the Gaussian state. - This is the inverse function of Qmat. - - Args: - Q (array): :math:`2N\times 2N` Husimi Q matrix - hbar (float): the value of :math:`\hbar` in the commutation - relation :math:`[\x,\p]=i\hbar`. - - Returns: - array: the :math:`xp`-ordered covariance matrix in the xp-ordering. - """ - # number of modes - n = len(Q) // 2 - I = np.identity(n) - N = Q[0:n, 0:n] - I - M = Q[n : 2 * n, 0:n] - mm11a = 2 * (N.real + M.real) + np.identity(n) - mm22a = 2 * (N.real - M.real) + np.identity(n) - mm12a = 2 * (M.imag + N.imag) - cov = np.block([[mm11a, mm12a], [mm12a.T, mm22a]]) - - return (hbar / 2) * cov - - -def Amat(cov, hbar=2, cov_is_qmat=False): - r"""Returns the :math:`A` matrix of the Gaussian state whose hafnian gives the photon number probabilities. - - Args: - cov (array): :math:`2N\times 2N` covariance matrix - hbar (float): the value of :math:`\hbar` in the commutation - relation :math:`[\x,\p]=i\hbar`. - cov_is_qmat (bool): if ``True``, it is assumed that ``cov`` is in fact the Q matrix. - - Returns: - array: the :math:`A` matrix. - """ - # number of modes - N = len(cov) // 2 - X = Xmat(N) - - # inverse Q matrix - if cov_is_qmat: - Q = cov - else: - Q = Qmat(cov, hbar=hbar) - - Qinv = np.linalg.inv(Q) - - # calculate Hamilton's A matrix: A = X.(I-Q^{-1})* - A = X @ (np.identity(2 * N) - Qinv).conj() - return A - - -def Beta(mu, hbar=2): - r"""Returns the vector of complex displacements and conjugate displacements. - - Args: - mu (array): length-:math:`2N` means vector - hbar (float): the value of :math:`\hbar` in the commutation - relation :math:`[\x,\p]=i\hbar`. - - Returns: - array: the expectation values - :math:`[\langle a_1\rangle, \langle a_2\rangle,\dots,\langle a_N\rangle, \langle a^\dagger_1\rangle, \dots, \langle a^\dagger_N\rangle]` - """ - N = len(mu) // 2 - # mean displacement of each mode - alpha = (mu[:N] + 1j * mu[N:]) / np.sqrt(2 * hbar) - # the expectation values (, ,...,, , ..., ) - return np.concatenate([alpha, alpha.conj()]) - - -def Means(beta, hbar=2): - r"""Returns the vector of real quadrature displacements. - - Args: - beta (array): length-:math:`2N` means bivector - hbar (float): the value of :math:`\hbar` in the commutation - relation :math:`[\x,\p]=i\hbar`. - - Returns: - array: the quadrature expectation values - :math:`[\langle q_1\rangle, \langle q_2\rangle,\dots,\langle q_N\rangle, \langle p_1\rangle, \dots, \langle p_N\rangle]` - """ - - N = len(beta) // 2 - alpha = beta[0:N] - return np.sqrt(2 * hbar) * np.concatenate([alpha.real, alpha.imag]) - - -def prefactor(mu, cov, hbar=2): - r"""Returns the prefactor. - - .. math:: prefactor = \frac{e^{-\beta Q^{-1}\beta^*/2}}{n_1!\cdots n_m! \sqrt{|Q|}} - - Args: - mu (array): length-:math:`2N` vector of mean values :math:`[\alpha,\alpha^*]` - cov (array): length-:math:`2N` `xp`-covariance matrix - - Returns: - float: the prefactor - """ - Q = Qmat(cov, hbar=hbar) - beta = Beta(mu, hbar=hbar) - Qinv = np.linalg.inv(Q) - return np.exp(-0.5 * beta @ Qinv @ beta.conj()) / np.sqrt(np.linalg.det(Q)) - - -def density_matrix_element(mu, cov, i, j, include_prefactor=True, tol=1e-10, hbar=2): - r"""Returns the :math:`\langle i | \rho | j \rangle` element of the density matrix - of a Gaussian state defined by covariance matrix cov. - - Args: - mu (array): length-:math:`2N` quadrature displacement vector - cov (array): length-:math:`2N` covariance matrix - i (list): list of density matrix rows - j (list): list of density matrix columns - include_prefactor (bool): if ``True``, the prefactor is automatically calculated - used to scale the result. - tol (float): tolerance for determining if displacement is negligible - hbar (float): the value of :math:`\hbar` in the commutation - relation :math:`[\x,\p]=i\hbar`. - - Returns: - complex: the density matrix element - """ - rpt = i + j - beta = Beta(mu, hbar=hbar) - A = Amat(cov, hbar=hbar) - if np.linalg.norm(beta) < tol: - # no displacement - if np.prod([k + 1 for k in rpt]) ** (1 / len(rpt)) < 3: - A_rpt = reduction(A, rpt) - haf = hafnian(A_rpt) - else: - haf = hafnian_repeated(A, rpt) - else: - # replace the diagonal of A with gamma - gamma = beta.conj() - A @ beta - if np.prod([k + 1 for k in rpt]) ** (1 / len(rpt)) < 3: - A_rpt = reduction(A, rpt) - np.fill_diagonal(A_rpt, reduction(gamma, rpt)) - haf = hafnian(A_rpt, loop=True) - else: - haf = hafnian_repeated(A, rpt, mu=gamma, loop=True) - - if include_prefactor: - haf *= prefactor(mu, cov, hbar=hbar) - - return haf / np.sqrt(np.prod(fac(rpt))) - - -def density_matrix(mu, cov, post_select=None, normalize=False, cutoff=5, hbar=2): - r"""Returns the density matrix of a (PNR post-selected) Gaussian state. - - The resulting density matrix will have shape - - .. math:: \underbrace{D\times D \times \cdots \times D}_{2M} - - where :math:`D` is the Fock space cutoff, and :math:`M` is the - number of *non* post-selected modes, i.e. ``M = len(mu)//2 - len(post_select)``. - - Note that we use the Strawberry Fields convention for indexing the density - matrix; the first two dimensions correspond to subsystem 1, the second two - dimensions correspond to subsystem 2, etc. - If post_select is None then the density matrix elements are calculated using - the multidimensional Hermite polynomials which provide a significantly faster - evaluation. - - Args: - mu (array): length-:math:`2N` means vector in xp-ordering - cov (array): :math:`2N\times 2N` covariance matrix in xp-ordering - post_select (dict): dictionary containing the post-selected modes, of - the form ``{mode: value}``. If post_select is None the whole non post-selected density matrix - is calculated directly using (multidimensional) Hermite polynomials, which is significantly faster - than calculating one hafnian at a time. - normalize (bool): If ``True``, a post-selected density matrix is re-normalized. - cutoff (dim): the final length (i.e., Hilbert space dimension) of each - mode in the density matrix. - hbar (float): the value of :math:`\hbar` in the commutation - relation :math:`[\x,\p]=i\hbar`. - - Returns: - np.array[complex]: the density matrix of the Gaussian state - """ - N = len(mu) // 2 - pref = prefactor(mu, cov, hbar=hbar) - - if post_select is None: - A = Amat(cov, hbar=hbar).conj() - sf_order = tuple(chain.from_iterable([[i, i + N] for i in range(N)])) - - if np.allclose(mu, np.zeros_like(mu)): - tensor = np.real_if_close(pref) * hermite_multidimensional( - -A, cutoff, renorm=True, modified=True - ) - return tensor.transpose(sf_order) - beta = Beta(mu, hbar=hbar) - y = beta - A @ beta.conj() - tensor = np.real_if_close(pref) * hermite_multidimensional( - -A, cutoff, y=y, renorm=True, modified=True - ) - return tensor.transpose(sf_order) - - M = N - len(post_select) - rho = np.zeros([cutoff] * (2 * M), dtype=np.complex128) - - for idx in product(range(cutoff), repeat=2 * M): - el = [] - - counter = count(0) - modes = (np.arange(2 * N) % N).tolist() - el = [post_select[i] if i in post_select else idx[next(counter)] for i in modes] - - el = np.array(el).reshape(2, -1) - el0 = el[0].tolist() - el1 = el[1].tolist() - - sf_idx = np.array(idx).reshape(2, -1) - sf_el = tuple(sf_idx[::-1].T.flatten()) - - rho[sf_el] = density_matrix_element(mu, cov, el0, el1, include_prefactor=False, hbar=hbar) - - rho *= pref - - if normalize: - # construct the standard 2D density matrix, and take the trace - new_ax = np.arange(2 * M).reshape([M, 2]).T.flatten() - tr = np.trace(rho.transpose(new_ax).reshape([cutoff ** M, cutoff ** M])).real - # renormalize - rho /= tr - - return rho - - -def pure_state_amplitude(mu, cov, i, include_prefactor=True, tol=1e-10, hbar=2, check_purity=True): - r"""Returns the :math:`\langle i | \psi\rangle` element of the state ket - of a Gaussian state defined by covariance matrix cov. - - - Args: - mu (array): length-:math:`2N` quadrature displacement vector - cov (array): length-:math:`2N` covariance matrix - i (list): list of amplitude elements - include_prefactor (bool): if ``True``, the prefactor is automatically calculated - used to scale the result. - tol (float): tolerance for determining if displacement is negligible - hbar (float): the value of :math:`\hbar` in the commutation - relation :math:`[\x,\p]=i\hbar`. - check_purity (bool): if ``True``, the purity of the Gaussian state is checked - before calculating the state vector. - - Returns: - complex: the pure state amplitude - """ - if check_purity: - if not is_pure_cov(cov, hbar=hbar, rtol=1e-05, atol=1e-08): - raise ValueError("The covariance matrix does not correspond to a pure state") - - rpt = i - beta = Beta(mu, hbar=hbar) - Q = Qmat(cov, hbar=hbar) - A = Amat(cov, hbar=hbar) - (n, _) = cov.shape - N = n // 2 - B = A[0:N, 0:N].conj() - alpha = beta[0:N] - - if np.linalg.norm(alpha) < tol: - # no displacement - if np.prod([k + 1 for k in rpt]) ** (1 / len(rpt)) < 3: - B_rpt = reduction(B, rpt) - haf = hafnian(B_rpt) - else: - haf = hafnian_repeated(B, rpt) - else: - 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(gamma, rpt)) - haf = hafnian(B_rpt, loop=True) - else: - haf = hafnian_repeated(B, rpt, mu=gamma, loop=True) - - if include_prefactor: - 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))) - - -def state_vector( - mu, cov, post_select=None, normalize=False, cutoff=5, hbar=2, check_purity=True, **kwargs -): - r"""Returns the state vector of a (PNR post-selected) Gaussian state. - - The resulting density matrix will have shape - - .. math:: \underbrace{D\times D \times \cdots \times D}_M - - where :math:`D` is the Fock space cutoff, and :math:`M` is the - number of *non* post-selected modes, i.e. ``M = len(mu)//2 - len(post_select)``. - - If post_select is None then the density matrix elements are calculated using - the multidimensional Hermite polynomials which provide a significantly faster - evaluation. - - - Args: - mu (array): length-:math:`2N` means vector in xp-ordering - cov (array): :math:`2N\times 2N` covariance matrix in xp-ordering - post_select (dict): dictionary containing the post-selected modes, of - the form ``{mode: value}``. - normalize (bool): If ``True``, a post-selected density matrix is re-normalized. - cutoff (dim): the final length (i.e., Hilbert space dimension) of each - mode in the density matrix. - hbar (float): the value of :math:`\hbar` in the commutation - relation :math:`[\x,\p]=i\hbar`. - check_purity (bool): if ``True``, the purity of the Gaussian state is checked - before calculating the state vector. - - Keyword Args: - choi_r (float or None): Value of the two-mode squeezing parameter used in Choi-Jamiolkoski - trick in :func:`~.fock_tensor`. This keyword argument should only be used when ``state_vector`` - is called by :func:`~.fock_tensor`. - - Returns: - np.array[complex]: the state vector of the Gaussian state - """ - if check_purity: - if not is_pure_cov(cov, hbar=hbar, rtol=1e-05, atol=1e-08): - raise ValueError("The covariance matrix does not correspond to a pure state") - - beta = Beta(mu, hbar=hbar) - A = Amat(cov, hbar=hbar) - Q = Qmat(cov, hbar=hbar) - - (n, _) = cov.shape - N = n // 2 - - B = A[0:N, 0:N] - alpha = beta[0:N] - gamma = np.conj(alpha) - B @ alpha - prefexp = -0.5 * (np.linalg.norm(alpha) ** 2 - alpha @ B @ alpha) - pref = np.exp(prefexp.conj()) - if post_select is None: - choi_r = kwargs.get("choi_r", None) - if choi_r is None: - denom = np.sqrt(np.sqrt(np.linalg.det(Q).real)) - else: - rescaling = np.concatenate( - [np.ones([N // 2]), (1.0 / np.tanh(choi_r)) * np.ones([N // 2])] - ) - B = np.diag(rescaling) @ B @ np.diag(rescaling) - gamma = rescaling * gamma - denom = np.sqrt(np.sqrt(np.linalg.det(Q / np.cosh(choi_r)).real)) - - psi = ( - np.real_if_close(pref) - * hafnian_batched(B.conj(), cutoff, mu=gamma.conj(), renorm=True) - / denom - ) - else: - M = N - len(post_select) - psi = np.zeros([cutoff] * (M), dtype=np.complex128) - - for idx in product(range(cutoff), repeat=M): - el = [] - - counter = count(0) - modes = (np.arange(N)).tolist() - el = [post_select[i] if i in post_select else idx[next(counter)] for i in modes] - psi[idx] = pure_state_amplitude( - mu, cov, el, check_purity=False, include_prefactor=False, hbar=hbar - ) - - psi = psi * pref - - if normalize: - norm = np.sqrt(np.sum(np.abs(psi) ** 2)) - psi = psi / norm - - return psi - -def mean_number_of_clicks(cov, hbar=2): - r""" Calculates the total mean number of clicks when a zero-mean gaussian state - is measured using threshold detectors. - - Args - cov (array): :math:`2N\times 2N` covariance matrix in xp-ordering - hbar (float): the value of :math:`\hbar` in the commutation relation :math:`[\x,\p]=i\hbar` - - Returns - float: mean number of clicks - """ - n, _ = cov.shape - nmodes = n // 2 - Q = Qmat(cov, hbar=hbar) - meanc = 1.0 * nmodes - - for i in range(nmodes): - det_val = np.real(Q[i, i] * Q[i + nmodes, i + nmodes] - Q[i + nmodes, i] * Q[i, i + nmodes]) - meanc -= 1.0 / np.sqrt(det_val) - return meanc - - -def variance_number_of_clicks(cov, hbar=2): - r""" Calculates the variance of the total number of clicks when a zero-mean gaussian state - is measured using threshold detectors. - - Args - cov (array): :math:`2N\times 2N` covariance matrix in xp-ordering - hbar (float): the value of :math:`\hbar` in the commutation relation :math:`[\x,\p]=i\hbar` - - Returns - float: variance in the total number of clicks - """ - n, _ = cov.shape - means = np.zeros([n]) - nmodes = n // 2 - Q = Qmat(cov, hbar=hbar) - vac_probs = np.array( - [ - np.real(Q[i, i] * Q[i + nmodes, i + nmodes] - Q[i + nmodes, i] * Q[i, i + nmodes]) - for i in range(nmodes) - ] - ) - vac_probs = np.sqrt(vac_probs) - vac_probs = 1 / vac_probs - term1 = np.sum(vac_probs * (1 - vac_probs)) - term2 = 0 - for i in range(nmodes): - for j in range(i): - _, Qij = reduced_gaussian(means, Q, [i, j]) - prob_vac_ij = np.linalg.det(Qij).real - prob_vac_ij = 1.0 / np.sqrt(prob_vac_ij) - term2 += prob_vac_ij - vac_probs[i] * vac_probs[j] - - return term1 + 2 * term2 - - -def mean_number_of_clicks_graph(A): - r""" Given an adjacency matrix this function calculates the mean number of clicks. - For this to make sense the user must provide a matrix with singular values - less than or equal to one. See Appendix A.3 of `_ - by Banchi et al. - - Args: - A (array): rescaled adjacency matrix - - Returns: - float: mean number of clicks - """ - n, _ = A.shape - idn = np.identity(n) - X = np.block([[0 * idn, idn], [idn, 0 * idn]]) - B = np.block([[A, 0 * A], [0 * A, np.conj(A)]]) - Q = np.linalg.inv(np.identity(2 * n) - X @ B) - return mean_number_of_clicks(Covmat(Q)) - - -def find_scaling_adjacency_matrix_torontonian(A, c_mean): - r""" Returns the scaling parameter by which the adjacency matrix A - should be rescaled so that the Gaussian state that encodes it has - give a mean number of clicks equal to ``c_mean`` when measured with - threshold detectors. - - Args: - A (array): adjacency matrix - c_mean (float): mean photon number of the Gaussian state - - Returns: - float: scaling parameter - """ - n, _ = A.shape - if c_mean < 0 or c_mean > n: - raise ValueError("The mean number of clicks should be smaller than the number of modes") - - vals = np.linalg.svd(A, compute_uv=False) - localA = A / vals[0] # rescale the matrix so that the singular values are between 0 and 1. - - def cost(x): - r""" Cost function giving the difference between the wanted number of clicks and the number - of clicks at a given scaling value. It assumes that the adjacency matrix has been rescaled - so that it has singular values between 0 and 1. - - Args: - x (float): scaling value - - Return: - float: difference between desired and obtained mean number of clicks - """ - if x >= 1.0: - return c_mean - n - if x <= 0: - return c_mean - return c_mean - mean_number_of_clicks_graph(x * localA) - - res = root_scalar(cost, x0=0.5, bracket=(0.0, 1.0)) # Do the optimization - - if not res.converged: - raise ValueError("The search for a scaling value failed") - return res.root / vals[0] - - -def find_scaling_adjacency_matrix(A, n_mean): - r""" Returns the scaling parameter by which the adjacency matrix A - should be rescaled so that the Gaussian state that endodes it has - a total mean photon number n_mean. - - Args: - A (array): Adjacency matrix - n_mean (float): Mean photon number of the Gaussian state - - Returns: - float: Scaling parameter - """ - eps = 1e-10 - ls = np.linalg.svd(A, compute_uv=False) - max_sv = ls[0] - a_lim = 0.0 - b_lim = 1.0 / (eps + max_sv) - x_init = 0.5 * b_lim - - if 1000 * eps >= max_sv: - raise ValueError("The singular values of the matrix A are too small.") - - def mean_photon_number(x, vals): - r""" Returns the mean number of photons in the Gaussian state that - encodes the adjacency matrix x*A where vals are the singular values of A - - Args: - x (float): Scaling parameter - vals (array): Singular values of the matrix A - - Returns: - n_mean: Mean photon number in the Gaussian state - """ - vals2 = (x * vals) ** 2 - n = np.sum(vals2 / (1.0 - vals2)) - return n - - # The following function is implicitly tested in test_find_scaling_adjacency_matrix - def grad_mean_photon_number(x, vals): # pragma: no cover - r""" Returns the gradient od the mean number of photons in the Gaussian state that - encodes the adjacency matrix x*A with respect to x. - vals are the singular values of A - - Args: - x (float): Scaling parameter - vals (array): Singular values of the matrix A - - Returns: - d_n_mean: Derivative of the mean photon number in the Gaussian state - with respect to x - """ - vals1 = vals * x - dn = (2.0 / x) * np.sum((vals1 / (1 - vals1 ** 2)) ** 2) - return dn - - f = lambda x: mean_photon_number(x, ls) - n_mean - df = lambda x: grad_mean_photon_number(x, ls) - res = root_scalar(f, fprime=df, x0=x_init, bracket=(a_lim, b_lim)) - - if not res.converged: - raise ValueError("The search for a scaling value failed") - - return res.root - - -def gen_Qmat_from_graph(A, n_mean): - r""" Returns the Qmat xp-covariance matrix associated to a graph with - adjacency matrix :math:`A` and with mean photon number :math:`n_{mean}`. - - Args: - A (array): a :math:`N\times N` ``np.float64`` (symmetric) adjacency matrix - n_mean (float): mean photon number of the Gaussian state - - Returns: - array: the :math:`2N\times 2N` Q matrix. - """ - n, m = A.shape - - if n != m: - raise ValueError("Matrix must be square.") - - sc = find_scaling_adjacency_matrix(A, n_mean) - Asc = sc * A - A = np.block([[Asc, 0 * Asc], [0 * Asc, Asc.conj()]]) - I = np.identity(2 * n) - X = Xmat(n) - Q = np.linalg.inv(I - X @ A) - return Q - -def photon_number_mean(mu, cov, j, hbar=2): - r""" Calculate the mean photon number of mode j of a Gaussian state. - - Args: - mu (array): vector of means of the Gaussian state using the ordering - :math:`[q_1, q_2, \dots, q_n, p_1, p_2, \dots, p_n]` - cov (array): the covariance matrix of the Gaussian state - j (int): the j :sup:`th` mode - hbar (float): the ``hbar`` convention used in the commutation - relation :math:`[q, p]=i\hbar` - - Returns: - float: the mean photon number in mode :math:`j`. - """ - num_modes = len(mu) // 2 - return ( - mu[j] ** 2 - + mu[j + num_modes] ** 2 - + cov[j, j] - + cov[j + num_modes, j + num_modes] - - hbar - ) / (2 * hbar) - -def photon_number_covar(mu, cov, j, k, hbar=2): - r""" Calculate the variance/covariance of the photon number distribution - of a Gaussian state. - - Implements the covariance matrix of the photon number distribution of a - Gaussian state according to the Last two eq. of Part II. in - `'Multidimensional Hermite polynomials and photon distribution for polymode - mixed light', Dodonov et al. `_ - - .. math:: - \sigma_{n_j n_j} &= \frac{1}{2}\left(T_j^2 - 2d_j - \frac{1}{2}\right) - + \left<\mathbf{Q}_j\right>\mathcal{M}_j\left<\mathbf{Q}_j\right>, \\ - \sigma_{n_j n_k} &= \frac{1}{2}\mathrm{Tr}\left(\Lambda_j \mathbf{M} \Lambda_k \mathbf{M}\right) - + \frac{1}{2}\left<\mathbf{Q}\right>\Lambda_j \mathbf{M} \Lambda_k\left<\mathbf{Q}\right>, - - where :math:`T_j` and :math:`d_j` are the trace and the determinant of - :math:`2 \times 2` matrix :math:`\mathcal{M}_j` whose elements coincide - with the nonzero elements of matrix :math:`\mathbf{M}_j = \Lambda_j \mathbf{M} \Lambda_k` - while the two-vector :math:`\mathbf{Q}_j` has the components :math:`(q_j, p_j)`. - :math:`2N \times 2N` projector matrix :math:`\Lambda_j` has only two nonzero - elements: :math:`\left(\Lambda_j\right)_{jj} = \left(\Lambda_j\right)_{j+N,j+N} = 1`. - Note that the convention for ``mu`` used here differs from the one used in Dodonov et al., - They both provide the same results in this particular case. - - Args: - mu (array): vector of means of the Gaussian state using the ordering - :math:`[q_1, q_2, \dots, q_n, p_1, p_2, \dots, p_n]` - cov (array): the covariance matrix of the Gaussian state - j (int): the j :sup:`th` mode - k (int): the k :sup:`th` mode - hbar (float): the ``hbar`` convention used in the commutation - relation :math:`[q, p]=i\hbar` - - Returns: - float: the covariance for the photon numbers at modes :math:`j` and :math:`k`. - """ - # renormalise the covariance matrix - cov = cov / hbar - - N = len(mu) // 2 - mu = np.array(mu) / np.sqrt(hbar) - - lambda_1 = np.zeros((2 * N, 2 * N)) - lambda_1[j, j] = lambda_1[j + N, j + N] = 1 - - lambda_2 = np.zeros((2 * N, 2 * N)) - lambda_2[k, k] = lambda_2[k + N, k + N] = 1 - - if j == k: - idxs = ((j, j, j + N, j + N), (j, j + N, j, j + N)) - M = (lambda_1 @ cov @ lambda_2)[idxs].reshape(2, 2) - - term_1 = (np.trace(M) ** 2 - 2 * np.linalg.det(M) - 0.5) / 2 - term_2 = mu[[j, j + N]] @ M @ mu[[j, j + N]] - else: - term_1 = np.trace(lambda_1 @ cov @ lambda_2 @ cov) / 2 - term_2 = (mu @ lambda_1 @ cov @ lambda_2 @ mu) / 2 - - return term_1 + term_2 - - -def photon_number_covmat(mu, cov, hbar=2): - r""" Calculate the covariance matrix of the photon number distribution of a - Gaussian state. - - Args: - mu (array): vector of means of the Gaussian state using the ordering - :math:`[q_1, q_2, \dots, q_n, p_1, p_2, \dots, p_n]` - cov (array): the covariance matrix of the Gaussian state - hbar (float): the ``hbar`` convention used in the commutation - relation :math:`[q, p]=i\hbar` - - Returns: - array: the covariance matrix of the photon number distribution - """ - N = len(mu) // 2 - pnd_cov = np.zeros((N, N)) - for i in range(N): - for j in range(i+1): - pnd_cov[i][j] = photon_number_covar(mu, cov, i, j, hbar=hbar) - pnd_cov[j][i] = pnd_cov[i][j] - return pnd_cov - - -def photon_number_mean_vector(mu, cov, hbar=2): - r""" Calculate the mean photon number of each of the modes in a Gaussian state - - Args: - mu (array): vector of means of the Gaussian state using the ordering - :math:`[q_1, q_2, \dots, q_n, p_1, p_2, \dots, p_n]` - cov (array): the covariance matrix of the Gaussian state - hbar (float): the ``hbar`` convention used in the commutation - relation :math:`[q, p]=i\hbar` - - Returns: - array: the vector of means of the photon number distribution - """ - - N = len(mu) // 2 - return np.array([photon_number_mean(mu, cov, j, hbar=hbar) for j in range(N)]) - -def is_valid_cov(cov, hbar=2, rtol=1e-05, atol=1e-08): - r""" Checks if the covariance matrix is a valid quantum covariance matrix. - - Args: - cov (array): a covariance matrix - hbar (float): value of hbar in the uncertainty relation - rtol (float): the relative tolerance parameter used in `np.allclose` - atol (float): the absolute tolerance parameter used in `np.allclose` - - Returns: - (bool): whether the given covariance matrix is a valid covariance matrix - """ - (n, m) = cov.shape - if n != m: - # raise ValueError("The input matrix must be square") - return False - if not np.allclose(cov, np.transpose(cov), rtol=rtol, atol=atol): - # raise ValueError("The input matrix is not symmetric") - return False - if n % 2 != 0: - # raise ValueError("The input matrix is of even dimension") - return False - - nmodes = n // 2 - vals = np.linalg.eigvalsh(cov + 0.5j * hbar * sympmat(nmodes)) - vals[np.abs(vals) < atol] = 0.0 - if np.all(vals >= 0): - # raise ValueError("The input matrix violates the uncertainty relation") - return True - - return False - - -def is_pure_cov(cov, hbar=2, rtol=1e-05, atol=1e-08): - r""" Checks if the covariance matrix is a valid quantum covariance matrix - that corresponds to a quantum pure state - - Args: - cov (array): a covariance matrix - hbar (float): value of hbar in the uncertainty relation - rtol (float): the relative tolerance parameter used in `np.allclose` - atol (float): the absolute tolerance parameter used in `np.allclose` - - Returns: - (bool): whether the given covariance matrix corresponds to a pure state - """ - if is_valid_cov(cov, hbar=hbar, rtol=rtol, atol=atol): - purity = 1 / np.sqrt(np.linalg.det(2 * cov / hbar)) - if np.allclose(purity, 1.0, rtol=rtol, atol=atol): - return True - - return False - - -def is_classical_cov(cov, hbar=2, atol=1e-08): - r""" Checks if the covariance matrix can be efficiently sampled. - - Args: - cov (array): a covariance matrix - hbar (float): value of hbar in the uncertainty relation - atol (float): the absolute tolerance parameter used in `np.allclose` - - Returns: - (bool): whether the given covariance matrix corresponds to a classical state - """ - - if is_valid_cov(cov, hbar=hbar, atol=atol): - (n, _) = cov.shape - vals = np.linalg.eigvalsh(cov - 0.5 * hbar * np.identity(n)) - vals[np.abs(vals) < atol] = 0.0 - - if np.all(vals >= 0): - return True - return False - - -def find_classical_subsystem(cov, hbar=2, atol=1e-08): - """Find the largest integer ``k`` so that subsystem in modes ``[0,1,...,k-1]`` is a classical state. - - - Args: - cov (array): a covariance matrix - hbar (float): value of hbar in the uncertainty relation - atol (float): the absolute tolerance parameter used when determining if the state is classical - - Returns: - int: the largest k so that modes ``[0,1,...,k-1]`` are in a classical state. - """ - n, _ = cov.shape - nmodes = n // 2 - if is_classical_cov(cov, hbar=hbar, atol=atol): - return nmodes - k = 0 - mu = np.zeros(n) - is_classical = True - while is_classical: - _, Vk = reduced_gaussian(mu, cov, list(range(k + 1))) - is_classical = is_classical_cov(Vk, hbar=hbar, atol=atol) - k += 1 - return k - 1 - -def gen_single_mode_dist(s, cutoff=50, N=1): - """Generate the photon number distribution of :math:`N` identical single mode squeezed states. - - Args: - s (float): squeezing parameter - cutoff (int): Fock cutoff - N (float): number of squeezed states - - Returns: - (array): Photon number distribution - """ - r = 0.5 * N - q = 1.0 - np.tanh(s) ** 2 - N = cutoff // 2 - ps_tot = np.zeros(cutoff) - if cutoff % 2 == 0: - ps = nbinom.pmf(np.arange(N), p=q, n=r) - ps_tot[0::2] = ps - else: - ps = nbinom.pmf(np.arange(N + 1), p=q, n=r) - ps_tot[0:-1][0::2] = ps[0:-1] - ps_tot[-1] = ps[-1] - - return ps_tot - - -def gen_multi_mode_dist(s, cutoff=50, padding_factor=2): - """Generates the total photon number distribution of single mode squeezed states with different squeezing values. - - Args: - s (array): array of squeezing parameters - cutoff (int): Fock cutoff - - Returns: - (array[int]): total photon number distribution - """ - scale = padding_factor - cutoff_sc = scale * cutoff - ps = np.zeros(cutoff_sc) - ps[0] = 1.0 - for s_val in s: - ps = np.convolve(ps, gen_single_mode_dist(s_val, cutoff_sc))[0:cutoff_sc] - return ps - - -def total_photon_num_dist_pure_state(cov, cutoff=50, hbar=2, padding_factor=2): - r""" Calculates the total photon number distribution of a pure state - with zero mean. - - Args: - cov (array): :math:`2N\times 2N` covariance matrix in xp-ordering - cutoff (int): Fock cutoff - tol (float): tolerance for determining if displacement is negligible - hbar (float): the value of :math:`\hbar` in the commutation - padding_factor (int): expanded size of the photon distribution to avoid accumulation of errors - - Returns: - (array): Total photon number distribution - """ - if is_pure_cov(cov): - A = Amat(cov, hbar=hbar) - (n, _) = A.shape - N = n // 2 - B = A[0:N, 0:N] - rs = np.arctanh(np.linalg.svd(B, compute_uv=False)) - return gen_multi_mode_dist(rs, cutoff=cutoff, padding_factor=padding_factor)[0:cutoff] - raise ValueError("The Gaussian state is not pure") - - -def fock_tensor( - S, - alpha, - cutoff, - choi_r=np.arcsinh(1.0), - check_symplectic=True, - sf_order=False, - rtol=1e-05, - atol=1e-08, -): - r""" - Calculates the Fock representation of a Gaussian unitary parametrized by - the symplectic matrix S and the displacements alpha up to cutoff in Fock space. - - Args: - S (array): symplectic matrix - alpha (array): complex vector of displacements - cutoff (int): cutoff in Fock space - choi_r (float): squeezing parameter used for the Choi expansion - check_symplectic (boolean): checks whether the input matrix is symplectic - sf_order (boolean): reshapes the tensor so that it follows the sf ordering of indices - rtol (float): the relative tolerance parameter used in `np.allclose` - atol (float): the absolute tolerance parameter used in `np.allclose` - - Return: - (array): Tensor containing the Fock representation of the Gaussian unitary - """ - # Check the matrix is symplectic - if check_symplectic: - if not is_symplectic(S, rtol=rtol, atol=atol): - raise ValueError("The matrix S is not symplectic") - - # And that S and alpha have compatible dimensions - m, _ = S.shape - l = m // 2 - if l != len(alpha): - raise ValueError( - "The matrix S and the vector alpha do not have compatible dimensions" - ) - # Check if S corresponds to an interferometer, if so use optimized routines - if np.allclose(S @ S.T, np.identity(m), rtol=rtol, atol=atol) and np.allclose( - alpha, 0, rtol=rtol, atol=atol - ): - reU = S[:l, :l] - imU = S[:l, l:] - if np.allclose(imU, 0, rtol=rtol, atol=atol): - Ub = np.block([[0 * reU, -reU], [-reU.T, 0 * reU]]) - tensor = interferometer_real(Ub, cutoff) - else: - U = reU - 1j * imU - Ub = np.block([[0 * U, -U], [-U.T, 0 * U]]) - tensor = interferometer(Ub, cutoff) - else: - # Construct the covariance matrix of l two-mode squeezed vacua pairing modes i and i+l - ch = np.cosh(choi_r) * np.identity(l) - sh = np.sinh(choi_r) * np.identity(l) - zh = np.zeros([l, l]) - Schoi = np.block( - [[ch, sh, zh, zh], [sh, ch, zh, zh], [zh, zh, ch, -sh], [zh, zh, -sh, ch]] - ) - # And then its Choi expanded symplectic - S_exp = expand(S, list(range(l)), 2 * l) @ Schoi - # And this is the corresponding covariance matrix - cov = S_exp @ S_exp.T - alphat = np.array(list(alpha) + ([0] * l)) - x = 2 * alphat.real - p = 2 * alphat.imag - mu = np.concatenate([x, p]) - - tensor = state_vector( - mu, - cov, - normalize=False, - cutoff=cutoff, - hbar=2, - check_purity=False, - choi_r=choi_r, - ) - - if sf_order: - sf_indexing = tuple(chain.from_iterable([[i, i + l] for i in range(l)])) - return tensor.transpose(sf_indexing) - - return tensor - - - -def probabilities(mu, cov, cutoff, parallel=False, hbar=2.0, rtol=1e-05, atol=1e-08): - r"""Generate the Fock space probabilities of a Gaussian state up to a Fock space cutoff. - - .. note:: - - Individual density matrix elements are computed using multithreading by OpenMP. - Setting ``parallel=True`` will further result in *multiple* density matrix elements - being computed in parallel. - - When setting ``parallel=True``, OpenMP will need to be turned off by setting the - environment variable ``OMP_NUM_THREADS=1`` (forcing single threaded use for individual - matrix elements). Remove the environment variable or set it to ``OMP_NUM_THREADS=''`` - to again use multithreading with OpenMP. - - Args: - mu (array): vector of means of length ``2*n_modes`` - cov (array): covariance matrix of shape ``[2*n_modes, 2*n_modes]`` - cutoff (int): cutoff in Fock space - parallel (bool): if ``True``, uses ``dask`` for parallelization instead of OpenMP - hbar (float): value of :math:`\hbar` in the commutation relation :math;`[\hat{x}, \hat{p}]=i\hbar` - rtol (float): the relative tolerance parameter used in ``np.allclose`` - atol (float): the absolute tolerance parameter used in ``np.allclose`` - - Returns: - (array): Fock space probabilities up to cutoff. The shape of this tensor is ``[cutoff]*num_modes``. - """ - if is_pure_cov(cov, hbar=hbar, rtol=rtol, atol=atol): # Check if the covariance matrix cov is pure - return np.abs(state_vector(mu, cov, cutoff=cutoff, hbar=hbar, check_purity=False)) ** 2 - num_modes = len(mu) // 2 - - if parallel: - compute_list = [] - # create a list of parallelizable computations - for i in product(range(cutoff), repeat=num_modes): - compute_list.append(dask.delayed(density_matrix_element)(mu, cov, i, i, hbar=hbar)) - - probs = np.maximum( - 0.0, np.real_if_close(dask.compute(*compute_list, scheduler="processes")) - ).reshape([cutoff] * num_modes) - # maximum is needed because sometimes a probability is very close to zero from below - else: - probs = np.zeros([cutoff] * num_modes) - for i in product(range(cutoff), repeat=num_modes): - probs[i] = np.maximum( - 0.0, np.real_if_close(density_matrix_element(mu, cov, i, i, hbar=hbar)) - ) - # maximum is needed because sometimes a probability is very close to zero from below - return probs - - - -@jit(nopython=True) -def loss_mat(eta, cutoff): # pragma: no cover - r"""Constructs a binomial loss matrix with transmission eta up to n photons. - - Args: - eta (float): Transmission coefficient. ``eta=0.0`` corresponds to complete loss and ``eta=1.0`` corresponds to no loss. - cutoff (int): cutoff in Fock space. - - Returns: - array: :math:`n\times n` matrix representing the loss. - """ - # If full transmission return the identity - - if eta < 0.0 or eta > 1.0: - raise ValueError("The transmission parameter eta should be a number between 0 and 1.") - - if eta == 1.0: - return np.identity(cutoff) - - # Otherwise construct the matrix elements recursively - lm = np.zeros((cutoff, cutoff)) - mu = 1.0 - eta - lm[:, 0] = mu ** (np.arange(cutoff)) - for i in range(cutoff): - for j in range(1, i + 1): - lm[i, j] = lm[i, j - 1] * (eta / mu) * (i - j + 1) / (j) - return lm - - -def update_probabilities_with_loss(etas, probs): - """Given a list of transmissivities a tensor of probabilitites, calculate - an updated tensor of probabilities after loss is applied. - - Args: - etas (list): List of transmissitivities describing the loss in each of the modes - probs (array): Array of probabilitites in the different modes - - Returns: - array: List of loss-updated probabilities with the same shape as probs. - """ - - probs_shape = probs.shape - if len(probs_shape) != len(etas): - raise ValueError("The list of transmission etas and the tensor of probabilities probs have incompatible dimensions.") - - alphabet = "abcdefghijklmnopqrstuvwxyz" - cutoff = probs_shape[0] - for i, eta in enumerate(etas): - einstrings = "ij,{}i...->{}j...".format(alphabet[:i], alphabet[:i]) - - qein = np.zeros_like(probs) - qein = np.einsum(einstrings, loss_mat(eta, cutoff), probs) - probs = np.copy(qein) - return qein - -@jit(nopython=True) -def _update_1d(probs, one_d, cutoff): # pragma: no cover - """ Performs a convolution of the two arrays. The first one does not need to be one dimensional, which is why we do not use ``np.convolve``. - - Args: - probs (array): (multidimensional) array - one_d (array): one dimensional array - cutoff (int): cutoff in Fock space for the first array - - Returns: - (array): the convolution of the two arrays, with the same shape as ``probs``. - """ - new_d = np.zeros_like(probs) - for i in range(cutoff): - for j in range(min(i + 1, len(one_d))): - new_d[i] += probs[i - j] * one_d[j] - return new_d - -def update_probabilities_with_noise(probs_noise, probs): - """Given a list of noise probability distributions for each of the modes and a tensor of - probabilitites, calculate an updated tensor of probabilities after noise is applied. - - Args: - probs_noise (list): List of probability distributions describing the noise in each of the modes - probs (array): Array of probabilitites in the different modes - - Returns: - array: List of noise-updated probabilities with the same shape as probs. - """ - probs_shape = probs.shape - num_modes = len(probs_shape) - cutoff = probs_shape[0] - if num_modes != len(probs_noise): - raise ValueError( - "The list of probability distributions probs_noise and the tensor of probabilities probs have incompatible dimensions." - ) - - for k in range(num_modes): #update one mode at a time - perm = np.arange(num_modes) - perm[0] = k - perm[k] = 0 - one_d = probs_noise[k] - probs_masked = np.transpose(probs, axes=perm) - probs_masked = _update_1d(probs_masked, one_d, cutoff) - probs = np.transpose(probs_masked, axes=perm) - return probs - - -def fidelity(mu1, cov1, mu2, cov2, hbar=2, rtol=1e-05, atol=1e-08): - """Calculates the fidelity between two Gaussian quantum states. - - Note that if the covariance matrices correspond to pure states this - function reduces to the modulus square of the overlap of their state vectors. - For the derivation see `'Quantum Fidelity for Arbitrary Gaussian States', Banchi et al. <10.1103/PhysRevLett.115.260501>`_. - - Args: - mu1 (array): vector of means of the first state - cov1 (array): covariance matrix of the first state - mu2 (array): vector of means of the second state - cov2 (array): covariance matrix of the second state - hbar (float): value of hbar in the uncertainty relation - rtol (float): the relative tolerance parameter used in `np.allclose` - atol (float): the absolute tolerance parameter used in `np.allclose` - - Returns: - (float): value of the fidelity between the two states - """ - - n0, n1 = cov1.shape - m0, m1 = cov2.shape - (l0,) = mu1.shape - (l1,) = mu1.shape - if not n0 == n1 == m0 == m1 == l0 == l1: - raise ValueError("The inputs have incompatible shapes") - - v1 = cov1 / hbar - v2 = cov2 / hbar - deltar = (mu1 - mu2) / np.sqrt(hbar / 2) - n0, n1 = cov1.shape - n = n0 // 2 - W = sympmat(n) - - si12 = np.linalg.inv(v1 + v2) - vaux = W.T @ si12 @ (0.25 * W + v2 @ W @ v1) - p1 = vaux @ W - p1 = p1 @ p1 - p1 = np.identity(2 * n) + 0.25 * np.linalg.inv(p1) - if np.allclose(p1, 0, rtol=rtol, atol=atol): - p1 = np.zeros_like(p1) - else: - p1 = sqrtm(p1) - p1 = 2 * (p1 + np.identity(2 * n)) - p1 = p1 @ vaux - f = np.sqrt(np.linalg.det(si12) * np.linalg.det(p1)) * np.exp( - -0.25 * deltar @ si12 @ deltar - ) - return f - -def normal_ordered_expectation(mu, cov, rpt, hbar=2): - r"""Calculates the expectation value of the normal ordered product - :math:`\prod_{i=0}^{N-1} a_i^{\dagger n_i} \prod_{j=0}^{N-1} a_j^{m_j}` with respect to an N-mode Gaussian state, - where :math:`\text{rpt}=(n_0, n_1, \ldots, n_{N-1}, m_0, m_1, \ldots, m_{N-1})`. - - Args: - mu (array): length-:math:`2N` means vector in xp-ordering. - cov (array): :math:`2N\times 2N` covariance matrix in xp-ordering. - rpt (list): integers specifying the terms to calculate. - hbar (float): value of hbar in the uncertainty relation. - - Returns: - (float): expectation value of the normal ordered product of operators - """ - alpha = Beta(mu, hbar=hbar) - n = len(cov) - V = (Qmat(cov, hbar=hbar) - np.identity(n)) @ Xmat(n // 2) - A = reduction(V, rpt) - if np.allclose(mu, 0): - res = np.conj(hafnian(A)) - else: - np.fill_diagonal(A, reduction(np.conj(alpha), rpt)) - res = np.conj(hafnian(A, loop=True)) - return np.conj(res) - -def photon_number_expectation(mu, cov, modes, hbar=2): - r"""Calculates the expectation value of the product of the number operator of the modes in a Gaussian state. - - Args: - mu (array): length-:math:`2N` means vector in xp-ordering. - cov (array): :math:`2N\times 2N` covariance matrix in xp-ordering. - modes (list): list of modes - hbar (float): value of hbar in the uncertainty relation. - - Returns: - (float): expectation value of the product of the number operators of the modes. - """ - n, _ = cov.shape - n_modes = n // 2 - rpt = np.zeros([n], dtype=int) - for i in modes: - rpt[i] = 1 - rpt[i + n_modes] = 1 - - return normal_ordered_expectation(mu, cov, rpt, hbar=hbar) - - -def photon_number_squared_expectation(mu, cov, modes, hbar=2): - r"""Calculates the expectation value of the square of the product of the number operator of the modes in - a Gaussian state. - - Args: - mu (array): length-:math:`2N` means vector in xp-ordering. - cov (array): :math:`2N\times 2N` covariance matrix in xp-ordering. - modes (list): list of modes - hbar (float): value of hbar in the uncertainty relation. - - Returns: - (float): expectation value of the square of the product of the number operator of the modes. - """ - n_modes = len(modes) - - mu_red, cov_red = reduced_gaussian(mu, cov, modes) - result = 0 - for item in product([1, 2], repeat=n_modes): - rpt = item + item - term = normal_ordered_expectation(mu_red, cov_red, rpt, hbar=hbar) - result += term - return result From c50d73c78271d6ac584a15f89b8d762c7d386ef3 Mon Sep 17 00:00:00 2001 From: Theodor Isacsson Date: Fri, 4 Sep 2020 18:10:39 -0400 Subject: [PATCH 04/17] Fix import --- thewalrus/quantum/means_and_variances.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/thewalrus/quantum/means_and_variances.py b/thewalrus/quantum/means_and_variances.py index 4ed183b3a..889f3ecf9 100644 --- a/thewalrus/quantum/means_and_variances.py +++ b/thewalrus/quantum/means_and_variances.py @@ -20,7 +20,7 @@ import numpy as np -from .covariance_matrices import normal_ordered_expectation +from .covariance_matrices import normal_ordered_expectation, reduced_gaussian ################################################################################ From fb25fd2b7230217b000a1beb6b8652a41c908c13 Mon Sep 17 00:00:00 2001 From: Theodor Isacsson Date: Fri, 4 Sep 2020 18:11:59 -0400 Subject: [PATCH 05/17] Fix pylint issues --- thewalrus/quantum/__init__.py | 2 +- thewalrus/quantum/fock_states_and_tensors.py | 1 + 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/thewalrus/quantum/__init__.py b/thewalrus/quantum/__init__.py index d3c9fa1e4..d3909152c 100644 --- a/thewalrus/quantum/__init__.py +++ b/thewalrus/quantum/__init__.py @@ -193,4 +193,4 @@ "update_probabilities_with_noise", "normal_ordered_expectation", "fidelity", -] \ No newline at end of file +] diff --git a/thewalrus/quantum/fock_states_and_tensors.py b/thewalrus/quantum/fock_states_and_tensors.py index cbaabf23b..08b67ad72 100644 --- a/thewalrus/quantum/fock_states_and_tensors.py +++ b/thewalrus/quantum/fock_states_and_tensors.py @@ -15,6 +15,7 @@ Set of functions for calculating various state representations, probabilites and fidelities of Gaussian states. """ +# pylint: disable=too-many-arguments from itertools import count, product, chain From ca8ba325d3714778b15ef9a013709e092a7b0016 Mon Sep 17 00:00:00 2001 From: Theodor Isacsson Date: Fri, 4 Sep 2020 18:26:38 -0400 Subject: [PATCH 06/17] Add import --- thewalrus/__init__.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/thewalrus/__init__.py b/thewalrus/__init__.py index 55dc957c6..6e655d603 100644 --- a/thewalrus/__init__.py +++ b/thewalrus/__init__.py @@ -89,6 +89,8 @@ import numpy as np +import thewalrus.quantum + from ._hafnian import ( haf_complex, haf_int, From 3e381f1b86de2531310c8d217495d2c7154ee8f1 Mon Sep 17 00:00:00 2001 From: Theodor Isacsson Date: Wed, 9 Sep 2020 10:16:39 -0400 Subject: [PATCH 07/17] Add find_packages to setup --- setup.py | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/setup.py b/setup.py index 8510d98a3..24f8688f9 100644 --- a/setup.py +++ b/setup.py @@ -17,7 +17,7 @@ import os import platform -import setuptools +from setuptools import find_packages with open("thewalrus/_version.py") as f: @@ -137,10 +137,7 @@ def cythonize(x, compile_time_env=None): 'maintainer_email': 'nicolas@xanadu.ai', 'url': 'https://github.com/XanaduAI/thewalrus', 'license': 'Apache License 2.0', - 'packages': [ - 'thewalrus', - 'thewalrus.tests' - ], + 'packages': find_packages(where="."), 'description': 'Open source library for hafnian calculation', 'long_description': open('README.rst').read(), 'provides': ["thewalrus"], From 3a0eebd6bc585da406190503a030693407347085 Mon Sep 17 00:00:00 2001 From: Theodor Date: Wed, 9 Sep 2020 15:50:59 -0400 Subject: [PATCH 08/17] Apply suggestions from code review Co-authored-by: Nicolas Quesada --- thewalrus/quantum/__init__.py | 1 - thewalrus/quantum/photon_number_distributions.py | 2 +- 2 files changed, 1 insertion(+), 2 deletions(-) diff --git a/thewalrus/quantum/__init__.py b/thewalrus/quantum/__init__.py index d3909152c..a4e36314a 100644 --- a/thewalrus/quantum/__init__.py +++ b/thewalrus/quantum/__init__.py @@ -119,7 +119,6 @@ total_photon_num_dist_pure_state gen_single_mode_dist gen_multi_mode_dist - normal_ordered_complex_cov Details diff --git a/thewalrus/quantum/photon_number_distributions.py b/thewalrus/quantum/photon_number_distributions.py index 526bd48bf..1b158277e 100644 --- a/thewalrus/quantum/photon_number_distributions.py +++ b/thewalrus/quantum/photon_number_distributions.py @@ -81,7 +81,7 @@ def gen_multi_mode_dist(s, cutoff=50, padding_factor=2): Args: s (array): array of squeezing parameters cutoff (int): Fock cutoff - + padding_factor (int): expanded size of the photon distribution to avoid accumulation of errors Returns: (array[int]): total photon number distribution """ From 40bcb4a593d60b7541a4f749eaa5a6834b97ebad Mon Sep 17 00:00:00 2001 From: Theodor Isacsson Date: Wed, 9 Sep 2020 17:27:26 -0400 Subject: [PATCH 09/17] Finalise refactor --- thewalrus/quantum/__init__.py | 57 ++++--- thewalrus/quantum/adjacency_matrices.py | 88 +---------- ...{covariance_matrices.py => conversions.py} | 125 +-------------- ..._states_and_tensors.py => fock_tensors.py} | 83 ++-------- thewalrus/quantum/gaussian_state_tests.py | 149 ++++++++++++++++++ thewalrus/quantum/means_and_variances.py | 117 +++++++++++++- .../quantum/photon_number_distributions.py | 17 +- thewalrus/tests/test_quantum.py | 64 ++++---- 8 files changed, 356 insertions(+), 344 deletions(-) rename thewalrus/quantum/{covariance_matrices.py => conversions.py} (56%) rename thewalrus/quantum/{fock_states_and_tensors.py => fock_tensors.py} (90%) create mode 100644 thewalrus/quantum/gaussian_state_tests.py diff --git a/thewalrus/quantum/__init__.py b/thewalrus/quantum/__init__.py index a4e36314a..439c41bd9 100644 --- a/thewalrus/quantum/__init__.py +++ b/thewalrus/quantum/__init__.py @@ -100,9 +100,8 @@ Qmat Covmat Amat - Beta - Means - prefactor + complex_to_real_displacements + real_to_complex_displacements find_scaling_adjacency_matrix mean_number_of_clicks mean_number_of_clicks_graph @@ -116,54 +115,49 @@ is_pure_cov is_classical_cov find_classical_subsystem - total_photon_num_dist_pure_state - gen_single_mode_dist - gen_multi_mode_dist + pure_state_distribution Details ^^^^^^^ """ -from .fock_states_and_tensors import ( +from .fock_tensors import ( pure_state_amplitude, state_vector, density_matrix_element, density_matrix, fock_tensor, probabilities, + loss_mat, update_probabilities_with_loss, update_probabilities_with_noise, - fidelity, - find_classical_subsystem, - prefactor, - loss_mat, + _prefactor, ) from .adjacency_matrices import ( - mean_number_of_clicks, - mean_number_of_clicks_graph, - variance_number_of_clicks, find_scaling_adjacency_matrix, find_scaling_adjacency_matrix_torontonian, gen_Qmat_from_graph, ) -from .covariance_matrices import ( +from .gaussian_state_tests import ( + is_valid_cov, + is_pure_cov, + is_classical_cov, + fidelity, +) + +from .conversions import ( reduced_gaussian, Xmat, Qmat, Covmat, Amat, - normal_ordered_expectation, - Beta, - Means, - - is_valid_cov, - is_pure_cov, - is_classical_cov + complex_to_real_displacements, + real_to_complex_displacements, ) from .means_and_variances import ( @@ -173,12 +167,16 @@ photon_number_covmat, photon_number_expectation, photon_number_squared_expectation, + normal_ordered_expectation, + mean_number_of_clicks, + variance_number_of_clicks, + mean_number_of_clicks_graph, ) from .photon_number_distributions import ( - total_photon_num_dist_pure_state, - gen_single_mode_dist, - gen_multi_mode_dist, + pure_state_distribution, + _squeezed_state_distribution, + _convolve_squeezed_state_distribution, ) __all__ = [ @@ -193,3 +191,12 @@ "normal_ordered_expectation", "fidelity", ] + + +# old names for functions; remove in due time +Means = real_to_complex_displacements +Beta = complex_to_real_displacements +prefactor = _prefactor +gen_multi_mode_dist = _convolve_squeezed_state_distribution +gen_single_mode_dist = _squeezed_state_distribution +total_photon_num_dist_pure_state = pure_state_distribution diff --git a/thewalrus/quantum/adjacency_matrices.py b/thewalrus/quantum/adjacency_matrices.py index 882895415..132bce2a3 100644 --- a/thewalrus/quantum/adjacency_matrices.py +++ b/thewalrus/quantum/adjacency_matrices.py @@ -12,95 +12,15 @@ # See the License for the specific language governing permissions and # limitations under the License. """ -Functions for rescaling adjacency matrices as well as for calculating the mean -number of clicks for an adjacency matrix. +Functions for rescaling adjacency matrices as well as calculating the Qmat +covariance matrix from an adjacency matrix. """ import numpy as np from scipy.optimize import root_scalar -from .covariance_matrices import Xmat, Covmat, Qmat, reduced_gaussian - -def mean_number_of_clicks(cov, hbar=2): - r""" Calculates the total mean number of clicks when a zero-mean gaussian state - is measured using threshold detectors. - - Args - cov (array): :math:`2N\times 2N` covariance matrix in xp-ordering - hbar (float): the value of :math:`\hbar` in the commutation relation :math:`[\x,\p]=i\hbar` - - Returns - float: mean number of clicks - """ - n, _ = cov.shape - nmodes = n // 2 - Q = Qmat(cov, hbar=hbar) - meanc = 1.0 * nmodes - - for i in range(nmodes): - det_val = np.real(Q[i, i] * Q[i + nmodes, i + nmodes] - Q[i + nmodes, i] * Q[i, i + nmodes]) - meanc -= 1.0 / np.sqrt(det_val) - return meanc - - -def variance_number_of_clicks(cov, hbar=2): - r""" Calculates the variance of the total number of clicks when a zero-mean gaussian state - is measured using threshold detectors. - - Args - cov (array): :math:`2N\times 2N` covariance matrix in xp-ordering - hbar (float): the value of :math:`\hbar` in the commutation relation :math:`[\x,\p]=i\hbar` - - Returns - float: variance in the total number of clicks - """ - n, _ = cov.shape - means = np.zeros([n]) - nmodes = n // 2 - Q = Qmat(cov, hbar=hbar) - vac_probs = np.array( - [ - np.real(Q[i, i] * Q[i + nmodes, i + nmodes] - Q[i + nmodes, i] * Q[i, i + nmodes]) - for i in range(nmodes) - ] - ) - vac_probs = np.sqrt(vac_probs) - vac_probs = 1 / vac_probs - term1 = np.sum(vac_probs * (1 - vac_probs)) - term2 = 0 - for i in range(nmodes): - for j in range(i): - _, Qij = reduced_gaussian(means, Q, [i, j]) - prob_vac_ij = np.linalg.det(Qij).real - prob_vac_ij = 1.0 / np.sqrt(prob_vac_ij) - term2 += prob_vac_ij - vac_probs[i] * vac_probs[j] - - return term1 + 2 * term2 - - -def mean_number_of_clicks_graph(A): - r""" Given an adjacency matrix this function calculates the mean number of clicks. - For this to make sense the user must provide a matrix with singular values - less than or equal to one. See Appendix A.3 of `_ - by Banchi et al. - - Args: - A (array): rescaled adjacency matrix - - Returns: - float: mean number of clicks - """ - n, _ = A.shape - idn = np.identity(n) - X = np.block([[0 * idn, idn], [idn, 0 * idn]]) - B = np.block([[A, 0 * A], [0 * A, np.conj(A)]]) - Q = np.linalg.inv(np.identity(2 * n) - X @ B) - return mean_number_of_clicks(Covmat(Q)) - - -################################################################################ -# Rescale adjacency matrices -################################################################################ +from .conversions import Xmat +from .means_and_variances import mean_number_of_clicks_graph def find_scaling_adjacency_matrix_torontonian(A, c_mean): diff --git a/thewalrus/quantum/covariance_matrices.py b/thewalrus/quantum/conversions.py similarity index 56% rename from thewalrus/quantum/covariance_matrices.py rename to thewalrus/quantum/conversions.py index 917b5d637..ca2d8e8a3 100644 --- a/thewalrus/quantum/covariance_matrices.py +++ b/thewalrus/quantum/conversions.py @@ -13,19 +13,11 @@ # limitations under the License. """ Functions for transforming one type of covariance-matrix-like object into -another as well as various property tests for covariance matrices. +another as well as reducing a gaussian state. """ import numpy as np -from ..symplectic import sympmat -from .._hafnian import hafnian, reduction - - -################################################################################ -# Construct the reduced means and cov of a Gaussian state -################################################################################ - def reduced_gaussian(mu, cov, modes): r""" Returns the vector of means and the covariance matrix of the specified modes. @@ -60,11 +52,6 @@ def reduced_gaussian(mu, cov, modes): return mu[ind], cov[rows, cols] -################################################################################ -# Transform one type of covariance-matrix-like object into another -################################################################################ - - def Xmat(N): r"""Returns the matrix :math:`X_n = \begin{bmatrix}0 & I_n\\ I_n & 0\end{bmatrix}` @@ -163,33 +150,7 @@ def Amat(cov, hbar=2, cov_is_qmat=False): return A -def normal_ordered_expectation(mu, cov, rpt, hbar=2): - r"""Calculates the expectation value of the normal ordered product - :math:`\prod_{i=0}^{N-1} a_i^{\dagger n_i} \prod_{j=0}^{N-1} a_j^{m_j}` with respect to an N-mode Gaussian state, - where :math:`\text{rpt}=(n_0, n_1, \ldots, n_{N-1}, m_0, m_1, \ldots, m_{N-1})`. - - Args: - mu (array): length-:math:`2N` means vector in xp-ordering. - cov (array): :math:`2N\times 2N` covariance matrix in xp-ordering. - rpt (list): integers specifying the terms to calculate. - hbar (float): value of hbar in the uncertainty relation. - - Returns: - (float): expectation value of the normal ordered product of operators - """ - alpha = Beta(mu, hbar=hbar) - n = len(cov) - V = (Qmat(cov, hbar=hbar) - np.identity(n)) @ Xmat(n // 2) - A = reduction(V, rpt) - if np.allclose(mu, 0): - res = np.conj(hafnian(A)) - else: - np.fill_diagonal(A, reduction(np.conj(alpha), rpt)) - res = np.conj(hafnian(A, loop=True)) - return np.conj(res) - - -def Beta(mu, hbar=2): +def complex_to_real_displacements(mu, hbar=2): r"""Returns the vector of complex displacements and conjugate displacements. Args: @@ -208,7 +169,7 @@ def Beta(mu, hbar=2): return np.concatenate([alpha, alpha.conj()]) -def Means(beta, hbar=2): +def real_to_complex_displacements(beta, hbar=2): r"""Returns the vector of real quadrature displacements. Args: @@ -224,83 +185,3 @@ def Means(beta, hbar=2): N = len(beta) // 2 alpha = beta[0:N] return np.sqrt(2 * hbar) * np.concatenate([alpha.real, alpha.imag]) - -################################################################################ -# Test properties of covariance matrices -################################################################################ - - -def is_valid_cov(cov, hbar=2, rtol=1e-05, atol=1e-08): - r""" Checks if the covariance matrix is a valid quantum covariance matrix. - - Args: - cov (array): a covariance matrix - hbar (float): value of hbar in the uncertainty relation - rtol (float): the relative tolerance parameter used in `np.allclose` - atol (float): the absolute tolerance parameter used in `np.allclose` - - Returns: - (bool): whether the given covariance matrix is a valid covariance matrix - """ - (n, m) = cov.shape - if n != m: - # raise ValueError("The input matrix must be square") - return False - if not np.allclose(cov, np.transpose(cov), rtol=rtol, atol=atol): - # raise ValueError("The input matrix is not symmetric") - return False - if n % 2 != 0: - # raise ValueError("The input matrix is of even dimension") - return False - - nmodes = n // 2 - vals = np.linalg.eigvalsh(cov + 0.5j * hbar * sympmat(nmodes)) - vals[np.abs(vals) < atol] = 0.0 - if np.all(vals >= 0): - # raise ValueError("The input matrix violates the uncertainty relation") - return True - - return False - - -def is_pure_cov(cov, hbar=2, rtol=1e-05, atol=1e-08): - r""" Checks if the covariance matrix is a valid quantum covariance matrix - that corresponds to a quantum pure state - - Args: - cov (array): a covariance matrix - hbar (float): value of hbar in the uncertainty relation - rtol (float): the relative tolerance parameter used in `np.allclose` - atol (float): the absolute tolerance parameter used in `np.allclose` - - Returns: - (bool): whether the given covariance matrix corresponds to a pure state - """ - if is_valid_cov(cov, hbar=hbar, rtol=rtol, atol=atol): - purity = 1 / np.sqrt(np.linalg.det(2 * cov / hbar)) - if np.allclose(purity, 1.0, rtol=rtol, atol=atol): - return True - - return False - - -def is_classical_cov(cov, hbar=2, atol=1e-08): - r""" Checks if the covariance matrix can be efficiently sampled. - - Args: - cov (array): a covariance matrix - hbar (float): value of hbar in the uncertainty relation - atol (float): the absolute tolerance parameter used in `np.allclose` - - Returns: - (bool): whether the given covariance matrix corresponds to a classical state - """ - - if is_valid_cov(cov, hbar=hbar, atol=atol): - (n, _) = cov.shape - vals = np.linalg.eigvalsh(cov - 0.5 * hbar * np.identity(n)) - vals[np.abs(vals) < atol] = 0.0 - - if np.all(vals >= 0): - return True - return False diff --git a/thewalrus/quantum/fock_states_and_tensors.py b/thewalrus/quantum/fock_tensors.py similarity index 90% rename from thewalrus/quantum/fock_states_and_tensors.py rename to thewalrus/quantum/fock_tensors.py index 08b67ad72..6925a0b32 100644 --- a/thewalrus/quantum/fock_states_and_tensors.py +++ b/thewalrus/quantum/fock_tensors.py @@ -13,7 +13,7 @@ # limitations under the License. """ Set of functions for calculating various state representations, probabilites and -fidelities of Gaussian states. +classical subsystems of Gaussian states. """ # pylint: disable=too-many-arguments @@ -23,22 +23,24 @@ import dask from scipy.special import factorial as fac -from scipy.linalg import sqrtm from numba import jit -from ..symplectic import expand, sympmat, is_symplectic +from ..symplectic import expand, is_symplectic from ..libwalrus import interferometer, interferometer_real from .._hafnian import hafnian, hafnian_repeated, reduction from .._hermite_multidimensional import hermite_multidimensional, hafnian_batched -from .covariance_matrices import ( - is_classical_cov, - is_pure_cov, +from .conversions import ( Amat, - Beta, Qmat, reduced_gaussian, + complex_to_real_displacements, +) + +from .gaussian_state_tests import ( + is_classical_cov, + is_pure_cov, ) @@ -67,7 +69,7 @@ def pure_state_amplitude(mu, cov, i, include_prefactor=True, tol=1e-10, hbar=2, raise ValueError("The covariance matrix does not correspond to a pure state") rpt = i - beta = Beta(mu, hbar=hbar) + beta = complex_to_real_displacements(mu, hbar=hbar) Q = Qmat(cov, hbar=hbar) A = Amat(cov, hbar=hbar) (n, _) = cov.shape @@ -140,7 +142,7 @@ def state_vector( if not is_pure_cov(cov, hbar=hbar, rtol=1e-05, atol=1e-08): raise ValueError("The covariance matrix does not correspond to a pure state") - beta = Beta(mu, hbar=hbar) + beta = complex_to_real_displacements(mu, hbar=hbar) A = Amat(cov, hbar=hbar) Q = Qmat(cov, hbar=hbar) @@ -211,7 +213,7 @@ def density_matrix_element(mu, cov, i, j, include_prefactor=True, tol=1e-10, hba complex: the density matrix element """ rpt = i + j - beta = Beta(mu, hbar=hbar) + beta = complex_to_real_displacements(mu, hbar=hbar) A = Amat(cov, hbar=hbar) if np.linalg.norm(beta) < tol: # no displacement @@ -231,7 +233,7 @@ def density_matrix_element(mu, cov, i, j, include_prefactor=True, tol=1e-10, hba haf = hafnian_repeated(A, rpt, mu=gamma, loop=True) if include_prefactor: - haf *= prefactor(mu, cov, hbar=hbar) + haf *= _prefactor(mu, cov, hbar=hbar) return haf / np.sqrt(np.prod(fac(rpt))) @@ -270,7 +272,7 @@ def density_matrix(mu, cov, post_select=None, normalize=False, cutoff=5, hbar=2) np.array[complex]: the density matrix of the Gaussian state """ N = len(mu) // 2 - pref = prefactor(mu, cov, hbar=hbar) + pref = _prefactor(mu, cov, hbar=hbar) if post_select is None: A = Amat(cov, hbar=hbar).conj() @@ -281,7 +283,7 @@ def density_matrix(mu, cov, post_select=None, normalize=False, cutoff=5, hbar=2) -A, cutoff, renorm=True, modified=True ) return tensor.transpose(sf_order) - beta = Beta(mu, hbar=hbar) + beta = complex_to_real_displacements(mu, hbar=hbar) y = beta - A @ beta.conj() tensor = np.real_if_close(pref) * hermite_multidimensional( -A, cutoff, y=y, renorm=True, modified=True @@ -556,57 +558,6 @@ def update_probabilities_with_noise(probs_noise, probs): return probs -def fidelity(mu1, cov1, mu2, cov2, hbar=2, rtol=1e-05, atol=1e-08): - """Calculates the fidelity between two Gaussian quantum states. - - Note that if the covariance matrices correspond to pure states this - function reduces to the modulus square of the overlap of their state vectors. - For the derivation see `'Quantum Fidelity for Arbitrary Gaussian States', Banchi et al. <10.1103/PhysRevLett.115.260501>`_. - - Args: - mu1 (array): vector of means of the first state - cov1 (array): covariance matrix of the first state - mu2 (array): vector of means of the second state - cov2 (array): covariance matrix of the second state - hbar (float): value of hbar in the uncertainty relation - rtol (float): the relative tolerance parameter used in `np.allclose` - atol (float): the absolute tolerance parameter used in `np.allclose` - - Returns: - (float): value of the fidelity between the two states - """ - - n0, n1 = cov1.shape - m0, m1 = cov2.shape - (l0,) = mu1.shape - (l1,) = mu1.shape - if not n0 == n1 == m0 == m1 == l0 == l1: - raise ValueError("The inputs have incompatible shapes") - - v1 = cov1 / hbar - v2 = cov2 / hbar - deltar = (mu1 - mu2) / np.sqrt(hbar / 2) - n0, n1 = cov1.shape - n = n0 // 2 - W = sympmat(n) - - si12 = np.linalg.inv(v1 + v2) - vaux = W.T @ si12 @ (0.25 * W + v2 @ W @ v1) - p1 = vaux @ W - p1 = p1 @ p1 - p1 = np.identity(2 * n) + 0.25 * np.linalg.inv(p1) - if np.allclose(p1, 0, rtol=rtol, atol=atol): - p1 = np.zeros_like(p1) - else: - p1 = sqrtm(p1) - p1 = 2 * (p1 + np.identity(2 * n)) - p1 = p1 @ vaux - f = np.sqrt(np.linalg.det(si12) * np.linalg.det(p1)) * np.exp( - -0.25 * deltar @ si12 @ deltar - ) - return f - - def find_classical_subsystem(cov, hbar=2, atol=1e-08): """Find the largest integer ``k`` so that subsystem in modes ``[0,1,...,k-1]`` is a classical state. @@ -633,7 +584,7 @@ def find_classical_subsystem(cov, hbar=2, atol=1e-08): return k - 1 -def prefactor(mu, cov, hbar=2): +def _prefactor(mu, cov, hbar=2): r"""Returns the prefactor. .. math:: prefactor = \frac{e^{-\beta Q^{-1}\beta^*/2}}{n_1!\cdots n_m! \sqrt{|Q|}} @@ -646,6 +597,6 @@ def prefactor(mu, cov, hbar=2): float: the prefactor """ Q = Qmat(cov, hbar=hbar) - beta = Beta(mu, hbar=hbar) + beta = complex_to_real_displacements(mu, hbar=hbar) Qinv = np.linalg.inv(Q) return np.exp(-0.5 * beta @ Qinv @ beta.conj()) / np.sqrt(np.linalg.det(Q)) diff --git a/thewalrus/quantum/gaussian_state_tests.py b/thewalrus/quantum/gaussian_state_tests.py new file mode 100644 index 000000000..9aedc00b6 --- /dev/null +++ b/thewalrus/quantum/gaussian_state_tests.py @@ -0,0 +1,149 @@ +# Copyright 2019-2020 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 various properties of covariance matrices as well as fidelity +calculations for gaussian states. +""" + +import numpy as np +from scipy.linalg import sqrtm + +from ..symplectic import sympmat + + +def is_valid_cov(cov, hbar=2, rtol=1e-05, atol=1e-08): + r""" Checks if the covariance matrix is a valid quantum covariance matrix. + + Args: + cov (array): a covariance matrix + hbar (float): value of hbar in the uncertainty relation + rtol (float): the relative tolerance parameter used in `np.allclose` + atol (float): the absolute tolerance parameter used in `np.allclose` + + Returns: + (bool): whether the given covariance matrix is a valid covariance matrix + """ + (n, m) = cov.shape + if n != m: + # raise ValueError("The input matrix must be square") + return False + if not np.allclose(cov, np.transpose(cov), rtol=rtol, atol=atol): + # raise ValueError("The input matrix is not symmetric") + return False + if n % 2 != 0: + # raise ValueError("The input matrix is of even dimension") + return False + + nmodes = n // 2 + vals = np.linalg.eigvalsh(cov + 0.5j * hbar * sympmat(nmodes)) + vals[np.abs(vals) < atol] = 0.0 + if np.all(vals >= 0): + # raise ValueError("The input matrix violates the uncertainty relation") + return True + + return False + + +def is_pure_cov(cov, hbar=2, rtol=1e-05, atol=1e-08): + r""" Checks if the covariance matrix is a valid quantum covariance matrix + that corresponds to a quantum pure state + + Args: + cov (array): a covariance matrix + hbar (float): value of hbar in the uncertainty relation + rtol (float): the relative tolerance parameter used in `np.allclose` + atol (float): the absolute tolerance parameter used in `np.allclose` + + Returns: + (bool): whether the given covariance matrix corresponds to a pure state + """ + if is_valid_cov(cov, hbar=hbar, rtol=rtol, atol=atol): + purity = 1 / np.sqrt(np.linalg.det(2 * cov / hbar)) + if np.allclose(purity, 1.0, rtol=rtol, atol=atol): + return True + + return False + + +def is_classical_cov(cov, hbar=2, atol=1e-08): + r""" Checks if the covariance matrix can be efficiently sampled. + + Args: + cov (array): a covariance matrix + hbar (float): value of hbar in the uncertainty relation + atol (float): the absolute tolerance parameter used in `np.allclose` + + Returns: + (bool): whether the given covariance matrix corresponds to a classical state + """ + + if is_valid_cov(cov, hbar=hbar, atol=atol): + (n, _) = cov.shape + vals = np.linalg.eigvalsh(cov - 0.5 * hbar * np.identity(n)) + vals[np.abs(vals) < atol] = 0.0 + + if np.all(vals >= 0): + return True + return False + + +def fidelity(mu1, cov1, mu2, cov2, hbar=2, rtol=1e-05, atol=1e-08): + """Calculates the fidelity between two Gaussian quantum states. + + Note that if the covariance matrices correspond to pure states this + function reduces to the modulus square of the overlap of their state vectors. + For the derivation see `'Quantum Fidelity for Arbitrary Gaussian States', Banchi et al. <10.1103/PhysRevLett.115.260501>`_. + + Args: + mu1 (array): vector of means of the first state + cov1 (array): covariance matrix of the first state + mu2 (array): vector of means of the second state + cov2 (array): covariance matrix of the second state + hbar (float): value of hbar in the uncertainty relation + rtol (float): the relative tolerance parameter used in `np.allclose` + atol (float): the absolute tolerance parameter used in `np.allclose` + + Returns: + (float): value of the fidelity between the two states + """ + + n0, n1 = cov1.shape + m0, m1 = cov2.shape + (l0,) = mu1.shape + (l1,) = mu1.shape + if not n0 == n1 == m0 == m1 == l0 == l1: + raise ValueError("The inputs have incompatible shapes") + + v1 = cov1 / hbar + v2 = cov2 / hbar + deltar = (mu1 - mu2) / np.sqrt(hbar / 2) + n0, n1 = cov1.shape + n = n0 // 2 + W = sympmat(n) + + si12 = np.linalg.inv(v1 + v2) + vaux = W.T @ si12 @ (0.25 * W + v2 @ W @ v1) + p1 = vaux @ W + p1 = p1 @ p1 + p1 = np.identity(2 * n) + 0.25 * np.linalg.inv(p1) + if np.allclose(p1, 0, rtol=rtol, atol=atol): + p1 = np.zeros_like(p1) + else: + p1 = sqrtm(p1) + p1 = 2 * (p1 + np.identity(2 * n)) + p1 = p1 @ vaux + f = np.sqrt(np.linalg.det(si12) * np.linalg.det(p1)) * np.exp( + -0.25 * deltar @ si12 @ deltar + ) + return f diff --git a/thewalrus/quantum/means_and_variances.py b/thewalrus/quantum/means_and_variances.py index 889f3ecf9..0854895c4 100644 --- a/thewalrus/quantum/means_and_variances.py +++ b/thewalrus/quantum/means_and_variances.py @@ -13,19 +13,23 @@ # limitations under the License. """ Functions for constructing/calculating the means, variances and covariances of -Gaussian states and photon number distributions of Gaussian states. +Gaussian states. """ from itertools import product import numpy as np -from .covariance_matrices import normal_ordered_expectation, reduced_gaussian +from .._hafnian import hafnian, reduction +from .conversions import ( + reduced_gaussian, + Covmat, + Qmat, + Xmat, + complex_to_real_displacements +) -################################################################################ -# Calculate means or variances of photon number -################################################################################ def photon_number_mean(mu, cov, j, hbar=2): r""" Calculate the mean photon number of mode j of a Gaussian state. @@ -197,3 +201,106 @@ def photon_number_squared_expectation(mu, cov, modes, hbar=2): term = normal_ordered_expectation(mu_red, cov_red, rpt, hbar=hbar) result += term return result + + +def normal_ordered_expectation(mu, cov, rpt, hbar=2): + r"""Calculates the expectation value of the normal ordered product + :math:`\prod_{i=0}^{N-1} a_i^{\dagger n_i} \prod_{j=0}^{N-1} a_j^{m_j}` with respect to an N-mode Gaussian state, + where :math:`\text{rpt}=(n_0, n_1, \ldots, n_{N-1}, m_0, m_1, \ldots, m_{N-1})`. + + Args: + mu (array): length-:math:`2N` means vector in xp-ordering. + cov (array): :math:`2N\times 2N` covariance matrix in xp-ordering. + rpt (list): integers specifying the terms to calculate. + hbar (float): value of hbar in the uncertainty relation. + + Returns: + (float): expectation value of the normal ordered product of operators + """ + alpha = complex_to_real_displacements(mu, hbar=hbar) + n = len(cov) + V = (Qmat(cov, hbar=hbar) - np.identity(n)) @ Xmat(n // 2) + A = reduction(V, rpt) + if np.allclose(mu, 0): + res = np.conj(hafnian(A)) + else: + np.fill_diagonal(A, reduction(np.conj(alpha), rpt)) + res = np.conj(hafnian(A, loop=True)) + return np.conj(res) + + +def mean_number_of_clicks(cov, hbar=2): + r""" Calculates the total mean number of clicks when a zero-mean gaussian state + is measured using threshold detectors. + + Args + cov (array): :math:`2N\times 2N` covariance matrix in xp-ordering + hbar (float): the value of :math:`\hbar` in the commutation relation :math:`[\x,\p]=i\hbar` + + Returns + float: mean number of clicks + """ + n, _ = cov.shape + nmodes = n // 2 + Q = Qmat(cov, hbar=hbar) + meanc = 1.0 * nmodes + + for i in range(nmodes): + det_val = np.real(Q[i, i] * Q[i + nmodes, i + nmodes] - Q[i + nmodes, i] * Q[i, i + nmodes]) + meanc -= 1.0 / np.sqrt(det_val) + return meanc + + +def variance_number_of_clicks(cov, hbar=2): + r""" Calculates the variance of the total number of clicks when a zero-mean gaussian state + is measured using threshold detectors. + + Args + cov (array): :math:`2N\times 2N` covariance matrix in xp-ordering + hbar (float): the value of :math:`\hbar` in the commutation relation :math:`[\x,\p]=i\hbar` + + Returns + float: variance in the total number of clicks + """ + n, _ = cov.shape + means = np.zeros([n]) + nmodes = n // 2 + Q = Qmat(cov, hbar=hbar) + vac_probs = np.array( + [ + np.real(Q[i, i] * Q[i + nmodes, i + nmodes] - Q[i + nmodes, i] * Q[i, i + nmodes]) + for i in range(nmodes) + ] + ) + vac_probs = np.sqrt(vac_probs) + vac_probs = 1 / vac_probs + term1 = np.sum(vac_probs * (1 - vac_probs)) + term2 = 0 + for i in range(nmodes): + for j in range(i): + _, Qij = reduced_gaussian(means, Q, [i, j]) + prob_vac_ij = np.linalg.det(Qij).real + prob_vac_ij = 1.0 / np.sqrt(prob_vac_ij) + term2 += prob_vac_ij - vac_probs[i] * vac_probs[j] + + return term1 + 2 * term2 + + +def mean_number_of_clicks_graph(A): + r""" Given an adjacency matrix this function calculates the mean number of clicks. + For this to make sense the user must provide a matrix with singular values + less than or equal to one. See Appendix A.3 of `_ + by Banchi et al. + + Args: + A (array): rescaled adjacency matrix + + Returns: + float: mean number of clicks + """ + n, _ = A.shape + idn = np.identity(n) + X = np.block([[0 * idn, idn], [idn, 0 * idn]]) + B = np.block([[A, 0 * A], [0 * A, np.conj(A)]]) + Q = np.linalg.inv(np.identity(2 * n) - X @ B) + return mean_number_of_clicks(Covmat(Q)) diff --git a/thewalrus/quantum/photon_number_distributions.py b/thewalrus/quantum/photon_number_distributions.py index 1b158277e..1939af67b 100644 --- a/thewalrus/quantum/photon_number_distributions.py +++ b/thewalrus/quantum/photon_number_distributions.py @@ -18,14 +18,11 @@ import numpy as np from scipy.stats import nbinom -from .covariance_matrices import Amat, is_pure_cov +from .conversions import Amat +from .gaussian_state_tests import is_pure_cov -################################################################################ -# Calculate the total photon number distribution -################################################################################ - -def total_photon_num_dist_pure_state(cov, cutoff=50, hbar=2, padding_factor=2): +def pure_state_distribution(cov, cutoff=50, hbar=2, padding_factor=2): r""" Calculates the total photon number distribution of a pure state with zero mean. @@ -45,11 +42,11 @@ def total_photon_num_dist_pure_state(cov, cutoff=50, hbar=2, padding_factor=2): N = n // 2 B = A[0:N, 0:N] rs = np.arctanh(np.linalg.svd(B, compute_uv=False)) - return gen_multi_mode_dist(rs, cutoff=cutoff, padding_factor=padding_factor)[0:cutoff] + return _convolve_squeezed_state_distribution(rs, cutoff=cutoff, padding_factor=padding_factor)[0:cutoff] raise ValueError("The Gaussian state is not pure") -def gen_single_mode_dist(s, cutoff=50, N=1): +def _squeezed_state_distribution(s, cutoff=50, N=1): """Generate the photon number distribution of :math:`N` identical single mode squeezed states. Args: @@ -75,7 +72,7 @@ def gen_single_mode_dist(s, cutoff=50, N=1): return ps_tot -def gen_multi_mode_dist(s, cutoff=50, padding_factor=2): +def _convolve_squeezed_state_distribution(s, cutoff=50, padding_factor=2): """Generates the total photon number distribution of single mode squeezed states with different squeezing values. Args: @@ -90,5 +87,5 @@ def gen_multi_mode_dist(s, cutoff=50, padding_factor=2): ps = np.zeros(cutoff_sc) ps[0] = 1.0 for s_val in s: - ps = np.convolve(ps, gen_single_mode_dist(s_val, cutoff_sc))[0:cutoff_sc] + ps = np.convolve(ps, _squeezed_state_distribution(s_val, cutoff_sc))[0:cutoff_sc] return ps diff --git a/thewalrus/tests/test_quantum.py b/thewalrus/tests/test_quantum.py index ca58a2da6..702992b7f 100644 --- a/thewalrus/tests/test_quantum.py +++ b/thewalrus/tests/test_quantum.py @@ -29,8 +29,8 @@ Xmat, Qmat, Amat, - Beta, - prefactor, + complex_to_real_displacements, + _prefactor, density_matrix_element, density_matrix, find_scaling_adjacency_matrix, @@ -38,7 +38,7 @@ mean_number_of_clicks_graph, Covmat, gen_Qmat_from_graph, - Means, + real_to_complex_displacements, photon_number_covmat, is_valid_cov, is_pure_cov, @@ -46,8 +46,8 @@ state_vector, is_classical_cov, find_classical_subsystem, - total_photon_num_dist_pure_state, - gen_single_mode_dist, + pure_state_distribution, + _squeezed_state_distribution, fock_tensor, photon_number_mean_vector, photon_number_mean, @@ -180,7 +180,7 @@ def test_Amat_TMS_using_Q(): def test_beta(): """test the correct beta is returned""" mu = np.arange(4) - res = Beta(mu) + res = complex_to_real_displacements(mu) alpha = (mu[:2] + 1j * mu[2:]) / np.sqrt(2 * 2) ex = np.concatenate([alpha, alpha.conj()]) @@ -190,8 +190,8 @@ def test_beta(): def test_Means(): """test the correct beta is returned""" res = np.arange(4) - mu = Beta(res) - ex = Means(mu) + mu = complex_to_real_displacements(res) + ex = real_to_complex_displacements(mu) assert np.allclose(res, ex) @@ -200,7 +200,7 @@ def test_prefactor_vacuum(): Q = np.identity(2) beta = np.zeros([2]) - res = prefactor(Means(beta), Covmat(Q)) + res = _prefactor(real_to_complex_displacements(beta), Covmat(Q)) ex = 1 assert np.allclose(res, ex) @@ -213,7 +213,7 @@ def test_prefactor_TMS(): beta = np.zeros([4]) - res = prefactor(Means(beta), Covmat(Q)) + res = _prefactor(real_to_complex_displacements(beta), Covmat(Q)) ex = 0.5 assert np.allclose(res, ex) @@ -228,7 +228,7 @@ def test_prefactor_with_displacement(): vect = 1.2 * np.ones([2]) + 1j * np.ones(2) beta = np.concatenate([vect, vect.conj()]) - res = prefactor(Means(beta), Covmat(Q), hbar=2) + res = _prefactor(real_to_complex_displacements(beta), Covmat(Q), hbar=2) ex = np.exp(-0.5 * beta @ Qinv @ beta.conj()) / np.sqrt(np.linalg.det(Q)) assert np.allclose(res, ex) @@ -240,18 +240,18 @@ def test_density_matrix_element_vacuum(): el = [[0], [0]] ex = 1 - res = density_matrix_element(Means(beta), Covmat(Q), el[0], el[1]) + res = density_matrix_element(real_to_complex_displacements(beta), Covmat(Q), el[0], el[1]) assert np.allclose(ex, res) el = [[1], [1]] # res = density_matrix_element(beta, A, Q, el[0], el[1]) - res = density_matrix_element(Means(beta), Covmat(Q), el[0], el[1]) + res = density_matrix_element(real_to_complex_displacements(beta), Covmat(Q), el[0], el[1]) assert np.allclose(0, res) el = [[1], [0]] # res = density_matrix_element(beta, A, Q, el[0], el[1]) - res = density_matrix_element(Means(beta), Covmat(Q), el[0], el[1]) + res = density_matrix_element(real_to_complex_displacements(beta), Covmat(Q), el[0], el[1]) assert np.allclose(0, res) @@ -282,12 +282,12 @@ def test_density_matrix_element_vacuum(): @pytest.mark.parametrize("t", [t0, t1, t2, t3, t4]) def test_density_matrix_element_disp(t): """Test density matrix elements for a state with displacement""" - beta = Beta(mu) + beta = complex_to_real_displacements(mu) Q = Qmat(V) el = t[0] ex = t[1] - res = density_matrix_element(Means(beta), Covmat(Q), el[0], el[1]) + res = density_matrix_element(real_to_complex_displacements(beta), Covmat(Q), el[0], el[1]) assert np.allclose(ex, res) @@ -302,12 +302,12 @@ def test_density_matrix_element_disp(t): @pytest.mark.parametrize("t", [t0, t1, t2, t3, t4]) def test_density_matrix_element_no_disp(t): """Test density matrix elements for a state with no displacement""" - beta = Beta(np.zeros([6])) + beta = complex_to_real_displacements(np.zeros([6])) Q = Qmat(V) el = t[0] ex = t[1] - res = density_matrix_element(Means(beta), Covmat(Q), el[0], el[1]) + res = density_matrix_element(real_to_complex_displacements(beta), Covmat(Q), el[0], el[1]) assert np.allclose(ex, res) @@ -565,7 +565,7 @@ def test_pure_state_amplitude_coherent(i): """ Tests pure state amplitude for a coherent state """ cov = np.identity(2) mu = np.array([1.0, 2.0]) - beta = Beta(mu) + beta = complex_to_real_displacements(mu) alpha = beta[0] exact = np.exp(-0.5 * np.abs(alpha) ** 2) * alpha ** i / np.sqrt(np.math.factorial(i)) num = pure_state_amplitude(mu, cov, [i]) @@ -628,7 +628,7 @@ def test_state_vector_coherent(): cutoff = 5 cov = np.identity(2) mu = np.array([1.0, 2.0]) - beta = Beta(mu) + beta = complex_to_real_displacements(mu) alpha = beta[0] exact = np.array([(np.exp(-0.5 * np.abs(alpha) ** 2) * alpha ** i / np.sqrt(np.math.factorial(i))) for i in range(cutoff)]) num = state_vector(mu, cov, cutoff=cutoff) @@ -669,15 +669,15 @@ def test_is_classical_cov_thermal(nbar): @pytest.mark.parametrize("cutoff", [50, 51, 52, 53]) -def test_total_photon_num_dist_pure_state(cutoff): +def test_pure_state_distribution(cutoff): """ Test the correct photon number distribution is obtained for n modes with nmean number of photons up to Fock cutoff nmax""" n = 3 nmean = 1.0 rs = np.arcsinh(np.sqrt(nmean)) * np.ones([n]) cov = np.diag(np.concatenate([np.exp(2 * rs), np.exp(-2 * rs)])) - p1 = total_photon_num_dist_pure_state(cov, cutoff=cutoff) - p2 = gen_single_mode_dist(np.arcsinh(np.sqrt(nmean)), N=n, cutoff=cutoff) + p1 = pure_state_distribution(cov, cutoff=cutoff) + p2 = _squeezed_state_distribution(np.arcsinh(np.sqrt(nmean)), N=n, cutoff=cutoff) assert np.allclose(p1, p2) @@ -1033,14 +1033,14 @@ def test_update_with_noise_coherent(num_modes, parallel, monkeypatch): noise_dists = np.array([poisson.pmf(np.arange(cutoff), nbar) for nbar in nbar_vals]) hbar = 2 beta = np.random.rand(num_modes) + 1j * np.random.rand(num_modes) - means = Means(np.concatenate((beta, beta.conj())), hbar=hbar) + means = real_to_complex_displacements(np.concatenate((beta, beta.conj())), hbar=hbar) cov = hbar * np.identity(2 * num_modes) / 2 cutoff = 10 probs = probabilities(means, cov, cutoff, parallel=parallel, hbar=2) updated_probs = update_probabilities_with_noise(noise_dists, probs) beta_expected = np.sqrt(nbar_vals + np.abs(beta) ** 2) - means_expected = Means( + means_expected = real_to_complex_displacements( np.concatenate((beta_expected, beta_expected.conj())), hbar=hbar ) expected = probabilities(means_expected, cov, cutoff, parallel=parallel, hbar=2) @@ -1055,7 +1055,7 @@ def test_update_with_noise_coherent_value_error(): noise_dists = np.array([poisson.pmf(np.arange(cutoff), nbar) for nbar in nbar_vals]) hbar = 2 beta = np.random.rand(num_modes) + 1j * np.random.rand(num_modes) - means = Means(np.concatenate((beta, beta.conj())), hbar=hbar) + means = real_to_complex_displacements(np.concatenate((beta, beta.conj())), hbar=hbar) cov = hbar * np.identity(2 * num_modes) / 2 cutoff = 10 probs = probabilities(means, cov, cutoff, hbar=2) @@ -1099,8 +1099,8 @@ def test_fidelity_coherent_state(num_modes, hbar): """Test the fidelity of two multimode coherent states""" beta1 = np.random.rand(num_modes) + 1j * np.random.rand(num_modes) beta2 = np.random.rand(num_modes) + 1j * np.random.rand(num_modes) - means1 = Means(np.concatenate([beta1, beta1.conj()]), hbar=hbar) - means2 = Means(np.concatenate([beta2, beta2.conj()]), hbar=hbar) + means1 = real_to_complex_displacements(np.concatenate([beta1, beta1.conj()]), hbar=hbar) + means2 = real_to_complex_displacements(np.concatenate([beta2, beta2.conj()]), hbar=hbar) cov1 = hbar * np.identity(2 * num_modes) / 2 cov2 = hbar * np.identity(2 * num_modes) / 2 fid = fidelity(means1, cov1, means2, cov2, hbar=hbar) @@ -1114,7 +1114,7 @@ def test_fidelity_coherent_state(num_modes, hbar): def test_fidelity_vac_to_displaced_squeezed(r, alpha, hbar): """Calculates the fidelity between a coherent squeezed state and vacuum""" cov1 = np.diag([np.exp(2 * r), np.exp(-2 * r)]) * hbar / 2 - means1 = Means(np.array([alpha, np.conj(alpha)]), hbar=hbar) + means1 = real_to_complex_displacements(np.array([alpha, np.conj(alpha)]), hbar=hbar) means2 = np.zeros([2]) cov2 = np.identity(2) * hbar / 2 expected = ( @@ -1152,7 +1152,7 @@ def test_fidelity_wrong_shape(): def test_expectation_normal_ordered_coherent(alpha, hbar): """Test the correct evaluation of the normal ordered expectation value for a product of coherent states""" beta = np.concatenate([alpha, np.conj(alpha)]) - means = Means(beta, hbar=hbar) + means = real_to_complex_displacements(beta, hbar=hbar) cov = np.identity(len(beta)) * hbar / 2 pattern = np.random.randint(low=0, high=3, size=len(beta)) result = normal_ordered_expectation(means, cov, pattern, hbar=hbar) @@ -1219,7 +1219,7 @@ def test_single_mode_displaced_squeezed(r, phi, x, y): cov = 0.5 * hbar * S @ S.T beta = np.array([x + 1j * y, x - 1j * y]) alpha = beta[0] - means = Means(beta, hbar=hbar) + means = real_to_complex_displacements(beta, hbar=hbar) a = alpha adxa = np.abs(alpha) ** 2 + np.sinh(r) ** 2 a2 = alpha ** 2 - np.exp(1j * phi) * 0.5 * np.sinh(2 * r) @@ -1267,7 +1267,7 @@ def test_expt_two_mode_squeezed(r, phi): def test_photon_number_expectation_displaced(alpha, hbar): """Tests the correct photon number expectation for coherent states""" beta = np.concatenate([alpha, np.conj(alpha)]) - means = Means(beta, hbar=hbar) + means = real_to_complex_displacements(beta, hbar=hbar) cov = np.identity(len(beta)) * hbar / 2 val = photon_number_expectation( means, cov, modes=list(range(len(alpha))), hbar=hbar From 45fa38646855521f6b7e6fc25f28e9d48654660a Mon Sep 17 00:00:00 2001 From: Theodor Isacsson Date: Wed, 9 Sep 2020 17:31:58 -0400 Subject: [PATCH 10/17] Fix spelling --- thewalrus/quantum/conversions.py | 2 +- thewalrus/quantum/fock_tensors.py | 2 +- thewalrus/quantum/gaussian_state_tests.py | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/thewalrus/quantum/conversions.py b/thewalrus/quantum/conversions.py index ca2d8e8a3..b445b9555 100644 --- a/thewalrus/quantum/conversions.py +++ b/thewalrus/quantum/conversions.py @@ -13,7 +13,7 @@ # limitations under the License. """ Functions for transforming one type of covariance-matrix-like object into -another as well as reducing a gaussian state. +another as well as reducing a Gaussian state. """ import numpy as np diff --git a/thewalrus/quantum/fock_tensors.py b/thewalrus/quantum/fock_tensors.py index 6925a0b32..dc694a7b5 100644 --- a/thewalrus/quantum/fock_tensors.py +++ b/thewalrus/quantum/fock_tensors.py @@ -12,7 +12,7 @@ # See the License for the specific language governing permissions and # limitations under the License. """ -Set of functions for calculating various state representations, probabilites and +Set of functions for calculating various state representations, probabilities and classical subsystems of Gaussian states. """ # pylint: disable=too-many-arguments diff --git a/thewalrus/quantum/gaussian_state_tests.py b/thewalrus/quantum/gaussian_state_tests.py index 9aedc00b6..12a0bc36a 100644 --- a/thewalrus/quantum/gaussian_state_tests.py +++ b/thewalrus/quantum/gaussian_state_tests.py @@ -13,7 +13,7 @@ # limitations under the License. """ Tests for various properties of covariance matrices as well as fidelity -calculations for gaussian states. +calculations for Gaussian states. """ import numpy as np From c1526a8e7ca0b5b004810d668fd76e6d522757c9 Mon Sep 17 00:00:00 2001 From: Theodor Isacsson Date: Thu, 10 Sep 2020 13:27:30 -0400 Subject: [PATCH 11/17] Updates from code-review --- thewalrus/quantum/__init__.py | 33 ++++++------- thewalrus/quantum/adjacency_matrices.py | 36 ++++++++++---- thewalrus/quantum/fock_tensors.py | 2 +- ...sian_state_tests.py => gaussian_checks.py} | 0 thewalrus/quantum/means_and_variances.py | 25 +--------- .../quantum/photon_number_distributions.py | 2 +- thewalrus/tests/test_quantum.py | 48 +++++++++---------- 7 files changed, 71 insertions(+), 75 deletions(-) rename thewalrus/quantum/{gaussian_state_tests.py => gaussian_checks.py} (100%) diff --git a/thewalrus/quantum/__init__.py b/thewalrus/quantum/__init__.py index 439c41bd9..bb755ed4c 100644 --- a/thewalrus/quantum/__init__.py +++ b/thewalrus/quantum/__init__.py @@ -102,11 +102,10 @@ Amat complex_to_real_displacements real_to_complex_displacements - find_scaling_adjacency_matrix - mean_number_of_clicks - mean_number_of_clicks_graph - find_scaling_adjacency_matrix_torontonian - gen_Qmat_from_graph + adj_scaling + mean_clicks + adj_scaling_torontonian + adj_to_qmat photon_number_mean photon_number_mean_vector photon_number_covar @@ -133,16 +132,15 @@ update_probabilities_with_loss, update_probabilities_with_noise, find_classical_subsystem, - _prefactor, ) from .adjacency_matrices import ( - find_scaling_adjacency_matrix, - find_scaling_adjacency_matrix_torontonian, - gen_Qmat_from_graph, + adj_scaling, + adj_scaling_torontonian, + adj_to_qmat, ) -from .gaussian_state_tests import ( +from .gaussian_checks import ( is_valid_cov, is_pure_cov, is_classical_cov, @@ -168,15 +166,12 @@ photon_number_expectation, photon_number_squared_expectation, normal_ordered_expectation, - mean_number_of_clicks, - variance_number_of_clicks, - mean_number_of_clicks_graph, + mean_clicks, + variance_clicks, ) from .photon_number_distributions import ( pure_state_distribution, - _squeezed_state_distribution, - _convolve_squeezed_state_distribution, ) __all__ = [ @@ -196,7 +191,9 @@ # old names for functions; remove in due time Means = real_to_complex_displacements Beta = complex_to_real_displacements -prefactor = _prefactor -gen_multi_mode_dist = _convolve_squeezed_state_distribution -gen_single_mode_dist = _squeezed_state_distribution total_photon_num_dist_pure_state = pure_state_distribution +gen_Qmat_from_graph = adj_to_qmat +find_scaling_adjacency_matrix = adj_scaling +find_scaling_adjacency_matrix_torontonian = adj_scaling_torontonian +mean_number_of_clicks = mean_clicks +variance_number_of_clicks = variance_clicks diff --git a/thewalrus/quantum/adjacency_matrices.py b/thewalrus/quantum/adjacency_matrices.py index 132bce2a3..75c2749e3 100644 --- a/thewalrus/quantum/adjacency_matrices.py +++ b/thewalrus/quantum/adjacency_matrices.py @@ -19,11 +19,11 @@ import numpy as np from scipy.optimize import root_scalar -from .conversions import Xmat -from .means_and_variances import mean_number_of_clicks_graph +from .conversions import Xmat, Covmat +from .means_and_variances import mean_clicks -def find_scaling_adjacency_matrix_torontonian(A, c_mean): +def adj_scaling_torontonian(A, c_mean): r""" Returns the scaling parameter by which the adjacency matrix A should be rescaled so that the Gaussian state that encodes it has give a mean number of clicks equal to ``c_mean`` when measured with @@ -58,7 +58,7 @@ def cost(x): return c_mean - n if x <= 0: return c_mean - return c_mean - mean_number_of_clicks_graph(x * localA) + return c_mean - _mean_clicks_adj(x * localA) res = root_scalar(cost, x0=0.5, bracket=(0.0, 1.0)) # Do the optimization @@ -67,7 +67,27 @@ def cost(x): return res.root / vals[0] -def find_scaling_adjacency_matrix(A, n_mean): +def _mean_clicks_adj(A): + r""" Given an adjacency matrix this function calculates the mean number of clicks. + For this to make sense the user must provide a matrix with singular values + less than or equal to one. See Appendix A.3 of `_ + by Banchi et al. + + Args: + A (array): rescaled adjacency matrix + + Returns: + float: mean number of clicks + """ + n, _ = A.shape + idn = np.identity(n) + X = np.block([[0 * idn, idn], [idn, 0 * idn]]) + B = np.block([[A, 0 * A], [0 * A, np.conj(A)]]) + Q = np.linalg.inv(np.identity(2 * n) - X @ B) + return mean_clicks(Covmat(Q)) + + +def adj_scaling(A, n_mean): r""" Returns the scaling parameter by which the adjacency matrix A should be rescaled so that the Gaussian state that endodes it has a total mean photon number n_mean. @@ -104,7 +124,7 @@ def mean_photon_number(x, vals): n = np.sum(vals2 / (1.0 - vals2)) return n - # The following function is implicitly tested in test_find_scaling_adjacency_matrix + # The following function is implicitly tested in test_adj_scaling def grad_mean_photon_number(x, vals): # pragma: no cover r""" Returns the gradient od the mean number of photons in the Gaussian state that encodes the adjacency matrix x*A with respect to x. @@ -132,7 +152,7 @@ def grad_mean_photon_number(x, vals): # pragma: no cover return res.root -def gen_Qmat_from_graph(A, n_mean): +def adj_to_qmat(A, n_mean): r""" Returns the Qmat xp-covariance matrix associated to a graph with adjacency matrix :math:`A` and with mean photon number :math:`n_{mean}`. @@ -148,7 +168,7 @@ def gen_Qmat_from_graph(A, n_mean): if n != m: raise ValueError("Matrix must be square.") - sc = find_scaling_adjacency_matrix(A, n_mean) + sc = adj_scaling(A, n_mean) Asc = sc * A A = np.block([[Asc, 0 * Asc], [0 * Asc, Asc.conj()]]) I = np.identity(2 * n) diff --git a/thewalrus/quantum/fock_tensors.py b/thewalrus/quantum/fock_tensors.py index dc694a7b5..fc1be6ca5 100644 --- a/thewalrus/quantum/fock_tensors.py +++ b/thewalrus/quantum/fock_tensors.py @@ -38,7 +38,7 @@ complex_to_real_displacements, ) -from .gaussian_state_tests import ( +from .gaussian_checks import ( is_classical_cov, is_pure_cov, ) diff --git a/thewalrus/quantum/gaussian_state_tests.py b/thewalrus/quantum/gaussian_checks.py similarity index 100% rename from thewalrus/quantum/gaussian_state_tests.py rename to thewalrus/quantum/gaussian_checks.py diff --git a/thewalrus/quantum/means_and_variances.py b/thewalrus/quantum/means_and_variances.py index 0854895c4..82fd6fa7f 100644 --- a/thewalrus/quantum/means_and_variances.py +++ b/thewalrus/quantum/means_and_variances.py @@ -24,7 +24,6 @@ from .conversions import ( reduced_gaussian, - Covmat, Qmat, Xmat, complex_to_real_displacements @@ -229,7 +228,7 @@ def normal_ordered_expectation(mu, cov, rpt, hbar=2): return np.conj(res) -def mean_number_of_clicks(cov, hbar=2): +def mean_clicks(cov, hbar=2): r""" Calculates the total mean number of clicks when a zero-mean gaussian state is measured using threshold detectors. @@ -251,7 +250,7 @@ def mean_number_of_clicks(cov, hbar=2): return meanc -def variance_number_of_clicks(cov, hbar=2): +def variance_clicks(cov, hbar=2): r""" Calculates the variance of the total number of clicks when a zero-mean gaussian state is measured using threshold detectors. @@ -284,23 +283,3 @@ def variance_number_of_clicks(cov, hbar=2): term2 += prob_vac_ij - vac_probs[i] * vac_probs[j] return term1 + 2 * term2 - - -def mean_number_of_clicks_graph(A): - r""" Given an adjacency matrix this function calculates the mean number of clicks. - For this to make sense the user must provide a matrix with singular values - less than or equal to one. See Appendix A.3 of `_ - by Banchi et al. - - Args: - A (array): rescaled adjacency matrix - - Returns: - float: mean number of clicks - """ - n, _ = A.shape - idn = np.identity(n) - X = np.block([[0 * idn, idn], [idn, 0 * idn]]) - B = np.block([[A, 0 * A], [0 * A, np.conj(A)]]) - Q = np.linalg.inv(np.identity(2 * n) - X @ B) - return mean_number_of_clicks(Covmat(Q)) diff --git a/thewalrus/quantum/photon_number_distributions.py b/thewalrus/quantum/photon_number_distributions.py index 1939af67b..da572167c 100644 --- a/thewalrus/quantum/photon_number_distributions.py +++ b/thewalrus/quantum/photon_number_distributions.py @@ -19,7 +19,7 @@ from scipy.stats import nbinom from .conversions import Amat -from .gaussian_state_tests import is_pure_cov +from .gaussian_checks import is_pure_cov def pure_state_distribution(cov, cutoff=50, hbar=2, padding_factor=2): diff --git a/thewalrus/tests/test_quantum.py b/thewalrus/tests/test_quantum.py index 702992b7f..116a659fa 100644 --- a/thewalrus/tests/test_quantum.py +++ b/thewalrus/tests/test_quantum.py @@ -24,20 +24,21 @@ from thewalrus.random import random_covariance, random_interferometer +from thewalrus.quantum.fock_tensors import _prefactor +from thewalrus.quantum.photon_number_distributions import _squeezed_state_distribution +from thewalrus.quantum.adjacency_matrices import _mean_clicks_adj from thewalrus.quantum import ( reduced_gaussian, Xmat, Qmat, Amat, complex_to_real_displacements, - _prefactor, density_matrix_element, density_matrix, - find_scaling_adjacency_matrix, - find_scaling_adjacency_matrix_torontonian, - mean_number_of_clicks_graph, + adj_scaling, + adj_scaling_torontonian, Covmat, - gen_Qmat_from_graph, + adj_to_qmat, real_to_complex_displacements, photon_number_covmat, is_valid_cov, @@ -47,7 +48,6 @@ is_classical_cov, find_classical_subsystem, pure_state_distribution, - _squeezed_state_distribution, fock_tensor, photon_number_mean_vector, photon_number_mean, @@ -59,8 +59,8 @@ normal_ordered_expectation, photon_number_expectation, photon_number_squared_expectation, - variance_number_of_clicks, - mean_number_of_clicks, + variance_clicks, + mean_clicks, ) @@ -393,54 +393,54 @@ def test_density_matrix_displaced_squeezed_postselect(): assert np.allclose(res, expected) -def test_find_scaling_adjacency_matrix(): +def test_adj_scaling(): """Test the find_scaling_adjacency matrix for a the one mode case""" r = 0.75 + 0.9j rabs = np.abs(r) n_mean = 10.0 A = r * np.identity(1) sc_exact = np.sqrt(n_mean / (1.0 + n_mean)) / rabs - sc_num = find_scaling_adjacency_matrix(A, n_mean) + sc_num = adj_scaling(A, n_mean) assert np.allclose(sc_exact, sc_num) -def test_find_scaling_adjacency_matrix_torontonian(): - """Test the find_scaling_adjacency_matrix_torontonian for a multimode problem""" +def test_adj_scaling_torontonian(): + """Test the adj_scaling_torontonian for a multimode problem""" n = 10 A = np.random.rand(n, n) + 1j * np.random.rand(n, n) A += A.T nc = 3.0 - x = find_scaling_adjacency_matrix_torontonian(A, nc) - assert np.allclose(mean_number_of_clicks_graph(x * A), nc) + x = adj_scaling_torontonian(A, nc) + assert np.allclose(_mean_clicks_adj(x * A), nc) -def test_mean_number_of_clicks_graph(): +def test_mean_clicks_adj(): """Test that a two mode squeezed vacuum with parameter r has mean number of clicks equal to 2*tanh(r)""" r = 3.0 tr = np.tanh(r) A = np.array([[0, tr], [tr, 0]]) - value = mean_number_of_clicks_graph(A) + value = _mean_clicks_adj(A) expected = 2 * tr**2 assert np.allclose(expected, value) @pytest.mark.parametrize("hbar", [1.0 / 137, 1, 2, 0.5]) @pytest.mark.parametrize("theta", [0, 1, 2, 3]) @pytest.mark.parametrize("r", [0, 1, 2, 3]) -def test_variance_number_of_clicks(r, theta, hbar): +def test_variance_clicks(r, theta, hbar): """Test one gets the correct variance of the number of clicks""" r = np.arcsinh(1) V = two_mode_squeezing(2 * r, theta) * hbar / 2 - var = variance_number_of_clicks(V, hbar=hbar) + var = variance_clicks(V, hbar=hbar) expected = (4 * np.tanh(r) ** 2) * (1 - np.tanh(r) ** 2) assert np.allclose(var, expected) @pytest.mark.parametrize("hbar", [1.0 / 137, 1, 2, 0.5]) @pytest.mark.parametrize("theta", [0, 1, 2, 3]) @pytest.mark.parametrize("r", [0, 1, 2, 3]) -def test_mean_number_of_clicks(r, theta, hbar): +def test_mean_clicks(r, theta, hbar): """Test one gets the correct mean of the number of clicks""" r = np.arcsinh(1) V = two_mode_squeezing(2 * r, theta) * hbar / 2 - mean = mean_number_of_clicks(V, hbar=hbar) + mean = mean_clicks(V, hbar=hbar) expected = 2 * np.tanh(r) ** 2 assert np.allclose(mean, expected) @@ -449,7 +449,7 @@ def test_Covmat(): n = 1 B = np.random.rand(n, n) + 1j * np.random.rand(n, n) B = B + B.T - sc = find_scaling_adjacency_matrix(B, 1) + sc = adj_scaling(B, 1) idm = np.identity(2 * n) X = Xmat(n) Bsc = sc * B @@ -460,11 +460,11 @@ def test_Covmat(): assert np.allclose(Q, Qrec) -def test_gen_Qmat_from_graph(): - """ Test the gen_Qmat_from_graph for the analytically solvable case of a single mode""" +def test_adj_to_qmat(): + """ Test the adj_to_qmat for the analytically solvable case of a single mode""" A = np.array([[10.0]]) n_mean = 1.0 - cov = Covmat(gen_Qmat_from_graph(A, n_mean)) + cov = Covmat(adj_to_qmat(A, n_mean)) r = np.arcsinh(np.sqrt(n_mean)) cov_e = np.diag([(np.exp(2 * r)), (np.exp(-2 * r))]) assert np.allclose(cov, cov_e) From 7df040815fac0d076780c6a68277241822317f04 Mon Sep 17 00:00:00 2001 From: Theodor Isacsson Date: Thu, 10 Sep 2020 14:44:19 -0400 Subject: [PATCH 12/17] Update doc --- thewalrus/quantum/__init__.py | 102 ++++++++++++++++++++-------------- 1 file changed, 59 insertions(+), 43 deletions(-) diff --git a/thewalrus/quantum/__init__.py b/thewalrus/quantum/__init__.py index bb755ed4c..0efef5a15 100644 --- a/thewalrus/quantum/__init__.py +++ b/thewalrus/quantum/__init__.py @@ -41,7 +41,7 @@ Fock states and tensors ------------------------ +^^^^^^^^^^^^^^^^^^^^^^^ .. autosummary:: @@ -51,47 +51,32 @@ density_matrix fock_tensor probabilities + loss_mat update_probabilities_with_loss update_probabilities_with_noise - normal_ordered_expectation - fidelity - -Details -^^^^^^^ - -.. autofunction:: - pure_state_amplitude - -.. autofunction:: - state_vector - -.. autofunction:: - density_matrix_element - -.. autofunction:: - density_matrix + find_classical_subsystem -.. autofunction:: - fock_tensor +Adjacency matrices +^^^^^^^^^^^^^^^^^^ -.. autofunction:: - probabilities +.. autosummary:: -.. autofunction:: - update_probabilities_with_loss + adj_scaling + adj_scaling_torontonian + adj_to_qmat -.. autofunction:: - update_probabilities_with_noise +Gaussian checks +^^^^^^^^^^^^^^^ -.. autofunction:: - normal_ordered_expectation +.. autosummary:: -.. autofunction:: + is_valid_cov + is_pure_cov + is_classical_cov fidelity - -Utility functions ------------------ +Conversions +^^^^^^^^^^^ .. autosummary:: @@ -102,20 +87,28 @@ Amat complex_to_real_displacements real_to_complex_displacements - adj_scaling - mean_clicks - adj_scaling_torontonian - adj_to_qmat + +Means and variances +^^^^^^^^^^^^^^^^^^^ + +.. autosummary:: + photon_number_mean photon_number_mean_vector photon_number_covar photon_number_covmat - is_valid_cov - is_pure_cov - is_classical_cov - find_classical_subsystem - pure_state_distribution + photon_number_expectation + photon_number_squared_expectation + normal_ordered_expectation + mean_clicks + variance_clicks + +Photon number distributions +^^^^^^^^^^^^^^^^^^^^^^^^^^^ +.. autosummary:: + + pure_state_distribution Details ^^^^^^^ @@ -149,7 +142,6 @@ from .conversions import ( reduced_gaussian, - Xmat, Qmat, Covmat, @@ -181,10 +173,34 @@ "density_matrix", "fock_tensor", "probabilities", + "loss_mat", "update_probabilities_with_loss", "update_probabilities_with_noise", - "normal_ordered_expectation", + "find_classical_subsystem", + "adj_scaling", + "adj_scaling_torontonian", + "adj_to_qmat", + "is_valid_cov", + "is_pure_cov", + "is_classical_cov", "fidelity", + "reduced_gaussian", + "Xmat", + "Qmat", + "Covmat", + "Amat", + "complex_to_real_displacements", + "real_to_complex_displacements", + "photon_number_mean", + "photon_number_mean_vector", + "photon_number_covar", + "photon_number_covmat", + "photon_number_expectation", + "photon_number_squared_expectation", + "normal_ordered_expectation", + "mean_clicks", + "variance_clicks", + "pure_state_distribution", ] From 44b33dec6b4c4c3b664a2b6bbecabdc25fa0c1c9 Mon Sep 17 00:00:00 2001 From: Theodor Isacsson Date: Thu, 10 Sep 2020 16:00:39 -0400 Subject: [PATCH 13/17] Add deprecation warning --- thewalrus/quantum/__init__.py | 27 +++++++++++++++++++-------- 1 file changed, 19 insertions(+), 8 deletions(-) diff --git a/thewalrus/quantum/__init__.py b/thewalrus/quantum/__init__.py index 0efef5a15..5e39c7bad 100644 --- a/thewalrus/quantum/__init__.py +++ b/thewalrus/quantum/__init__.py @@ -113,6 +113,7 @@ Details ^^^^^^^ """ +import warnings from .fock_tensors import ( pure_state_amplitude, @@ -203,13 +204,23 @@ "pure_state_distribution", ] +def deprecate(new_func): + """Wrapper for deprecated functions to raise warning""" + def wrapper(*args, **kwargs): + warnings.warn( + f"Use {new_func.__name__} instead.", + DeprecationWarning + ) + return new_func(*args, **kwargs) + return wrapper # old names for functions; remove in due time -Means = real_to_complex_displacements -Beta = complex_to_real_displacements -total_photon_num_dist_pure_state = pure_state_distribution -gen_Qmat_from_graph = adj_to_qmat -find_scaling_adjacency_matrix = adj_scaling -find_scaling_adjacency_matrix_torontonian = adj_scaling_torontonian -mean_number_of_clicks = mean_clicks -variance_number_of_clicks = variance_clicks +Means = deprecate(real_to_complex_displacements) +Beta = deprecate(complex_to_real_displacements) +total_photon_num_dist_pure_state = deprecate(pure_state_distribution) +gen_Qmat_from_graph = deprecate(adj_to_qmat) +find_scaling_adjacency_matrix = deprecate(adj_scaling) +find_scaling_adjacency_matrix_torontonian = deprecate(adj_scaling_torontonian) +mean_number_of_clicks = deprecate(mean_clicks) +variance_number_of_clicks = deprecate(variance_clicks) + From e94a99f4b1482ef433d04c021bc8ac70454e4dbe Mon Sep 17 00:00:00 2001 From: Theodor Isacsson Date: Thu, 10 Sep 2020 16:05:10 -0400 Subject: [PATCH 14/17] Fix pylint --- thewalrus/quantum/__init__.py | 1 - thewalrus/quantum/gaussian_checks.py | 1 + 2 files changed, 1 insertion(+), 1 deletion(-) diff --git a/thewalrus/quantum/__init__.py b/thewalrus/quantum/__init__.py index 5e39c7bad..00b906c61 100644 --- a/thewalrus/quantum/__init__.py +++ b/thewalrus/quantum/__init__.py @@ -223,4 +223,3 @@ def wrapper(*args, **kwargs): find_scaling_adjacency_matrix_torontonian = deprecate(adj_scaling_torontonian) mean_number_of_clicks = deprecate(mean_clicks) variance_number_of_clicks = deprecate(variance_clicks) - diff --git a/thewalrus/quantum/gaussian_checks.py b/thewalrus/quantum/gaussian_checks.py index 12a0bc36a..af546af3f 100644 --- a/thewalrus/quantum/gaussian_checks.py +++ b/thewalrus/quantum/gaussian_checks.py @@ -15,6 +15,7 @@ Tests for various properties of covariance matrices as well as fidelity calculations for Gaussian states. """ +# pylint: disable=too-many-arguments import numpy as np from scipy.linalg import sqrtm From cbc6416c9ee4ba08f701003a2e17270eede54563 Mon Sep 17 00:00:00 2001 From: Theodor Date: Mon, 14 Sep 2020 12:35:31 -0400 Subject: [PATCH 15/17] Use `functools.wrap` Co-authored-by: Josh Izaac --- thewalrus/quantum/__init__.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/thewalrus/quantum/__init__.py b/thewalrus/quantum/__init__.py index 00b906c61..58fd204e2 100644 --- a/thewalrus/quantum/__init__.py +++ b/thewalrus/quantum/__init__.py @@ -206,6 +206,8 @@ def deprecate(new_func): """Wrapper for deprecated functions to raise warning""" + + @functools.wraps(new_func) def wrapper(*args, **kwargs): warnings.warn( f"Use {new_func.__name__} instead.", From f1ee2c87f4a13059316272ca70d609ace7ad6827 Mon Sep 17 00:00:00 2001 From: Theodor Date: Mon, 14 Sep 2020 12:36:11 -0400 Subject: [PATCH 16/17] Update thewalrus/quantum/__init__.py Co-authored-by: Josh Izaac --- thewalrus/quantum/__init__.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/thewalrus/quantum/__init__.py b/thewalrus/quantum/__init__.py index 58fd204e2..46018adf3 100644 --- a/thewalrus/quantum/__init__.py +++ b/thewalrus/quantum/__init__.py @@ -210,7 +210,7 @@ def deprecate(new_func): @functools.wraps(new_func) def wrapper(*args, **kwargs): warnings.warn( - f"Use {new_func.__name__} instead.", + f"This function is deprecated and will be removed. Use {new_func.__name__} instead.", DeprecationWarning ) return new_func(*args, **kwargs) From 82cd63f5facb16c2c1caa23d7ca2fcb96153a7c1 Mon Sep 17 00:00:00 2001 From: Theodor Isacsson Date: Mon, 14 Sep 2020 12:37:16 -0400 Subject: [PATCH 17/17] Import functools --- thewalrus/quantum/__init__.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/thewalrus/quantum/__init__.py b/thewalrus/quantum/__init__.py index 46018adf3..79ee698f5 100644 --- a/thewalrus/quantum/__init__.py +++ b/thewalrus/quantum/__init__.py @@ -114,6 +114,7 @@ ^^^^^^^ """ import warnings +import functools from .fock_tensors import ( pure_state_amplitude, @@ -206,7 +207,7 @@ def deprecate(new_func): """Wrapper for deprecated functions to raise warning""" - + @functools.wraps(new_func) def wrapper(*args, **kwargs): warnings.warn(