From ae0d3a496d87669410eb10af883d03948113a455 Mon Sep 17 00:00:00 2001 From: ziofil Date: Sun, 30 May 2021 08:22:53 -0400 Subject: [PATCH 01/29] adds sparse hafnian --- thewalrus/_hafnian.py | 33 +++++++++++++++++++++++++++++++++ 1 file changed, 33 insertions(+) diff --git a/thewalrus/_hafnian.py b/thewalrus/_hafnian.py index 0f95d4ab5..7d1e3296e 100644 --- a/thewalrus/_hafnian.py +++ b/thewalrus/_hafnian.py @@ -178,6 +178,39 @@ def hafnian( ) +def hafnian_sparse(A, D: set = None): + r"""Returns the hafnian of a sparse symmetric matrix. + This pure python implementation is very slow on full matrices, but faster the sparser a matrix is. + As a rule of thumb, the crossover in runtime with respect to hafnian() happens around 50% sparsity. + + Args: + A (array): the symmetric matrix of which we want to compute the hafnian + D (set): set of indices that identify a submatrix. If None (default) it computes + the hafnian of the whole matrix. + + Returns: + (float) hafnian of A or of the submatrix of A defined by the set of indices D. + """ + if D is None: + D = frozenset([i for i in range(len(A)) if not np.isclose(A[i, i], 0)]) + else: + D = frozenset(D) + + @lru_cache(maxsize=2 ** len(D)) + def indices(d, k): + return d.intersection(set(np.nonzero(A[k])[0])) + + @lru_cache(maxsize=2 ** len(D)) + def lhaf(d: frozenset) -> float: + if not d: + return 1 + d_ = set(d) + k = d_.pop() + return sum(A[i, k] * lhaf(frozenset(d_).difference({i})) for i in indices(d, k)) + + return lhaf(D) + + def hafnian_repeated(A, rpt, mu=None, loop=False, rtol=1e-05, atol=1e-08): r"""Returns the hafnian of matrix with repeated rows/columns. From 276e1f0fb7d98761b59feaa4075bf6c8becf4523 Mon Sep 17 00:00:00 2001 From: ziofil Date: Sun, 30 May 2021 08:57:23 -0400 Subject: [PATCH 02/29] adds sparse hafnian --- thewalrus/__init__.py | 2 ++ thewalrus/_hafnian.py | 16 +++++++++++----- thewalrus/tests/conftest.py | 4 +++- thewalrus/tests/test_hafnian.py | 14 +++++++++++++- 4 files changed, 29 insertions(+), 7 deletions(-) diff --git a/thewalrus/__init__.py b/thewalrus/__init__.py index b445f56f8..047be4ff4 100644 --- a/thewalrus/__init__.py +++ b/thewalrus/__init__.py @@ -100,6 +100,7 @@ hafnian, hafnian_repeated, reduction, + hafnian_sparse, ) from ._low_rank_haf import low_rank_hafnian from ._hermite_multidimensional import hafnian_batched, hermite_multidimensional @@ -112,6 +113,7 @@ "hafnian", "hafnian_repeated", "hafnian_batched", + "hafnian_sparse", "tor", "perm", "permanent_repeated", diff --git a/thewalrus/_hafnian.py b/thewalrus/_hafnian.py index 7d1e3296e..a14436695 100644 --- a/thewalrus/_hafnian.py +++ b/thewalrus/_hafnian.py @@ -15,9 +15,9 @@ Hafnian Python interface """ import numpy as np - +from functools import lru_cache from .libwalrus import haf_complex, haf_int, haf_real, haf_rpt_complex, haf_rpt_real - +from collections import Counter def input_validation(A, rtol=1e-05, atol=1e-08): """Checks that the matrix A satisfies the requirements for Hafnian calculation. @@ -192,15 +192,21 @@ def hafnian_sparse(A, D: set = None): (float) hafnian of A or of the submatrix of A defined by the set of indices D. """ if D is None: - D = frozenset([i for i in range(len(A)) if not np.isclose(A[i, i], 0)]) + D = frozenset([i for i in range(len(A))]) else: D = frozenset(D) - @lru_cache(maxsize=2 ** len(D)) + if np.allclose(A, np.zeros_like(A)): + return 0.0 + + r, _ = np.nonzero(A) + m = max(Counter(r).values()) # max nonzero values per row/column + + @lru_cache(maxsize=2 ** m) def indices(d, k): return d.intersection(set(np.nonzero(A[k])[0])) - @lru_cache(maxsize=2 ** len(D)) + @lru_cache(maxsize=2 ** m) def lhaf(d: frozenset) -> float: if not d: return 1 diff --git a/thewalrus/tests/conftest.py b/thewalrus/tests/conftest.py index 4af8feffa..fd0c3928d 100644 --- a/thewalrus/tests/conftest.py +++ b/thewalrus/tests/conftest.py @@ -41,10 +41,12 @@ def random_matrix(dtype): """Returns a random symmetric matrix of type dtype and of size n x n""" - def _wrapper(n): + def _wrapper(n, fill_factor=1.0): """wrapper function""" A = np.complex128(np.random.random([n, n])) A += 1j * np.random.random([n, n]) + mask = np.random.binomial(1, p=fill_factor, size=(n, n)) + A *= mask A += A.T if not np.issubdtype(dtype, np.complexfloating): diff --git a/thewalrus/tests/test_hafnian.py b/thewalrus/tests/test_hafnian.py index d055f369b..7c9c6aa95 100644 --- a/thewalrus/tests/test_hafnian.py +++ b/thewalrus/tests/test_hafnian.py @@ -19,7 +19,7 @@ from scipy.special import factorial as fac import thewalrus as hf -from thewalrus import hafnian, reduction +from thewalrus import hafnian, reduction, hafnian_sparse from thewalrus.libwalrus import haf_complex, haf_real, haf_int @@ -277,3 +277,15 @@ def test_diag(self, n): haf = hafnian(A, loop=True) expected = np.prod(v) assert np.allclose(haf, expected) + + +@pytest.mark.parametrize("n", [2, 3, 5, 10, 15, 20]) +@pytest.mark.parametrize("fill", [0.5, 0.2, 0.1, 0.05]) +class TestHafnianSparse: + """Tests that the loop hafnian sparse code gives + the same results as the full loop hafnian + """ + + def test_valid_output(self, random_matrix, n, fill): + A = random_matrix(n, fill_factor=fill) + assert np.allclose(hafnian(A, loop=True), hafnian_sparse(A)) From 1bc870a5e401febae899320512ebe2e4c3657807 Mon Sep 17 00:00:00 2001 From: ziofil Date: Sun, 30 May 2021 09:02:10 -0400 Subject: [PATCH 03/29] fixed codecov issues --- thewalrus/_hafnian.py | 8 +++++--- thewalrus/tests/test_hafnian.py | 1 + 2 files changed, 6 insertions(+), 3 deletions(-) diff --git a/thewalrus/_hafnian.py b/thewalrus/_hafnian.py index a14436695..ed1da4b21 100644 --- a/thewalrus/_hafnian.py +++ b/thewalrus/_hafnian.py @@ -14,10 +14,12 @@ """ Hafnian Python interface """ -import numpy as np from functools import lru_cache -from .libwalrus import haf_complex, haf_int, haf_real, haf_rpt_complex, haf_rpt_real from collections import Counter +import numpy as np + +from .libwalrus import haf_complex, haf_int, haf_real, haf_rpt_complex, haf_rpt_real + def input_validation(A, rtol=1e-05, atol=1e-08): """Checks that the matrix A satisfies the requirements for Hafnian calculation. @@ -192,7 +194,7 @@ def hafnian_sparse(A, D: set = None): (float) hafnian of A or of the submatrix of A defined by the set of indices D. """ if D is None: - D = frozenset([i for i in range(len(A))]) + D = frozenset(range(len(A))) else: D = frozenset(D) diff --git a/thewalrus/tests/test_hafnian.py b/thewalrus/tests/test_hafnian.py index 7c9c6aa95..12306e9df 100644 --- a/thewalrus/tests/test_hafnian.py +++ b/thewalrus/tests/test_hafnian.py @@ -287,5 +287,6 @@ class TestHafnianSparse: """ def test_valid_output(self, random_matrix, n, fill): + "tests that sparse loop hafnian matches full implementation" A = random_matrix(n, fill_factor=fill) assert np.allclose(hafnian(A, loop=True), hafnian_sparse(A)) From 9935fcba81fa1ddf2d6b993b9831ebba8072dc0d Mon Sep 17 00:00:00 2001 From: ziofil Date: Sun, 30 May 2021 09:03:11 -0400 Subject: [PATCH 04/29] fixed CodeFactor issue --- thewalrus/tests/test_hafnian.py | 13 ++++--------- 1 file changed, 4 insertions(+), 9 deletions(-) diff --git a/thewalrus/tests/test_hafnian.py b/thewalrus/tests/test_hafnian.py index 12306e9df..b1fb9d30a 100644 --- a/thewalrus/tests/test_hafnian.py +++ b/thewalrus/tests/test_hafnian.py @@ -281,12 +281,7 @@ def test_diag(self, n): @pytest.mark.parametrize("n", [2, 3, 5, 10, 15, 20]) @pytest.mark.parametrize("fill", [0.5, 0.2, 0.1, 0.05]) -class TestHafnianSparse: - """Tests that the loop hafnian sparse code gives - the same results as the full loop hafnian - """ - - def test_valid_output(self, random_matrix, n, fill): - "tests that sparse loop hafnian matches full implementation" - A = random_matrix(n, fill_factor=fill) - assert np.allclose(hafnian(A, loop=True), hafnian_sparse(A)) +def test_valid_output(random_matrix, n, fill): + "tests that sparse loop hafnian matches full implementation" + A = random_matrix(n, fill_factor=fill) + assert np.allclose(hafnian(A, loop=True), hafnian_sparse(A)) From 3438904cf9f6f0f9e18535b9529a8bf914d3839e Mon Sep 17 00:00:00 2001 From: ziofil Date: Sun, 30 May 2021 11:19:10 -0400 Subject: [PATCH 05/29] complete coverage --- thewalrus/tests/test_hafnian.py | 1 + 1 file changed, 1 insertion(+) diff --git a/thewalrus/tests/test_hafnian.py b/thewalrus/tests/test_hafnian.py index b1fb9d30a..fb9b24e1f 100644 --- a/thewalrus/tests/test_hafnian.py +++ b/thewalrus/tests/test_hafnian.py @@ -284,4 +284,5 @@ def test_diag(self, n): def test_valid_output(random_matrix, n, fill): "tests that sparse loop hafnian matches full implementation" A = random_matrix(n, fill_factor=fill) + assert np.allclose(hafnian_sparse(A), hafnian_sparse(A, D=set(range(len(A))))) assert np.allclose(hafnian(A, loop=True), hafnian_sparse(A)) From db91ea979dea83ec386c50bca31b163d3101e3ee Mon Sep 17 00:00:00 2001 From: ziofil Date: Sun, 30 May 2021 13:29:56 -0400 Subject: [PATCH 06/29] updated changelog --- .github/CHANGELOG.md | 2 ++ 1 file changed, 2 insertions(+) diff --git a/.github/CHANGELOG.md b/.github/CHANGELOG.md index 54aca53fb..c1f4f20f4 100644 --- a/.github/CHANGELOG.md +++ b/.github/CHANGELOG.md @@ -2,6 +2,7 @@ ### New features +* Adds the function `hafnian_sparse` to compute sparse loop hafnians (pure python implementation). ### Improvements * Speeds up the calculation of photon number variances/covariances [#224](https://github.com/XanaduAI/thewalrus/pull/224) @@ -14,6 +15,7 @@ This release contains contributions from (in alphabetical order): +Filippo Miatto Nicolas Quesada --- From 05d315b707a4319d2b300d1bee3b9c776908e1e3 Mon Sep 17 00:00:00 2001 From: ziofil Date: Sun, 30 May 2021 14:43:47 -0400 Subject: [PATCH 07/29] Update thewalrus/_hafnian.py Co-authored-by: Nicolas Quesada --- thewalrus/_hafnian.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/thewalrus/_hafnian.py b/thewalrus/_hafnian.py index ed1da4b21..112b82437 100644 --- a/thewalrus/_hafnian.py +++ b/thewalrus/_hafnian.py @@ -180,7 +180,7 @@ def hafnian( ) -def hafnian_sparse(A, D: set = None): +def hafnian_sparse(A, D = None): r"""Returns the hafnian of a sparse symmetric matrix. This pure python implementation is very slow on full matrices, but faster the sparser a matrix is. As a rule of thumb, the crossover in runtime with respect to hafnian() happens around 50% sparsity. From 263a521edc141e6f7b88a24873372ee79cc32698 Mon Sep 17 00:00:00 2001 From: ziofil Date: Sun, 30 May 2021 17:53:48 -0400 Subject: [PATCH 08/29] Update .github/CHANGELOG.md Co-authored-by: Nicolas Quesada --- .github/CHANGELOG.md | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/.github/CHANGELOG.md b/.github/CHANGELOG.md index c1f4f20f4..a291e5631 100644 --- a/.github/CHANGELOG.md +++ b/.github/CHANGELOG.md @@ -15,8 +15,7 @@ This release contains contributions from (in alphabetical order): -Filippo Miatto -Nicolas Quesada +Filippo Miatto, Nicolas Quesada --- From 33dbb8c59d30db4d2fd3a4ce2f973140e2f6174e Mon Sep 17 00:00:00 2001 From: ziofil Date: Tue, 1 Jun 2021 10:20:46 -0400 Subject: [PATCH 09/29] Update .github/CHANGELOG.md Co-authored-by: Nicolas Quesada --- .github/CHANGELOG.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/CHANGELOG.md b/.github/CHANGELOG.md index a291e5631..3e52a5777 100644 --- a/.github/CHANGELOG.md +++ b/.github/CHANGELOG.md @@ -5,7 +5,7 @@ * Adds the function `hafnian_sparse` to compute sparse loop hafnians (pure python implementation). ### Improvements -* Speeds up the calculation of photon number variances/covariances [#224](https://github.com/XanaduAI/thewalrus/pull/224) +* Speeds up the calculation of photon number variances/covariances [#245](https://github.com/XanaduAI/thewalrus/pull/245) ### Bug fixes From 21ddadea84b410112b2d3b07f529ade19e016961 Mon Sep 17 00:00:00 2001 From: ziofil Date: Wed, 9 Jun 2021 14:59:57 -0400 Subject: [PATCH 10/29] Update thewalrus/tests/test_hafnian.py Co-authored-by: Nicolas Quesada --- thewalrus/tests/test_hafnian.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/thewalrus/tests/test_hafnian.py b/thewalrus/tests/test_hafnian.py index fb9b24e1f..8763679cb 100644 --- a/thewalrus/tests/test_hafnian.py +++ b/thewalrus/tests/test_hafnian.py @@ -282,7 +282,7 @@ def test_diag(self, n): @pytest.mark.parametrize("n", [2, 3, 5, 10, 15, 20]) @pytest.mark.parametrize("fill", [0.5, 0.2, 0.1, 0.05]) def test_valid_output(random_matrix, n, fill): - "tests that sparse loop hafnian matches full implementation" + """tests that sparse loop hafnian matches full implementation""" A = random_matrix(n, fill_factor=fill) assert np.allclose(hafnian_sparse(A), hafnian_sparse(A, D=set(range(len(A))))) assert np.allclose(hafnian(A, loop=True), hafnian_sparse(A)) From 20cd25212529e34ab9250ce94a036b7738a2c243 Mon Sep 17 00:00:00 2001 From: ziofil Date: Thu, 10 Jun 2021 20:59:20 -0400 Subject: [PATCH 11/29] Update thewalrus/_hafnian.py Co-authored-by: antalszava --- thewalrus/_hafnian.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/thewalrus/_hafnian.py b/thewalrus/_hafnian.py index 112b82437..c947aac6b 100644 --- a/thewalrus/_hafnian.py +++ b/thewalrus/_hafnian.py @@ -188,7 +188,7 @@ def hafnian_sparse(A, D = None): Args: A (array): the symmetric matrix of which we want to compute the hafnian D (set): set of indices that identify a submatrix. If None (default) it computes - the hafnian of the whole matrix. + the hafnian of the whole matrix. Returns: (float) hafnian of A or of the submatrix of A defined by the set of indices D. From e0d6c0fec85a074626e7fa32e1ed9d4957be0e26 Mon Sep 17 00:00:00 2001 From: ziofil Date: Thu, 10 Jun 2021 20:59:35 -0400 Subject: [PATCH 12/29] Update thewalrus/_hafnian.py Co-authored-by: antalszava --- thewalrus/_hafnian.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/thewalrus/_hafnian.py b/thewalrus/_hafnian.py index c947aac6b..29a175671 100644 --- a/thewalrus/_hafnian.py +++ b/thewalrus/_hafnian.py @@ -187,7 +187,7 @@ def hafnian_sparse(A, D = None): Args: A (array): the symmetric matrix of which we want to compute the hafnian - D (set): set of indices that identify a submatrix. If None (default) it computes + D (set): set of indices that identify a submatrix. If ``None`` (default) it computes the hafnian of the whole matrix. Returns: From 24ba72b60c24a8d70afea794548e10bdd830ad9d Mon Sep 17 00:00:00 2001 From: ziofil Date: Thu, 10 Jun 2021 20:59:51 -0400 Subject: [PATCH 13/29] Update thewalrus/_hafnian.py Co-authored-by: antalszava --- thewalrus/_hafnian.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/thewalrus/_hafnian.py b/thewalrus/_hafnian.py index 29a175671..c4d81849c 100644 --- a/thewalrus/_hafnian.py +++ b/thewalrus/_hafnian.py @@ -183,7 +183,7 @@ def hafnian( def hafnian_sparse(A, D = None): r"""Returns the hafnian of a sparse symmetric matrix. This pure python implementation is very slow on full matrices, but faster the sparser a matrix is. - As a rule of thumb, the crossover in runtime with respect to hafnian() happens around 50% sparsity. + As a rule of thumb, the crossover in runtime with respect to :func:`~.hafnian` happens around 50% sparsity. Args: A (array): the symmetric matrix of which we want to compute the hafnian From 881250491a787598c280a80a41d2f36d5db40b42 Mon Sep 17 00:00:00 2001 From: ziofil Date: Thu, 10 Jun 2021 21:00:13 -0400 Subject: [PATCH 14/29] Update thewalrus/_hafnian.py Co-authored-by: antalszava --- thewalrus/_hafnian.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/thewalrus/_hafnian.py b/thewalrus/_hafnian.py index c4d81849c..745d8a076 100644 --- a/thewalrus/_hafnian.py +++ b/thewalrus/_hafnian.py @@ -191,7 +191,7 @@ def hafnian_sparse(A, D = None): the hafnian of the whole matrix. Returns: - (float) hafnian of A or of the submatrix of A defined by the set of indices D. + (float) hafnian of ``A`` or of the submatrix of ``A`` defined by the set of indices ``D``. """ if D is None: D = frozenset(range(len(A))) From 397350c2f8bed80e86cd59015561f3a0cfeaec9c Mon Sep 17 00:00:00 2001 From: ziofil Date: Thu, 10 Jun 2021 21:02:09 -0400 Subject: [PATCH 15/29] Update thewalrus/tests/test_hafnian.py Co-authored-by: antalszava --- thewalrus/tests/test_hafnian.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/thewalrus/tests/test_hafnian.py b/thewalrus/tests/test_hafnian.py index 8763679cb..efc30b5aa 100644 --- a/thewalrus/tests/test_hafnian.py +++ b/thewalrus/tests/test_hafnian.py @@ -282,7 +282,7 @@ def test_diag(self, n): @pytest.mark.parametrize("n", [2, 3, 5, 10, 15, 20]) @pytest.mark.parametrize("fill", [0.5, 0.2, 0.1, 0.05]) def test_valid_output(random_matrix, n, fill): - """tests that sparse loop hafnian matches full implementation""" + """Tests that sparse loop hafnian matches full implementation""" A = random_matrix(n, fill_factor=fill) assert np.allclose(hafnian_sparse(A), hafnian_sparse(A, D=set(range(len(A))))) assert np.allclose(hafnian(A, loop=True), hafnian_sparse(A)) From d660550ed58261137a95314272abd629e317de13 Mon Sep 17 00:00:00 2001 From: ziofil Date: Thu, 10 Jun 2021 21:13:55 -0400 Subject: [PATCH 16/29] Update .github/CHANGELOG.md Co-authored-by: antalszava --- .github/CHANGELOG.md | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/.github/CHANGELOG.md b/.github/CHANGELOG.md index 3e52a5777..37070a811 100644 --- a/.github/CHANGELOG.md +++ b/.github/CHANGELOG.md @@ -2,7 +2,8 @@ ### New features -* Adds the function `hafnian_sparse` to compute sparse loop hafnians (pure python implementation). +* Adds the function `hafnian_sparse` to compute sparse loop hafnians (pure python implementation). [#245](https://github.com/XanaduAI/thewalrus/pull/245) + ### Improvements * Speeds up the calculation of photon number variances/covariances [#245](https://github.com/XanaduAI/thewalrus/pull/245) From f53e75d0446df3326b57887a441e0ce6bfffdbec Mon Sep 17 00:00:00 2001 From: ziofil Date: Thu, 10 Jun 2021 21:14:03 -0400 Subject: [PATCH 17/29] Update .github/CHANGELOG.md Co-authored-by: antalszava --- .github/CHANGELOG.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/CHANGELOG.md b/.github/CHANGELOG.md index 37070a811..79e35e157 100644 --- a/.github/CHANGELOG.md +++ b/.github/CHANGELOG.md @@ -6,7 +6,7 @@ ### Improvements -* Speeds up the calculation of photon number variances/covariances [#245](https://github.com/XanaduAI/thewalrus/pull/245) +* Speeds up the calculation of photon number variances/covariances [#244](https://github.com/XanaduAI/thewalrus/pull/244) ### Bug fixes From 2922b3a3cf043f8f960159339904cdc9e2e276df Mon Sep 17 00:00:00 2001 From: ziofil Date: Thu, 10 Jun 2021 21:14:29 -0400 Subject: [PATCH 18/29] Update thewalrus/_hafnian.py Co-authored-by: antalszava --- thewalrus/_hafnian.py | 1 + 1 file changed, 1 insertion(+) diff --git a/thewalrus/_hafnian.py b/thewalrus/_hafnian.py index 745d8a076..687fa57f5 100644 --- a/thewalrus/_hafnian.py +++ b/thewalrus/_hafnian.py @@ -182,6 +182,7 @@ def hafnian( def hafnian_sparse(A, D = None): r"""Returns the hafnian of a sparse symmetric matrix. + This pure python implementation is very slow on full matrices, but faster the sparser a matrix is. As a rule of thumb, the crossover in runtime with respect to :func:`~.hafnian` happens around 50% sparsity. From a8bf92db3e15d83d23ec4dd9e7a058663ac6a6ff Mon Sep 17 00:00:00 2001 From: ziofil Date: Thu, 10 Jun 2021 21:19:19 -0400 Subject: [PATCH 19/29] Update test_hafnian.py --- thewalrus/tests/test_hafnian.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/thewalrus/tests/test_hafnian.py b/thewalrus/tests/test_hafnian.py index efc30b5aa..52b5505d6 100644 --- a/thewalrus/tests/test_hafnian.py +++ b/thewalrus/tests/test_hafnian.py @@ -279,10 +279,10 @@ def test_diag(self, n): assert np.allclose(haf, expected) -@pytest.mark.parametrize("n", [2, 3, 5, 10, 15, 20]) -@pytest.mark.parametrize("fill", [0.5, 0.2, 0.1, 0.05]) -def test_valid_output(random_matrix, n, fill): - """Tests that sparse loop hafnian matches full implementation""" - A = random_matrix(n, fill_factor=fill) - assert np.allclose(hafnian_sparse(A), hafnian_sparse(A, D=set(range(len(A))))) - assert np.allclose(hafnian(A, loop=True), hafnian_sparse(A)) + @pytest.mark.parametrize("n", [2, 3, 5, 10, 15, 20]) + @pytest.mark.parametrize("fill", [0.5, 0.2, 0.1, 0.05]) + def test_valid_output(self, random_matrix, n, fill): + """Tests that sparse loop hafnian matches full implementation""" + A = random_matrix(n, fill_factor=fill) + assert np.allclose(hafnian_sparse(A), hafnian_sparse(A, D=set(range(len(A))))) + assert np.allclose(hafnian(A, loop=True), hafnian_sparse(A)) From 3ec17da4ab940232358074d97b063433f3809c42 Mon Sep 17 00:00:00 2001 From: ziofil Date: Thu, 10 Jun 2021 21:41:03 -0400 Subject: [PATCH 20/29] adds loop bool --- thewalrus/_hafnian.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/thewalrus/_hafnian.py b/thewalrus/_hafnian.py index 687fa57f5..db42329f0 100644 --- a/thewalrus/_hafnian.py +++ b/thewalrus/_hafnian.py @@ -180,7 +180,7 @@ def hafnian( ) -def hafnian_sparse(A, D = None): +def hafnian_sparse(A, D=None, loop=False): r"""Returns the hafnian of a sparse symmetric matrix. This pure python implementation is very slow on full matrices, but faster the sparser a matrix is. @@ -190,6 +190,7 @@ def hafnian_sparse(A, D = None): A (array): the symmetric matrix of which we want to compute the hafnian D (set): set of indices that identify a submatrix. If ``None`` (default) it computes the hafnian of the whole matrix. + loop (bool): If ``True``, the loop hafnian is returned. Default is ``False``. Returns: (float) hafnian of ``A`` or of the submatrix of ``A`` defined by the set of indices ``D``. From 2e84a88fd68de74d6d31b711823b0cf31629f311 Mon Sep 17 00:00:00 2001 From: ziofil Date: Thu, 10 Jun 2021 22:02:29 -0400 Subject: [PATCH 21/29] adds tests --- thewalrus/tests/test_hafnian.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/thewalrus/tests/test_hafnian.py b/thewalrus/tests/test_hafnian.py index 52b5505d6..bae85a674 100644 --- a/thewalrus/tests/test_hafnian.py +++ b/thewalrus/tests/test_hafnian.py @@ -278,11 +278,12 @@ def test_diag(self, n): expected = np.prod(v) assert np.allclose(haf, expected) - @pytest.mark.parametrize("n", [2, 3, 5, 10, 15, 20]) @pytest.mark.parametrize("fill", [0.5, 0.2, 0.1, 0.05]) def test_valid_output(self, random_matrix, n, fill): """Tests that sparse loop hafnian matches full implementation""" A = random_matrix(n, fill_factor=fill) - assert np.allclose(hafnian_sparse(A), hafnian_sparse(A, D=set(range(len(A))))) - assert np.allclose(hafnian(A, loop=True), hafnian_sparse(A)) + assert np.allclose(hafnian_sparse(A, loop=True), hafnian_sparse(A, D=set(range(len(A))), loop=True)) + assert np.allclose(hafnian_sparse(A, loop=False), hafnian_sparse(A, D=set(range(len(A))), loop=False)) + assert np.allclose(hafnian(A, loop=True), hafnian_sparse(A, loop=True)) + assert np.allclose(hafnian(A, loop=False), hafnian_sparse(A, loop=False)) From 8ef103a6a598729034fa39eb8f84a62634b5cd1b Mon Sep 17 00:00:00 2001 From: ziofil Date: Thu, 10 Jun 2021 22:02:48 -0400 Subject: [PATCH 22/29] adds loop logic --- .DS_Store | Bin 0 -> 8196 bytes .vscode/settings.json | 3 + .../.ipynb_checkpoints/_hafnian-checkpoint.py | 296 ++++++++++++++++++ .../_hermite_multidimensional-checkpoint.py | 146 +++++++++ thewalrus/_hafnian.py | 3 + 5 files changed, 448 insertions(+) create mode 100644 .DS_Store create mode 100644 .vscode/settings.json create mode 100644 thewalrus/.ipynb_checkpoints/_hafnian-checkpoint.py create mode 100644 thewalrus/.ipynb_checkpoints/_hermite_multidimensional-checkpoint.py diff --git a/.DS_Store b/.DS_Store new file mode 100644 index 0000000000000000000000000000000000000000..a8be244a784cb8c110c82e0064cf849911c7ef08 GIT binary patch literal 8196 zcmeHLy>1gh5S|4GE+7N~AyQDBsAv);G>D2Lm=X~wH4?el*hW6i!gioM@dQ+%%^Ohg z01@I1=t2P$@CZnB@XhYp?Cjo2ra(euSGpZ*=bM}FzWq7R5s})Qb~cDsh$z9u_QF|o zHI4mz1?_|#S%+1KC)%YUbHyQRE-{-A!u1|8Q&hKB*T;{Pw z=`slYBJa11^PZmNr9Mq57uo5_!w}(dSbvvj4&Hq{CJiXThYn_Z_GpY*Dyc(V$%pP+ zJ`Td;b7<#>`3RIyrU51S=u>uya_BO^H=mB+^5NEXaB;kn2i>r`2!zX_4k7X%Rjfw~ zDwI%TEVAGJ{=oHlx%aIG^H@0>6!=BnyE<=zZ(J{EhYA|$?_DxKU7aF-fRQ|^*{s*# zTqT-uf){NeCD!YJp6We34c3R+DMCa8C+|K|p*bTDCHSZ>%;%Jzoe7lM^H*zfM$00v7`&S^u?9?Dce)DH#8JYWF@El4Eamr3T z9%;q@m9e1B@4n zbD!bk&#wn||1B8%`>FmkRy-#*YHWc5C(XzT*Z=pLzyAkuaUd!X75IM?P_6cMdlS!N zYwO_2xYkZ_J;KF}ah0(;1RXta3#OP31Q_27J4*^bH6W{;w K{m*OhIJRF5_w7^w literal 0 HcmV?d00001 diff --git a/.vscode/settings.json b/.vscode/settings.json new file mode 100644 index 000000000..f5e8a594e --- /dev/null +++ b/.vscode/settings.json @@ -0,0 +1,3 @@ +{ + "python.pythonPath": "/Users/filippo/miniforge3/envs/tf24/bin/python" +} \ No newline at end of file diff --git a/thewalrus/.ipynb_checkpoints/_hafnian-checkpoint.py b/thewalrus/.ipynb_checkpoints/_hafnian-checkpoint.py new file mode 100644 index 000000000..db42329f0 --- /dev/null +++ b/thewalrus/.ipynb_checkpoints/_hafnian-checkpoint.py @@ -0,0 +1,296 @@ +# 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. +""" +Hafnian Python interface +""" +from functools import lru_cache +from collections import Counter +import numpy as np + +from .libwalrus import haf_complex, haf_int, haf_real, haf_rpt_complex, haf_rpt_real + + +def input_validation(A, rtol=1e-05, atol=1e-08): + """Checks that the matrix A satisfies the requirements for Hafnian calculation. + + These include: + + * That the ``A`` is a NumPy array + * That ``A`` is square + * That ``A`` does not contain any NaNs + * That ``A`` is symmetric + + Args: + A (array): a NumPy array. + rtol (float): the relative tolerance parameter used in ``np.allclose``. + atol (float): the absolute tolerance parameter used in ``np.allclose``. + + Returns: + bool: returns True if the matrix satisfies all requirements. + """ + + if not isinstance(A, np.ndarray): + raise TypeError("Input matrix must be a NumPy array.") + + n = A.shape + + if n[0] != n[1]: + raise ValueError("Input matrix must be square.") + + if np.isnan(A).any(): + raise ValueError("Input matrix must not contain NaNs.") + + if not np.allclose(A, A.T, rtol=rtol, atol=atol): + raise ValueError("Input matrix must be symmetric.") + + return True + + +def reduction(A, rpt): + r"""Calculates the reduction of an array by a vector of indices. + + This is equivalent to repeating the ith row/column of :math:`A`, :math:`rpt_i` times. + + Args: + A (array): matrix of size [N, N] + rpt (Sequence): sequence of N positive integers indicating the corresponding rows/columns + of A to be repeated. + Returns: + array: the reduction of A by the index vector rpt + """ + rows = [i for sublist in [[idx] * j for idx, j in enumerate(rpt)] for i in sublist] + + if A.ndim == 1: + return A[rows] + + return A[:, rows][rows] + +# pylint: disable=too-many-arguments +def hafnian( + A, loop=False, recursive=True, rtol=1e-05, atol=1e-08, quad=True, approx=False, num_samples=1000 +): # pylint: disable=too-many-arguments + """Returns the hafnian of a matrix. + + For more direct control, you may wish to call :func:`haf_real`, + :func:`haf_complex`, or :func:`haf_int` directly. + + Args: + A (array): a square, symmetric array of even dimensions. + loop (bool): If ``True``, the loop hafnian is returned. Default is ``False``. + recursive (bool): If ``True``, the recursive algorithm is used. Note: + the recursive algorithm does not currently support the loop hafnian. + If ``loop=True``, then this keyword argument is ignored. + rtol (float): the relative tolerance parameter used in ``np.allclose``. + atol (float): the absolute tolerance parameter used in ``np.allclose``. + quad (bool): If ``True``, the hafnian algorithm is performed with quadruple precision. + approx (bool): If ``True``, an approximation algorithm is used to estimate the hafnian. Note that + the approximation algorithm can only be applied to matrices ``A`` that only have non-negative entries. + num_samples (int): If ``approx=True``, the approximation algorithm performs ``num_samples`` iterations + for estimation of the hafnian of the non-negative matrix ``A``. + + Returns: + np.int64 or np.float64 or np.complex128: the hafnian of matrix A. + """ + # pylint: disable=too-many-return-statements,too-many-branches + input_validation(A, rtol=rtol, atol=atol) + + matshape = A.shape + + if matshape == (0, 0): + return 1 + + if matshape[0] % 2 != 0 and not loop: + return 0.0 + + if np.allclose(np.diag(np.diag(A)), A, rtol=rtol, atol=atol): + if loop: + return np.prod(np.diag(A)) + return 0 + + if matshape[0] % 2 != 0 and loop: + A = np.pad(A, pad_width=((0, 1), (0, 1)), mode="constant") + A[-1, -1] = 1.0 + + matshape = A.shape + + if matshape[0] == 2: + if loop: + return A[0, 1] + A[0, 0] * A[1, 1] + return A[0][1] + + if matshape[0] == 4: + if loop: + result = ( + A[0, 1] * A[2, 3] + + A[0, 2] * A[1, 3] + + A[0, 3] * A[1, 2] + + A[0, 0] * A[1, 1] * A[2, 3] + + A[0, 1] * A[2, 2] * A[3, 3] + + A[0, 2] * A[1, 1] * A[3, 3] + + A[0, 0] * A[2, 2] * A[1, 3] + + A[0, 0] * A[3, 3] * A[1, 2] + + A[0, 3] * A[1, 1] * A[2, 2] + + A[0, 0] * A[1, 1] * A[2, 2] * A[3, 3] + ) + return result + + return A[0, 1] * A[2, 3] + A[0, 2] * A[1, 3] + A[0, 3] * A[1, 2] + + if approx: + if np.any(np.iscomplex(A)): + raise ValueError("Input matrix must be real") + + if np.any(A < 0): + raise ValueError("Input matrix must not have negative entries") + + if A.dtype == np.complex: + # array data is complex type + if np.any(np.iscomplex(A)): + # array values contain non-zero imaginary parts + return haf_complex(A, loop=loop, recursive=recursive, quad=quad) + + # all array values have zero imaginary parts + return haf_real(np.float64(A.real), loop=loop, recursive=recursive, quad=quad) + + if np.issubdtype(A.dtype, np.integer) and not loop: + # array data is an integer type, and the user is not + # requesting the loop hafnian + return haf_int(np.int64(A)) + + if np.issubdtype(A.dtype, np.integer) and loop: + # array data is an integer type, and the user is + # requesting the loop hafnian. Currently no + # integer function for loop hafnians, have to instead + # convert to float and use haf_real + A = np.float64(A) + + return haf_real( + A, loop=loop, recursive=recursive, quad=quad, approx=approx, nsamples=num_samples + ) + + +def hafnian_sparse(A, D=None, loop=False): + r"""Returns the hafnian of a sparse symmetric matrix. + + This pure python implementation is very slow on full matrices, but faster the sparser a matrix is. + As a rule of thumb, the crossover in runtime with respect to :func:`~.hafnian` happens around 50% sparsity. + + Args: + A (array): the symmetric matrix of which we want to compute the hafnian + D (set): set of indices that identify a submatrix. If ``None`` (default) it computes + the hafnian of the whole matrix. + loop (bool): If ``True``, the loop hafnian is returned. Default is ``False``. + + Returns: + (float) hafnian of ``A`` or of the submatrix of ``A`` defined by the set of indices ``D``. + """ + if D is None: + D = frozenset(range(len(A))) + else: + D = frozenset(D) + + if np.allclose(A, np.zeros_like(A)): + return 0.0 + + r, _ = np.nonzero(A) + m = max(Counter(r).values()) # max nonzero values per row/column + + @lru_cache(maxsize=2 ** m) + def indices(d, k): + return d.intersection(set(np.nonzero(A[k])[0])) + + @lru_cache(maxsize=2 ** m) + def lhaf(d: frozenset) -> float: + if not d: + return 1 + d_ = set(d) + k = d_.pop() + return sum(A[i, k] * lhaf(frozenset(d_).difference({i})) for i in indices(d, k)) + + return lhaf(D) + + +def hafnian_repeated(A, rpt, mu=None, loop=False, rtol=1e-05, atol=1e-08): + r"""Returns the hafnian of matrix with repeated rows/columns. + + The :func:`reduction` function may be used to show the resulting matrix + with repeated rows and columns as per ``rpt``. + + As a result, the following are identical: + + >>> hafnian_repeated(A, rpt) + >>> hafnian(reduction(A, rpt)) + + However, using ``hafnian_repeated`` in the case where there are a large number + of repeated rows and columns (:math:`\sum_{i}rpt_i \gg N`) can be + significantly faster. + + .. note:: + + If :math:`rpt=(1, 1, \dots, 1)`, then + + >>> hafnian_repeated(A, rpt) == hafnian(A) + + For more direct control, you may wish to call :func:`haf_rpt_real` or + :func:`haf_rpt_complex` directly. + + Args: + A (array): a square, symmetric :math:`N\times N` array. + rpt (Sequence): a length-:math:`N` positive integer sequence, corresponding + to the number of times each row/column of matrix :math:`A` is repeated. + mu (array): a vector of length :math:`N` representing the vector of means/displacement. + If not provided, ``mu`` is set to the diagonal of matrix ``A``. Note that this + only affects the loop hafnian. + loop (bool): If ``True``, the loop hafnian is returned. Default is ``False``. + use_eigen (bool): if True (default), the Eigen linear algebra library + is used for matrix multiplication. If the hafnian library was compiled + with BLAS/Lapack support, then BLAS will be used for matrix multiplication. + rtol (float): the relative tolerance parameter used in ``np.allclose``. + atol (float): the absolute tolerance parameter used in ``np.allclose``. + + Returns: + np.int64 or np.float64 or np.complex128: the hafnian of matrix A. + """ + # pylint: disable=too-many-return-statements,too-many-branches + input_validation(A, atol=atol, rtol=rtol) + + if len(rpt) != len(A): + raise ValueError("the rpt argument must be 1-dimensional sequence of length len(A).") + + nud = np.array(rpt, dtype=np.int32) + + if not np.all(np.mod(nud, 1) == 0) or np.any(nud < 0): + raise ValueError("the rpt argument must contain non-negative integers.") + + if np.all(nud == 0): + return 1.0 + + if np.sum(nud) % 2 != 0 and not loop: + return 0.0 + + if mu is None: + mu = A.diagonal().copy() + + if np.allclose(A, 0, rtol=rtol, atol=atol): + if loop: + return np.prod(mu ** rpt) + return 0 + + if len(mu) != len(A): + raise ValueError("Length of means vector must be the same length as the matrix A.") + + if A.dtype == np.complex or mu.dtype == np.complex: + return haf_rpt_complex(A, nud, mu=mu, loop=loop) + + return haf_rpt_real(A, nud, mu=mu, loop=loop) diff --git a/thewalrus/.ipynb_checkpoints/_hermite_multidimensional-checkpoint.py b/thewalrus/.ipynb_checkpoints/_hermite_multidimensional-checkpoint.py new file mode 100644 index 000000000..6190e66a1 --- /dev/null +++ b/thewalrus/.ipynb_checkpoints/_hermite_multidimensional-checkpoint.py @@ -0,0 +1,146 @@ +# 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. +""" +Hermite Multidimensional Python interface +""" +import numpy as np + +from .libwalrus import hermite_multidimensional as hm +from .libwalrus import hermite_multidimensional_real as hmr + +from .libwalrus import renorm_hermite_multidimensional as rhm +from .libwalrus import renorm_hermite_multidimensional_real as rhmr + +from ._hafnian import input_validation + +# pylint: disable=too-many-arguments +def hermite_multidimensional(R, cutoff, y=None, renorm=False, make_tensor=True, modified=False, rtol=1e-05, atol=1e-08): + r"""Returns the multidimensional Hermite polynomials :math:`H_k^{(R)}(y)`. + + Here :math:`R` is an :math:`n \times n` square matrix, and + :math:`y` is an :math:`n` dimensional vector. The polynomials are + parametrized by the multi-index :math:`k=(k_0,k_1,\ldots,k_{n-1})`, + and are calculated for all values :math:`0 \leq k_j < \text{cutoff}`, + thus a tensor of dimensions :math:`\text{cutoff}^n` is returned. + + This tensor can either be flattened into a vector or returned as an actual + tensor with :math:`n` indices. + + .. note:: + + Note that if :math:`R = (1)` then :math:`H_k^{(R)}(y)` + are precisely the well known **probabilists' Hermite polynomials** :math:`He_k(y)`, + and if :math:`R = (2)` then :math:`H_k^{(R)}(y)` are precisely the well known + **physicists' Hermite polynomials** :math:`H_k(y)`. + + Args: + R (array): square matrix parametrizing the Hermite polynomial family + cutoff (int): maximum size of the subindices in the Hermite polynomial + y (array): vector argument of the Hermite polynomial + renorm (bool): If ``True``, normalizes the returned multidimensional Hermite + polynomials such that :math:`H_k^{(R)}(y)/\prod_i k_i!` + make_tensor (bool): If ``False``, returns a flattened one dimensional array + containing the values of the polynomial + modified (bool): whether to return the modified multidimensional Hermite polynomials or the standard ones + rtol (float): the relative tolerance parameter used in ``np.allclose`` + atol (float): the absolute tolerance parameter used in ``np.allclose`` + Returns: + (array): the multidimensional Hermite polynomials + """ + + input_validation(R, atol=atol, rtol=rtol) + n, _ = R.shape + + if (modified is False) and (y is not None): + m = y.shape[0] + if m == n: + ym = R @ y + return hermite_multidimensional( + R, cutoff, y=ym, renorm=renorm, make_tensor=make_tensor, modified=True + ) + + if y is None: + y = np.zeros([n], dtype=complex) + + m = y.shape[0] + if m != n: + raise ValueError("The matrix R and vector y have incompatible dimensions") + + + Rt = np.real_if_close(R) + yt = np.real_if_close(y) + + if Rt.dtype == np.float and yt.dtype == np.float: + if renorm: + values = np.array(rhmr(Rt, yt, cutoff)) + else: + values = np.array(hmr(Rt, yt, cutoff)) + else: + if renorm: + values = np.array(rhm(np.complex128(R), np.complex128(y), cutoff)) + else: + values = np.array(hm(np.complex128(R), np.complex128(y), cutoff)) + + if make_tensor: + shape = cutoff * np.ones([n], dtype=int) + values = np.reshape(values, shape) + + return values + + +# pylint: disable=too-many-arguments +def hafnian_batched(A, cutoff, mu=None, rtol=1e-05, atol=1e-08, renorm=False, make_tensor=True): + r"""Calculates the hafnian of :func:`reduction(A, k) ` + for all possible values of vector ``k`` below the specified cutoff. + + Here, + + * :math:`A` is am :math:`n\times n` square matrix + * :math:`k` is a vector of (non-negative) integers with the same dimensions as :math:`A`, + i.e., :math:`k = (k_0,k_1,\ldots,k_{n-1})`, and where :math:`0 \leq k_j < \texttt{cutoff}`. + + The function :func:`~.hafnian_repeated` can be used to calculate the reduced hafnian + for a *specific* value of :math:`k`; see the documentation for more information. + + .. note:: + + If ``mu`` is not ``None``, this function instead returns + ``hafnian(np.fill_diagonal(reduction(A, k), reduction(mu, k)), loop=True)``. + This calculation can only be performed if the matrix :math:`A` is invertible. + + Args: + A (array): a square, symmetric :math:`N\times N` array. + cutoff (int): maximum size of the subindices in the Hermite polynomial + mu (array): a vector of length :math:`N` representing the vector of means/displacement + renorm (bool): If ``True``, the returned hafnians are *normalized*, that is, + :math:`haf(reduction(A, k))/\ \sqrt{prod_i k_i!}` + (or :math:`lhaf(fill\_diagonal(reduction(A, k),reduction(mu, k)))` if + ``mu`` is not None) + make_tensor: If ``False``, returns a flattened one dimensional array instead + of a tensor with the values of the hafnians. + rtol (float): the relative tolerance parameter used in ``np.allclose``. + atol (float): the absolute tolerance parameter used in ``np.allclose``. + Returns: + (array): the values of the hafnians for each value of :math:`k` up to the cutoff + """ + input_validation(A, atol=atol, rtol=rtol) + n, _ = A.shape + + if mu is None: + mu = np.zeros([n], dtype=complex) + + return hermite_multidimensional( + -A, cutoff, y=mu, renorm=renorm, make_tensor=make_tensor, modified=True + ) +# Note the minus signs in the arguments. Those are intentional and are due to the fact that Dodonov et al. in PRA 50, 813 (1994) use (p,q) ordering instead of (q,p) ordering diff --git a/thewalrus/_hafnian.py b/thewalrus/_hafnian.py index db42329f0..a461b9915 100644 --- a/thewalrus/_hafnian.py +++ b/thewalrus/_hafnian.py @@ -200,6 +200,9 @@ def hafnian_sparse(A, D=None, loop=False): else: D = frozenset(D) + if not loop: + A = A - np.diag(np.diag(A)) + if np.allclose(A, np.zeros_like(A)): return 0.0 From 240893fdb6e7547f975df1262ba8132b3fd992b3 Mon Sep 17 00:00:00 2001 From: Nicolas Quesada Date: Tue, 22 Jun 2021 15:29:43 -0400 Subject: [PATCH 23/29] removes unused files --- .vscode/settings.json | 3 - .../.ipynb_checkpoints/_hafnian-checkpoint.py | 296 ------------------ .../_hermite_multidimensional-checkpoint.py | 146 --------- 3 files changed, 445 deletions(-) delete mode 100644 .vscode/settings.json delete mode 100644 thewalrus/.ipynb_checkpoints/_hafnian-checkpoint.py delete mode 100644 thewalrus/.ipynb_checkpoints/_hermite_multidimensional-checkpoint.py diff --git a/.vscode/settings.json b/.vscode/settings.json deleted file mode 100644 index f5e8a594e..000000000 --- a/.vscode/settings.json +++ /dev/null @@ -1,3 +0,0 @@ -{ - "python.pythonPath": "/Users/filippo/miniforge3/envs/tf24/bin/python" -} \ No newline at end of file diff --git a/thewalrus/.ipynb_checkpoints/_hafnian-checkpoint.py b/thewalrus/.ipynb_checkpoints/_hafnian-checkpoint.py deleted file mode 100644 index db42329f0..000000000 --- a/thewalrus/.ipynb_checkpoints/_hafnian-checkpoint.py +++ /dev/null @@ -1,296 +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. -""" -Hafnian Python interface -""" -from functools import lru_cache -from collections import Counter -import numpy as np - -from .libwalrus import haf_complex, haf_int, haf_real, haf_rpt_complex, haf_rpt_real - - -def input_validation(A, rtol=1e-05, atol=1e-08): - """Checks that the matrix A satisfies the requirements for Hafnian calculation. - - These include: - - * That the ``A`` is a NumPy array - * That ``A`` is square - * That ``A`` does not contain any NaNs - * That ``A`` is symmetric - - Args: - A (array): a NumPy array. - rtol (float): the relative tolerance parameter used in ``np.allclose``. - atol (float): the absolute tolerance parameter used in ``np.allclose``. - - Returns: - bool: returns True if the matrix satisfies all requirements. - """ - - if not isinstance(A, np.ndarray): - raise TypeError("Input matrix must be a NumPy array.") - - n = A.shape - - if n[0] != n[1]: - raise ValueError("Input matrix must be square.") - - if np.isnan(A).any(): - raise ValueError("Input matrix must not contain NaNs.") - - if not np.allclose(A, A.T, rtol=rtol, atol=atol): - raise ValueError("Input matrix must be symmetric.") - - return True - - -def reduction(A, rpt): - r"""Calculates the reduction of an array by a vector of indices. - - This is equivalent to repeating the ith row/column of :math:`A`, :math:`rpt_i` times. - - Args: - A (array): matrix of size [N, N] - rpt (Sequence): sequence of N positive integers indicating the corresponding rows/columns - of A to be repeated. - Returns: - array: the reduction of A by the index vector rpt - """ - rows = [i for sublist in [[idx] * j for idx, j in enumerate(rpt)] for i in sublist] - - if A.ndim == 1: - return A[rows] - - return A[:, rows][rows] - -# pylint: disable=too-many-arguments -def hafnian( - A, loop=False, recursive=True, rtol=1e-05, atol=1e-08, quad=True, approx=False, num_samples=1000 -): # pylint: disable=too-many-arguments - """Returns the hafnian of a matrix. - - For more direct control, you may wish to call :func:`haf_real`, - :func:`haf_complex`, or :func:`haf_int` directly. - - Args: - A (array): a square, symmetric array of even dimensions. - loop (bool): If ``True``, the loop hafnian is returned. Default is ``False``. - recursive (bool): If ``True``, the recursive algorithm is used. Note: - the recursive algorithm does not currently support the loop hafnian. - If ``loop=True``, then this keyword argument is ignored. - rtol (float): the relative tolerance parameter used in ``np.allclose``. - atol (float): the absolute tolerance parameter used in ``np.allclose``. - quad (bool): If ``True``, the hafnian algorithm is performed with quadruple precision. - approx (bool): If ``True``, an approximation algorithm is used to estimate the hafnian. Note that - the approximation algorithm can only be applied to matrices ``A`` that only have non-negative entries. - num_samples (int): If ``approx=True``, the approximation algorithm performs ``num_samples`` iterations - for estimation of the hafnian of the non-negative matrix ``A``. - - Returns: - np.int64 or np.float64 or np.complex128: the hafnian of matrix A. - """ - # pylint: disable=too-many-return-statements,too-many-branches - input_validation(A, rtol=rtol, atol=atol) - - matshape = A.shape - - if matshape == (0, 0): - return 1 - - if matshape[0] % 2 != 0 and not loop: - return 0.0 - - if np.allclose(np.diag(np.diag(A)), A, rtol=rtol, atol=atol): - if loop: - return np.prod(np.diag(A)) - return 0 - - if matshape[0] % 2 != 0 and loop: - A = np.pad(A, pad_width=((0, 1), (0, 1)), mode="constant") - A[-1, -1] = 1.0 - - matshape = A.shape - - if matshape[0] == 2: - if loop: - return A[0, 1] + A[0, 0] * A[1, 1] - return A[0][1] - - if matshape[0] == 4: - if loop: - result = ( - A[0, 1] * A[2, 3] - + A[0, 2] * A[1, 3] - + A[0, 3] * A[1, 2] - + A[0, 0] * A[1, 1] * A[2, 3] - + A[0, 1] * A[2, 2] * A[3, 3] - + A[0, 2] * A[1, 1] * A[3, 3] - + A[0, 0] * A[2, 2] * A[1, 3] - + A[0, 0] * A[3, 3] * A[1, 2] - + A[0, 3] * A[1, 1] * A[2, 2] - + A[0, 0] * A[1, 1] * A[2, 2] * A[3, 3] - ) - return result - - return A[0, 1] * A[2, 3] + A[0, 2] * A[1, 3] + A[0, 3] * A[1, 2] - - if approx: - if np.any(np.iscomplex(A)): - raise ValueError("Input matrix must be real") - - if np.any(A < 0): - raise ValueError("Input matrix must not have negative entries") - - if A.dtype == np.complex: - # array data is complex type - if np.any(np.iscomplex(A)): - # array values contain non-zero imaginary parts - return haf_complex(A, loop=loop, recursive=recursive, quad=quad) - - # all array values have zero imaginary parts - return haf_real(np.float64(A.real), loop=loop, recursive=recursive, quad=quad) - - if np.issubdtype(A.dtype, np.integer) and not loop: - # array data is an integer type, and the user is not - # requesting the loop hafnian - return haf_int(np.int64(A)) - - if np.issubdtype(A.dtype, np.integer) and loop: - # array data is an integer type, and the user is - # requesting the loop hafnian. Currently no - # integer function for loop hafnians, have to instead - # convert to float and use haf_real - A = np.float64(A) - - return haf_real( - A, loop=loop, recursive=recursive, quad=quad, approx=approx, nsamples=num_samples - ) - - -def hafnian_sparse(A, D=None, loop=False): - r"""Returns the hafnian of a sparse symmetric matrix. - - This pure python implementation is very slow on full matrices, but faster the sparser a matrix is. - As a rule of thumb, the crossover in runtime with respect to :func:`~.hafnian` happens around 50% sparsity. - - Args: - A (array): the symmetric matrix of which we want to compute the hafnian - D (set): set of indices that identify a submatrix. If ``None`` (default) it computes - the hafnian of the whole matrix. - loop (bool): If ``True``, the loop hafnian is returned. Default is ``False``. - - Returns: - (float) hafnian of ``A`` or of the submatrix of ``A`` defined by the set of indices ``D``. - """ - if D is None: - D = frozenset(range(len(A))) - else: - D = frozenset(D) - - if np.allclose(A, np.zeros_like(A)): - return 0.0 - - r, _ = np.nonzero(A) - m = max(Counter(r).values()) # max nonzero values per row/column - - @lru_cache(maxsize=2 ** m) - def indices(d, k): - return d.intersection(set(np.nonzero(A[k])[0])) - - @lru_cache(maxsize=2 ** m) - def lhaf(d: frozenset) -> float: - if not d: - return 1 - d_ = set(d) - k = d_.pop() - return sum(A[i, k] * lhaf(frozenset(d_).difference({i})) for i in indices(d, k)) - - return lhaf(D) - - -def hafnian_repeated(A, rpt, mu=None, loop=False, rtol=1e-05, atol=1e-08): - r"""Returns the hafnian of matrix with repeated rows/columns. - - The :func:`reduction` function may be used to show the resulting matrix - with repeated rows and columns as per ``rpt``. - - As a result, the following are identical: - - >>> hafnian_repeated(A, rpt) - >>> hafnian(reduction(A, rpt)) - - However, using ``hafnian_repeated`` in the case where there are a large number - of repeated rows and columns (:math:`\sum_{i}rpt_i \gg N`) can be - significantly faster. - - .. note:: - - If :math:`rpt=(1, 1, \dots, 1)`, then - - >>> hafnian_repeated(A, rpt) == hafnian(A) - - For more direct control, you may wish to call :func:`haf_rpt_real` or - :func:`haf_rpt_complex` directly. - - Args: - A (array): a square, symmetric :math:`N\times N` array. - rpt (Sequence): a length-:math:`N` positive integer sequence, corresponding - to the number of times each row/column of matrix :math:`A` is repeated. - mu (array): a vector of length :math:`N` representing the vector of means/displacement. - If not provided, ``mu`` is set to the diagonal of matrix ``A``. Note that this - only affects the loop hafnian. - loop (bool): If ``True``, the loop hafnian is returned. Default is ``False``. - use_eigen (bool): if True (default), the Eigen linear algebra library - is used for matrix multiplication. If the hafnian library was compiled - with BLAS/Lapack support, then BLAS will be used for matrix multiplication. - rtol (float): the relative tolerance parameter used in ``np.allclose``. - atol (float): the absolute tolerance parameter used in ``np.allclose``. - - Returns: - np.int64 or np.float64 or np.complex128: the hafnian of matrix A. - """ - # pylint: disable=too-many-return-statements,too-many-branches - input_validation(A, atol=atol, rtol=rtol) - - if len(rpt) != len(A): - raise ValueError("the rpt argument must be 1-dimensional sequence of length len(A).") - - nud = np.array(rpt, dtype=np.int32) - - if not np.all(np.mod(nud, 1) == 0) or np.any(nud < 0): - raise ValueError("the rpt argument must contain non-negative integers.") - - if np.all(nud == 0): - return 1.0 - - if np.sum(nud) % 2 != 0 and not loop: - return 0.0 - - if mu is None: - mu = A.diagonal().copy() - - if np.allclose(A, 0, rtol=rtol, atol=atol): - if loop: - return np.prod(mu ** rpt) - return 0 - - if len(mu) != len(A): - raise ValueError("Length of means vector must be the same length as the matrix A.") - - if A.dtype == np.complex or mu.dtype == np.complex: - return haf_rpt_complex(A, nud, mu=mu, loop=loop) - - return haf_rpt_real(A, nud, mu=mu, loop=loop) diff --git a/thewalrus/.ipynb_checkpoints/_hermite_multidimensional-checkpoint.py b/thewalrus/.ipynb_checkpoints/_hermite_multidimensional-checkpoint.py deleted file mode 100644 index 6190e66a1..000000000 --- a/thewalrus/.ipynb_checkpoints/_hermite_multidimensional-checkpoint.py +++ /dev/null @@ -1,146 +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. -""" -Hermite Multidimensional Python interface -""" -import numpy as np - -from .libwalrus import hermite_multidimensional as hm -from .libwalrus import hermite_multidimensional_real as hmr - -from .libwalrus import renorm_hermite_multidimensional as rhm -from .libwalrus import renorm_hermite_multidimensional_real as rhmr - -from ._hafnian import input_validation - -# pylint: disable=too-many-arguments -def hermite_multidimensional(R, cutoff, y=None, renorm=False, make_tensor=True, modified=False, rtol=1e-05, atol=1e-08): - r"""Returns the multidimensional Hermite polynomials :math:`H_k^{(R)}(y)`. - - Here :math:`R` is an :math:`n \times n` square matrix, and - :math:`y` is an :math:`n` dimensional vector. The polynomials are - parametrized by the multi-index :math:`k=(k_0,k_1,\ldots,k_{n-1})`, - and are calculated for all values :math:`0 \leq k_j < \text{cutoff}`, - thus a tensor of dimensions :math:`\text{cutoff}^n` is returned. - - This tensor can either be flattened into a vector or returned as an actual - tensor with :math:`n` indices. - - .. note:: - - Note that if :math:`R = (1)` then :math:`H_k^{(R)}(y)` - are precisely the well known **probabilists' Hermite polynomials** :math:`He_k(y)`, - and if :math:`R = (2)` then :math:`H_k^{(R)}(y)` are precisely the well known - **physicists' Hermite polynomials** :math:`H_k(y)`. - - Args: - R (array): square matrix parametrizing the Hermite polynomial family - cutoff (int): maximum size of the subindices in the Hermite polynomial - y (array): vector argument of the Hermite polynomial - renorm (bool): If ``True``, normalizes the returned multidimensional Hermite - polynomials such that :math:`H_k^{(R)}(y)/\prod_i k_i!` - make_tensor (bool): If ``False``, returns a flattened one dimensional array - containing the values of the polynomial - modified (bool): whether to return the modified multidimensional Hermite polynomials or the standard ones - rtol (float): the relative tolerance parameter used in ``np.allclose`` - atol (float): the absolute tolerance parameter used in ``np.allclose`` - Returns: - (array): the multidimensional Hermite polynomials - """ - - input_validation(R, atol=atol, rtol=rtol) - n, _ = R.shape - - if (modified is False) and (y is not None): - m = y.shape[0] - if m == n: - ym = R @ y - return hermite_multidimensional( - R, cutoff, y=ym, renorm=renorm, make_tensor=make_tensor, modified=True - ) - - if y is None: - y = np.zeros([n], dtype=complex) - - m = y.shape[0] - if m != n: - raise ValueError("The matrix R and vector y have incompatible dimensions") - - - Rt = np.real_if_close(R) - yt = np.real_if_close(y) - - if Rt.dtype == np.float and yt.dtype == np.float: - if renorm: - values = np.array(rhmr(Rt, yt, cutoff)) - else: - values = np.array(hmr(Rt, yt, cutoff)) - else: - if renorm: - values = np.array(rhm(np.complex128(R), np.complex128(y), cutoff)) - else: - values = np.array(hm(np.complex128(R), np.complex128(y), cutoff)) - - if make_tensor: - shape = cutoff * np.ones([n], dtype=int) - values = np.reshape(values, shape) - - return values - - -# pylint: disable=too-many-arguments -def hafnian_batched(A, cutoff, mu=None, rtol=1e-05, atol=1e-08, renorm=False, make_tensor=True): - r"""Calculates the hafnian of :func:`reduction(A, k) ` - for all possible values of vector ``k`` below the specified cutoff. - - Here, - - * :math:`A` is am :math:`n\times n` square matrix - * :math:`k` is a vector of (non-negative) integers with the same dimensions as :math:`A`, - i.e., :math:`k = (k_0,k_1,\ldots,k_{n-1})`, and where :math:`0 \leq k_j < \texttt{cutoff}`. - - The function :func:`~.hafnian_repeated` can be used to calculate the reduced hafnian - for a *specific* value of :math:`k`; see the documentation for more information. - - .. note:: - - If ``mu`` is not ``None``, this function instead returns - ``hafnian(np.fill_diagonal(reduction(A, k), reduction(mu, k)), loop=True)``. - This calculation can only be performed if the matrix :math:`A` is invertible. - - Args: - A (array): a square, symmetric :math:`N\times N` array. - cutoff (int): maximum size of the subindices in the Hermite polynomial - mu (array): a vector of length :math:`N` representing the vector of means/displacement - renorm (bool): If ``True``, the returned hafnians are *normalized*, that is, - :math:`haf(reduction(A, k))/\ \sqrt{prod_i k_i!}` - (or :math:`lhaf(fill\_diagonal(reduction(A, k),reduction(mu, k)))` if - ``mu`` is not None) - make_tensor: If ``False``, returns a flattened one dimensional array instead - of a tensor with the values of the hafnians. - rtol (float): the relative tolerance parameter used in ``np.allclose``. - atol (float): the absolute tolerance parameter used in ``np.allclose``. - Returns: - (array): the values of the hafnians for each value of :math:`k` up to the cutoff - """ - input_validation(A, atol=atol, rtol=rtol) - n, _ = A.shape - - if mu is None: - mu = np.zeros([n], dtype=complex) - - return hermite_multidimensional( - -A, cutoff, y=mu, renorm=renorm, make_tensor=make_tensor, modified=True - ) -# Note the minus signs in the arguments. Those are intentional and are due to the fact that Dodonov et al. in PRA 50, 813 (1994) use (p,q) ordering instead of (q,p) ordering From abb44319c98086b0d87bb93b5f20e1cd1f4877a0 Mon Sep 17 00:00:00 2001 From: Nicolas Quesada Date: Tue, 22 Jun 2021 15:30:48 -0400 Subject: [PATCH 24/29] Apply suggestions from code review --- thewalrus/__init__.py | 1 - 1 file changed, 1 deletion(-) diff --git a/thewalrus/__init__.py b/thewalrus/__init__.py index a28548e18..f7186f266 100644 --- a/thewalrus/__init__.py +++ b/thewalrus/__init__.py @@ -108,7 +108,6 @@ hafnian, hafnian_repeated, reduction, - sparse_hafnian hafnian_sparse, hafnian_banded, ) From 5b997d5d43f1921c53be11cabe853f40023c4409 Mon Sep 17 00:00:00 2001 From: ziofil Date: Tue, 22 Jun 2021 20:58:46 -0400 Subject: [PATCH 25/29] Update thewalrus/_hafnian.py Co-authored-by: Nicolas Quesada --- thewalrus/_hafnian.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/thewalrus/_hafnian.py b/thewalrus/_hafnian.py index 85d1ac3bf..7beb0c8a2 100644 --- a/thewalrus/_hafnian.py +++ b/thewalrus/_hafnian.py @@ -235,7 +235,7 @@ def hafnian_sparse(A, D=None, loop=False): if not loop: A = A - np.diag(np.diag(A)) - if np.allclose(A, np.zeros_like(A)): + if np.allclose(A, 0): return 0.0 r, _ = np.nonzero(A) From aa442e7fccf19f58ad8114e3487c459528c2b217 Mon Sep 17 00:00:00 2001 From: Nicolas Quesada Date: Tue, 22 Jun 2021 21:06:23 -0400 Subject: [PATCH 26/29] Apply suggestions from code review --- .github/CHANGELOG.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/CHANGELOG.md b/.github/CHANGELOG.md index a48d56dd6..dfb448326 100644 --- a/.github/CHANGELOG.md +++ b/.github/CHANGELOG.md @@ -24,7 +24,7 @@ This release contains contributions from (in alphabetical order): -Jake Bulmer, Timjan Kalajdzievski, Filippo Miatto,, Nicolas Quesada +Jake Bulmer, Timjan Kalajdzievski, Filippo Miatto, Nicolas Quesada --- From 2f7aa5f6444e5ddb9f0a54a89daf4b475a82778a Mon Sep 17 00:00:00 2001 From: Nicolas Quesada Date: Thu, 24 Jun 2021 14:06:24 -0400 Subject: [PATCH 27/29] Adds the new algorithms into the docs --- docs/algorithms.rst | 15 ++++++++++++++- thewalrus/__init__.py | 7 ++++++- thewalrus/_hafnian.py | 2 +- 3 files changed, 21 insertions(+), 3 deletions(-) diff --git a/docs/algorithms.rst b/docs/algorithms.rst index e2bb3d486..97a1f084a 100644 --- a/docs/algorithms.rst +++ b/docs/algorithms.rst @@ -202,4 +202,17 @@ Consider now the :math:`r`-partitions of the integer :math:`n`, there are :math: .. tip:: - *Implemented as* :func:`thewalrus.low_rank_hafnian`. This function takes as argument the matrix :math:`\bm{G} \in \mathbb{C}^{n \times r}` and returns the value of the hafnian of the matrix :math:`\bm{A} = \bm{G} \bm{G}^T` \ No newline at end of file + *Implemented as* :func:`thewalrus.low_rank_hafnian`. This function takes as argument the matrix :math:`\bm{G} \in \mathbb{C}^{n \times r}` and returns the value of the hafnian of the matrix :math:`\bm{A} = \bm{G} \bm{G}^T` + +Sparse and banded algorithms +---------------------------- +These algorithms take advantage of the Laplace expansion of the hafnian to calculate hafnians of matrices with many zeros. These matrices can be either banded or sparse. The Laplace expansion of the hafnian is given by + +.. math:: + \text{haf}(\bm{B}) = \sum_{j \neq i} B_{i,j} \text{haf}(\bm{B}_{-i-j}), + +where :math:`j` is a fixed index and :math:`\bm{B}_{-i-j}` is the matrix obtained from :math:`\bm{B}` by removing rows and columns :math:`i` and :math:`j`. The banded algorithm was introduced in Sec. V of Qi *et al* :cite:`qi2020efficient`. + +.. tip:: + + *Implemented as* :func:`thewalrus.hafnian_sparse` and :func:`thewalrus.hafnian_banded`. The loop hafnian calculation can be done by setting the option ``loop=True``. \ No newline at end of file diff --git a/thewalrus/__init__.py b/thewalrus/__init__.py index f7186f266..e1eb3fc56 100644 --- a/thewalrus/__init__.py +++ b/thewalrus/__init__.py @@ -1,4 +1,4 @@ -# Copyright 2019 Xanadu Quantum Technologies Inc. +# 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. @@ -69,6 +69,11 @@ *Efficient sampling from shallow Gaussian quantum-optical circuits with local interactions*, :cite:`qi2020efficient`. +Sparse hafnian algorithm + An algorithm that calculates the hafnian of a sparse matrix by taking advantage of the Laplace expansion and memoization, to store + only the relevant paths that contribute non-zero values to the final calculation. + + Python wrappers --------------- diff --git a/thewalrus/_hafnian.py b/thewalrus/_hafnian.py index 7beb0c8a2..a1fbb0584 100644 --- a/thewalrus/_hafnian.py +++ b/thewalrus/_hafnian.py @@ -1,4 +1,4 @@ -# Copyright 2019 Xanadu Quantum Technologies Inc. +# 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. From 66f2b222a2c147816ec74816240da04b15d5878e Mon Sep 17 00:00:00 2001 From: Nicolas Quesada Date: Thu, 24 Jun 2021 14:09:28 -0400 Subject: [PATCH 28/29] Adds period --- docs/algorithms.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/algorithms.rst b/docs/algorithms.rst index 97a1f084a..8c57ac483 100644 --- a/docs/algorithms.rst +++ b/docs/algorithms.rst @@ -202,7 +202,7 @@ Consider now the :math:`r`-partitions of the integer :math:`n`, there are :math: .. tip:: - *Implemented as* :func:`thewalrus.low_rank_hafnian`. This function takes as argument the matrix :math:`\bm{G} \in \mathbb{C}^{n \times r}` and returns the value of the hafnian of the matrix :math:`\bm{A} = \bm{G} \bm{G}^T` + *Implemented as* :func:`thewalrus.low_rank_hafnian`. This function takes as argument the matrix :math:`\bm{G} \in \mathbb{C}^{n \times r}` and returns the value of the hafnian of the matrix :math:`\bm{A} = \bm{G} \bm{G}^T`. Sparse and banded algorithms ---------------------------- From 3f6fb5cdeda8090f04dbec09304e77212a94b11b Mon Sep 17 00:00:00 2001 From: Nicolas Quesada Date: Thu, 24 Jun 2021 14:20:35 -0400 Subject: [PATCH 29/29] corrects typos --- thewalrus/__init__.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/thewalrus/__init__.py b/thewalrus/__init__.py index e1eb3fc56..13afb8f36 100644 --- a/thewalrus/__init__.py +++ b/thewalrus/__init__.py @@ -36,9 +36,9 @@ This algorithm scales like :math:`\mathcal{O}(n^4 2^{n/2})`. This algorithm does not currently support the loop hafnian. -Repeating hafnian algorithm +Repeated hafnian algorithm The algorithm described in *From moments of sum to moments of product*, :cite:`kan2008moments`. - This method is more efficient for matrices with repeated rows and columns, and supports caclulation of + This method is more efficient for matrices with repeated rows and columns, and supports calculation of the loop hafnian. Approximate hafnian algorithm