feat(gaussian): Wigner function

The calculation of the Wigner function has been implemented under
`GaussianState`. Since the implementation got too large, the core of the
calcutation got moved to `functions.py` in the same folder.

The `gaussian_wigner_function` accepts not only a single canonical
coordinate vector, but a list of vectors.

**Note**

The `.flake8` file got extended with exclusions, because `flake8` ran on
those folders as well.

.flake8

@@ -1,5 +1,6 @@
 [flake8]
 # Longer line-length is more convenient for math
 max-line-length = 88
+exclude = .tox,pennylane
 extend-ignore = W605,  # Ignoring latex syntax in docstrings
                 E201   # Ignoring whitespace after '(' or '[' for matrices

piquasso/functions.py

@@ -0,0 +1,29 @@
+#
+# Copyright (C) 2020 by TODO - All rights reserved.
+#
+
+import numpy as np
+
+
+def gaussian_wigner_function(quadrature_matrix, *, d, mean, cov):
+    assert len(quadrature_matrix[0]) == len(mean), (
+        "'quadrature_matrix' elements should have the same dimension as 'mean': "
+        f"dim(quadrature_matrix[0])={len(quadrature_matrix[0])}, dim(mean)={len(mean)}."
+    )
+
+    return [
+        _gaussian_wigner_function_for_scalar(quadrature_array, d=d, mean=mean, cov=cov)
+        for quadrature_array in quadrature_matrix
+    ]
+
+
+def _gaussian_wigner_function_for_scalar(X, *, d, mean, cov):
+    return (
+        (1 / (np.pi ** d))
+        * np.sqrt((1 / np.linalg.det(cov)))
+        * np.exp(
+            - (X - mean)
+            @ np.linalg.inv(cov)
+            @ (X - mean)
+        )
+    ).real

piquasso/gaussian/state.py

@@ -4,7 +4,7 @@
 
 import numpy as np
 
-from piquasso import constants
+from piquasso import constants, functions
 from piquasso.context import Context
 from piquasso.state import State
 
@@ -360,6 +360,43 @@ def reduced_rotated_mean_and_cov(self, modes, phi):
 
         return transformed_state.mu, transformed_state.cov
 
+    def wigner_function(self, quadrature_matrix, modes=None):
+        r"""
+        Calculates the Wigner function values at the specified `quadrature_matrix`,
+        according to the equation
+
+        .. math::
+            W(r) = \frac{1}{\pi^d \sqrt{\mathrm{det} \sigma}}
+                \exp \big (
+                    - (r - \mu)^T
+                    \sigma^{-1}
+                    (r - \mu)
+                \big ).
+
+        Args:
+            quadrature_matrix (list): list of canonical coordinates vectors.
+            modes (tuple or None): modes where Wigner function should be calculcated.
+
+        Returns:
+            tuple: The Wigner function values in the shape of `quadrature_matrix`.
+        """
+
+        if modes:
+            reduced_state = self.reduced(modes)
+            return functions.gaussian_wigner_function(
+                quadrature_matrix,
+                d=reduced_state.d,
+                mean=reduced_state.mu,
+                cov=reduced_state.cov
+            )
+
+        return functions.gaussian_wigner_function(
+            quadrature_matrix,
+            d=self.d,
+            mean=self.mu,
+            cov=self.cov,
+        )
+
     def apply_passive(self, T, modes):
         r"""Applies a passive transformation to the quantum state.
 

piquasso/gaussian/tests/test_state.py

@@ -133,6 +133,25 @@ def test_quad_representation(
             self.state.cov,
         )
 
+    def test_wigner_function(self, mu, cov):
+        quadrature_array = np.array(
+            [
+                [1, 2, 3, 4],
+                [5, 6, 7, 8],
+            ]
+        )
+
+        expected_result = np.array(
+            [
+                0.00040656676635938727,
+                0.01037619639200025,
+            ]
+        )
+
+        actual_result = self.state.wigner_function(quadrature_array)
+
+        assert np.allclose(expected_result, actual_result)
+
 
 class TestGaussianStateOperations:
     @pytest.fixture

piquasso/tests/test_functions.py

@@ -0,0 +1,64 @@
+#
+# Copyright (C) 2020 by TODO - All rights reserved.
+#
+
+import pytest
+import numpy as np
+
+from piquasso import functions
+
+
+@pytest.fixture
+def d():
+    return 1
+
+
+@pytest.fixture
+def mean():
+    return np.array([1, 2])
+
+
+@pytest.fixture
+def cov():
+    return np.array(
+        [
+            [1, 2],
+            [-3, 4],
+        ]
+    )
+
+
+def test_wigner_function_at_scalar(d, mean, cov):
+    quadrature_array = np.array([1, 2])
+
+    expected = 0.10065842420897406
+
+    actual = functions._gaussian_wigner_function_for_scalar(
+        X=quadrature_array, d=d, mean=mean, cov=cov
+    )
+
+    assert np.allclose(expected, actual)
+
+
+def test_gaussian_wigner_function_handles_vectors(d, mean, cov):
+    quadrature_matrix = np.array(
+        [
+            [1, 2],
+            [3, 4],
+            [5, 6],
+        ]
+    )
+
+    expected = np.array(
+        [
+            0.10065842420897406,
+            0.009131526225575573,
+            6.81746788883418e-06,
+        ],
+    )
+
+    actual = functions.gaussian_wigner_function(
+        quadrature_matrix, d=d, mean=mean, cov=cov
+    )
+
+    assert np.allclose(expected, actual)