From 0f99b20f90473f30d42a30c4851dd6653dc5ce67 Mon Sep 17 00:00:00 2001 From: rachelchadwick <33873960+rachelchadwick@users.noreply.github.com> Date: Mon, 29 Apr 2024 19:52:06 -0400 Subject: [PATCH] Squashed commit of the following: MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit commit 7c63905c9623684237091dc5a2779fabea1e79e9 Merge: 127a1a7 f935053 Author: Eli Bourassa <53090166+elib20@users.noreply.github.com> Date: Mon Apr 29 14:16:22 2024 -0400 Merge branch 'master' into internal_modes commit f935053ad9065ad505781298bcd0734534705ab9 Author: Alberto Fumagalli Date: Mon Apr 29 14:14:51 2024 -0400 Add codecov token to test workflow (#390) commit fe9e1124ca902050408e426f42b642ebb2930982 Author: Nicolas Quesada <991946+nquesada@users.noreply.github.com> Date: Wed Apr 10 15:31:34 2024 -0400 Updates references (#384) * updates references * Updates biblio --------- Co-authored-by: Nicolas Quesada commit 25b5f0820520865a311de2e7d44779e7338adc1d Author: Nicolas Quesada <991946+nquesada@users.noreply.github.com> Date: Wed Feb 28 20:35:28 2024 -0500 Implements the pre-Iwasawa and Iwasawa decompositions. (#382) * First working version * First working version * First working version * better tests * better tests * Finally correct * Finally correct * putting the traingles in the right places * Removes unnecesary T * More tests * More tests and spell correction * More tests and spell correction * Removes unnecesary complex number usage * Apply suggestions from code review Co-authored-by: Sebastián Duque Mesa <675763+sduquemesa@users.noreply.github.com> * Apply suggestions from code review Co-authored-by: Sebastián Duque Mesa <675763+sduquemesa@users.noreply.github.com> --------- Co-authored-by: Nicolas Quesada Co-authored-by: Sebastián Duque Mesa <675763+sduquemesa@users.noreply.github.com> commit 127a1a77c91710c52254f52efdb571a1869fa361 Merge: e6de063 8759363 Author: Nicolas Quesada <991946+nquesada@users.noreply.github.com> Date: Tue Feb 13 16:43:20 2024 -0500 Merge branch 'master' into internal_modes commit 875936316fb2e34b9a4e9e778c8a7b9a705367a3 Author: Nicolas Quesada <991946+nquesada@users.noreply.github.com> Date: Thu Feb 1 16:57:23 2024 -0500 Better blochmessiah (#381) * Nicer BM * Updates changelog * Unnecesary import * More simplifications * More simplifications * one more * Updates docstring --------- Co-authored-by: Nicolas Quesada commit 6247fc885d3902a7a5aab7e7dbb63687614b1c09 Author: Nicolas Quesada <991946+nquesada@users.noreply.github.com> Date: Wed Jan 24 10:43:57 2024 -0500 Updates docstring and simplifies Williamson (#380) Co-authored-by: Nicolas Quesada commit 0b1af633a0d745ebb896c4f01c6499873ff16bb0 Author: Sebastián Duque Mesa <675763+sduquemesa@users.noreply.github.com> Date: Thu Jan 11 17:09:32 2024 -0500 Increment version number to `0.22.0-dev` (#379) commit 544c4570905d165b5f3194017ea42bd8b10f424e Author: Sebastián Duque Mesa <675763+sduquemesa@users.noreply.github.com> Date: Wed Jan 10 16:52:38 2024 -0500 bump version to `0.21.0` (#378) commit 2268bf04c8e1abf4c9089a63cddd4d334eb9cf5c Author: Nicolas Quesada <991946+nquesada@users.noreply.github.com> Date: Thu Jan 4 11:12:27 2024 -0500 Adds extra degenerate tests with nonvac nullspace (#377) * Adds extra degenerate tests with nonvac nullspace * updates changelog * updates changelog * updates changelog * Removes unnecesary comment * Still breaking something * Simplifications still work * Simplifications still work * Fixes bug --------- Co-authored-by: Nicolas Quesada commit 8ec2172e5bdb48a7d32886742fe3e7c829e48fd9 Author: Nicolas Quesada <991946+nquesada@users.noreply.github.com> Date: Wed Nov 29 15:53:11 2023 -0500 Implements the Montrealer (#374) * Saving changes * Passes black * Passes black * Black * updates acknowledgements and changelog * updates acknowledgements and changelog * Adds montrealer * Adds Montrealer tests (#375) * Yanic the pirate * mtl test * mtl first tests * mtl test * test mtl * test montrealer * test montrealer * test Montrealer * test montrealer * test montrealer * test montrealer * test montrealer * test montrealer * tests montrealer * montrealer test * montrealer test * montrealer test * montrealer test * montrealer test * test montrealer * test montrealer * montrealer test * montrealer test * montrealer test * montrealer tests * montrealer tests * montrealer tests * montrealer tests * functions header update * test montrealer * test montrealer * test montrealer * test montrealer * test montrealer * test montrealer * test montrealer * test montrealer ready for review * test montrealer * test montrealer * test montrealer * test montrealer * test montrealer * test montrealer - failed tests * test montrealer * test montrealer * test montrealer * test montrealer * test montrealer * the rspms are now ordered * Passes black * garbage test file to be deleted later * minor changes * Done --------- Co-authored-by: Nicolas Quesada * making bots happy * Apply suggestions from code review Co-authored-by: Sebastián Duque Mesa <675763+sduquemesa@users.noreply.github.com> * Apply suggestions from code review Co-authored-by: Sebastián Duque Mesa <675763+sduquemesa@users.noreply.github.com> * Apply suggestions from code review * Apply suggestions from code review Co-authored-by: Sebastián Duque Mesa <675763+sduquemesa@users.noreply.github.com> * Apply suggestions from code review Co-authored-by: Sebastián Duque Mesa <675763+sduquemesa@users.noreply.github.com> * Adds extra tests --------- Co-authored-by: Nicolas Quesada Co-authored-by: Yanic Cardin Co-authored-by: Sebastián Duque Mesa <675763+sduquemesa@users.noreply.github.com> commit 0e163bfb3c4bb92b691d35b2e25ea817719132e4 Author: Nicolas Quesada <991946+nquesada@users.noreply.github.com> Date: Tue Aug 22 10:19:08 2023 -0400 Update symplectic.py (#372) * Update decompositions.py * New takagi * New takagi * Passes black * Simplifies blochmessiah * Simplifies blochmessiah * Found a case that breaks Takagi * Fixes all the tests * Fixes issues found by the linter * Adds extra test * Adds extra test * dummy * dummy * Update symplectic.py Removes `autonne` * relocates test * trying to make pylint happy --------- Co-authored-by: Nicolas Quesada commit e999e4937c75977b5bee7ae33fbf5afb718aec35 Author: Nicolas Quesada <991946+nquesada@users.noreply.github.com> Date: Thu Aug 17 12:55:45 2023 -0400 Better way to determine if a matrix is a phase times a real matrix (#373) * Add a better way to determine if a matrix is a phase times a real matrix and tests * CHANGELOG updates --------- Co-authored-by: Nicolas Quesada --- .github/ACKNOWLEDGMENTS.md | 2 + .github/CHANGELOG.md | 42 ++++- .github/workflows/tests.yml | 2 + docs/references.bib | 14 ++ thewalrus/__init__.py | 8 + thewalrus/_montrealer.py | 134 ++++++++++++++ thewalrus/_torontonian.py | 49 +++--- thewalrus/_version.py | 2 +- thewalrus/decompositions.py | 232 +++++++++++++++---------- thewalrus/reference.py | 124 ++++++++++++- thewalrus/tests/test_decompositions.py | 187 +++++++++++++++++++- thewalrus/tests/test_montrealer.py | 181 +++++++++++++++++++ 12 files changed, 848 insertions(+), 129 deletions(-) create mode 100644 thewalrus/_montrealer.py create mode 100644 thewalrus/tests/test_montrealer.py diff --git a/.github/ACKNOWLEDGMENTS.md b/.github/ACKNOWLEDGMENTS.md index 6eec27778..833bfc3f1 100644 --- a/.github/ACKNOWLEDGMENTS.md +++ b/.github/ACKNOWLEDGMENTS.md @@ -67,3 +67,5 @@ * [Martin Houde](https://github.com/MHoude2) (Polytechnique Montréal) - 🙃 Minister of amplification * Will McCutcheon (Heriot-Watt University) - 🧅 Gaussian Onion Merchant + +* [Yanic Cardin](https://github.com/yaniccd) (Polytechnique Montréal) - 🦜 Pirate of the permutations diff --git a/.github/CHANGELOG.md b/.github/CHANGELOG.md index 4c0b38751..66dc415ad 100644 --- a/.github/CHANGELOG.md +++ b/.github/CHANGELOG.md @@ -1,19 +1,17 @@ -# Release 0.21.0-dev +# Release 0.22.0-dev ### New features -* Adds the Takagi decomposition [(#363)](https://github.com/XanaduAI/thewalrus/pull/338) +* Implements the pre-Iwasawa and Iwasawa decompositions for symplectic matrices [(#382)](https://github.com/XanaduAI/thewalrus/pull/382). ### Breaking changes ### Improvements +* Further simplifies the implementation of `decompositions.williamson` and corrects its docstring [(#380)](https://github.com/XanaduAI/thewalrus/pull/380). -* Tighten power-trace bound of odd loop Hafnian. [(#362)](https://github.com/XanaduAI/thewalrus/pull/362) - -* Simplifies the internal working of Bloch-Messiah decomposition [(#363)](https://github.com/XanaduAI/thewalrus/pull/338). +* Further simplifies the implementation of `decompositions.blochmessiah` [(#381)](https://github.com/XanaduAI/thewalrus/pull/381). -* Simplifies the internal working of Williamson decomposition [(#366)](https://github.com/XanaduAI/thewalrus/pull/338). ### Bug fixes @@ -21,9 +19,37 @@ ### Contributors -This release contains contributions from (in alphabetical order): +This release contains contributions from (in alphabetical order): + +Will McCutcheon, Nicolas Quesada + +--- + +# Release 0.21.0 + +### New features + +* Adds the Takagi decomposition [(#363)](https://github.com/XanaduAI/thewalrus/pull/363) + +* Adds the Montrealer and Loop Montrealer functions [(#363)](https://github.com/XanaduAI/thewalrus/pull/374). + +### Improvements + +* Tighten power-trace bound of odd loop Hafnian. [(#362)](https://github.com/XanaduAI/thewalrus/pull/362) + +* Simplifies the internal working of Bloch-Messiah decomposition [(#363)](https://github.com/XanaduAI/thewalrus/pull/338). + +* Simplifies the internal working of Williamson decomposition [(#366)](https://github.com/XanaduAI/thewalrus/pull/338). + +* Improves the handling of an edge case in Takagi [(#373)](https://github.com/XanaduAI/thewalrus/pull/373). + +* Adds extra tests for the Takagi decomposition [(#377)](https://github.com/XanaduAI/thewalrus/pull/377) + +### Contributors + +This release contains contributions from (in alphabetical order): -Gregory Morse, Nicolas Quesada +Yanic Cardin, Gregory Morse, Nicolas Quesada --- diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml index 0b3ab8198..cea2cb083 100644 --- a/.github/workflows/tests.yml +++ b/.github/workflows/tests.yml @@ -42,3 +42,5 @@ jobs: with: files: ./coverage.xml fail_ci_if_error: true + env: + CODECOV_TOKEN: ${{ secrets.CODECOV_TOKEN }} diff --git a/docs/references.bib b/docs/references.bib index cfc2a4be3..e5dd69776 100644 --- a/docs/references.bib +++ b/docs/references.bib @@ -642,3 +642,17 @@ @article{bulmer2022threshold journal={arXiv preprint arXiv:2202.04600}, year={2022} } + +@article{houde2024matrix, + title={Matrix decompositions in Quantum Optics: Takagi/Autonne, Bloch-Messiah/Euler, Iwasawa, and Williamson}, + author={Houde, Martin and McCutcheon, Will and Quesada, Nicol{\'a}s}, + journal={arXiv preprint arXiv:2403.04596}, + year={2024} +} + +@article{cardin2022photon, + title={Photon-number moments and cumulants of Gaussian states}, + author={Cardin, Yanic and Quesada, Nicol{\'a}s}, + journal={arXiv preprint arXiv:2212.06067}, + year={2022} +} diff --git a/thewalrus/__init__.py b/thewalrus/__init__.py index 87ec881bd..3e0c84c7c 100644 --- a/thewalrus/__init__.py +++ b/thewalrus/__init__.py @@ -132,6 +132,12 @@ rec_torontonian, rec_ltorontonian, ) + +from ._montrealer import ( + mtl, + lmtl, +) + from ._version import __version__ @@ -152,6 +158,8 @@ "reduction", "hermite_multidimensional", "grad_hermite_multidimensional", + "mtl", + "lmtl", "version", ] diff --git a/thewalrus/_montrealer.py b/thewalrus/_montrealer.py new file mode 100644 index 000000000..296d3056a --- /dev/null +++ b/thewalrus/_montrealer.py @@ -0,0 +1,134 @@ +""" +Montrealer Python interface + +* See Yanic Cardin and Nicolás Quesada. "Photon-number moments and cumulants of Gaussian states" + `arxiv:12212.06067 (2023) `_ :cite:`cardin2022photon` +""" +import numpy as np +import numba +from thewalrus.quantum.conversions import Xmat +from thewalrus.charpoly import powertrace +from ._hafnian import nb_ix +from ._torontonian import tor_input_checks + + +@numba.jit(nopython=True, cache=True) +def dec2bin(num, n): # pragma: no cover + """Helper function to convert an integer into an element of the power-set of ``n`` objects. + + Args: + num (int): label to convert + n (int): number of elements in the set + + Returns: + (array): array containing the labels of the elements to be selected + """ + digits = np.zeros((n), dtype=type(num)) + nn = num + counter = -1 + while nn >= 1: + digits[counter] = nn % 2 + counter -= 1 + nn //= 2 + return np.nonzero(digits)[0] + + +@numba.jit(nopython=True) +def montrealer(Sigma): # pragma: no cover + """Calculates the loop-montrealer of the zero-displacement Gaussian state with the given complex covariance matrix. + + Args: + Sigma (array): adjacency matrix of the Gaussian state + + Returns: + (np.complex128): the montrealer of ``Sigma`` + """ + n = len(Sigma) // 2 + tot_num = 2**n + val = np.complex128(0) + for p in numba.prange(tot_num): + pos = dec2bin(p, n) + lenpos = len(pos) + pos = np.concatenate((pos, n + pos)) + submat = nb_ix(Sigma, pos, pos) + sign = (-1) ** (lenpos + 1) + val += (sign) * powertrace(submat, n + 1)[-1] + return (-1) ** (n + 1) * val / (2 * n) + + +@numba.jit(nopython=True) +def power_loop(Sigma, zeta, n): # pragma: no cover + """Auxiliary function to calculate the product ``np.conj(zeta) @ Sigma^{n-1} @ zeta``. + + Args: + Sigma (array): square complex matrix + zeta (array): complex vector + n (int): sought after power + + Returns: + (np.complex128 or np.float64): the product np.conj(zeta) @ Sigma^{n-1} @ zeta + """ + vec = zeta + for _ in range(n - 1): + vec = Sigma @ vec + return np.conj(zeta) @ vec + + +@numba.jit(nopython=True, cache=True) +def lmontrealer(Sigma, zeta): # pragma: no cover + """Calculates the loop-montrealer of the displaced Gaussian state with the given complex covariance matrix and vector of displacements. + + Args: + Sigma (array): complex Glauber covariance matrix of the Gaussian state + zeta (array): vector of displacements + + Returns: + (np.complex128): the montrealer of ``Sigma`` + """ + n = len(Sigma) // 2 + tot_num = 2**n + val = np.complex128(0) + val_loops = np.complex128(0) + for p in numba.prange(tot_num): + pos = dec2bin(p, n) + lenpos = len(pos) + pos = np.concatenate((pos, n + pos)) + subvec = zeta[pos] + submat = nb_ix(Sigma, pos, pos) + sign = (-1) ** (lenpos + 1) + val_loops += sign * power_loop(submat, subvec, n) + val += sign * powertrace(submat, n + 1)[-1] + return (-1) ** (n + 1) * (val / (2 * n) + val_loops / 2) + + +def lmtl(A, zeta): + """Returns the montrealer of an NxN matrix and an N-length vector. + + Args: + A (array): an NxN array of even dimensions + zeta (array): an N-length vector of even dimensions + + Returns: + np.float64 or np.complex128: the loop montrealer of matrix A, vector zeta + """ + + tor_input_checks(A, zeta) + n = len(A) // 2 + Sigma = Xmat(n) @ A + return lmontrealer(Sigma, zeta) + + +def mtl(A): + """Returns the montrealer of an NxN matrix. + + Args: + A (array): an NxN array of even dimensions. + + Returns: + np.float64 or np.complex128: the montrealer of matrix ``A`` + """ + + tor_input_checks(A) + n = len(A) // 2 + Sigma = Xmat(n) @ A + return montrealer(Sigma) diff --git a/thewalrus/_torontonian.py b/thewalrus/_torontonian.py index 8626c46ba..1679a5130 100644 --- a/thewalrus/_torontonian.py +++ b/thewalrus/_torontonian.py @@ -20,19 +20,15 @@ from ._hafnian import reduction, find_kept_edges, nb_ix -def tor(A, recursive=True): - """Returns the Torontonian of a matrix. +def tor_input_checks(A, loops=None): + """Checks the correctness of the inputs for the torontonian/montrealer. Args: - A (array): a square array of even dimensions. - recursive: use the faster recursive implementation. - - Returns: - np.float64 or np.complex128: the torontonian of matrix A. + A (array): an NxN array of even dimensions + loops (array): optional argument, an N-length vector of even dimensions """ if not isinstance(A, np.ndarray): raise TypeError("Input matrix must be a NumPy array.") - matshape = A.shape if matshape[0] != matshape[1]: @@ -40,6 +36,25 @@ def tor(A, recursive=True): if matshape[0] % 2 != 0: raise ValueError("matrix dimension must be even") + + if loops is not None: + if not isinstance(loops, np.ndarray): + raise TypeError("Input matrix must be a NumPy array.") + if matshape[0] != len(loops): + raise ValueError("gamma must be a vector matching the dimension of A") + + +def tor(A, recursive=True): + """Returns the Torontonian of a matrix. + + Args: + A (array): a square array of even dimensions. + recursive: use the faster recursive implementation. + + Returns: + np.float64 or np.complex128: the torontonian of matrix A. + """ + tor_input_checks(A) return rec_torontonian(A) if recursive else numba_tor(A) @@ -54,23 +69,7 @@ def ltor(A, gamma, recursive=True): Returns: np.float64 or np.complex128: the loop torontonian of matrix A, vector gamma """ - - if not isinstance(A, np.ndarray): - raise TypeError("Input matrix must be a NumPy array.") - - if not isinstance(gamma, np.ndarray): - raise TypeError("Input matrix must be a NumPy array.") - - matshape = A.shape - - if matshape[0] != matshape[1]: - raise ValueError("Input matrix must be square.") - - if matshape[0] != len(gamma): - raise ValueError("gamma must be a vector matching the dimension of A") - - if matshape[0] % 2 != 0: - raise ValueError("matrix dimension must be even") + tor_input_checks(A, gamma) return rec_ltorontonian(A, gamma) if recursive else numba_ltor(A, gamma) diff --git a/thewalrus/_version.py b/thewalrus/_version.py index 8c079ffc5..f839ff5d4 100644 --- a/thewalrus/_version.py +++ b/thewalrus/_version.py @@ -16,4 +16,4 @@ Version number (major.minor.patch[-label]) """ -__version__ = "0.21.0-dev" +__version__ = "0.22.0-dev" diff --git a/thewalrus/decompositions.py b/thewalrus/decompositions.py index fdcd49c73..78f2575cf 100644 --- a/thewalrus/decompositions.py +++ b/thewalrus/decompositions.py @@ -22,6 +22,10 @@ This module implements common shared matrix decompositions that are used to perform gate decompositions. +For mathematical details of these decompositions see + +`Houde et al. Matrix decompositions in Quantum Optics: Takagi/Autonne, Bloch-Messiah/Euler, Iwasawa, and Williamson `_ +:cite:`houde2024matrix` Summary ------- @@ -30,13 +34,15 @@ symplectic_eigenvals takagi williamson + pre_iwasawa + iwasawa Code details ------------ """ import numpy as np -from scipy.linalg import block_diag, sqrtm, schur +from scipy.linalg import sqrtm, schur, polar, qr from thewalrus.symplectic import sympmat from thewalrus.quantum.gaussian_checks import is_symplectic @@ -44,8 +50,7 @@ def williamson(V, rtol=1e-05, atol=1e-08): r"""Williamson decomposition of positive-definite (real) symmetric matrix. - See `this thread `_ - and the `Williamson decomposition documentation `_ + See `Houde et al. Matrix decompositions in Quantum Optics: Takagi/Autonne, Bloch-Messiah/Euler, Iwasawa, and Williamson `_ Args: V (array[float]): positive definite symmetric (real) matrix @@ -54,7 +59,7 @@ def williamson(V, rtol=1e-05, atol=1e-08): Returns: tuple[array,array]: ``(Db, S)`` where ``Db`` is a diagonal matrix - and ``S`` is a symplectic matrix such that :math:`V = S^T Db S` + and ``S`` is a symplectic matrix such that :math:`V = S Db S^T` """ (n, m) = V.shape @@ -74,28 +79,27 @@ def williamson(V, rtol=1e-05, atol=1e-08): if not np.all(vals > 0): raise ValueError("Input matrix is not positive definite") - Mm12 = sqrtm(np.linalg.inv(V)).real - r1 = Mm12 @ omega @ Mm12 - s1, K = schur(r1) - # In what follows a permutation matrix perm1 is constructed so that the Schur matrix has + M12 = np.real_if_close(sqrtm(V)) + Mm12 = np.linalg.inv(M12) + Gamma = Mm12 @ omega @ Mm12 + a, O = schur(Gamma) + # In what follows a permutation matrix perm is constructed so that the Schur matrix has # only positive elements above the diagonal - # Also the Schur matrix uses the x_1,p_1, ..., x_n,p_n ordering thus a permutation perm2 is used + # Also the Schur matrix uses the x_1,p_1, ..., x_n,p_n ordering thus the permutation perm is updated # to go to the ordering x_1, ..., x_n, p_1, ... , p_n - perm1 = np.arange(2 * n) + perm = np.arange(2 * n) for i in range(n): - if s1[2 * i, 2 * i + 1] <= 0: - (perm1[2 * i], perm1[2 * i + 1]) = (perm1[2 * i + 1], perm1[2 * i]) - - perm2 = np.array([perm1[2 * i] for i in range(n)] + [perm1[2 * i + 1] for i in range(n)]) + if a[2 * i, 2 * i + 1] <= 0: + (perm[2 * i], perm[2 * i + 1]) = (perm[2 * i + 1], perm[2 * i]) - Ktt = K[:, perm2] - s1t = s1[:, perm1][perm1] + perm = np.array([perm[2 * i] for i in range(n)] + [perm[2 * i + 1] for i in range(n)]) - dd = np.array([1 / s1t[2 * i, 2 * i + 1] for i in range(n)]) - dd = np.concatenate([dd, dd]) - ddsqrt = np.sqrt(dd) - S = Mm12 @ Ktt * ddsqrt - return np.diag(dd), np.linalg.inv(S).T + O = O[:, perm] + phi = np.abs(np.diag(a, k=1)[::2]) + dd = np.concatenate([1 / phi, 1 / phi]) + ddsqrt = 1 / np.sqrt(dd) + S = M12 @ O * ddsqrt + return np.diag(dd), S def symplectic_eigenvals(cov): @@ -107,68 +111,53 @@ def symplectic_eigenvals(cov): Returns: (array): symplectic eigenvalues """ - M = int(len(cov) / 2) + M = len(cov) // 2 Omega = sympmat(M) return np.real_if_close(-1j * np.linalg.eigvals(Omega @ cov))[::2] def blochmessiah(S): - """Returns the Bloch-Messiah decomposition of a symplectic matrix S = uff @ dff @ vff - where uff and vff are orthogonal symplectic matrices and dff is a diagonal matrix - of the form diag(d1,d2,...,dn,d1^-1, d2^-1,...,dn^-1), + """Returns the Bloch-Messiah decomposition of a symplectic matrix S = O @ D @ Q + where O and Q are orthogonal symplectic matrices and D is a positive-definite diagonal matrix + of the form diag(d1,d2,...,dn,d1^-1, d2^-1,...,dn^-1). + + See `Houde et al. Matrix decompositions in Quantum Optics: Takagi/Autonne, Bloch-Messiah/Euler, Iwasawa, and Williamson `_ + Args: S (array[float]): 2N x 2N real symplectic matrix Returns: - tuple(array[float], : orthogonal symplectic matrix uff - array[float], : diagonal matrix dff - array[float]) : orthogonal symplectic matrix vff + tuple(array[float], : orthogonal symplectic matrix O + array[float], : diagonal matrix D + array[float]) : orthogonal symplectic matrix Q """ N, _ = S.shape if not is_symplectic(S): raise ValueError("Input matrix is not symplectic.") - - # Changing Basis - R = (1 / np.sqrt(2)) * np.block( - [[np.eye(N // 2), 1j * np.eye(N // 2)], [np.eye(N // 2), -1j * np.eye(N // 2)]] - ) - Sc = R @ S @ np.conjugate(R).T - # Polar Decomposition - u1, d1, v1 = np.linalg.svd(Sc) - Sig = u1 @ np.diag(d1) @ np.conjugate(u1).T - Unitary = u1 @ v1 - # Blocks of Unitary and Hermitian symplectics - alpha = Unitary[0 : N // 2, 0 : N // 2] - beta = Sig[0 : N // 2, N // 2 : N] - # Bloch-Messiah in this Basis - d2, takagibeta = takagi(beta) - sval = np.arcsinh(d2) - uf = block_diag(takagibeta, takagibeta.conj()) - blc = np.conjugate(takagibeta).T @ alpha - vf = block_diag(blc, blc.conj()) - df = np.block( - [ - [np.diag(np.cosh(sval)), np.diag(np.sinh(sval))], - [np.diag(np.sinh(sval)), np.diag(np.cosh(sval))], - ] - ) - # Rotating Back to Original Basis - uff = np.conjugate(R).T @ uf @ R - vff = np.conjugate(R).T @ vf @ R - dff = np.conjugate(R).T @ df @ R - dff = np.real_if_close(dff) - vff = np.real_if_close(vff) - uff = np.real_if_close(uff) - return uff, dff, vff + N = N // 2 + V, P = polar(S, side="left") + A = P[:N, :N] + B = P[:N, N:] + C = P[N:, N:] + M = A - C + 1j * (B + B.T) + Lam, W = takagi(M) + Lam = 0.5 * Lam + O = np.block([[W.real, -W.imag], [W.imag, W.real]]) + Q = O.T @ V + sqrt1pLam2 = np.sqrt(1 + Lam**2) + D = np.diag(np.concatenate([sqrt1pLam2 + Lam, sqrt1pLam2 - Lam])) + return O, D, Q def takagi(A, svd_order=True): r"""Autonne-Takagi decomposition of a complex symmetric (not Hermitian!) matrix. - Note that the input matrix is internally symmetrized. If the input matrix is indeed symmetric this leaves it unchanged. - See `Carl Caves note. `_ + Note that the input matrix is internally symmetrized by taking its upper triangular part. + If the input matrix is indeed symmetric this leaves it unchanged. + + See `Houde et al. Matrix decompositions in Quantum Optics: Takagi/Autonne, Bloch-Messiah/Euler, Iwasawa, and Williamson `_ Args: A (array): square, symmetric matrix @@ -182,8 +171,8 @@ def takagi(A, svd_order=True): n, m = A.shape if n != m: raise ValueError("The input matrix is not square") - # Here we force symmetrize the matrix - A = 0.5 * (A + A.T) + # Here we build a Symmetric matrix from the top right triangular part + A = np.triu(A) + np.triu(A, k=1).T A = np.real_if_close(A) @@ -193,26 +182,21 @@ def takagi(A, svd_order=True): if np.isrealobj(A): # If the matrix A is real one can be more clever and use its eigendecomposition ls, U = np.linalg.eigh(A) - U = U / np.exp(1j * np.angle(U)[0]) vals = np.abs(ls) # These are the Takagi eigenvalues - phases = -np.ones(vals.shape[0], dtype=np.complex128) - for j, l in enumerate(ls): - if np.allclose(l, 0) or l > 0: - phases[j] = 1 - phases = np.sqrt(phases) - Uc = U @ np.diag(phases) # One needs to readjust the phases - signs = np.sign(Uc.real)[0] - for k, s in enumerate(signs): - if np.allclose(s, 0): - signs[k] = 1 - Uc = np.real_if_close(Uc / signs) - list_vals = [(vals[i], i) for i in range(len(vals))] - # And also rearrange the unitary and values so that they are decreasingly ordered - list_vals.sort(reverse=svd_order) - sorted_ls, permutation = zip(*list_vals) - return np.array(sorted_ls), Uc[:, np.array(permutation)] - - phi = np.angle(A[0, 0]) + signs = (-1) ** (1 + np.heaviside(ls, 1)) + phases = np.sqrt(np.complex128(signs)) + Uc = U * phases # One needs to readjust the phases + # Find the permutation to sort in decreasing order + perm = np.argsort(vals) + # if svd_order reverse it + if svd_order: + perm = perm[::-1] + return vals[perm], Uc[:, perm] + + # Find the element with the largest absolute value + pos = np.unravel_index(np.argmax(np.abs(A)), (n, n)) + # Use it to find whether the input is a global phase times a real matrix + phi = np.angle(A[pos]) Amr = np.real_if_close(np.exp(-1j * phi) * A) if np.isrealobj(Amr): vals, U = takagi(Amr, svd_order=svd_order) @@ -220,10 +204,82 @@ def takagi(A, svd_order=True): u, d, v = np.linalg.svd(A) U = u @ sqrtm((v @ np.conjugate(u)).T) - # The line above could be simplifed to the line below if the product v @ np.conjugate(u) is diagonal - # Which it should be according to Caves http://info.phys.unm.edu/~caves/courses/qinfo-s17/lectures/polarsingularAutonne.pdf - # U = u * np.sqrt(0j + np.diag(np.conjugate(u) @ v)) - # This however breaks test_degenerate if svd_order is False: return d[::-1], U[:, ::-1] return d, U + + +def pre_iwasawa(S): + """Pre-Iwasawa decomposition of a symplectic matrix. + See `Arvind et al. The Real Symplectic Groups in Quantum Mechanics and Optics `_ + and `Houde et al. Matrix decompositions in Quantum Optics: Takagi/Autonne, Bloch-Messiah/Euler, Iwasawa, and Williamson `_ + + + Args: + S (array): the symplectic matrix + + Returns: + tuple[array, array, array]: (E,D,F) symplectic matrices such that E @ D @ F = S and, + E = np.block([[np.eye(N), np.zeros(N,N)],[X, np.eye(N)]]) with X == X.T, + D is block diagonal with the top left block being the inverse of the bottom right block, + F is symplectic orthogonal. + """ + + if not is_symplectic(S): + raise ValueError("Input matrix is not symplectic.") + + N, _ = S.shape + N = N // 2 + zerom = np.zeros([N, N]) + idm = np.eye(N) + A = S[:N, :N] + B = S[:N, N:] + C = S[N:, :N] + D = S[N:, N:] + A0 = sqrtm(A @ A.T + B @ B.T) + A0inv = np.linalg.inv(A0) + X = A0inv @ A + Y = A0inv @ B + C0 = (C @ A.T + D @ B.T) @ A0inv + E = np.block([[idm, zerom], [C0 @ A0inv, idm]]) + D = np.block([[A0, zerom], [zerom, A0inv]]) + F = np.block([[X, Y], [-Y, X]]) + return E, D, F + + +def iwasawa(S): + """Iwasawa decomposition of a symplectic matrix. + See `Arvind et al. The Real Symplectic Groups in Quantum Mechanics and Optics `_ + and `Houde et al. Matrix decompositions in Quantum Optics: Takagi/Autonne, Bloch-Messiah/Euler, Iwasawa, and Williamson `_ + + + Args: + S (array): the symplectic matrix + + Returns: + tuple[array, array, array]: (E,D,F) symplectic matrices such that E @ D @ F = S, + EE = np.block([[AA, np.zeros(N,N)],[CC, np.linalg.inv(A.T)]]) with A.T @ C == C.T @ A, and AA upper trinagular with ones in the diagonal + DD is diagonal and symplectic, + FF is symplectic orthogonal. + """ + + E, D, F = pre_iwasawa(S) + N, _ = S.shape + N = N // 2 + DNN = D[:N, :N] + Q, R = qr(DNN) + R = R.T + Q = Q.T + dR = np.diag(R) + dd = np.abs(dR) + ds = np.sign(dR) + R = R * (1 / dR) + RinvT = np.linalg.inv(R).T + DD = np.diag(np.concatenate([dd, 1 / dd])) + zerom = np.zeros([N, N]) + OO = np.block([[R, zerom], [zerom, RinvT]]) + Q = ds[:, None] * Q + AA = np.block([[Q, zerom], [zerom, Q]]) + EE = E @ OO + FF = AA @ F + return EE, DD, FF diff --git a/thewalrus/reference.py b/thewalrus/reference.py index 00d7934eb..9bd2e0417 100644 --- a/thewalrus/reference.py +++ b/thewalrus/reference.py @@ -34,12 +34,14 @@ .. autosummary:: hafnian + montrealer Code details ------------ .. autofunction:: hafnian + montrealer Auxiliary functions ------------------- @@ -49,6 +51,8 @@ partitions spm pmp + rspm + rpmp T Code details @@ -58,7 +62,7 @@ # pylint: disable=too-many-arguments from collections import OrderedDict -from itertools import tee +from itertools import tee, product, permutations, chain from types import GeneratorType MAXSIZE = 1000 @@ -278,3 +282,121 @@ def hafnian(M, loop=False): tot_sum = tot_sum + result return tot_sum + + +def mapper(x, objects): + """Helper function to turn a permutation and bistring into an element of rpmp. + + Args: + x (tuple): tuple containing a permutation and a bistring + objects (list): list objects to permute + + Returns: + tuple: permuted objects + """ + (perm, bit) = x + m = len(bit) + Blist = [list(range(m)), list(range(m, 2 * m))] + for i, j in enumerate(bit): + if int(j): + (Blist[0][i], Blist[1][i]) = (Blist[1][i], Blist[0][i]) + Blist = [Blist[0][i] for i in tuple((0,) + perm)] + [Blist[1][i] for i in tuple((0,) + perm)] + dico_list = {j: i + 1 for i, j in enumerate(Blist)} + new_mapping_list = { + objects[dico_list[i] - 1]: objects[dico_list[j] - 1] + for i, j in zip(list(range(0, m - 1)) + [m], list(range(m + 1, 2 * m)) + [m - 1]) + } + return tuple(new_mapping_list.items()) + + +def bitstrings(n): + """Returns the bistrings from 0 to n/2 + + Args: + n (int): Twice the highest bitstring value. + + Returns: + (iterator): An iterable of all bistrings. + """ + for binary in map("".join, product("01", repeat=n - 1)): + yield "0" + binary + + +def rpmp(s): + """Generates the restricted set of perfect matchings matching permutations. + + Args: + s (tuple): tuple of labels to be used + + Returns: + generator: the set of restricted perfect matching permutations of the tuple ``s`` + """ + m = len(s) // 2 + + def local_mapper(x): + """Helper function to define a local mapper based on the symbols s + Args: + x (iterable): object to be mapped + """ + return mapper(x, s) + + for i in product(permutations(range(1, m)), bitstrings(m)): + yield local_mapper(i) + + +def splitter(elem): + """Takes an element from the restricted perfect matching permutations (rpmp) and returns all the associated elements in the restricted single pair matchings (rspm) + + Args: + elem (tuple): tuple representing an element of rpmp + + Returns: + (iterator): all the associated elements in rspm + """ + num_elem = len(elem) + net = [elem] + for i in range(num_elem): + left = (elem[j] for j in range(i)) + middle_left = ((elem[i][0], elem[i][0]),) + middle_right = ((elem[i][1], elem[i][1]),) + right = (elem[j] for j in range(i + 1, num_elem)) + net.append(tuple(middle_right) + tuple(right) + tuple(left) + tuple(middle_left)) + for i in net: + yield i + + +def rspm(s): + """Generates the restricted set of single-pair matchings. + + Args: + s (tuple): tuple of labels to be used + + Returns: + generator: the set of restricted perfect matching permutations of the tuple s + """ + gen = rpmp(s) + return chain(*(splitter(i) for i in gen)) + + +def mtl(A, loop=False): + """Returns the Montrealer of an NxN matrix and an N-length vector. + + Args: + A (array): an NxN array of even dimensions. Can be symbolic. + loop (boolean): if set to ``True``, the loop montrealer is returned + + Returns: + np.float64, np.complex128 or sympy.core.add.Add: the Montrealer of matrix ``A``. + """ + n, _ = A.shape + net_sum = 0 + + perm = rspm(range(n)) if loop else rpmp(range(n)) + for s in perm: + net_prod = 1 + for a in s: + net_prod *= A[a[0], a[1]] + + net_sum += net_prod + + return net_sum diff --git a/thewalrus/tests/test_decompositions.py b/thewalrus/tests/test_decompositions.py index e6684df18..e4b75e234 100644 --- a/thewalrus/tests/test_decompositions.py +++ b/thewalrus/tests/test_decompositions.py @@ -19,7 +19,7 @@ from thewalrus.random import random_interferometer as haar_measure from thewalrus.random import random_symplectic -from thewalrus.decompositions import williamson, blochmessiah, takagi +from thewalrus.decompositions import williamson, blochmessiah, takagi, pre_iwasawa, iwasawa from thewalrus.symplectic import sympmat as omega from thewalrus.quantum.gaussian_checks import is_symplectic @@ -274,23 +274,30 @@ def test_takagi(n, datatype, svd_order): assert np.all(np.diff(r) >= 0) +# pylint: disable=too-many-arguments @pytest.mark.parametrize("n", [5, 10, 50]) @pytest.mark.parametrize("datatype", [np.complex128, np.float64]) @pytest.mark.parametrize("svd_order", [True, False]) @pytest.mark.parametrize("half_rank", [0, 1]) @pytest.mark.parametrize("phase", [0, 1]) -def test_degenerate(n, datatype, svd_order, half_rank, phase): +@pytest.mark.parametrize("null_space", [0, 5, 10]) +@pytest.mark.parametrize("offset", [0, 0.5]) +def test_degenerate(n, datatype, svd_order, half_rank, phase, null_space, offset): """Tests Takagi produces the correct result for very degenerate cases""" nhalf = n // 2 - diags = [half_rank * np.random.rand()] * nhalf + [np.random.rand()] * (n - nhalf) + diags = ( + [half_rank * np.random.rand()] * nhalf + + [np.random.rand() - offset] * (n - nhalf) + + [0] * null_space + ) if datatype is np.complex128: - U = haar_measure(n) + U = haar_measure(n + null_space) if datatype is np.float64: - U = np.exp(1j * phase) * haar_measure(n, real=True) + U = np.exp(1j * phase) * haar_measure(n + null_space, real=True) A = U @ np.diag(diags) @ U.T r, U = takagi(A, svd_order=svd_order) assert np.allclose(A, U @ np.diag(r) @ U.T) - assert np.allclose(U @ U.T.conj(), np.eye(n)) + assert np.allclose(U @ U.T.conj(), np.eye(n + null_space)) assert np.all(r >= 0) if svd_order is True: assert np.all(np.diff(r) <= 0) @@ -394,3 +401,171 @@ def test_real_degenerate(): rl, U = takagi(mat) assert np.allclose(U @ U.conj().T, np.eye(len(mat))) assert np.allclose(U @ np.diag(rl) @ U.T, mat) + + +@pytest.mark.parametrize("n", [5, 10, 50]) +@pytest.mark.parametrize("datatype", [np.complex128, np.float64]) +@pytest.mark.parametrize("svd_order", [True, False]) +def test_autonne_takagi(n, datatype, svd_order): + """Checks the correctness of the Autonne decomposition function""" + if datatype is np.complex128: + A = np.random.rand(n, n) + 1j * np.random.rand(n, n) + if datatype is np.float64: + A = np.random.rand(n, n) + A += A.T + r, U = takagi(A, svd_order=svd_order) + assert np.allclose(A, U @ np.diag(r) @ U.T) + assert np.all(r >= 0) + if svd_order is True: + assert np.all(np.diff(r) <= 0) + else: + assert np.all(np.diff(r) >= 0) + + +@pytest.mark.parametrize("size", [10, 20, 100]) +def test_flat_phase(size): + """Test that the correct decomposition is obtained even if the first entry is 0""" + A = np.random.rand(size, size) + 1j * np.random.rand(size, size) + A += A.T + A[0, 0] = 0 + l, u = takagi(A) + assert np.allclose(A, u * l @ u.T) + + +def test_real_input_edge(): + """Adapted from https://math.stackexchange.com/questions/4418925/why-does-this-algorithm-for-the-takagi-factorization-fail-here""" + rng = np.random.default_rng(0) # Important for reproducibility + A = (rng.random((100, 100)) - 0.5) * 114 + A = A * A.T # make A symmetric + l, u = takagi(A) + # Now, reconstruct A, see + Ar = u * l @ u.T + assert np.allclose(A, Ar) + + +@pytest.mark.parametrize("rank1", [2, 4, 5]) +@pytest.mark.parametrize("rank2", [2, 4, 5]) +@pytest.mark.parametrize("rankrand", [2, 4, 5]) +@pytest.mark.parametrize("rankzero", [2, 4, 5]) +@pytest.mark.parametrize("symmetric", [True, False]) +@pytest.mark.parametrize("unitary", [True, False]) +def test_pre_iwasawa(rank1, rank2, rankrand, rankzero, symmetric, unitary): + """Tests the pre_iwasawa decomposition""" + vals = np.array( + [np.random.rand(1)[0]] * rank1 + + [np.random.rand(1)[0]] * rank2 + + list(np.random.rand(rankrand)) + + [1] * rankzero + ) + if unitary is True: + vals = np.ones_like(vals) + dd = np.concatenate([vals, 1 / vals]) + dim = len(vals) + U = haar_measure(dim) + O = np.block([[U.real, -U.imag], [U.imag, U.real]]) + if symmetric is False: + V = haar_measure(dim) + P = np.block([[V.real, -V.imag], [V.imag, V.real]]) + else: + P = O.T + + S = (O * dd) @ P + EE, DD, FF = pre_iwasawa(S) + assert np.allclose(EE @ DD @ FF, S) + assert is_symplectic(EE) + assert is_symplectic(FF) + assert is_symplectic(FF) + assert np.allclose(FF @ FF.T, np.identity(2 * dim)) + assert np.allclose(DD[:dim, :dim] @ DD[dim:, dim:], np.identity(dim)) + assert np.allclose(DD[:dim, dim:], 0) + assert np.allclose(DD[dim:, :dim], 0) + A = EE[:dim, :dim] + B = EE[:dim, dim:] + C = EE[dim:, :dim] + D = EE[dim:, dim:] + assert np.allclose(A, np.eye(dim)) + assert np.allclose(B, 0) + assert np.allclose(C, C.T) + assert np.allclose(D, np.eye(dim)) + + +@pytest.mark.parametrize("rank1", [2, 4, 5]) +@pytest.mark.parametrize("rank2", [2, 4, 5]) +@pytest.mark.parametrize("rankrand", [2, 4, 5]) +@pytest.mark.parametrize("rankzero", [2, 4, 5]) +@pytest.mark.parametrize("symmetric", [True, False]) +@pytest.mark.parametrize("unitary", [True, False]) +def test_iwasawa(rank1, rank2, rankrand, rankzero, symmetric, unitary): + """Tests the Iwasawa decomposition""" + vals = np.array( + [np.random.rand(1)[0]] * rank1 + + [np.random.rand(1)[0]] * rank2 + + list(np.random.rand(rankrand)) + + [1] * rankzero + ) + if unitary is True: + vals = np.ones_like(vals) + dd = np.concatenate([vals, 1 / vals]) + dim = len(vals) + U = haar_measure(dim) + O = np.block([[U.real, -U.imag], [U.imag, U.real]]) + if symmetric is False: + V = haar_measure(dim) + P = np.block([[V.real, -V.imag], [V.imag, V.real]]) + else: + P = O.T + S = (O * dd) @ P + EE, DD, FF = iwasawa(S) + assert np.allclose(EE @ DD @ FF, S) + assert is_symplectic(EE) + assert is_symplectic(FF) + assert is_symplectic(FF) + assert np.allclose(FF @ FF.T, np.identity(2 * dim)) + assert np.allclose(DD, np.diag(np.diag(DD))) + assert np.allclose(DD[:dim, :dim] @ DD[dim:, dim:], np.identity(dim)) + A = EE[:dim, :dim] + B = EE[:dim, dim:] + C = EE[dim:, :dim] + D = EE[dim:, dim:] + assert np.allclose(B, 0) + XX = A.T @ C + assert np.allclose(XX, XX.T) + assert np.allclose(A @ D.T, np.eye(dim)) + assert np.allclose(np.diag(EE), 1) + assert np.allclose(np.tril(A), A) + assert np.allclose(np.triu(D), D) + + +def test_pre_iwasawa_error(): + """Tests error is raised when input not symplectic""" + M = np.random.rand(4, 5) + with pytest.raises(ValueError, match="Input matrix is not symplectic."): + pre_iwasawa(M) + + +def test_iwasawa_error(): + """Tests error is raised when input not symplectic""" + M = np.random.rand(4, 5) + with pytest.raises(ValueError, match="Input matrix is not symplectic."): + iwasawa(M) + + +def test_iwasawa2x2(): + """Compares numerics against exact result for 2x2 matrices in Arvind 1995""" + num_tests = 100 + for _ in range(num_tests): + S = random_symplectic(1) + A, N, K = iwasawa(S) + a = S[0, 0] + b = S[0, 1] + c = S[1, 0] + d = S[1, 1] + eta = a**2 + b**2 + xi = (a * c + b * d) / eta + eta = np.sqrt(eta) + AA = np.array([[1, 0], [xi, 1]]) + NN = np.diag([eta, 1 / eta]) + KK = np.array([[a, b], [-b, a]]) / eta + assert np.allclose(A, AA) + assert np.allclose(K, KK) + assert np.allclose(N, NN) diff --git a/thewalrus/tests/test_montrealer.py b/thewalrus/tests/test_montrealer.py new file mode 100644 index 000000000..f511d0fc3 --- /dev/null +++ b/thewalrus/tests/test_montrealer.py @@ -0,0 +1,181 @@ +# Copyright 2021 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 adj 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. +"""Montrealer tests +Yanic Cardin and Nicolás Quesada. "Photon-number moments and cumulants of Gaussian states" +`arxiv:12212.06067 (2023) `_ +""" + +import pytest +import numpy as np +from thewalrus import mtl, lmtl +from thewalrus.reference import mapper +from thewalrus.quantum import Qmat, Xmat +from thewalrus.reference import rspm, rpmp, mtl as mtl_symb +from thewalrus.random import random_covariance +from scipy.special import factorial2 +from scipy.stats import unitary_group + + +@pytest.mark.parametrize("n", range(1, 8)) +def test_montrealer_all_ones(n): + """Test that the Montrealer of a matrix of ones gives (2n-2)!!""" + adj = np.ones([2 * n, 2 * n]) + mtl_val = mtl(adj) + mtl_expect = factorial2(2 * n - 2) + assert np.allclose(mtl_val, mtl_expect) + + +@pytest.mark.parametrize("n", range(1, 8)) +def test_loop_montrealer_all_ones(n): + """Test that the loop Montrealer of a matrix of ones gives (n+1)(2n-2)!!""" + adj = np.ones([2 * n, 2 * n]) + lmtl_val = lmtl(adj, zeta=np.diag(adj)) + lmtl_expect = (n + 1) * factorial2(2 * n - 2) + assert np.allclose(lmtl_val, lmtl_expect) + + +@pytest.mark.parametrize("n", range(1, 8)) +def test_size_of_rpmp(n): + """rpmp(2n) should have (2n-2)!! elements""" + terms_rpmp = len(list(rpmp(range(2 * n)))) + terms_theo = factorial2(2 * n - 2) + assert terms_rpmp == terms_theo + + +@pytest.mark.parametrize("n", range(1, 8)) +def test_size_of_rspm(n): + """rspm(2n) should have (n+1)(2n-2)!! elements""" + terms_rspm = sum(1 for _ in rspm(range(2 * n))) + terms_theo = (n + 1) * factorial2(2 * n - 2) + assert terms_rspm == terms_theo + + +@pytest.mark.parametrize("n", range(2, 8)) +def test_rpmp_alternating_walk(n): + """The rpmp must form a Y-alternating walk without loops""" + test = True + for perfect in rpmp(range(1, 2 * n + 1)): + last = perfect[0][1] # starting point + reduced_last = last - n if last > n else last + # different mode in every tuple + if reduced_last == 1: + test = False + + for i in perfect[1:]: + reduced = i[0] - n if i[0] > n else i[0], i[1] - n if i[1] > n else i[1] + # different mode in every tuple + if reduced[0] == reduced[1]: + test = False + # consecutive tuple contain the same mode + if reduced_last not in reduced: + test = False + + last = i[0] if reduced[1] == reduced_last else i[1] + reduced_last = last - n if last > n else last + + # last mode most coincide with the first one + if reduced_last != 1: + test = False + + assert test + + +@pytest.mark.parametrize("n", range(1, 8)) +def test_mtl_functions_agree(n): + """Make sure both mtl functions agree with one another""" + V = random_covariance(n) + Aad = Xmat(n) @ (Qmat(V) - np.identity(2 * n)) + assert np.allclose(mtl_symb(Aad), mtl(Aad)) + + +@pytest.mark.parametrize("n", range(1, 8)) +def test_lmtl_functions_agree(n): + """Make sure both lmtl functions agree with one another""" + V = random_covariance(n) + Aad = Xmat(n) @ (Qmat(V) - np.identity(2 * n)) + zeta = np.diag(Aad).conj() + assert np.allclose(lmtl(Aad, zeta), mtl_symb(Aad, loop=True)) + + +@pytest.mark.parametrize("n", range(1, 8)) +def test_mtl_lmtl_agree(n): + """Make sure mtl and lmtl give the same result if zeta = 0""" + V = random_covariance(n) + Aad = Xmat(n) @ (Qmat(V) - np.identity(2 * n)) + zeta = np.zeros(2 * n, dtype=np.complex128) + assert np.allclose(lmtl(Aad, zeta), lmtl(Aad, zeta)) + + +@pytest.mark.parametrize("n", range(1, 8)) +def test_mtl_lmtl_reference_agree(n): + """Make sure mtl and lmtl from .reference give the same result if zeta = 0""" + V = random_covariance(n) + Aad = Xmat(n) @ (Qmat(V) - np.identity(2 * n)) + zeta = np.zeros(2 * n, dtype=np.complex128) + np.fill_diagonal(Aad, zeta) + assert np.allclose(mtl_symb(Aad, loop=True), mtl_symb(Aad)) + + +@pytest.mark.parametrize("n", range(1, 8)) +def test_mtl_permutation(n): + """Make sure the mtl is invariant under permutation + cf. Eq. 44 of `arxiv:12212.06067 (2023) `_""" + V = random_covariance(n) + Aad = Xmat(n) @ (Qmat(V) - np.identity(2 * n)) + perm = np.random.permutation(n) + perm = np.concatenate((perm, [i + n for i in perm])) + assert np.allclose(mtl(Aad), mtl(Aad[perm][:, perm])) + + +@pytest.mark.parametrize("n", range(2, 5)) +def test_mtl_associated_adjacency(n): + """Make sure the mtl of a matrix in which each block is block diaognal is zero. + cf. Eq. 45 of `arxiv:12212.06067 (2023) `_""" + u_zero = np.zeros((n, n), dtype=np.complex128) + + u_n1 = unitary_group.rvs(n) + u_n2 = unitary_group.rvs(n) + u_n = np.block([[u_n1, u_zero], [u_zero, u_n2]]) + u_n = u_n + u_n.conj().T + + u_m1 = unitary_group.rvs(n) + u_m2 = unitary_group.rvs(n) + u_m = np.block([[u_m1, u_zero], [u_zero, u_m2]]) + u_m_r = u_m + u_m.T + + u_m3 = unitary_group.rvs(n) + u_m4 = unitary_group.rvs(n) + u_m = np.block([[u_m3, u_zero], [u_zero, u_m4]]) + u_m_l = u_m + u_m.T + + adj = np.block([[u_m_r, u_n], [u_n.T, u_m_l]]) + + assert np.allclose(mtl(adj), 0) + + +@pytest.mark.parametrize("n", range(1, 8)) +def test_mtl_diagonal_trace(n): + """Make sure the mtl of A times a diagonal matrix gives the product of the norms of the diagonal matrix times the mtl of A + cf. Eq. 41 of `arxiv:12212.06067 (2023) `_""" + gamma = np.random.uniform(-1, 1, n) + 1.0j * np.random.uniform(-1, 1, n) + product = np.prod([abs(i) ** 2 for i in gamma]) + gamma = np.diag(np.concatenate((gamma, gamma.conj()))) + V = random_covariance(n) + Aad = Xmat(n) @ (Qmat(V) - np.identity(2 * n)) + assert np.allclose(mtl(gamma @ Aad @ gamma), product * mtl(Aad)) + + +def test_mapper_hard_coded(): + """Tests the the mapper function for a particular hardcoded value""" + assert mapper(((1, 2, 3), "0000"), (0, 1, 2, 3, 4, 5, 6, 7)) == ((0, 5), (1, 6), (2, 7), (4, 3))