Fix period. in Fourier gen, update tests, examples

examples/00_misc/06_fourier.py

@@ -21,16 +21,12 @@
 model = gs.Gaussian(dim=2, var=1, len_scale=200)
 
 # Next, we hand the cov. model to the spatial random field class
-# and set the generator to `Fourier`. We will let the class figure out the
-# modes internally, by handing over `period` and `mode_rel_cutoff` which is the cutoff
-# value of the spectral density, relative to the maximum spectral density at
-# the origin. Simply put, we will use `mode_rel_cutoff`% of the spectral
-# density for the calculations. The argument `period` is set to the domain
-# size.
+# and set the generator to `Fourier`. The `mode_no` argument sets the number of
+# Fourier modes per dimension. The argument `period` is set to the domain size.
 srf = gs.SRF(
     model,
     generator="Fourier",
-    mode_rel_cutoff=0.99,
+    mode_no=[32, 32],
     period=L,
     seed=1681903,
 )
@@ -40,24 +36,3 @@
 
 # GSTools has a few simple visualization methods built in.
 srf.plot()
-
-# Alternatively, we could calculate the modes ourselves and hand them over to
-# GSTools. Therefore, we set the cutoff values to absolut values in Fourier
-# space. But always check, if you cover enough of the spectral density to not
-# run into numerical problems.
-modes_cutoff = [1.0, 1.0]
-
-# Next, we have to compute the numerical step size in Fourier space. This choice
-# influences the periodicity, which we want to set to the domain size by
-modes_delta = 2 * np.pi / L
-
-# Now, we calculate the modes with
-modes = [np.arange(0, modes_cutoff[d], modes_delta[d]) for d in range(2)]
-
-# And we can create a new instance of the SRF class with our own modes.
-srf_modes = gs.SRF(
-    model,
-    generator="Fourier",
-    modes=modes,
-    seed=494754,
-)

examples/00_misc/07_fourier_trans.py

@@ -24,7 +24,7 @@
 srf = gs.SRF(
     model,
     generator="Fourier",
-    mode_rel_cutoff=0.999,
+    mode_no=[32, 32],
     period=L,
     seed=1681903,
 )

src/gstools/field/generator.py

@@ -599,49 +599,35 @@ class Fourier(Generator):
     def __init__(
         self,
         model,
-        mode_rel_cutoff=None,
-        period=None,
-        modes=None,
+        mode_no,
+        period,
         seed=None,
         verbose=False,
         **kwargs,
     ):
         if kwargs:
             warnings.warn("gstools.Fourier: **kwargs are ignored")
-        if (
-            mode_rel_cutoff is not None
-            and period is not None
-            and modes is not None
-        ):
-            warnings.warn(
-                "gstools.Fourier: mode_rel_cutoff & period are ignored, as "
-                "modes is provided."
-            )
-        if mode_rel_cutoff is None and period is None and modes is None:
-            raise ValueError("Fourier: No mode information provided.")
-
         dim = model.dim
-        if (
-            modes is None
-            and period is not None
-            and mode_rel_cutoff is not None
-        ):
-            if len(period) != dim:
-                raise ValueError("Fourier: Dimension mismatch.")
-            self._mode_rel_cutoff = mode_rel_cutoff
-            self._period = np.array(period)
-            self._delta_k = 2.0 * np.pi / self._period
-            modes_cutoff = self.calc_modes_cutoff(model, self._mode_rel_cutoff)
-            self._modes_cutoff = self._fill_to_dim(dim, modes_cutoff)
-            modes = [
-                np.arange(-self._modes_cutoff[d], self._modes_cutoff[d], self._delta_k[d])
-                for d in range(dim)
-            ]
-        elif modes is not None:
-            self._delta_k = np.array(
-                [modes[d][1] - modes[d][0] for d in range(dim)]
+        if len(mode_no) != dim:
+            raise ValueError(
+                "Fourier: Dimension mismatch in argument mode_no."
+            )
+        if len(period) != dim:
+            raise ValueError("Fourier: Dimension mismatch in argument period.")
+        if (np.asarray([m % 2 for m in mode_no]) != 0).any():
+            raise ValueError("Fourier: Odd mode_no not supported.")
+
+        self._period = np.array(period)
+        self._delta_k = 2.0 * np.pi / self._period
+        anis = np.insert(model.anis.copy(), 0, 1.0)
+        modes = [
+            np.arange(
+                -mode_no[d] / 2.0 * self._delta_k[d] / anis[d],
+                mode_no[d] / 2.0 * self._delta_k[d] / anis[d],
+                self._delta_k[d],
             )
-            self._period = 2.0 * np.pi / self._delta_k
+            for d in range(dim)
+        ]
 
         # initialize attributes
         self._modes = generate_grid(modes)
@@ -813,41 +799,6 @@ def _fill_to_dim(
             r = np.pad(r, (0, dim - len(r)), "edge")
         return r
 
-    def calc_modes_cutoff(
-        self, model, mode_rel_cutoff
-    ):  # pylint: disable=R6301
-        """Find the cutoff value so that `mode_rel_cutoff`% of the spectrum is kept.
-
-        This helper function uses a least squares algorithm to determine the
-        cutoff value so that `mode_rel_cutoff`% of the spectrum is kept.
-
-        Parameters
-        ----------
-        model : :any:`CovModel`
-            Covariance model
-        mode_rel_cutoff : :class:`float`
-
-        Returns
-        -------
-        :class:`float`
-            the cutoff value
-        """
-        norm = model.spectral_density(0)
-        # the first len_scale is a good enough first guess
-        try:
-            len_scale = model.len_scale[0]
-        except IndexError:
-            len_scale = model.len_scale
-        k_cutoff0 = np.sqrt(1.0 / len_scale)
-        mode_rel_cutoff_r = 1.0 - mode_rel_cutoff
-        res = least_squares(
-            lambda k, mode_rel_cutoff_r: model.spectral_density(k)
-            - mode_rel_cutoff_r * norm,
-            k_cutoff0,
-            args=(mode_rel_cutoff_r,),
-        )
-        return res.x[0]
-
     @property
     def seed(self):
         """:class:`int`: Seed of the master RNG.
@@ -873,6 +824,16 @@ def model(self):
     def model(self, model):
         self.update(model)
 
+    @property
+    def modes(self):
+        """:class:`numpy.ndarray`: Modes on which the spectrum is calculated."""
+        return self._modes
+
+    @property
+    def period(self):
+        """:class:`numpy.ndarray`: Period length of the spatial random field."""
+        return self._period
+
     @property
     def verbose(self):
         """:class:`bool`: Verbosity of the generator."""

tests/test_fouriergen.py

@@ -2,13 +2,11 @@
 This is the unittest of the Fourier class.
 """
 
-import copy
 import unittest
 
 import numpy as np
 
 import gstools as gs
-from gstools.field.generator import Fourier
 
 
 class TestFourier(unittest.TestCase):
@@ -17,62 +15,83 @@ def setUp(self):
         self.cov_model_1d = gs.Gaussian(dim=1, var=0.5, len_scale=10.0)
         self.cov_model_2d = gs.Gaussian(dim=2, var=2.0, len_scale=30.0)
         self.cov_model_3d = gs.Gaussian(dim=3, var=2.1, len_scale=21.0)
-        L = [80, 30, 91]
-        self.x = np.linspace(0, L[0], 11)
-        self.y = np.linspace(0, L[1], 31)
-        self.z = np.linspace(0, L[2], 13)
+        self.L = [80, 30, 91]
+        self.x = np.linspace(0, self.L[0], 11)
+        self.y = np.linspace(0, self.L[1], 31)
+        self.z = np.linspace(0, self.L[2], 13)
 
-        cutoff_rel = 0.999
-        cutoff_abs = 1
-        dk = [2 * np.pi / l for l in L]
-
-        self.modes_1d = [np.arange(0, cutoff_abs, dk[0])]
-        self.modes_2d = self.modes_1d + [np.arange(0, cutoff_abs, dk[1])]
-        self.modes_3d = self.modes_2d + [np.arange(0, cutoff_abs, dk[2])]
+        self.mode_no = [12, 6, 14]
 
         self.srf_1d = gs.SRF(
             self.cov_model_1d,
             generator="Fourier",
-            modes=self.modes_1d,
+            mode_no=[self.mode_no[0]],
+            period=[self.L[0]],
             seed=self.seed,
         )
         self.srf_2d = gs.SRF(
             self.cov_model_2d,
             generator="Fourier",
-            modes=self.modes_2d,
+            mode_no=self.mode_no[:2],
+            period=self.L[:2],
             seed=self.seed,
         )
         self.srf_3d = gs.SRF(
             self.cov_model_3d,
             generator="Fourier",
-            mode_rel_cutoff=cutoff_rel,
-            period=L,
+            mode_no=self.mode_no,
+            period=self.L,
             seed=self.seed,
         )
 
     def test_1d(self):
         field = self.srf_1d((self.x,), mesh_type="structured")
-        self.assertAlmostEqual(field[0], 0.40939877176695477)
+        self.assertAlmostEqual(field[0], 0.6236929351309081)
 
     def test_2d(self):
         field = self.srf_2d((self.x, self.y), mesh_type="structured")
-        self.assertAlmostEqual(field[0, 0], 1.6338790313270515)
+        self.assertAlmostEqual(field[0, 0], -0.1431996611581266)
 
     def test_3d(self):
         field = self.srf_3d((self.x, self.y, self.z), mesh_type="structured")
-        self.assertAlmostEqual(field[0, 0, 0], 0.3866689424599251)
+        self.assertAlmostEqual(field[0, 0, 0], -1.0433325279452803)
+
+    def test_periodicity_1d(self):
+        field = self.srf_1d((self.x,), mesh_type="structured")
+        self.assertAlmostEqual(field[0], field[-1])
 
-    def test_periodicity(self):
+    def test_periodicity_2d(self):
         field = self.srf_2d((self.x, self.y), mesh_type="structured")
         self.assertAlmostEqual(
             field[0, len(self.y) // 2], field[-1, len(self.y) // 2]
         )
+        self.assertAlmostEqual(
+            field[len(self.x) // 2, 0], field[len(self.x) // 2, -1]
+        )
+
+    def test_periodicity_3d(self):
+        field = self.srf_3d((self.x, self.y, self.z), mesh_type="structured")
+        self.assertAlmostEqual(
+            field[0, len(self.y) // 2, 0], field[-1, len(self.y) // 2, 0]
+        )
+        self.assertAlmostEqual(field[0, 0, 0], field[0, -1, 0])
+        self.assertAlmostEqual(
+            field[len(self.x) // 2, len(self.y) // 2, 0],
+            field[len(self.x) // 2, len(self.y) // 2, -1],
+        )
 
     def test_assertions(self):
         # unstructured grids not supported
         self.assertRaises(ValueError, self.srf_2d, (self.x, self.y))
         self.assertRaises(
             ValueError, self.srf_2d, (self.x, self.y), mesh_type="unstructured"
         )
-        with self.assertRaises(ValueError):
-            gs.SRF(self.cov_model_2d, generator="Fourier")
+        self.assertRaises(
+            ValueError,
+            gs.SRF,
+            self.cov_model_2d,
+            generator="Fourier",
+            mode_no=[13, 50],
+            period=self.L[:2],
+            seed=self.seed,
+        )