From 58a5bf43a6e069cdb361427bc138700ebb51a0b8 Mon Sep 17 00:00:00 2001 From: sylviemonet Date: Thu, 21 Mar 2024 15:02:08 +0800 Subject: [PATCH 01/30] add real_gaussian_integral --- mrmustard/physics/utils.py | 78 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 78 insertions(+) create mode 100644 mrmustard/physics/utils.py diff --git a/mrmustard/physics/utils.py b/mrmustard/physics/utils.py new file mode 100644 index 000000000..36aed4176 --- /dev/null +++ b/mrmustard/physics/utils.py @@ -0,0 +1,78 @@ +# Copyright 2024 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. + +""" +This module contains helper functions for physics. +""" +import numpy as np + +from mrmustard import math + + +def real_gaussian_integral( + Abc: tuple, + idx: tuple[int, ...], +): + r"""Computes the Gaussian integral of the exponential of a real quadratic form. + The integral is defined as (note that in general we integrate over a subset of m dimensions): + + :math:`\int_{R^m} F(x) dx` + + where + + :math:`F(x) = \textrm{exp}(0.5 x^T A x + b^T x)` + + Here, ``x`` is an ``n``-dim real vector, ``A`` is an ``n x n`` real matrix, + ``b`` is an ``n``-dim real vector, ``c`` is a real scalar. The integral indices + are specified by ``idx``. + + Arguments: + A,b,c: the ``(A,b,c)`` triple + idx: the tuple of indices of the x variables + + Returns: + The ``(A,b,c)`` triple of the result of the integral. + """ + A, b, c = Abc + + n = len(idx) + + if not idx: + return A, b, c + not_idx = tuple(i for i in range(A.shape[-1]) if i not in idx) + + I = math.eye(n, dtype=A.dtype) + Z = math.zeros((n, n), dtype=A.dtype) + X = math.block([[Z, I], [I, Z]]) + M = math.gather(math.gather(A, idx, axis=-1), idx, axis=-2) + X + bM = math.gather(b, idx, axis=-1) + + not_idx = tuple(i for i in range(A.shape[-1]) if i not in idx) + if math.asnumpy(not_idx).shape != (0,): + D = math.gather(math.gather(A, idx, axis=-1), not_idx, axis=-2) + R = math.gather(math.gather(A, not_idx, axis=-1), not_idx, axis=-2) + bR = math.gather(b, not_idx, axis=-1) + A_post = R - math.matmul(D, math.inv(M), math.transpose(D)) + b_post = bR - math.matvec(D, math.solve(M, bM)) + else: + A_post = math.astensor([]) + b_post = math.astensor([]) + + c_post = ( + c + * math.sqrt((2 * np.pi) ** n / math.det(M)) + * math.exp(-0.5 * math.sum(bM * math.solve(M, bM))) + ) + + return A_post, b_post, c_post From 750fabc7d0e343b910455b56e9f927014fc674fe Mon Sep 17 00:00:00 2001 From: Yuan Date: Tue, 26 Mar 2024 07:17:48 +0000 Subject: [PATCH 02/30] add tests for real gaussian integral --- mrmustard/physics/utils.py | 8 +++----- tests/test_physics/test_utils.py | 22 ++++++++++++++++++++++ 2 files changed, 25 insertions(+), 5 deletions(-) create mode 100644 tests/test_physics/test_utils.py diff --git a/mrmustard/physics/utils.py b/mrmustard/physics/utils.py index 36aed4176..d643a2c3d 100644 --- a/mrmustard/physics/utils.py +++ b/mrmustard/physics/utils.py @@ -27,7 +27,8 @@ def real_gaussian_integral( r"""Computes the Gaussian integral of the exponential of a real quadratic form. The integral is defined as (note that in general we integrate over a subset of m dimensions): - :math:`\int_{R^m} F(x) dx` + .. :math:: + \int_{R^m} F(x) dx, where @@ -52,10 +53,7 @@ def real_gaussian_integral( return A, b, c not_idx = tuple(i for i in range(A.shape[-1]) if i not in idx) - I = math.eye(n, dtype=A.dtype) - Z = math.zeros((n, n), dtype=A.dtype) - X = math.block([[Z, I], [I, Z]]) - M = math.gather(math.gather(A, idx, axis=-1), idx, axis=-2) + X + M = math.gather(math.gather(A, idx, axis=-1), idx, axis=-2) bM = math.gather(b, idx, axis=-1) not_idx = tuple(i for i in range(A.shape[-1]) if i not in idx) diff --git a/tests/test_physics/test_utils.py b/tests/test_physics/test_utils.py new file mode 100644 index 000000000..d92d8f717 --- /dev/null +++ b/tests/test_physics/test_utils.py @@ -0,0 +1,22 @@ +import numpy as np + +from mrmustard import math +from mrmustard.physics.utils import real_gaussian_integral + + +def test_real_gaussian_integral(): + """Tests the ``real_gaussian_integral`` method.""" + A0 = np.random.random((3, 3)) + A = (A0 + A0.T) / 2 + b = np.arange(3) + c = 1.0 + res = real_gaussian_integral((A, b, c), idx=[0, 1]) + assert np.allclose(res[0], A[2, 2] - A[2:, :2] @ np.linalg.inv(A[:2, :2]) @ A[:2, 2:]) + assert np.allclose(res[1], b[2] - A[2:, :2] @ np.linalg.inv(A[:2, :2]) @ b[:2]) + assert np.allclose( + res[2], + c + * (2 * np.pi) + / np.sqrt(np.linalg.det(A[:2, :2])) + * np.exp(-0.5 * b[:2] @ np.linalg.inv(A[:2, :2]) @ b[:2]), + ) From 9e74514817aabb8d7d240b8e29f0740a30f0f0d9 Mon Sep 17 00:00:00 2001 From: Yuan Date: Tue, 26 Mar 2024 07:39:12 +0000 Subject: [PATCH 03/30] fix bugs and add changelog --- .github/CHANGELOG.md | 3 +++ tests/test_physics/test_utils.py | 18 ++++++++++++++++-- 2 files changed, 19 insertions(+), 2 deletions(-) diff --git a/.github/CHANGELOG.md b/.github/CHANGELOG.md index bffd66c27..fb6b79867 100644 --- a/.github/CHANGELOG.md +++ b/.github/CHANGELOG.md @@ -4,6 +4,9 @@ * Added a function ``to_fock`` to map different representations into Fock representation. [(#355)](https://github.com/XanaduAI/MrMustard/pull/355) +* Added a function ``real_gaussian_integral`` as helper function to map between different representations. + [(#371)](https://github.com/XanaduAI/MrMustard/pull/371) + ### Breaking changes ### Improvements diff --git a/tests/test_physics/test_utils.py b/tests/test_physics/test_utils.py index d92d8f717..6b9a174f0 100644 --- a/tests/test_physics/test_utils.py +++ b/tests/test_physics/test_utils.py @@ -1,6 +1,20 @@ -import numpy as np +# Copyright 2024 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 -from mrmustard import math +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +"""Tests for the helper functions in physics.""" + +import numpy as np from mrmustard.physics.utils import real_gaussian_integral From 84c3723d1f7c1445d1f66500d409f83150e81fc5 Mon Sep 17 00:00:00 2001 From: Yuan Date: Tue, 26 Mar 2024 07:42:46 +0000 Subject: [PATCH 04/30] add hardcoded A matrix in the tests --- tests/test_physics/test_utils.py | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/tests/test_physics/test_utils.py b/tests/test_physics/test_utils.py index 6b9a174f0..533d888ab 100644 --- a/tests/test_physics/test_utils.py +++ b/tests/test_physics/test_utils.py @@ -19,9 +19,14 @@ def test_real_gaussian_integral(): - """Tests the ``real_gaussian_integral`` method.""" - A0 = np.random.random((3, 3)) - A = (A0 + A0.T) / 2 + """Tests the ``real_gaussian_integral`` method with a hard-coded A matric from a Gaussian(3) state.""" + A = np.array( + [ + [0.35307718 - 0.09738001j, -0.01297994 + 0.26050244j, 0.05349344 - 0.13728068j], + [-0.01297994 + 0.26050244j, 0.05696707 - 0.2351408j, 0.18954838 - 0.42959383j], + [0.05349344 - 0.13728068j, 0.18954838 - 0.42959383j, -0.16931712 - 0.09205837j], + ] + ) b = np.arange(3) c = 1.0 res = real_gaussian_integral((A, b, c), idx=[0, 1]) From a8b59322b4b1c2da851f0bad2dcb5e35731e43bf Mon Sep 17 00:00:00 2001 From: Yuan Date: Tue, 26 Mar 2024 08:50:17 +0000 Subject: [PATCH 05/30] fix bugs in tests --- tests/test_physics/test_utils.py | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/tests/test_physics/test_utils.py b/tests/test_physics/test_utils.py index 533d888ab..d9f2b2236 100644 --- a/tests/test_physics/test_utils.py +++ b/tests/test_physics/test_utils.py @@ -15,27 +15,28 @@ """Tests for the helper functions in physics.""" import numpy as np +from mrmustard import math from mrmustard.physics.utils import real_gaussian_integral def test_real_gaussian_integral(): """Tests the ``real_gaussian_integral`` method with a hard-coded A matric from a Gaussian(3) state.""" - A = np.array( + A = math.astensor(np.array( [ [0.35307718 - 0.09738001j, -0.01297994 + 0.26050244j, 0.05349344 - 0.13728068j], [-0.01297994 + 0.26050244j, 0.05696707 - 0.2351408j, 0.18954838 - 0.42959383j], [0.05349344 - 0.13728068j, 0.18954838 - 0.42959383j, -0.16931712 - 0.09205837j], ] - ) - b = np.arange(3) - c = 1.0 + )) + b = math.astensor(np.arange(3)+0j) + c = 1.0 + 0j res = real_gaussian_integral((A, b, c), idx=[0, 1]) - assert np.allclose(res[0], A[2, 2] - A[2:, :2] @ np.linalg.inv(A[:2, :2]) @ A[:2, 2:]) - assert np.allclose(res[1], b[2] - A[2:, :2] @ np.linalg.inv(A[:2, :2]) @ b[:2]) + assert np.allclose(res[0], A[2, 2] - A[2:, :2] @ math.inv(A[:2, :2]) @ A[:2, 2:]) + assert np.allclose(res[1], b[2] - math.sum(A[2:, :2] * math.matvec(math.inv(A[:2, :2]), b[:2]))) assert np.allclose( res[2], c * (2 * np.pi) - / np.sqrt(np.linalg.det(A[:2, :2])) - * np.exp(-0.5 * b[:2] @ np.linalg.inv(A[:2, :2]) @ b[:2]), + / math.sqrt(math.det(A[:2, :2])) + * math.exp(-0.5 * math.sum(b[:2] * math.matvec(math.inv(A[:2, :2]), b[:2]))), ) From 6eb0c3c374e42a8c4ff9d3e936267c47763f4d29 Mon Sep 17 00:00:00 2001 From: Yuan Date: Tue, 26 Mar 2024 08:51:07 +0000 Subject: [PATCH 06/30] blacked --- tests/test_physics/test_utils.py | 18 ++++++++++-------- 1 file changed, 10 insertions(+), 8 deletions(-) diff --git a/tests/test_physics/test_utils.py b/tests/test_physics/test_utils.py index d9f2b2236..ae3f31973 100644 --- a/tests/test_physics/test_utils.py +++ b/tests/test_physics/test_utils.py @@ -21,14 +21,16 @@ def test_real_gaussian_integral(): """Tests the ``real_gaussian_integral`` method with a hard-coded A matric from a Gaussian(3) state.""" - A = math.astensor(np.array( - [ - [0.35307718 - 0.09738001j, -0.01297994 + 0.26050244j, 0.05349344 - 0.13728068j], - [-0.01297994 + 0.26050244j, 0.05696707 - 0.2351408j, 0.18954838 - 0.42959383j], - [0.05349344 - 0.13728068j, 0.18954838 - 0.42959383j, -0.16931712 - 0.09205837j], - ] - )) - b = math.astensor(np.arange(3)+0j) + A = math.astensor( + np.array( + [ + [0.35307718 - 0.09738001j, -0.01297994 + 0.26050244j, 0.05349344 - 0.13728068j], + [-0.01297994 + 0.26050244j, 0.05696707 - 0.2351408j, 0.18954838 - 0.42959383j], + [0.05349344 - 0.13728068j, 0.18954838 - 0.42959383j, -0.16931712 - 0.09205837j], + ] + ) + ) + b = math.astensor(np.arange(3) + 0j) c = 1.0 + 0j res = real_gaussian_integral((A, b, c), idx=[0, 1]) assert np.allclose(res[0], A[2, 2] - A[2:, :2] @ math.inv(A[:2, :2]) @ A[:2, 2:]) From 547a0ac3de16c359c408502c59886afb40323106 Mon Sep 17 00:00:00 2001 From: sylviemonet Date: Wed, 27 Mar 2024 17:47:29 +0800 Subject: [PATCH 07/30] ad more tests --- tests/test_physics/test_utils.py | 25 +++++++++++++++++++++++++ 1 file changed, 25 insertions(+) diff --git a/tests/test_physics/test_utils.py b/tests/test_physics/test_utils.py index ae3f31973..af4544e86 100644 --- a/tests/test_physics/test_utils.py +++ b/tests/test_physics/test_utils.py @@ -42,3 +42,28 @@ def test_real_gaussian_integral(): / math.sqrt(math.det(A[:2, :2])) * math.exp(-0.5 * math.sum(b[:2] * math.matvec(math.inv(A[:2, :2]), b[:2]))), ) + res2 = real_gaussian_integral((A, b, c), idx=[]) + assert np.allclose(res2[0], A) + assert np.allclose(res2[1], b) + assert np.allclose(res2[2], c) + + A2 = math.astensor( + np.array( + [ + [0.35307718 - 0.09738001j, -0.01297994 + 0.26050244j], + [-0.01297994 + 0.26050244j, 0.05696707 - 0.2351408j], + ] + ) + ) + b2 = math.astensor(np.arange(2) + 0j) + c2 = 1.0 + 0j + res3 = real_gaussian_integral((A2, b2, c2), idx=[0, 1]) + assert np.allclose(res3[0], math.astensor([])) + assert np.allclose(res3[1], math.astensor([])) + assert np.allclose( + res3[2], + c2 + * (2 * np.pi) + / math.sqrt(math.det(A2[:2, :2])) + * math.exp(-0.5 * math.sum(b2[:2] * math.matvec(math.inv(A2[:2, :2]), b2[:2]))), + ) From 73911983387450c78494cafd07f13703777a0613 Mon Sep 17 00:00:00 2001 From: Yuan YAO Date: Thu, 11 Apr 2024 19:14:44 +0000 Subject: [PATCH 08/30] remove the redundency line --- mrmustard/physics/utils.py | 1 - 1 file changed, 1 deletion(-) diff --git a/mrmustard/physics/utils.py b/mrmustard/physics/utils.py index d643a2c3d..9b8497353 100644 --- a/mrmustard/physics/utils.py +++ b/mrmustard/physics/utils.py @@ -56,7 +56,6 @@ def real_gaussian_integral( M = math.gather(math.gather(A, idx, axis=-1), idx, axis=-2) bM = math.gather(b, idx, axis=-1) - not_idx = tuple(i for i in range(A.shape[-1]) if i not in idx) if math.asnumpy(not_idx).shape != (0,): D = math.gather(math.gather(A, idx, axis=-1), not_idx, axis=-2) R = math.gather(math.gather(A, not_idx, axis=-1), not_idx, axis=-2) From e4d2f163867089dced6f17229ccb8acd785d1784 Mon Sep 17 00:00:00 2001 From: Yuan YAO Date: Thu, 11 Apr 2024 19:19:49 +0000 Subject: [PATCH 09/30] changed the long math as Fillippo suggested --- mrmustard/physics/utils.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/mrmustard/physics/utils.py b/mrmustard/physics/utils.py index 9b8497353..cf3970f04 100644 --- a/mrmustard/physics/utils.py +++ b/mrmustard/physics/utils.py @@ -60,8 +60,10 @@ def real_gaussian_integral( D = math.gather(math.gather(A, idx, axis=-1), not_idx, axis=-2) R = math.gather(math.gather(A, not_idx, axis=-1), not_idx, axis=-2) bR = math.gather(b, not_idx, axis=-1) - A_post = R - math.matmul(D, math.inv(M), math.transpose(D)) - b_post = bR - math.matvec(D, math.solve(M, bM)) + T = math.transpose + L = T(math.solve(T(M), T(D))) + A_post = R - math.matmul(L, T(D)) + b_post = bR - math.matvec(L, bM) else: A_post = math.astensor([]) b_post = math.astensor([]) From 93a49f2bdb20ea3ed772dadc4c126705c99b2ffd Mon Sep 17 00:00:00 2001 From: Yuan YAO Date: Thu, 11 Apr 2024 20:22:33 +0000 Subject: [PATCH 10/30] change the file name and move all related functions together --- mrmustard/physics/gaussian_integrals.py | 218 ++++++++++++++++++ tests/test_physics/test_gaussian_integrals.py | 169 ++++++++++++++ 2 files changed, 387 insertions(+) create mode 100644 mrmustard/physics/gaussian_integrals.py create mode 100644 tests/test_physics/test_gaussian_integrals.py diff --git a/mrmustard/physics/gaussian_integrals.py b/mrmustard/physics/gaussian_integrals.py new file mode 100644 index 000000000..2b274b2b8 --- /dev/null +++ b/mrmustard/physics/gaussian_integrals.py @@ -0,0 +1,218 @@ +# Copyright 2024 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. + +""" +This module contains real and comple gaussian integral functions and related helper functions. +""" +from typing import Sequence, Tuple + +import numpy as np + +from mrmustard import math +from mrmustard.utils.typing import ComplexMatrix, ComplexVector + + +def real_gaussian_integral( + Abc: tuple, + idx: tuple[int, ...], +): + r"""Computes the Gaussian integral of the exponential of a real quadratic form. + The integral is defined as (note that in general we integrate over a subset of m dimensions): + + .. :math:: + \int_{R^m} F(x) dx, + + where + + :math:`F(x) = \textrm{exp}(0.5 x^T A x + b^T x)` + + Here, ``x`` is an ``n``-dim real vector, ``A`` is an ``n x n`` real matrix, + ``b`` is an ``n``-dim real vector, ``c`` is a real scalar. The integral indices + are specified by ``idx``. + + Arguments: + A,b,c: the ``(A,b,c)`` triple + idx: the tuple of indices of the x variables + + Returns: + The ``(A,b,c)`` triple of the result of the integral. + """ + A, b, c = Abc + + n = len(idx) + + if not idx: + return A, b, c + not_idx = tuple(i for i in range(A.shape[-1]) if i not in idx) + + M = math.gather(math.gather(A, idx, axis=-1), idx, axis=-2) + bM = math.gather(b, idx, axis=-1) + + if math.asnumpy(not_idx).shape != (0,): + D = math.gather(math.gather(A, idx, axis=-1), not_idx, axis=-2) + R = math.gather(math.gather(A, not_idx, axis=-1), not_idx, axis=-2) + bR = math.gather(b, not_idx, axis=-1) + T = math.transpose + L = T(math.solve(T(M), T(D))) + A_post = R - math.matmul(L, T(D)) + b_post = bR - math.matvec(L, bM) + else: + A_post = math.astensor([]) + b_post = math.astensor([]) + + c_post = ( + c + * math.sqrt((2 * np.pi) ** n / math.det(M)) + * math.exp(-0.5 * math.sum(bM * math.solve(M, bM))) + ) + + return A_post, b_post, c_post + +def complex_gaussian_integral( + Abc: tuple, idx_z: tuple[int, ...], idx_zconj: tuple[int, ...], measure: float = -1 +): + r"""Computes the Gaussian integral of the exponential of a complex quadratic form. + The integral is defined as (note that in general we integrate over a subset of 2m dimensions): + + :math:`\int_{C^m} F(z) d\mu(z)` + + where + + :math:`F(z) = \textrm{exp}(-0.5 z^T A z + b^T z)` + + Here, ``z`` is an ``n``-dim complex vector, ``A`` is an ``n x n`` complex matrix, + ``b`` is an ``n``-dim complex vector, ``c`` is a complex scalar, and :math:`d\mu(z)` + is a non-holomorphic complex measure over a subset of m pairs of z,z* variables. These + are specified by the indices ``idx_z`` and ``idx_zconj``. The ``measure`` parameter is + the exponent of the measure: + + :math: `dmu(z) = \textrm{exp}(m * |z|^2) \frac{d^{2n}z}{\pi^n} = \frac{1}{\pi^n}\textrm{exp}(m * |z|^2) d\textrm{Re}(z) d\textrm{Im}(z)` + + Note that the indices must be a complex variable pairs with each other (idx_z, idx_zconj) to make this contraction meaningful. + Please make sure the corresponding complex variable with respect to your Abc triples. + For examples, if the indices of Abc denotes the variables ``(\alpha, \beta, \alpha^*, \beta^*, \gamma, \eta)``, the contraction only works + with the indices between ``(\alpha, \alpha^*)`` pairs and ``(\beta, \beta^*)`` pairs. + + Arguments: + A,b,c: the ``(A,b,c)`` triple + idx_z: the tuple of indices of the z variables + idx_zconj: the tuple of indices of the z* variables + measure: the exponent of the measure (default is -1: Bargmann measure) + + Returns: + The ``(A,b,c)`` triple of the result of the integral. + + Raises: + ValueError: If ``idx_z`` and ``idx_zconj`` have different lengths. + """ + A, b, c = Abc + if len(idx_z) != len(idx_zconj): + raise ValueError( + f"idx_z and idx_zconj must have the same length, got {len(idx_z)} and {len(idx_zconj)}" + ) + n = len(idx_z) + idx = tuple(idx_z) + tuple(idx_zconj) + if not idx: + return A, b, c + not_idx = tuple(i for i in range(A.shape[-1]) if i not in idx) + + I = math.eye(n, dtype=A.dtype) + Z = math.zeros((n, n), dtype=A.dtype) + X = math.block([[Z, I], [I, Z]]) + M = math.gather(math.gather(A, idx, axis=-1), idx, axis=-2) + X * measure + bM = math.gather(b, idx, axis=-1) + + not_idx = tuple(i for i in range(A.shape[-1]) if i not in idx) + if math.asnumpy(not_idx).shape != (0,): + D = math.gather(math.gather(A, idx, axis=-1), not_idx, axis=-2) + R = math.gather(math.gather(A, not_idx, axis=-1), not_idx, axis=-2) + bR = math.gather(b, not_idx, axis=-1) + A_post = R - math.matmul(D, math.inv(M), math.transpose(D)) + b_post = bR - math.matvec(D, math.solve(M, bM)) + else: + A_post = math.astensor([]) + b_post = math.astensor([]) + + c_post = ( + c * math.sqrt((-1) ** n / math.det(M)) * math.exp(-0.5 * math.sum(bM * math.solve(M, bM))) + ) + + return A_post, b_post, c_post + + +def join_Abc(Abc1, Abc2): + r"""Joins two ``(A,b,c)`` triples into a single ``(A,b,c)`` triple by block addition of the ``A`` + matrices and concatenating the ``b`` vectors. + + Arguments: + Abc1: the first ``(A,b,c)`` triple + Abc2: the second ``(A,b,c)`` triple + + Returns: + The joined ``(A,b,c)`` triple + """ + A1, b1, c1 = Abc1 + A2, b2, c2 = Abc2 + A12 = math.block_diag(A1, A2) + b12 = math.concat([b1, b2], axis=-1) + c12 = math.outer(c1, c2) + return A12, b12, c12 + + +def reorder_abc(Abc: tuple, order: Sequence[int]): + r""" + Reorders the indices of the A matrix and b vector of an (A,b,c) triple. + + Arguments: + Abc: the ``(A,b,c)`` triple + order: the new order of the indices + + Returns: + The reordered ``(A,b,c)`` triple + """ + A, b, c = Abc + A = math.gather(math.gather(A, order, axis=-1), order, axis=-2) + b = math.gather(b, order, axis=-1) + if len(c.shape) == len(order): + c = math.transpose(c, order) + return A, b, c + + +def contract_two_Abc( + Abc1: Tuple[ComplexMatrix, ComplexVector, complex], + Abc2: Tuple[ComplexMatrix, ComplexVector, complex], + idx1: Sequence[int], + idx2: Sequence[int], +): + r""" + Returns the contraction of two ``(A,b,c)`` triples with given indices. + + Note that the indices must be a complex variable pairs with each other to make this contraction meaningful. Please make sure + the corresponding complex variable with respect to your Abc triples. + For examples, if the indices of Abc1 denotes the variables ``(\alpha, \beta)``, the indices of Abc2 denotes the variables + ``(\alpha^*,\gamma)``, the contraction only works with ``idx1 = [0], idx2 = [0]``. + + Arguments: + Abc1: the first ``(A,b,c)`` triple + Abc2: the second ``(A,b,c)`` triple + idx1: the indices of the first ``(A,b,c)`` triple to contract + idx2: the indices of the second ``(A,b,c)`` triple to contract + + Returns: + The contracted ``(A,b,c)`` triple + """ + Abc = join_Abc(Abc1, Abc2) + return complex_gaussian_integral( + Abc, idx1, tuple(n + Abc1[0].shape[-1] for n in idx2), measure=-1.0 + ) \ No newline at end of file diff --git a/tests/test_physics/test_gaussian_integrals.py b/tests/test_physics/test_gaussian_integrals.py new file mode 100644 index 000000000..c46f29b4d --- /dev/null +++ b/tests/test_physics/test_gaussian_integrals.py @@ -0,0 +1,169 @@ +# Copyright 2024 Xanadu Quantum Technologies Inc. + +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at + +# http://www.apache.org/licenses/LICENSE-2.0 + +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +"""Tests for real and comple gaussian integral functions and related helper functions.""" + +import numpy as np +import pytest +from mrmustard import math +from mrmustard.physics import triples +from mrmustard.physics.gaussian_integrals import ( + real_gaussian_integral, + complex_gaussian_integral, + join_Abc, + contract_two_Abc, + reorder_abc, +) + + +def test_real_gaussian_integral(): + """Tests the ``real_gaussian_integral`` method with a hard-coded A matric from a Gaussian(3) state.""" + A = math.astensor( + np.array( + [ + [0.35307718 - 0.09738001j, -0.01297994 + 0.26050244j, 0.05349344 - 0.13728068j], + [-0.01297994 + 0.26050244j, 0.05696707 - 0.2351408j, 0.18954838 - 0.42959383j], + [0.05349344 - 0.13728068j, 0.18954838 - 0.42959383j, -0.16931712 - 0.09205837j], + ] + ) + ) + b = math.astensor(np.arange(3) + 0j) + c = 1.0 + 0j + res = real_gaussian_integral((A, b, c), idx=[0, 1]) + assert np.allclose(res[0], A[2, 2] - A[2:, :2] @ math.inv(A[:2, :2]) @ A[:2, 2:]) + assert np.allclose(res[1], b[2] - math.sum(A[2:, :2] * math.matvec(math.inv(A[:2, :2]), b[:2]))) + assert np.allclose( + res[2], + c + * (2 * np.pi) + / math.sqrt(math.det(A[:2, :2])) + * math.exp(-0.5 * math.sum(b[:2] * math.matvec(math.inv(A[:2, :2]), b[:2]))), + ) + res2 = real_gaussian_integral((A, b, c), idx=[]) + assert np.allclose(res2[0], A) + assert np.allclose(res2[1], b) + assert np.allclose(res2[2], c) + + A2 = math.astensor( + np.array( + [ + [0.35307718 - 0.09738001j, -0.01297994 + 0.26050244j], + [-0.01297994 + 0.26050244j, 0.05696707 - 0.2351408j], + ] + ) + ) + b2 = math.astensor(np.arange(2) + 0j) + c2 = 1.0 + 0j + res3 = real_gaussian_integral((A2, b2, c2), idx=[0, 1]) + assert np.allclose(res3[0], math.astensor([])) + assert np.allclose(res3[1], math.astensor([])) + assert np.allclose( + res3[2], + c2 + * (2 * np.pi) + / math.sqrt(math.det(A2[:2, :2])) + * math.exp(-0.5 * math.sum(b2[:2] * math.matvec(math.inv(A2[:2, :2]), b2[:2]))), + ) + + +def test_join_Abc(): + """Tests the ``join_Abc`` method.""" + A1, b1, c1 = triples.vacuum_state_Abc(2) + A2, b2, c2 = triples.displacement_gate_Abc(x=[0.1, 0.2], y=0.3) + + joined_Abc = join_Abc((A1, b1, c1), (A2, b2, c2)) + assert np.allclose(joined_Abc[0], math.block_diag(A1, A2)) + assert np.allclose(joined_Abc[1], math.concat([b1, b2], axis=-1)) + assert np.allclose(joined_Abc[2], math.outer(c1, c2)) + + A12 = math.block_diag(A1, A2) + b12 = math.concat([b1, b2], axis=-1) + c12 = math.outer(c1, c2) + return A12, b12, c12 + + +def test_complex_gaussian_integral(): + """Tests the ``complex_gaussian_integral`` method.""" + A1, b1, c1 = triples.vacuum_state_Abc(2) + A2, b2, c2 = triples.displacement_gate_Abc(x=[0.1, 0.2], y=0.3) + A3, b3, c3 = triples.displaced_squeezed_vacuum_state_Abc(x=[0.1, 0.2], y=0.3) + + joined_Abc = join_Abc((A1, b1, c1), (A2, b2, c2)) + + res1 = complex_gaussian_integral(joined_Abc, [], []) + assert np.allclose(res1[0], joined_Abc[0]) + assert np.allclose(res1[1], joined_Abc[1]) + assert np.allclose(res1[2], joined_Abc[2]) + + res2 = complex_gaussian_integral(joined_Abc, [0, 1], [4, 5]) + assert np.allclose(res2[0], A3) + assert np.allclose(res2[1], b3) + assert np.allclose(res2[2], c3) + + res3 = complex_gaussian_integral(join_Abc((A1, b1, c1), (A1, b1, c1)), [0, 1], [2, 3]) + assert np.allclose(res3[0], 0) + assert np.allclose(res3[1], 0) + assert np.allclose(res3[2], 1) + + +def test_complex_gaussian_integral_error(): + """Tests the error of the ``complex_gaussian_integral`` method.""" + A1, b1, c1 = triples.vacuum_state_Abc(2) + A2, b2, c2 = triples.displacement_gate_Abc(x=[0.1, 0.2], y=0.3) + + with pytest.raises(ValueError): + complex_gaussian_integral( + join_Abc((A1, b1, c1), (A2, b2, c2)), + [0, 1], + [ + 4, + ], + ) + + +def test_contract_two_Abc(): + """Tests the error of the ``contract_two_Abc`` method.""" + A1, b1, c1 = triples.vacuum_state_Abc(2) + A2, b2, c2 = triples.displacement_gate_Abc(x=[0.1, 0.2], y=0.3) + + res1 = contract_two_Abc((A1, b1, c1), (A2, b2, c2), [], []) + assert np.allclose(res1[0], math.block_diag(A1, A2)) + assert np.allclose(res1[1], [0, 0, 0.1 + 0.3j, 0.2 + 0.3j, -0.1 + 0.3j, -0.2 + 0.3j]) + assert np.allclose(res1[2], c1 * c2) + + res2 = contract_two_Abc((A1, b1, c1), (A2, b2, c2), [0, 1], [2, 3]) + assert np.allclose(res2[0], math.zeros((2, 2))) + assert np.allclose(res2[1], [0.1 + 0.3j, 0.2 + 0.3j]) + assert np.allclose(res2[2], c1 * c2) + + res3 = contract_two_Abc((A1, b1, c1), (A2, b2, c2), [0, 1], [0, 1]) + assert np.allclose(res3[0], math.zeros((2, 2))) + assert np.allclose(res3[1], [-0.1 + 0.3j, -0.2 + 0.3j]) + assert np.allclose(res3[2], c1 * c2) + + +def test_reorder_abc(): + """Test that the reorder_abc function works correctly""" + A = np.array([[1, 2], [2, 3]]) + b = np.array([4, 5]) + c = np.array(6) + same = reorder_abc((A, b, c), (0, 1)) + assert all(np.allclose(x, y) for x, y in zip(same, (A, b, c))) + flipped = reorder_abc((A, b, c), (1, 0)) + assert all(np.allclose(x, y) for x, y in zip(flipped, (A[[1, 0], :][:, [1, 0]], b[[1, 0]], c))) + c = np.array([[6, 7], [8, 9]]) + flipped = reorder_abc((A, b, c), (1, 0)) # test transposition of c + assert all( + np.allclose(x, y) for x, y in zip(flipped, (A[[1, 0], :][:, [1, 0]], b[[1, 0]], c.T)) + ) From 09805e6244fdcb8d259ce06bc259428c57a12bcc Mon Sep 17 00:00:00 2001 From: Yuan YAO Date: Thu, 11 Apr 2024 20:23:26 +0000 Subject: [PATCH 11/30] change things accordingly inside bargmann --- mrmustard/physics/bargmann.py | 142 ---------------------------- mrmustard/physics/utils.py | 77 --------------- tests/test_physics/test_bargmann.py | 99 ------------------- tests/test_physics/test_utils.py | 69 -------------- 4 files changed, 387 deletions(-) delete mode 100644 mrmustard/physics/utils.py delete mode 100644 tests/test_physics/test_utils.py diff --git a/mrmustard/physics/bargmann.py b/mrmustard/physics/bargmann.py index 922330082..a71e51fc0 100644 --- a/mrmustard/physics/bargmann.py +++ b/mrmustard/physics/bargmann.py @@ -15,13 +15,10 @@ """ This module contains functions for performing calculations on objects in the Bargmann representations. """ -from typing import Sequence, Tuple - import numpy as np from mrmustard import math, settings from mrmustard.physics.husimi import pq_to_aadag, wigner_to_husimi -from mrmustard.utils.typing import ComplexMatrix, ComplexVector def cayley(X, c): @@ -104,142 +101,3 @@ def wigner_to_bargmann_U(X, d): N = X.shape[-1] // 2 A, B, C = wigner_to_bargmann_Choi(X, math.zeros_like(X), d) return A[2 * N :, 2 * N :], B[2 * N :], math.sqrt(C) - - -def complex_gaussian_integral( - Abc: tuple, idx_z: tuple[int, ...], idx_zconj: tuple[int, ...], measure: float = -1 -): - r"""Computes the Gaussian integral of the exponential of a complex quadratic form. - The integral is defined as (note that in general we integrate over a subset of 2m dimensions): - - :math:`\int_{C^m} F(z) d\mu(z)` - - where - - :math:`F(z) = \textrm{exp}(-0.5 z^T A z + b^T z)` - - Here, ``z`` is an ``n``-dim complex vector, ``A`` is an ``n x n`` complex matrix, - ``b`` is an ``n``-dim complex vector, ``c`` is a complex scalar, and :math:`d\mu(z)` - is a non-holomorphic complex measure over a subset of m pairs of z,z* variables. These - are specified by the indices ``idx_z`` and ``idx_zconj``. The ``measure`` parameter is - the exponent of the measure: - - :math: `dmu(z) = \textrm{exp}(m * |z|^2) \frac{d^{2n}z}{\pi^n} = \frac{1}{\pi^n}\textrm{exp}(m * |z|^2) d\textrm{Re}(z) d\textrm{Im}(z)` - - Note that the indices must be a complex variable pairs with each other (idx_z, idx_zconj) to make this contraction meaningful. - Please make sure the corresponding complex variable with respect to your Abc triples. - For examples, if the indices of Abc denotes the variables ``(\alpha, \beta, \alpha^*, \beta^*, \gamma, \eta)``, the contraction only works - with the indices between ``(\alpha, \alpha^*)`` pairs and ``(\beta, \beta^*)`` pairs. - - Arguments: - A,b,c: the ``(A,b,c)`` triple - idx_z: the tuple of indices of the z variables - idx_zconj: the tuple of indices of the z* variables - measure: the exponent of the measure (default is -1: Bargmann measure) - - Returns: - The ``(A,b,c)`` triple of the result of the integral. - - Raises: - ValueError: If ``idx_z`` and ``idx_zconj`` have different lengths. - """ - A, b, c = Abc - if len(idx_z) != len(idx_zconj): - raise ValueError( - f"idx_z and idx_zconj must have the same length, got {len(idx_z)} and {len(idx_zconj)}" - ) - n = len(idx_z) - idx = tuple(idx_z) + tuple(idx_zconj) - if not idx: - return A, b, c - not_idx = tuple(i for i in range(A.shape[-1]) if i not in idx) - - I = math.eye(n, dtype=A.dtype) - Z = math.zeros((n, n), dtype=A.dtype) - X = math.block([[Z, I], [I, Z]]) - M = math.gather(math.gather(A, idx, axis=-1), idx, axis=-2) + X * measure - bM = math.gather(b, idx, axis=-1) - - not_idx = tuple(i for i in range(A.shape[-1]) if i not in idx) - if math.asnumpy(not_idx).shape != (0,): - D = math.gather(math.gather(A, idx, axis=-1), not_idx, axis=-2) - R = math.gather(math.gather(A, not_idx, axis=-1), not_idx, axis=-2) - bR = math.gather(b, not_idx, axis=-1) - A_post = R - math.matmul(D, math.inv(M), math.transpose(D)) - b_post = bR - math.matvec(D, math.solve(M, bM)) - else: - A_post = math.astensor([]) - b_post = math.astensor([]) - - c_post = ( - c * math.sqrt((-1) ** n / math.det(M)) * math.exp(-0.5 * math.sum(bM * math.solve(M, bM))) - ) - - return A_post, b_post, c_post - - -def join_Abc(Abc1, Abc2): - r"""Joins two ``(A,b,c)`` triples into a single ``(A,b,c)`` triple by block addition of the ``A`` - matrices and concatenating the ``b`` vectors. - - Arguments: - Abc1: the first ``(A,b,c)`` triple - Abc2: the second ``(A,b,c)`` triple - - Returns: - The joined ``(A,b,c)`` triple - """ - A1, b1, c1 = Abc1 - A2, b2, c2 = Abc2 - A12 = math.block_diag(A1, A2) - b12 = math.concat([b1, b2], axis=-1) - c12 = math.outer(c1, c2) - return A12, b12, c12 - - -def reorder_abc(Abc: tuple, order: Sequence[int]): - r""" - Reorders the indices of the A matrix and b vector of an (A,b,c) triple. - - Arguments: - Abc: the ``(A,b,c)`` triple - order: the new order of the indices - - Returns: - The reordered ``(A,b,c)`` triple - """ - A, b, c = Abc - A = math.gather(math.gather(A, order, axis=-1), order, axis=-2) - b = math.gather(b, order, axis=-1) - if len(c.shape) == len(order): - c = math.transpose(c, order) - return A, b, c - - -def contract_two_Abc( - Abc1: Tuple[ComplexMatrix, ComplexVector, complex], - Abc2: Tuple[ComplexMatrix, ComplexVector, complex], - idx1: Sequence[int], - idx2: Sequence[int], -): - r""" - Returns the contraction of two ``(A,b,c)`` triples with given indices. - - Note that the indices must be a complex variable pairs with each other to make this contraction meaningful. Please make sure - the corresponding complex variable with respect to your Abc triples. - For examples, if the indices of Abc1 denotes the variables ``(\alpha, \beta)``, the indices of Abc2 denotes the variables - ``(\alpha^*,\gamma)``, the contraction only works with ``idx1 = [0], idx2 = [0]``. - - Arguments: - Abc1: the first ``(A,b,c)`` triple - Abc2: the second ``(A,b,c)`` triple - idx1: the indices of the first ``(A,b,c)`` triple to contract - idx2: the indices of the second ``(A,b,c)`` triple to contract - - Returns: - The contracted ``(A,b,c)`` triple - """ - Abc = join_Abc(Abc1, Abc2) - return complex_gaussian_integral( - Abc, idx1, tuple(n + Abc1[0].shape[-1] for n in idx2), measure=-1.0 - ) diff --git a/mrmustard/physics/utils.py b/mrmustard/physics/utils.py deleted file mode 100644 index cf3970f04..000000000 --- a/mrmustard/physics/utils.py +++ /dev/null @@ -1,77 +0,0 @@ -# Copyright 2024 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. - -""" -This module contains helper functions for physics. -""" -import numpy as np - -from mrmustard import math - - -def real_gaussian_integral( - Abc: tuple, - idx: tuple[int, ...], -): - r"""Computes the Gaussian integral of the exponential of a real quadratic form. - The integral is defined as (note that in general we integrate over a subset of m dimensions): - - .. :math:: - \int_{R^m} F(x) dx, - - where - - :math:`F(x) = \textrm{exp}(0.5 x^T A x + b^T x)` - - Here, ``x`` is an ``n``-dim real vector, ``A`` is an ``n x n`` real matrix, - ``b`` is an ``n``-dim real vector, ``c`` is a real scalar. The integral indices - are specified by ``idx``. - - Arguments: - A,b,c: the ``(A,b,c)`` triple - idx: the tuple of indices of the x variables - - Returns: - The ``(A,b,c)`` triple of the result of the integral. - """ - A, b, c = Abc - - n = len(idx) - - if not idx: - return A, b, c - not_idx = tuple(i for i in range(A.shape[-1]) if i not in idx) - - M = math.gather(math.gather(A, idx, axis=-1), idx, axis=-2) - bM = math.gather(b, idx, axis=-1) - - if math.asnumpy(not_idx).shape != (0,): - D = math.gather(math.gather(A, idx, axis=-1), not_idx, axis=-2) - R = math.gather(math.gather(A, not_idx, axis=-1), not_idx, axis=-2) - bR = math.gather(b, not_idx, axis=-1) - T = math.transpose - L = T(math.solve(T(M), T(D))) - A_post = R - math.matmul(L, T(D)) - b_post = bR - math.matvec(L, bM) - else: - A_post = math.astensor([]) - b_post = math.astensor([]) - - c_post = ( - c - * math.sqrt((2 * np.pi) ** n / math.det(M)) - * math.exp(-0.5 * math.sum(bM * math.solve(M, bM))) - ) - - return A_post, b_post, c_post diff --git a/tests/test_physics/test_bargmann.py b/tests/test_physics/test_bargmann.py index 645871fc8..229536b96 100644 --- a/tests/test_physics/test_bargmann.py +++ b/tests/test_physics/test_bargmann.py @@ -1,35 +1,12 @@ import numpy as np -import pytest -from mrmustard import math from mrmustard.lab import Attenuator, Dgate, Gaussian, Ggate from mrmustard.physics.bargmann import ( - complex_gaussian_integral, - contract_two_Abc, - join_Abc, wigner_to_bargmann_Choi, wigner_to_bargmann_psi, wigner_to_bargmann_rho, wigner_to_bargmann_U, - reorder_abc, ) -from mrmustard.physics import triples - - -def test_reorder_abc(): - """Test that the reorder_abc function works correctly""" - A = np.array([[1, 2], [2, 3]]) - b = np.array([4, 5]) - c = np.array(6) - same = reorder_abc((A, b, c), (0, 1)) - assert all(np.allclose(x, y) for x, y in zip(same, (A, b, c))) - flipped = reorder_abc((A, b, c), (1, 0)) - assert all(np.allclose(x, y) for x, y in zip(flipped, (A[[1, 0], :][:, [1, 0]], b[[1, 0]], c))) - c = np.array([[6, 7], [8, 9]]) - flipped = reorder_abc((A, b, c), (1, 0)) # test transposition of c - assert all( - np.allclose(x, y) for x, y in zip(flipped, (A[[1, 0], :][:, [1, 0]], b[[1, 0]], c.T)) - ) def test_wigner_to_bargmann_psi(): @@ -74,79 +51,3 @@ def test_bargmann_numpy_transformation(): """Tests that the numpy option of the bargmann method of State works correctly""" transformation = Ggate(1) assert all(isinstance(thing, np.ndarray) for thing in transformation.bargmann(numpy=True)) - - -def test_join_Abc(): - """Tests the ``join_Abc`` method.""" - A1, b1, c1 = triples.vacuum_state_Abc(2) - A2, b2, c2 = triples.displacement_gate_Abc(x=[0.1, 0.2], y=0.3) - - joined_Abc = join_Abc((A1, b1, c1), (A2, b2, c2)) - assert np.allclose(joined_Abc[0], math.block_diag(A1, A2)) - assert np.allclose(joined_Abc[1], math.concat([b1, b2], axis=-1)) - assert np.allclose(joined_Abc[2], math.outer(c1, c2)) - - A12 = math.block_diag(A1, A2) - b12 = math.concat([b1, b2], axis=-1) - c12 = math.outer(c1, c2) - return A12, b12, c12 - - -def test_complex_gaussian_integral(): - """Tests the ``complex_gaussian_integral`` method.""" - A1, b1, c1 = triples.vacuum_state_Abc(2) - A2, b2, c2 = triples.displacement_gate_Abc(x=[0.1, 0.2], y=0.3) - A3, b3, c3 = triples.displaced_squeezed_vacuum_state_Abc(x=[0.1, 0.2], y=0.3) - - joined_Abc = join_Abc((A1, b1, c1), (A2, b2, c2)) - - res1 = complex_gaussian_integral(joined_Abc, [], []) - assert np.allclose(res1[0], joined_Abc[0]) - assert np.allclose(res1[1], joined_Abc[1]) - assert np.allclose(res1[2], joined_Abc[2]) - - res2 = complex_gaussian_integral(joined_Abc, [0, 1], [4, 5]) - assert np.allclose(res2[0], A3) - assert np.allclose(res2[1], b3) - assert np.allclose(res2[2], c3) - - res3 = complex_gaussian_integral(join_Abc((A1, b1, c1), (A1, b1, c1)), [0, 1], [2, 3]) - assert np.allclose(res3[0], 0) - assert np.allclose(res3[1], 0) - assert np.allclose(res3[2], 1) - - -def test_complex_gaussian_integral_error(): - """Tests the error of the ``complex_gaussian_integral`` method.""" - A1, b1, c1 = triples.vacuum_state_Abc(2) - A2, b2, c2 = triples.displacement_gate_Abc(x=[0.1, 0.2], y=0.3) - - with pytest.raises(ValueError): - complex_gaussian_integral( - join_Abc((A1, b1, c1), (A2, b2, c2)), - [0, 1], - [ - 4, - ], - ) - - -def test_contract_two_Abc(): - """Tests the error of the ``contract_two_Abc`` method.""" - A1, b1, c1 = triples.vacuum_state_Abc(2) - A2, b2, c2 = triples.displacement_gate_Abc(x=[0.1, 0.2], y=0.3) - - res1 = contract_two_Abc((A1, b1, c1), (A2, b2, c2), [], []) - assert np.allclose(res1[0], math.block_diag(A1, A2)) - assert np.allclose(res1[1], [0, 0, 0.1 + 0.3j, 0.2 + 0.3j, -0.1 + 0.3j, -0.2 + 0.3j]) - assert np.allclose(res1[2], c1 * c2) - - res2 = contract_two_Abc((A1, b1, c1), (A2, b2, c2), [0, 1], [2, 3]) - assert np.allclose(res2[0], math.zeros((2, 2))) - assert np.allclose(res2[1], [0.1 + 0.3j, 0.2 + 0.3j]) - assert np.allclose(res2[2], c1 * c2) - - res3 = contract_two_Abc((A1, b1, c1), (A2, b2, c2), [0, 1], [0, 1]) - assert np.allclose(res3[0], math.zeros((2, 2))) - assert np.allclose(res3[1], [-0.1 + 0.3j, -0.2 + 0.3j]) - assert np.allclose(res3[2], c1 * c2) diff --git a/tests/test_physics/test_utils.py b/tests/test_physics/test_utils.py deleted file mode 100644 index af4544e86..000000000 --- a/tests/test_physics/test_utils.py +++ /dev/null @@ -1,69 +0,0 @@ -# Copyright 2024 Xanadu Quantum Technologies Inc. - -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at - -# http://www.apache.org/licenses/LICENSE-2.0 - -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. - -"""Tests for the helper functions in physics.""" - -import numpy as np -from mrmustard import math -from mrmustard.physics.utils import real_gaussian_integral - - -def test_real_gaussian_integral(): - """Tests the ``real_gaussian_integral`` method with a hard-coded A matric from a Gaussian(3) state.""" - A = math.astensor( - np.array( - [ - [0.35307718 - 0.09738001j, -0.01297994 + 0.26050244j, 0.05349344 - 0.13728068j], - [-0.01297994 + 0.26050244j, 0.05696707 - 0.2351408j, 0.18954838 - 0.42959383j], - [0.05349344 - 0.13728068j, 0.18954838 - 0.42959383j, -0.16931712 - 0.09205837j], - ] - ) - ) - b = math.astensor(np.arange(3) + 0j) - c = 1.0 + 0j - res = real_gaussian_integral((A, b, c), idx=[0, 1]) - assert np.allclose(res[0], A[2, 2] - A[2:, :2] @ math.inv(A[:2, :2]) @ A[:2, 2:]) - assert np.allclose(res[1], b[2] - math.sum(A[2:, :2] * math.matvec(math.inv(A[:2, :2]), b[:2]))) - assert np.allclose( - res[2], - c - * (2 * np.pi) - / math.sqrt(math.det(A[:2, :2])) - * math.exp(-0.5 * math.sum(b[:2] * math.matvec(math.inv(A[:2, :2]), b[:2]))), - ) - res2 = real_gaussian_integral((A, b, c), idx=[]) - assert np.allclose(res2[0], A) - assert np.allclose(res2[1], b) - assert np.allclose(res2[2], c) - - A2 = math.astensor( - np.array( - [ - [0.35307718 - 0.09738001j, -0.01297994 + 0.26050244j], - [-0.01297994 + 0.26050244j, 0.05696707 - 0.2351408j], - ] - ) - ) - b2 = math.astensor(np.arange(2) + 0j) - c2 = 1.0 + 0j - res3 = real_gaussian_integral((A2, b2, c2), idx=[0, 1]) - assert np.allclose(res3[0], math.astensor([])) - assert np.allclose(res3[1], math.astensor([])) - assert np.allclose( - res3[2], - c2 - * (2 * np.pi) - / math.sqrt(math.det(A2[:2, :2])) - * math.exp(-0.5 * math.sum(b2[:2] * math.matvec(math.inv(A2[:2, :2]), b2[:2]))), - ) From af7ba6c36951e7cc13c803e502316ba7d80aa43c Mon Sep 17 00:00:00 2001 From: Yuan YAO Date: Thu, 11 Apr 2024 20:24:13 +0000 Subject: [PATCH 12/30] fix code factor issue --- mrmustard/physics/gaussian_integrals.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/mrmustard/physics/gaussian_integrals.py b/mrmustard/physics/gaussian_integrals.py index 2b274b2b8..e7ffd68cb 100644 --- a/mrmustard/physics/gaussian_integrals.py +++ b/mrmustard/physics/gaussian_integrals.py @@ -215,4 +215,4 @@ def contract_two_Abc( Abc = join_Abc(Abc1, Abc2) return complex_gaussian_integral( Abc, idx1, tuple(n + Abc1[0].shape[-1] for n in idx2), measure=-1.0 - ) \ No newline at end of file + ) From 554bedbe7e0a957d5b2fdba17475f598fe6c81ed Mon Sep 17 00:00:00 2001 From: Yuan YAO Date: Thu, 11 Apr 2024 20:26:10 +0000 Subject: [PATCH 13/30] blacked --- mrmustard/physics/gaussian_integrals.py | 1 + 1 file changed, 1 insertion(+) diff --git a/mrmustard/physics/gaussian_integrals.py b/mrmustard/physics/gaussian_integrals.py index e7ffd68cb..866eb9375 100644 --- a/mrmustard/physics/gaussian_integrals.py +++ b/mrmustard/physics/gaussian_integrals.py @@ -79,6 +79,7 @@ def real_gaussian_integral( return A_post, b_post, c_post + def complex_gaussian_integral( Abc: tuple, idx_z: tuple[int, ...], idx_zconj: tuple[int, ...], measure: float = -1 ): From 3e23674f9e4995ac07715a349a10534bc6d8c4de Mon Sep 17 00:00:00 2001 From: Yuan YAO <16817699+sylviemonet@users.noreply.github.com> Date: Thu, 11 Apr 2024 16:27:20 -0400 Subject: [PATCH 14/30] Update mrmustard/physics/gaussian_integrals.py --- mrmustard/physics/gaussian_integrals.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/mrmustard/physics/gaussian_integrals.py b/mrmustard/physics/gaussian_integrals.py index 866eb9375..1241b5dec 100644 --- a/mrmustard/physics/gaussian_integrals.py +++ b/mrmustard/physics/gaussian_integrals.py @@ -13,7 +13,7 @@ # limitations under the License. """ -This module contains real and comple gaussian integral functions and related helper functions. +This module contains gaussian integral functions and related helper functions. """ from typing import Sequence, Tuple From 79e4ae547ab84d977c62e05b72a8f42db9b5b18e Mon Sep 17 00:00:00 2001 From: Yuan YAO Date: Fri, 12 Apr 2024 18:34:59 +0000 Subject: [PATCH 15/30] add tag for contract_two_abc --- mrmustard/physics/gaussian_integrals.py | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/mrmustard/physics/gaussian_integrals.py b/mrmustard/physics/gaussian_integrals.py index 1241b5dec..d2f0f36b6 100644 --- a/mrmustard/physics/gaussian_integrals.py +++ b/mrmustard/physics/gaussian_integrals.py @@ -195,6 +195,7 @@ def contract_two_Abc( Abc2: Tuple[ComplexMatrix, ComplexVector, complex], idx1: Sequence[int], idx2: Sequence[int], + real: bool, ): r""" Returns the contraction of two ``(A,b,c)`` triples with given indices. @@ -209,11 +210,15 @@ def contract_two_Abc( Abc2: the second ``(A,b,c)`` triple idx1: the indices of the first ``(A,b,c)`` triple to contract idx2: the indices of the second ``(A,b,c)`` triple to contract + real: the tag to indicate whether the contraction using real or complex gaussian integral Returns: The contracted ``(A,b,c)`` triple """ Abc = join_Abc(Abc1, Abc2) - return complex_gaussian_integral( - Abc, idx1, tuple(n + Abc1[0].shape[-1] for n in idx2), measure=-1.0 - ) + if real: + return real_gaussian_integral(Abc, idx1 + list(n + Abc1[0].shape[-1] for n in idx2)) + else: + return complex_gaussian_integral( + Abc, idx1, tuple(n + Abc1[0].shape[-1] for n in idx2), measure=-1.0 + ) From 56a08f6f69e9d7dfa4912c76228d6630940d544a Mon Sep 17 00:00:00 2001 From: Yuan YAO Date: Fri, 12 Apr 2024 18:41:12 +0000 Subject: [PATCH 16/30] add tests --- tests/test_physics/test_gaussian_integrals.py | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/tests/test_physics/test_gaussian_integrals.py b/tests/test_physics/test_gaussian_integrals.py index c46f29b4d..cdac627cf 100644 --- a/tests/test_physics/test_gaussian_integrals.py +++ b/tests/test_physics/test_gaussian_integrals.py @@ -137,21 +137,26 @@ def test_contract_two_Abc(): A1, b1, c1 = triples.vacuum_state_Abc(2) A2, b2, c2 = triples.displacement_gate_Abc(x=[0.1, 0.2], y=0.3) - res1 = contract_two_Abc((A1, b1, c1), (A2, b2, c2), [], []) + res1 = contract_two_Abc((A1, b1, c1), (A2, b2, c2), [], [], real=False) assert np.allclose(res1[0], math.block_diag(A1, A2)) assert np.allclose(res1[1], [0, 0, 0.1 + 0.3j, 0.2 + 0.3j, -0.1 + 0.3j, -0.2 + 0.3j]) assert np.allclose(res1[2], c1 * c2) - res2 = contract_two_Abc((A1, b1, c1), (A2, b2, c2), [0, 1], [2, 3]) + res2 = contract_two_Abc((A1, b1, c1), (A2, b2, c2), [0, 1], [2, 3], real=False) assert np.allclose(res2[0], math.zeros((2, 2))) assert np.allclose(res2[1], [0.1 + 0.3j, 0.2 + 0.3j]) assert np.allclose(res2[2], c1 * c2) - res3 = contract_two_Abc((A1, b1, c1), (A2, b2, c2), [0, 1], [0, 1]) + res3 = contract_two_Abc((A1, b1, c1), (A2, b2, c2), [0, 1], [0, 1], real=False) assert np.allclose(res3[0], math.zeros((2, 2))) assert np.allclose(res3[1], [-0.1 + 0.3j, -0.2 + 0.3j]) assert np.allclose(res3[2], c1 * c2) + res4 = contract_two_Abc((A1, b1, c1), (A2, b2, c2), [], [], real=True) + assert np.allclose(res4[0], res1[0]) + assert np.allclose(res4[1], res1[1]) + assert np.allclose(res4[2], res1[2]) + def test_reorder_abc(): """Test that the reorder_abc function works correctly""" From 8dbd06308566b63e2170430a21fcf878880c98b6 Mon Sep 17 00:00:00 2001 From: Yuan YAO Date: Fri, 12 Apr 2024 18:47:39 +0000 Subject: [PATCH 17/30] fix the importation errors --- mrmustard/physics/representations.py | 12 +++--------- tests/test_physics/test_representations.py | 4 ++-- 2 files changed, 5 insertions(+), 11 deletions(-) diff --git a/mrmustard/physics/representations.py b/mrmustard/physics/representations.py index 5d99ac36a..602fb7f10 100644 --- a/mrmustard/physics/representations.py +++ b/mrmustard/physics/representations.py @@ -27,6 +27,7 @@ from mrmustard import math, settings from mrmustard.physics import bargmann +from mrmustard.physics.gaussian_integrals import contract_two_Abc, reorder_abc from mrmustard.physics.ansatze import Ansatz, PolyExpAnsatz, ArrayAnsatz from mrmustard.utils.typing import ( Batch, @@ -316,7 +317,7 @@ def reorder(self, order: tuple[int, ...] | list[int]) -> Bargmann: Returns: The reordered Bargmann object. """ - A, b, c = bargmann.reorder_abc((self.A, self.b, self.c), order) + A, b, c = reorder_abc((self.A, self.b, self.c), order) return self.__class__(A, b, c) def plot( @@ -440,14 +441,7 @@ def __matmul__(self, other: Union[Bargmann, Fock]) -> Union[Bargmann, Fock]: Abc = [] for A1, b1, c1 in zip(self.A, self.b, self.c): for A2, b2, c2 in zip(other.A, other.b, other.c): - Abc.append( - bargmann.contract_two_Abc( - (A1, b1, c1), - (A2, b2, c2), - idx_s, - idx_o, - ) - ) + Abc.append(contract_two_Abc((A1, b1, c1), (A2, b2, c2), idx_s, idx_o, real=False)) A, b, c = zip(*Abc) return Bargmann(A, b, c) diff --git a/tests/test_physics/test_representations.py b/tests/test_physics/test_representations.py index bb56b8840..cd9ce4b9f 100644 --- a/tests/test_physics/test_representations.py +++ b/tests/test_physics/test_representations.py @@ -20,7 +20,7 @@ from mrmustard import math, settings from mrmustard.physics.converters import to_fock from mrmustard.physics.triples import displacement_gate_Abc, attenuator_Abc -from mrmustard.physics.bargmann import contract_two_Abc, complex_gaussian_integral +from mrmustard.physics.gaussian_integrals import contract_two_Abc, complex_gaussian_integral from mrmustard.physics.representations import Bargmann, Fock from ..random import Abc_triple @@ -181,7 +181,7 @@ def test_matmul_barg_barg(self): triple2 = Abc_triple(3) res1 = Bargmann(*triple1) @ Bargmann(*triple2) - exp1 = contract_two_Abc(triple1, triple2, [], []) + exp1 = contract_two_Abc(triple1, triple2, [], [], real=False) assert np.allclose(res1.A, exp1[0]) assert np.allclose(res1.b, exp1[1]) assert np.allclose(res1.c, exp1[2]) From ccfebe0306901d823e1dc9de05f191f95db8d569 Mon Sep 17 00:00:00 2001 From: Yuan YAO <16817699+sylviemonet@users.noreply.github.com> Date: Fri, 12 Apr 2024 14:48:45 -0400 Subject: [PATCH 18/30] Update mrmustard/physics/gaussian_integrals.py Co-authored-by: SamFerracin --- mrmustard/physics/gaussian_integrals.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/mrmustard/physics/gaussian_integrals.py b/mrmustard/physics/gaussian_integrals.py index d2f0f36b6..b7ae53676 100644 --- a/mrmustard/physics/gaussian_integrals.py +++ b/mrmustard/physics/gaussian_integrals.py @@ -42,7 +42,7 @@ def real_gaussian_integral( are specified by ``idx``. Arguments: - A,b,c: the ``(A,b,c)`` triple + Abc: the ``(A,b,c)`` triple idx: the tuple of indices of the x variables Returns: From fca8dbd77474d319e7a21f5bbfb3e763c0be0814 Mon Sep 17 00:00:00 2001 From: Yuan YAO <16817699+sylviemonet@users.noreply.github.com> Date: Fri, 12 Apr 2024 14:48:53 -0400 Subject: [PATCH 19/30] Update mrmustard/physics/gaussian_integrals.py Co-authored-by: SamFerracin --- mrmustard/physics/gaussian_integrals.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/mrmustard/physics/gaussian_integrals.py b/mrmustard/physics/gaussian_integrals.py index b7ae53676..426fc7a42 100644 --- a/mrmustard/physics/gaussian_integrals.py +++ b/mrmustard/physics/gaussian_integrals.py @@ -28,7 +28,7 @@ def real_gaussian_integral( idx: tuple[int, ...], ): r"""Computes the Gaussian integral of the exponential of a real quadratic form. - The integral is defined as (note that in general we integrate over a subset of m dimensions): + The integral is defined as (note that in general we integrate over a subset of `m` dimensions): .. :math:: \int_{R^m} F(x) dx, From 89073aa14f2c1901cb7e0b7927f978c5367ffff7 Mon Sep 17 00:00:00 2001 From: Yuan YAO Date: Fri, 12 Apr 2024 18:50:10 +0000 Subject: [PATCH 20/30] fix small bugs --- mrmustard/physics/gaussian_integrals.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/mrmustard/physics/gaussian_integrals.py b/mrmustard/physics/gaussian_integrals.py index 426fc7a42..ac10fe98b 100644 --- a/mrmustard/physics/gaussian_integrals.py +++ b/mrmustard/physics/gaussian_integrals.py @@ -50,8 +50,6 @@ def real_gaussian_integral( """ A, b, c = Abc - n = len(idx) - if not idx: return A, b, c not_idx = tuple(i for i in range(A.shape[-1]) if i not in idx) @@ -73,7 +71,7 @@ def real_gaussian_integral( c_post = ( c - * math.sqrt((2 * np.pi) ** n / math.det(M)) + * math.sqrt((2 * np.pi) ** len(idx) / math.det(M)) * math.exp(-0.5 * math.sum(bM * math.solve(M, bM))) ) From bf2edcb6a341107332c53619e89b13d14ec96034 Mon Sep 17 00:00:00 2001 From: Yuan YAO Date: Fri, 12 Apr 2024 18:51:12 +0000 Subject: [PATCH 21/30] add type hints --- mrmustard/physics/gaussian_integrals.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/mrmustard/physics/gaussian_integrals.py b/mrmustard/physics/gaussian_integrals.py index ac10fe98b..bc2528384 100644 --- a/mrmustard/physics/gaussian_integrals.py +++ b/mrmustard/physics/gaussian_integrals.py @@ -150,7 +150,10 @@ def complex_gaussian_integral( return A_post, b_post, c_post -def join_Abc(Abc1, Abc2): +def join_Abc( + Abc1: Tuple[ComplexMatrix, ComplexVector, complex], + Abc2: Tuple[ComplexMatrix, ComplexVector, complex], +): r"""Joins two ``(A,b,c)`` triples into a single ``(A,b,c)`` triple by block addition of the ``A`` matrices and concatenating the ``b`` vectors. From 26cc481aa9790b8d93ce969ea24184ad54eb5462 Mon Sep 17 00:00:00 2001 From: Yuan YAO Date: Fri, 12 Apr 2024 18:52:27 +0000 Subject: [PATCH 22/30] fix the docs --- mrmustard/physics/gaussian_integrals.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/mrmustard/physics/gaussian_integrals.py b/mrmustard/physics/gaussian_integrals.py index bc2528384..1b61b7d1a 100644 --- a/mrmustard/physics/gaussian_integrals.py +++ b/mrmustard/physics/gaussian_integrals.py @@ -201,7 +201,7 @@ def contract_two_Abc( r""" Returns the contraction of two ``(A,b,c)`` triples with given indices. - Note that the indices must be a complex variable pairs with each other to make this contraction meaningful. Please make sure + Note that the indices of the complex_gaussian_integral (real = False) must be a complex variable pairs with each other to make this contraction meaningful. Please make sure the corresponding complex variable with respect to your Abc triples. For examples, if the indices of Abc1 denotes the variables ``(\alpha, \beta)``, the indices of Abc2 denotes the variables ``(\alpha^*,\gamma)``, the contraction only works with ``idx1 = [0], idx2 = [0]``. From bd69cbe46c6eff7583f1d0fb4e0613f8144e4387 Mon Sep 17 00:00:00 2001 From: Yuan YAO Date: Fri, 12 Apr 2024 18:55:15 +0000 Subject: [PATCH 23/30] fix the docs bug --- mrmustard/physics/representations.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/mrmustard/physics/representations.py b/mrmustard/physics/representations.py index 602fb7f10..e9668c74b 100644 --- a/mrmustard/physics/representations.py +++ b/mrmustard/physics/representations.py @@ -27,7 +27,7 @@ from mrmustard import math, settings from mrmustard.physics import bargmann -from mrmustard.physics.gaussian_integrals import contract_two_Abc, reorder_abc +from mrmustard.physics.gaussian_integrals import contract_two_Abc, reorder_abc, complex_gaussian_integral from mrmustard.physics.ansatze import Ansatz, PolyExpAnsatz, ArrayAnsatz from mrmustard.utils.typing import ( Batch, @@ -290,7 +290,7 @@ def trace(self, idx_z: tuple[int, ...], idx_zconj: tuple[int, ...]) -> Bargmann: ) A, b, c = [], [], [] for Abci in zip(self.A, self.b, self.c): - Aij, bij, cij = bargmann.complex_gaussian_integral(Abci, idx_z, idx_zconj, measure=-1.0) + Aij, bij, cij = complex_gaussian_integral(Abci, idx_z, idx_zconj, measure=-1.0) A.append(Aij) b.append(bij) c.append(cij) From cf383ae001dafc64fd76898843c9e572c583fc46 Mon Sep 17 00:00:00 2001 From: Yuan YAO Date: Fri, 12 Apr 2024 18:57:00 +0000 Subject: [PATCH 24/30] blacked --- mrmustard/physics/representations.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/mrmustard/physics/representations.py b/mrmustard/physics/representations.py index e9668c74b..3fe34a6eb 100644 --- a/mrmustard/physics/representations.py +++ b/mrmustard/physics/representations.py @@ -26,8 +26,11 @@ import numpy as np from mrmustard import math, settings -from mrmustard.physics import bargmann -from mrmustard.physics.gaussian_integrals import contract_two_Abc, reorder_abc, complex_gaussian_integral +from mrmustard.physics.gaussian_integrals import ( + contract_two_Abc, + reorder_abc, + complex_gaussian_integral, +) from mrmustard.physics.ansatze import Ansatz, PolyExpAnsatz, ArrayAnsatz from mrmustard.utils.typing import ( Batch, From 73b09e005eda9adf78a84f52d14fa3655299539b Mon Sep 17 00:00:00 2001 From: Yuan YAO Date: Fri, 12 Apr 2024 21:05:30 +0000 Subject: [PATCH 25/30] rerun --- mrmustard/physics/gaussian_integrals.py | 1 - 1 file changed, 1 deletion(-) diff --git a/mrmustard/physics/gaussian_integrals.py b/mrmustard/physics/gaussian_integrals.py index 1b61b7d1a..13732f065 100644 --- a/mrmustard/physics/gaussian_integrals.py +++ b/mrmustard/physics/gaussian_integrals.py @@ -16,7 +16,6 @@ This module contains gaussian integral functions and related helper functions. """ from typing import Sequence, Tuple - import numpy as np from mrmustard import math From 4b888c8ce84c38d2a62ff8398eb051fdce495f6b Mon Sep 17 00:00:00 2001 From: Yuan YAO Date: Mon, 15 Apr 2024 14:31:29 +0000 Subject: [PATCH 26/30] add the comments from Sam --- doc/code/physics/bargmann.rst | 1 + doc/code/physics/utils/gaussian_integral.rst | 8 ++++++++ mrmustard/physics/gaussian_integrals.py | 8 ++++---- 3 files changed, 13 insertions(+), 4 deletions(-) create mode 100644 doc/code/physics/utils/gaussian_integral.rst diff --git a/doc/code/physics/bargmann.rst b/doc/code/physics/bargmann.rst index 2562ef566..28499bcf1 100644 --- a/doc/code/physics/bargmann.rst +++ b/doc/code/physics/bargmann.rst @@ -6,3 +6,4 @@ Calculations on Bargmann objects utils/triples utils/bargmann_calculations + utils/gaussian_integral diff --git a/doc/code/physics/utils/gaussian_integral.rst b/doc/code/physics/utils/gaussian_integral.rst new file mode 100644 index 000000000..4071d4316 --- /dev/null +++ b/doc/code/physics/utils/gaussian_integral.rst @@ -0,0 +1,8 @@ +Calculations on Gaussian Integrals +================================ + +.. currentmodule:: mrmustard.physics.bargmann + +.. automodapi:: mrmustard.physics.bargmann + :no-heading: + :include-all-objects: \ No newline at end of file diff --git a/mrmustard/physics/gaussian_integrals.py b/mrmustard/physics/gaussian_integrals.py index 13732f065..7a276abb8 100644 --- a/mrmustard/physics/gaussian_integrals.py +++ b/mrmustard/physics/gaussian_integrals.py @@ -218,7 +218,7 @@ def contract_two_Abc( Abc = join_Abc(Abc1, Abc2) if real: return real_gaussian_integral(Abc, idx1 + list(n + Abc1[0].shape[-1] for n in idx2)) - else: - return complex_gaussian_integral( - Abc, idx1, tuple(n + Abc1[0].shape[-1] for n in idx2), measure=-1.0 - ) + + return complex_gaussian_integral( + Abc, idx1, tuple(n + Abc1[0].shape[-1] for n in idx2), measure=-1.0 + ) From 50ade055d7e42ba7da269554c378ea74ea5fbd46 Mon Sep 17 00:00:00 2001 From: SamFerracin Date: Mon, 15 Apr 2024 11:35:51 -0400 Subject: [PATCH 27/30] empty From 7c9a9d75acdc539edde8e8420b61adb01d0c164a Mon Sep 17 00:00:00 2001 From: Yuan YAO Date: Mon, 15 Apr 2024 19:04:08 +0000 Subject: [PATCH 28/30] change back contract_two_Abc --- mrmustard/physics/gaussian_integrals.py | 7 ------- mrmustard/physics/representations.py | 2 +- tests/test_physics/test_gaussian_integrals.py | 8 ++++---- tests/test_physics/test_representations.py | 2 +- 4 files changed, 6 insertions(+), 13 deletions(-) diff --git a/mrmustard/physics/gaussian_integrals.py b/mrmustard/physics/gaussian_integrals.py index 7a276abb8..e1794757b 100644 --- a/mrmustard/physics/gaussian_integrals.py +++ b/mrmustard/physics/gaussian_integrals.py @@ -195,13 +195,10 @@ def contract_two_Abc( Abc2: Tuple[ComplexMatrix, ComplexVector, complex], idx1: Sequence[int], idx2: Sequence[int], - real: bool, ): r""" Returns the contraction of two ``(A,b,c)`` triples with given indices. - Note that the indices of the complex_gaussian_integral (real = False) must be a complex variable pairs with each other to make this contraction meaningful. Please make sure - the corresponding complex variable with respect to your Abc triples. For examples, if the indices of Abc1 denotes the variables ``(\alpha, \beta)``, the indices of Abc2 denotes the variables ``(\alpha^*,\gamma)``, the contraction only works with ``idx1 = [0], idx2 = [0]``. @@ -210,15 +207,11 @@ def contract_two_Abc( Abc2: the second ``(A,b,c)`` triple idx1: the indices of the first ``(A,b,c)`` triple to contract idx2: the indices of the second ``(A,b,c)`` triple to contract - real: the tag to indicate whether the contraction using real or complex gaussian integral Returns: The contracted ``(A,b,c)`` triple """ Abc = join_Abc(Abc1, Abc2) - if real: - return real_gaussian_integral(Abc, idx1 + list(n + Abc1[0].shape[-1] for n in idx2)) - return complex_gaussian_integral( Abc, idx1, tuple(n + Abc1[0].shape[-1] for n in idx2), measure=-1.0 ) diff --git a/mrmustard/physics/representations.py b/mrmustard/physics/representations.py index 3fe34a6eb..94c899d37 100644 --- a/mrmustard/physics/representations.py +++ b/mrmustard/physics/representations.py @@ -444,7 +444,7 @@ def __matmul__(self, other: Union[Bargmann, Fock]) -> Union[Bargmann, Fock]: Abc = [] for A1, b1, c1 in zip(self.A, self.b, self.c): for A2, b2, c2 in zip(other.A, other.b, other.c): - Abc.append(contract_two_Abc((A1, b1, c1), (A2, b2, c2), idx_s, idx_o, real=False)) + Abc.append(contract_two_Abc((A1, b1, c1), (A2, b2, c2), idx_s, idx_o)) A, b, c = zip(*Abc) return Bargmann(A, b, c) diff --git a/tests/test_physics/test_gaussian_integrals.py b/tests/test_physics/test_gaussian_integrals.py index cdac627cf..3ad9906c3 100644 --- a/tests/test_physics/test_gaussian_integrals.py +++ b/tests/test_physics/test_gaussian_integrals.py @@ -137,22 +137,22 @@ def test_contract_two_Abc(): A1, b1, c1 = triples.vacuum_state_Abc(2) A2, b2, c2 = triples.displacement_gate_Abc(x=[0.1, 0.2], y=0.3) - res1 = contract_two_Abc((A1, b1, c1), (A2, b2, c2), [], [], real=False) + res1 = contract_two_Abc((A1, b1, c1), (A2, b2, c2), [], []) assert np.allclose(res1[0], math.block_diag(A1, A2)) assert np.allclose(res1[1], [0, 0, 0.1 + 0.3j, 0.2 + 0.3j, -0.1 + 0.3j, -0.2 + 0.3j]) assert np.allclose(res1[2], c1 * c2) - res2 = contract_two_Abc((A1, b1, c1), (A2, b2, c2), [0, 1], [2, 3], real=False) + res2 = contract_two_Abc((A1, b1, c1), (A2, b2, c2), [0, 1], [2, 3]) assert np.allclose(res2[0], math.zeros((2, 2))) assert np.allclose(res2[1], [0.1 + 0.3j, 0.2 + 0.3j]) assert np.allclose(res2[2], c1 * c2) - res3 = contract_two_Abc((A1, b1, c1), (A2, b2, c2), [0, 1], [0, 1], real=False) + res3 = contract_two_Abc((A1, b1, c1), (A2, b2, c2), [0, 1], [0, 1]) assert np.allclose(res3[0], math.zeros((2, 2))) assert np.allclose(res3[1], [-0.1 + 0.3j, -0.2 + 0.3j]) assert np.allclose(res3[2], c1 * c2) - res4 = contract_two_Abc((A1, b1, c1), (A2, b2, c2), [], [], real=True) + res4 = contract_two_Abc((A1, b1, c1), (A2, b2, c2), [], []) assert np.allclose(res4[0], res1[0]) assert np.allclose(res4[1], res1[1]) assert np.allclose(res4[2], res1[2]) diff --git a/tests/test_physics/test_representations.py b/tests/test_physics/test_representations.py index cd9ce4b9f..1ffba8d7e 100644 --- a/tests/test_physics/test_representations.py +++ b/tests/test_physics/test_representations.py @@ -181,7 +181,7 @@ def test_matmul_barg_barg(self): triple2 = Abc_triple(3) res1 = Bargmann(*triple1) @ Bargmann(*triple2) - exp1 = contract_two_Abc(triple1, triple2, [], [], real=False) + exp1 = contract_two_Abc(triple1, triple2, [], []) assert np.allclose(res1.A, exp1[0]) assert np.allclose(res1.b, exp1[1]) assert np.allclose(res1.c, exp1[2]) From 0260104be0bb71bf071fee31676375e2f5e34a02 Mon Sep 17 00:00:00 2001 From: Yuan YAO Date: Mon, 15 Apr 2024 19:05:10 +0000 Subject: [PATCH 29/30] change the text --- mrmustard/physics/gaussian_integrals.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/mrmustard/physics/gaussian_integrals.py b/mrmustard/physics/gaussian_integrals.py index e1794757b..a0baca448 100644 --- a/mrmustard/physics/gaussian_integrals.py +++ b/mrmustard/physics/gaussian_integrals.py @@ -198,7 +198,11 @@ def contract_two_Abc( ): r""" Returns the contraction of two ``(A,b,c)`` triples with given indices. - + + Note that the indices must be a complex variable pairs with each other to make this contraction meaningful. Please make sure + the corresponding complex variable with respect to your Abc triples. + For examples, if the indices of Abc1 denotes the variables ``(\alpha, \beta)``, the indices of Abc2 denotes the variables + ``(\alpha^*,\gamma)``, the contraction only works with ``idx1 = [0], idx2 = [0]``. For examples, if the indices of Abc1 denotes the variables ``(\alpha, \beta)``, the indices of Abc2 denotes the variables ``(\alpha^*,\gamma)``, the contraction only works with ``idx1 = [0], idx2 = [0]``. From c095dfe85e398b7265367df7d20126f53aec6534 Mon Sep 17 00:00:00 2001 From: Yuan YAO Date: Mon, 15 Apr 2024 19:05:23 +0000 Subject: [PATCH 30/30] fix docs --- mrmustard/physics/gaussian_integrals.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/mrmustard/physics/gaussian_integrals.py b/mrmustard/physics/gaussian_integrals.py index a0baca448..e16fa657d 100644 --- a/mrmustard/physics/gaussian_integrals.py +++ b/mrmustard/physics/gaussian_integrals.py @@ -198,13 +198,11 @@ def contract_two_Abc( ): r""" Returns the contraction of two ``(A,b,c)`` triples with given indices. - + Note that the indices must be a complex variable pairs with each other to make this contraction meaningful. Please make sure the corresponding complex variable with respect to your Abc triples. For examples, if the indices of Abc1 denotes the variables ``(\alpha, \beta)``, the indices of Abc2 denotes the variables ``(\alpha^*,\gamma)``, the contraction only works with ``idx1 = [0], idx2 = [0]``. - For examples, if the indices of Abc1 denotes the variables ``(\alpha, \beta)``, the indices of Abc2 denotes the variables - ``(\alpha^*,\gamma)``, the contraction only works with ``idx1 = [0], idx2 = [0]``. Arguments: Abc1: the first ``(A,b,c)`` triple