diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml index 7d1c7d01..a3031e6a 100644 --- a/.github/workflows/main.yml +++ b/.github/workflows/main.yml @@ -48,9 +48,9 @@ jobs: run: | python -m pylint src/gstools/ - - name: cython-lint check - run: | - cython-lint src/gstools/ + #- name: cython-lint check + #run: | + #cython-lint src/gstools/ build_wheels: name: wheels for ${{ matrix.os }} diff --git a/README.md b/README.md index fd5baf63..b90ad129 100644 --- a/README.md +++ b/README.md @@ -34,7 +34,7 @@ GeoStatTools provides geostatistical tools for various purposes: -- random field generation +- random field generation, including periodic boundaries - simple, ordinary, universal and external drift kriging - conditioned field generation - incompressible random vector field generation diff --git a/docs/source/index.rst b/docs/source/index.rst index 0f5e181f..7abb1e9e 100644 --- a/docs/source/index.rst +++ b/docs/source/index.rst @@ -24,7 +24,7 @@ Purpose GeoStatTools provides geostatistical tools for various purposes: -- random field generation +- random field generation, including periodic boundaries - simple, ordinary, universal and external drift kriging - conditioned field generation - incompressible random vector field generation diff --git a/examples/01_random_field/08_fourier.py b/examples/01_random_field/08_fourier.py new file mode 100644 index 00000000..8b700ae6 --- /dev/null +++ b/examples/01_random_field/08_fourier.py @@ -0,0 +1,44 @@ +""" +Generating a Simple Periodic Random Field +----------------------------------------- + +In this simple example we are going to learn how to generate periodic spatial +random fields. The Fourier method comes naturally with the property of +periodicity, so we'll use it to create the random field. +""" + +import numpy as np + +import gstools as gs + +# We start off by defining the spatial grid. For the sake of simplicity, we +# use a square domain. We set the optional argument `endpoint` to `False`, to +# not make the domain in each dimension one grid cell larger than the +# periodicity. +L = 500.0 +x = np.linspace(0, L, 256, endpoint=False) +y = np.linspace(0, L, 128, endpoint=False) + +# Now, we create a Gaussian covariance model with a correlation length which is +# roughly half the size of the grid. +model = gs.Gaussian(dim=2, var=1, len_scale=200) + +# Next, we hand the cov. model to the spatial random field class `SRF` +# and set the generator to `"Fourier"`. The argument `period` is set to the +# domain size. If only a single number is given, the same periodicity is +# applied in each dimension, as shown in this example. The `mode_no` argument +# sets the number of Fourier modes. If only an integer is given, that number +# of modes is used for all dimensions. +srf = gs.SRF( + model, + generator="Fourier", + period=L, + mode_no=32, + seed=1681903, +) + +# Now, we can calculate the field with the given parameters. +srf((x, y), mesh_type="structured") + +# GSTools has a few simple visualization methods built in. +srf.plot() diff --git a/examples/01_random_field/09_fourier_trans.py b/examples/01_random_field/09_fourier_trans.py new file mode 100644 index 00000000..81c303f3 --- /dev/null +++ b/examples/01_random_field/09_fourier_trans.py @@ -0,0 +1,49 @@ +""" +Generating a Transformed Periodic Random Field +---------------------------------------------- + +Building on the precious example, we are now going to generate periodic +spatial random fields with a transformation applied, resulting in a level set. +""" + +import numpy as np + +import gstools as gs + +# We start off by defining the spatial grid. As in the previous example, we do +# not want to include the endpoints. +L = np.array((500, 400)) +x = np.linspace(0, L[0], 300, endpoint=False) +y = np.linspace(0, L[1], 200, endpoint=False) + +# Instead of using a Gaussian covariance model, we will use the much rougher +# exponential model and we will introduce an anisotropy by using two different +# length scales in the x- and y-directions +model = gs.Exponential(dim=2, var=2, len_scale=[80, 20]) + +# Same as before, we set up the spatial random field. But this time, we will +# use a periodicity which is equal to the domain size in x-direction, but +# half the domain size in y-direction. And we will use different `mode_no` for +# the different dimensions. +srf = gs.SRF( + model, + generator="Fourier", + period=[L[0], L[1] / 2], + mode_no=[30, 20], + seed=1681903, +) +# and compute it on our spatial domain +srf((x, y), mesh_type="structured") + +# With the field generated, we can now apply transformations starting with a +# discretization of the field into 4 different values +thresholds = np.linspace(np.min(srf.field), np.max(srf.field), 4) +srf.transform("discrete", store="transform_discrete", values=thresholds) +srf.plot("transform_discrete") + +# This is already a nice result, but we want to pronounce the peaks of the +# field. We can do this by applying a log-normal transformation on top +srf.transform( + "lognormal", field="transform_discrete", store="transform_lognormal" +) +srf.plot("transform_lognormal") diff --git a/examples/01_random_field/README.rst b/examples/01_random_field/README.rst index 6b226b2f..ab723bae 100644 --- a/examples/01_random_field/README.rst +++ b/examples/01_random_field/README.rst @@ -11,6 +11,12 @@ semi-variogram. This is done by using the so-called randomization method. The spatial random field is represented by a stochastic Fourier integral and its discretised modes are evaluated at random frequencies. +In case you want to generate spatial random fields with periodic boundaries, +you can use the so-called Fourier method. See the corresponding examples for +how to do that. The spatial random field is represented by a stochastic +Fourier integral and its discretised modes are evaluated at equidistant +frequencies. + GSTools supports arbitrary and non-isotropic covariance models. Examples diff --git a/src/gstools/field/generator.py b/src/gstools/field/generator.py old mode 100644 new mode 100755 index d19bd6d3..d4cb0dea --- a/src/gstools/field/generator.py +++ b/src/gstools/field/generator.py @@ -11,6 +11,7 @@ Generator RandMeth IncomprRandMeth + Fourier """ # pylint: disable=C0103, W0222, C0412, W0231 @@ -23,15 +24,18 @@ from gstools import config from gstools.covmodel.base import CovModel from gstools.field.summator import summate as summate_c +from gstools.field.summator import summate_fourier as summate_fourier_c from gstools.field.summator import summate_incompr as summate_incompr_c from gstools.random.rng import RNG +from gstools.tools.geometric import generate_grid if config._GSTOOLS_CORE_AVAIL: # pylint: disable=W0212; # pragma: no cover # pylint: disable=E0401 from gstools_core import summate as summate_gsc + from gstools_core import summate_fourier as summate_fourier_gsc from gstools_core import summate_incompr as summate_incompr_gsc -__all__ = ["Generator", "RandMeth", "IncomprRandMeth"] +__all__ = ["Generator", "RandMeth", "IncomprRandMeth", "Fourier"] SAMPLING = ["auto", "inversion", "mcmc"] @@ -68,6 +72,20 @@ def _summate_incompr( return summate_incompr_fct(cov_samples, z_1, z_2, pos, num_threads) +def _summate_fourier(spectrum_factor, modes, z_1, z_2, pos, num_threads=None): + """A wrapper function for calling the Fourier algorithms.""" + if ( + config.USE_GSTOOLS_CORE + and config._GSTOOLS_CORE_AVAIL # pylint: disable=W0212 + ): + summate_fourier_fct = summate_fourier_gsc # pylint: disable=E0606 + else: + summate_fourier_fct = summate_fourier_c + return summate_fourier_fct( + spectrum_factor, modes, z_1, z_2, pos, num_threads + ) + + class Generator(ABC): """ Abstract generator class. @@ -303,7 +321,7 @@ def update(self, model=None, seed=np.nan): raise ValueError( "gstools.field.generator.RandMeth: no 'model' given" ) - # if the user tries to trick us, we beat him! + # if the user tries to trick us, we beat them! elif model is None and np.isnan(seed): if not ( isinstance(self._model, CovModel) @@ -564,3 +582,326 @@ def _create_unit_vector(self, broadcast_shape, axis=0): e1 = np.zeros(shape) e1[axis] = 1.0 return e1 + + +class Fourier(Generator): + r"""Fourier method for calculating periodic, isotropic random fields. + + Parameters + ---------- + model : :any:`CovModel` + Covariance model + period : :class:`list` or :class:`float` + The spatial periodicity of the field, is often the domain size. + mode_no : :class:`list` or :class:`float`, optional + Number of Fourier modes per dimension. + seed : :class:`int`, optional + The seed of the random number generator. + If "None", a random seed is used. Default: :any:`None` + + **kwargs + Placeholder for keyword-args + + Notes + ----- + The Fourier method is used to generate periodic isotropic spatial random + fields characterized by a given covariance model. + The calculation looks like: + + .. math:: + u\left(x\right)= + \sum_{i=1}^{N}\sqrt{2S(k_{i})\Delta k}\left( + Z_{1,i}\cdot\cos\left(\left\langle k_{i},x\right\rangle \right)+ + Z_{2,i}\cdot\sin\left(\left\langle k_{i},x\right\rangle \right) + \right) + + where: + + * :math:`S` : spectrum of the covariance model + * :math:`Z_{j,i}` : random samples from a normal distribution + * :math:`k_i` : the equidistant Fourier grid + """ + + def __init__( + self, + model, + period, + mode_no=32, + seed=None, + **kwargs, + ): + if kwargs: + warnings.warn("gstools.Fourier: **kwargs are ignored") + + # initialize private attributes + self._modes = None + self._period = None + self._mode_no = None + self._delta_k = None + self._model = None + self._seed = None + self._rng = None + self._z_1 = None + self._z_2 = None + self._spectrum_factor = None + self._value_type = "scalar" + # set model and seed + self.update(model, seed, period, mode_no) + + def __call__(self, pos, add_nugget=True): + """Calculate the modes for the Fourier method. + + This method calls the `summate_fourier` Cython method, which is the + heart of the Fourier method. + + Parameters + ---------- + pos : (d, n), :class:`numpy.ndarray` + the position tuple with d dimensions and n points. + add_nugget : :class:`bool` + Whether to add nugget noise to the field. + + Returns + ------- + :class:`numpy.ndarray` + the random modes + """ + pos = np.asarray(pos, dtype=np.double) + + summed_modes = _summate_fourier( + self._spectrum_factor, + self._modes, + self._z_1, + self._z_2, + pos, + config.NUM_THREADS, + ) + nugget = self.get_nugget(summed_modes.shape) if add_nugget else 0.0 + return summed_modes + nugget + + def get_nugget(self, shape): + """ + Generate normal distributed values for the nugget simulation. + + Parameters + ---------- + shape : :class:`tuple` + the shape of the summed modes + + Returns + ------- + nugget : :class:`numpy.ndarray` + the nugget in the same shape as the summed modes + """ + if self.model.nugget > 0: + nugget = np.sqrt(self.model.nugget) * self._rng.random.normal( + size=shape + ) + else: + nugget = 0.0 + return nugget + + def update(self, model=None, seed=np.nan, period=None, mode_no=None): + """Update the model and the seed. + + If model and seed are not different, nothing will be done. + + Parameters + ---------- + model : :any:`CovModel` or :any:`None`, optional + covariance model. Default: :any:`None` + seed : :class:`int` or :any:`None` or :any:`numpy.nan`, optional + the seed of the random number generator. + If :any:`None`, a random seed is used. If :any:`numpy.nan`, + the actual seed will be kept. Default: :any:`numpy.nan` + period : :class:`list` or :any:`None, optional + The spatial periodicity of the field, is often the domain size. + mode_no : :class:`list` or :any:`None`, optional + Number of Fourier modes per dimension. + """ + tmp_model = model if model is not None else self._model + dim = tmp_model.dim + if period is not None: + self._period = self._fill_to_dim(period, dim) + anis = np.insert(tmp_model.anis.copy(), 0, 1.0) + self._delta_k = 2.0 * np.pi / self._period * anis + if mode_no is None: + self._set_modes(self._mode_no, dim) + if mode_no is not None: + mode_no = self._fill_to_dim(mode_no, dim) + if (np.asarray([m % 2 for m in mode_no]) != 0).any(): + raise ValueError("Fourier: Odd mode_no not supported.") + self._set_modes(mode_no, dim) + + # check if a new model is given + if isinstance(model, CovModel): + if self.model != model: + self._model = dcp(model) + if seed is None or not np.isnan(seed): + self.reset_seed(seed) + else: + self.reset_seed(self._seed) + # just update the seed, if its a new one + elif seed is None or not np.isnan(seed): + self.seed = seed + # or just update the seed, when no model is given + elif model is None and (seed is None or not np.isnan(seed)): + if isinstance(self._model, CovModel): + self.seed = seed + else: + raise ValueError( + "gstools.field.generator.Fourier: no 'model' given" + ) + # but also update when mode mesh was modified + elif mode_no is not None or period is not None: + if seed is None or not np.isnan(seed): + self.reset_seed(seed) + else: + self.reset_seed(self._seed) + # if the user tries to trick us, we beat them! + elif model is None and np.isnan(seed): + if ( + isinstance(self._model, CovModel) + and self._z_1 is not None + and self._z_2 is not None + and self._spectrum_factor is not None + ): + raise ValueError( + "gstools.field.generator.Fourier: " + "neither 'model' nor 'seed' given!" + ) + # wrong model type + else: + raise ValueError( + "gstools.field.generator.Fourier: 'model' is not an " + "instance of 'gstools.CovModel'" + ) + + def reset_seed(self, seed=np.nan): + """ + Recalculate the random values with the given seed. + + Parameters + ---------- + seed : :class:`int` or :any:`None` or :any:`numpy.nan`, optional + the seed of the random number generator. + If :any:`None`, a random seed is used. If :any:`numpy.nan`, + the actual seed will be kept. Default: :any:`numpy.nan` + + Notes + ----- + Even if the given seed is the present one, modes will be recalculated. + """ + if seed is None or not np.isnan(seed): + self._seed = seed + self._rng = RNG(self._seed) + # normal distributed samples for randmeth + self._z_1 = self._rng.random.normal(size=np.prod(self._mode_no)) + self._z_2 = self._rng.random.normal(size=np.prod(self._mode_no)) + # pre calc. the spectrum for all wave numbers they are handed over to + # Cython, which doesn't have access to the CovModel + k_norm = np.linalg.norm(self._modes, axis=0) + self._spectrum_factor = np.sqrt( + self._model.spectrum(k_norm) * np.prod(self._delta_k) + ) + + def _fill_to_dim( + self, values, dim, dtype=float, default_value=None + ): # pylint: disable=R6301 + """Fill an array with last element up to len(dim).""" + r = np.atleast_1d(values) + if values is None: + if default_value is None: + raise ValueError("Fourier: Value has to be provided") + r = default_value + r = np.array(r, dtype=dtype) + r = np.atleast_1d(r)[:dim] + if len(r) > dim: + raise ValueError(f"Fourier: len(values) <= {dim=} not fulfilled") + # fill up values with values[-1], such that len()==dim + if len(r) < dim: + r = np.pad(r, (0, dim - len(r)), "edge") + return r + + def _set_modes(self, mode_no, dim): + """Calculate the mode mesh. + + Parameters + ---------- + mode_no : :class:`list` + Number of Fourier modes per dimension. + dim : :class:`int` + dimension of the model. + + Notes + ----- + `self._reset_seed` *has* to be called after this method! + """ + modes = [ + np.arange( + -mode_no[d] / 2.0 * self._delta_k[d], + mode_no[d] / 2.0 * self._delta_k[d], + self._delta_k[d], + ) + for d in range(dim) + ] + # initialize attributes + self._modes = generate_grid(modes) + self._mode_no = [len(m) for m in modes] + + @property + def seed(self): + """:class:`int`: Seed of the master RNG. + + Notes + ----- + If a new seed is given, the setter property not only saves the + new seed, but also creates new random modes with the new seed. + """ + return self._seed + + @seed.setter + def seed(self, new_seed): + if new_seed is not self._seed: + self.reset_seed(new_seed) + + @property + def model(self): + """:any:`CovModel`: Covariance model of the spatial random field.""" + return self._model + + @model.setter + 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 mode_no(self): + """:class:`numpy.ndarray`: Number of modes per dimension.""" + return self._mode_no + + @mode_no.setter + def mode_no(self, mode_no): + self.update(mode_no=mode_no) + + @property + def period(self): + """:class:`numpy.ndarray`: Period length of the spatial random field.""" + return self._period + + @period.setter + def period(self, period): + self.update(period=period) + + @property + def value_type(self): + """:class:`str`: Type of the field values (scalar, vector).""" + return self._value_type + + def __repr__(self): + """Return String representation.""" + return f"{self.name}(model={self.model}, seed={self.seed})" diff --git a/src/gstools/field/srf.py b/src/gstools/field/srf.py old mode 100644 new mode 100755 index d88e46c0..c1d5081e --- a/src/gstools/field/srf.py +++ b/src/gstools/field/srf.py @@ -14,7 +14,12 @@ import numpy as np from gstools.field.base import Field -from gstools.field.generator import Generator, IncomprRandMeth, RandMeth +from gstools.field.generator import ( + Fourier, + Generator, + IncomprRandMeth, + RandMeth, +) from gstools.field.upscaling import var_coarse_graining, var_no_scaling __all__ = ["SRF"] @@ -24,6 +29,7 @@ "IncomprRandMeth": IncomprRandMeth, "VectorField": IncomprRandMeth, "VelocityField": IncomprRandMeth, + "Fourier": Fourier, } """dict: Standard generators for spatial random fields.""" @@ -75,6 +81,7 @@ class SRF(Field): See: :any:`IncomprRandMeth` * "VectorField" : an alias for "IncomprRandMeth" * "VelocityField" : an alias for "IncomprRandMeth" + * "Fourier" : the periodic Fourier method Default: "RandMeth" **generator_kwargs diff --git a/src/gstools/field/summator.pyx b/src/gstools/field/summator.pyx old mode 100644 new mode 100755 index 5af23b91..416024f1 --- a/src/gstools/field/summator.pyx +++ b/src/gstools/field/summator.pyx @@ -10,7 +10,7 @@ if OPENMP: cimport openmp cimport numpy as np -from libc.math cimport cos, sin +from libc.math cimport cos, pi, sin, sqrt def set_num_threads(num_threads): @@ -97,3 +97,34 @@ def summate_incompr( proj[d] * (z_1[j] * cos(phase) + z_2[j] * sin(phase)) ) return np.asarray(summed_modes) + + +def summate_fourier( + const double[:] spectrum_factor, + const double[:, :] modes, + const double[:] z_1, + const double[:] z_2, + const double[:, :] pos, + num_threads=None, +): + cdef int i, j, d + cdef double phase + cdef int dim = pos.shape[0] + + cdef int X_len = pos.shape[1] + cdef int N = modes.shape[1] + + cdef double[:] summed_modes = np.zeros(X_len, dtype=float) + + cdef int num_threads_c = set_num_threads(num_threads) + + for i in prange(X_len, nogil=True, num_threads=num_threads_c): + for j in range(N): + phase = 0. + for d in range(dim): + phase += modes[d, j] * pos[d, i] + summed_modes[i] += ( + spectrum_factor[j] * (z_1[j] * cos(phase) + z_2[j] * sin(phase)) + ) + + return np.asarray(summed_modes) diff --git a/tests/test_fouriergen.py b/tests/test_fouriergen.py new file mode 100755 index 00000000..093b47d5 --- /dev/null +++ b/tests/test_fouriergen.py @@ -0,0 +1,111 @@ +""" +This is the unittest of the Fourier class. +""" + +import unittest + +import numpy as np + +import gstools as gs + + +class TestFourier(unittest.TestCase): + def setUp(self): + self.seed = 19900408 + 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) + 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) + + self.mode_no = [12, 6, 14] + + self.srf_1d = gs.SRF( + self.cov_model_1d, + generator="Fourier", + mode_no=[self.mode_no[0]], + period=[self.L[0]], + seed=self.seed, + ) + self.srf_2d = gs.SRF( + self.cov_model_2d, + generator="Fourier", + 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_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.6236929351309081) + + def test_2d(self): + field = self.srf_2d((self.x, self.y), mesh_type="structured") + 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], -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_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_setters(self): + new_period = [5, 10] + self.srf_2d.generator.period = new_period + np.testing.assert_almost_equal( + self.srf_2d.generator.period, + np.array(new_period), + ) + new_mode_no = [6, 6] + self.srf_2d.generator.mode_no = new_mode_no + np.testing.assert_almost_equal( + self.srf_2d.generator.mode_no, + np.array(new_mode_no), + ) + + 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" + ) + self.assertRaises( + ValueError, + gs.SRF, + self.cov_model_2d, + generator="Fourier", + mode_no=[13, 50], + period=self.L[:2], + seed=self.seed, + )