Skip to content

Commit 8ec2172

Browse files
nquesadaNicolas Quesadayaniccdsduquemesa
authored
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 <[email protected]> * making bots happy * Apply suggestions from code review Co-authored-by: Sebastián Duque Mesa <[email protected]> * Apply suggestions from code review Co-authored-by: Sebastián Duque Mesa <[email protected]> * Apply suggestions from code review * Apply suggestions from code review Co-authored-by: Sebastián Duque Mesa <[email protected]> * Apply suggestions from code review Co-authored-by: Sebastián Duque Mesa <[email protected]> * Adds extra tests --------- Co-authored-by: Nicolas Quesada <[email protected]> Co-authored-by: Yanic Cardin <[email protected]> Co-authored-by: Sebastián Duque Mesa <[email protected]>
1 parent 0e163bf commit 8ec2172

File tree

7 files changed

+475
-28
lines changed

7 files changed

+475
-28
lines changed

‎.github/ACKNOWLEDGMENTS.md‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -67,3 +67,5 @@
6767
* [Martin Houde](https://github.com/MHoude2) (Polytechnique Montréal) - 🙃 Minister of amplification
6868

6969
* Will McCutcheon (Heriot-Watt University) - 🧅 Gaussian Onion Merchant
70+
71+
* [Yanic Cardin](https://github.com/yaniccd) (Polytechnique Montréal) - 🦜 Pirate of the permutations

‎.github/CHANGELOG.md‎

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -2,8 +2,9 @@
22

33
### New features
44

5-
* Adds the Takagi decomposition [(#363)](https://github.com/XanaduAI/thewalrus/pull/338)
5+
* Adds the Takagi decomposition [(#363)](https://github.com/XanaduAI/thewalrus/pull/363)
66

7+
* Adds the Montrealer and Loop Montrealer functions [(#363)](https://github.com/XanaduAI/thewalrus/pull/374).
78
### Breaking changes
89

910
### Improvements
@@ -25,7 +26,7 @@
2526

2627
This release contains contributions from (in alphabetical order):
2728

28-
Gregory Morse, Nicolas Quesada
29+
Yanic Cardin, Gregory Morse, Nicolas Quesada
2930

3031
---
3132

‎thewalrus/__init__.py‎

Lines changed: 8 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -132,6 +132,12 @@
132132
rec_torontonian,
133133
rec_ltorontonian,
134134
)
135+
136+
from ._montrealer import (
137+
mtl,
138+
lmtl,
139+
)
140+
135141
from ._version import __version__
136142

137143

@@ -152,6 +158,8 @@
152158
"reduction",
153159
"hermite_multidimensional",
154160
"grad_hermite_multidimensional",
161+
"mtl",
162+
"lmtl",
155163
"version",
156164
]
157165

‎thewalrus/_montrealer.py‎

Lines changed: 134 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,134 @@
1+
"""
2+
Montrealer Python interface
3+
4+
* Yanic Cardin and Nicolás Quesada. "Photon-number moments and cumulants of Gaussian states"
5+
`arxiv:12212.06067 (2023) <https://arxiv.org/abs/2212.06067>`_
6+
"""
7+
import numpy as np
8+
import numba
9+
from thewalrus.quantum.conversions import Xmat
10+
from thewalrus.charpoly import powertrace
11+
from ._hafnian import nb_ix
12+
from ._torontonian import tor_input_checks
13+
14+
15+
@numba.jit(nopython=True, cache=True)
16+
def dec2bin(num, n): # pragma: no cover
17+
"""Helper function to convert an integer into an element of the power-set of ``n`` objects.
18+
19+
Args:
20+
num (int): label to convert
21+
n (int): number of elements in the set
22+
23+
Returns:
24+
(array): array containing the labels of the elements to be selected
25+
"""
26+
digits = np.zeros((n), dtype=type(num))
27+
nn = num
28+
counter = -1
29+
while nn >= 1:
30+
digits[counter] = nn % 2
31+
counter -= 1
32+
nn //= 2
33+
return np.nonzero(digits)[0]
34+
35+
36+
@numba.jit(nopython=True)
37+
def montrealer(Sigma): # pragma: no cover
38+
"""Calculates the loop-montrealer of the zero-displacement Gaussian state with the given complex covariance matrix.
39+
40+
Args:
41+
Sigma (array): adjacency matrix of the Gaussian state
42+
43+
Returns:
44+
(np.complex128): the montrealer of ``Sigma``
45+
"""
46+
n = len(Sigma) // 2
47+
tot_num = 2**n
48+
val = np.complex128(0)
49+
for p in numba.prange(tot_num):
50+
pos = dec2bin(p, n)
51+
lenpos = len(pos)
52+
pos = np.concatenate((pos, n + pos))
53+
submat = nb_ix(Sigma, pos, pos)
54+
sign = (-1) ** (lenpos + 1)
55+
val += (sign) * powertrace(submat, n + 1)[-1]
56+
return (-1) ** (n + 1) * val / (2 * n)
57+
58+
59+
@numba.jit(nopython=True)
60+
def power_loop(Sigma, zeta, n): # pragma: no cover
61+
"""Auxiliary function to calculate the product ``np.conj(zeta) @ Sigma^{n-1} @ zeta``.
62+
63+
Args:
64+
Sigma (array): square complex matrix
65+
zeta (array): complex vector
66+
n (int): sought after power
67+
68+
Returns:
69+
(np.complex128 or np.float64): the product np.conj(zeta) @ Sigma^{n-1} @ zeta
70+
"""
71+
vec = zeta
72+
for _ in range(n - 1):
73+
vec = Sigma @ vec
74+
return np.conj(zeta) @ vec
75+
76+
77+
@numba.jit(nopython=True, cache=True)
78+
def lmontrealer(Sigma, zeta): # pragma: no cover
79+
"""Calculates the loop-montrealer of the displaced Gaussian state with the given complex covariance matrix and vector of displacements.
80+
81+
Args:
82+
Sigma (array): complex Glauber covariance matrix of the Gaussian state
83+
zeta (array): vector of displacements
84+
85+
Returns:
86+
(np.complex128): the montrealer of ``Sigma``
87+
"""
88+
n = len(Sigma) // 2
89+
tot_num = 2**n
90+
val = np.complex128(0)
91+
val_loops = np.complex128(0)
92+
for p in numba.prange(tot_num):
93+
pos = dec2bin(p, n)
94+
lenpos = len(pos)
95+
pos = np.concatenate((pos, n + pos))
96+
subvec = zeta[pos]
97+
submat = nb_ix(Sigma, pos, pos)
98+
sign = (-1) ** (lenpos + 1)
99+
val_loops += sign * power_loop(submat, subvec, n)
100+
val += sign * powertrace(submat, n + 1)[-1]
101+
return (-1) ** (n + 1) * (val / (2 * n) + val_loops / 2)
102+
103+
104+
def lmtl(A, zeta):
105+
"""Returns the montrealer of an NxN matrix and an N-length vector.
106+
107+
Args:
108+
A (array): an NxN array of even dimensions
109+
zeta (array): an N-length vector of even dimensions
110+
111+
Returns:
112+
np.float64 or np.complex128: the loop montrealer of matrix A, vector zeta
113+
"""
114+
115+
tor_input_checks(A, zeta)
116+
n = len(A) // 2
117+
Sigma = Xmat(n) @ A
118+
return lmontrealer(Sigma, zeta)
119+
120+
121+
def mtl(A):
122+
"""Returns the montrealer of an NxN matrix.
123+
124+
Args:
125+
A (array): an NxN array of even dimensions.
126+
127+
Returns:
128+
np.float64 or np.complex128: the montrealer of matrix ``A``
129+
"""
130+
131+
tor_input_checks(A)
132+
n = len(A) // 2
133+
Sigma = Xmat(n) @ A
134+
return montrealer(Sigma)

‎thewalrus/_torontonian.py‎

Lines changed: 24 additions & 25 deletions
Original file line numberDiff line numberDiff line change
@@ -20,26 +20,41 @@
2020
from ._hafnian import reduction, find_kept_edges, nb_ix
2121

2222

23-
def tor(A, recursive=True):
24-
"""Returns the Torontonian of a matrix.
23+
def tor_input_checks(A, loops=None):
24+
"""Checks the correctness of the inputs for the torontonian/montrealer.
2525
2626
Args:
27-
A (array): a square array of even dimensions.
28-
recursive: use the faster recursive implementation.
29-
30-
Returns:
31-
np.float64 or np.complex128: the torontonian of matrix A.
27+
A (array): an NxN array of even dimensions
28+
loops (array): optional argument, an N-length vector of even dimensions
3229
"""
3330
if not isinstance(A, np.ndarray):
3431
raise TypeError("Input matrix must be a NumPy array.")
35-
3632
matshape = A.shape
3733

3834
if matshape[0] != matshape[1]:
3935
raise ValueError("Input matrix must be square.")
4036

4137
if matshape[0] % 2 != 0:
4238
raise ValueError("matrix dimension must be even")
39+
40+
if loops is not None:
41+
if not isinstance(loops, np.ndarray):
42+
raise TypeError("Input matrix must be a NumPy array.")
43+
if matshape[0] != len(loops):
44+
raise ValueError("gamma must be a vector matching the dimension of A")
45+
46+
47+
def tor(A, recursive=True):
48+
"""Returns the Torontonian of a matrix.
49+
50+
Args:
51+
A (array): a square array of even dimensions.
52+
recursive: use the faster recursive implementation.
53+
54+
Returns:
55+
np.float64 or np.complex128: the torontonian of matrix A.
56+
"""
57+
tor_input_checks(A)
4358
return rec_torontonian(A) if recursive else numba_tor(A)
4459

4560

@@ -54,23 +69,7 @@ def ltor(A, gamma, recursive=True):
5469
Returns:
5570
np.float64 or np.complex128: the loop torontonian of matrix A, vector gamma
5671
"""
57-
58-
if not isinstance(A, np.ndarray):
59-
raise TypeError("Input matrix must be a NumPy array.")
60-
61-
if not isinstance(gamma, np.ndarray):
62-
raise TypeError("Input matrix must be a NumPy array.")
63-
64-
matshape = A.shape
65-
66-
if matshape[0] != matshape[1]:
67-
raise ValueError("Input matrix must be square.")
68-
69-
if matshape[0] != len(gamma):
70-
raise ValueError("gamma must be a vector matching the dimension of A")
71-
72-
if matshape[0] % 2 != 0:
73-
raise ValueError("matrix dimension must be even")
72+
tor_input_checks(A, gamma)
7473

7574
return rec_ltorontonian(A, gamma) if recursive else numba_ltor(A, gamma)
7675

0 commit comments

Comments
 (0)