Squashed commit of the following:

commit 7c63905 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 f935053 Author: Alberto Fumagalli <alberto@xanadu.ai> Date: Mon Apr 29 14:14:51 2024 -0400 Add codecov token to test workflow (#390) commit fe9e112 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 <nquesada@pop-os.localdomain> commit 25b5f08 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 <nquesada@pop-os.localdomain> Co-authored-by: Sebastián Duque Mesa <675763+sduquemesa@users.noreply.github.com> commit 127a1a7 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 8759363 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 <nquesada@pop-os.localdomain> commit 6247fc8 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 <nquesada@pop-os.localdomain> commit 0b1af63 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 544c457 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 2268bf0 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 <nquesada@pop-os.localdomain> commit 8ec2172 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 <nquesada@pop-os.localdomain> * 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 <nquesada@pop-os.localdomain> Co-authored-by: Yanic Cardin <yanic.cardin.1@gmail.com> Co-authored-by: Sebastián Duque Mesa <675763+sduquemesa@users.noreply.github.com> commit 0e163bf 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 <nquesada@pop-os.localdomain> commit e999e49 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 <nquesada@pop-os.localdomain>
XanaduAI · Apr 29, 2024 · 0f99b20 · 0f99b20
1 parent 1ed51d9
commit 0f99b20
Show file tree

Hide file tree

Showing 12 changed files with 848 additions and 129 deletions.
diff --git 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
@@ -1,29 +1,55 @@
-# 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
 
 ### Documentation
 
 ### 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
@@ -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
@@ -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
@@ -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
@@ -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) <https://arxiv.org/abs/2212.06067>`_ :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
@@ -20,26 +20,41 @@
 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]:
         raise ValueError("Input matrix must be square.")
 
     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
@@ -16,4 +16,4 @@
    Version number (major.minor.patch[-label])
 """
 
-__version__ = "0.21.0-dev"
+__version__ = "0.22.0-dev"