483994410

Build Type

github

Committed by GitHub

Commit Message

<a href="https://github.com/GeoStat-Framework/GSTools/commit/<a class=hub.com/GeoStat-Framework/GSTools/commit/29998a6c42b2abe4c796663709777921e9cf1346">29998a6c4<a href="https://github.com/GeoStat-Framework/GSTools/commit/29998a6c42b2abe4c796663709777921e9cf1346">">Merge </a><a class="double-link" href="https://github.com/GeoStat-Framework/GSTools/commit/<a class="double-link" href="https://github.com/GeoStat-Framework/GSTools/commit/079f5cef4abaca3e147c36704653d24d56320f74">079f5cef4</a>">079f5cef4</a><a href="https://github.com/GeoStat-Framework/GSTools/commit/29998a6c42b2abe4c796663709777921e9cf1346"> into c1b56dadc">c1b56dadc</a>

Pull Request Pull Request #130: New Class for Conditioned Spatial Random Fields

Run Details

108 of 108 new or added lines in 7 files covered. (100.0%)

2111 of 2329 relevant lines covered (90.64%)

13.59 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

84.96

/gstools/field/generator.py

# -*- coding: utf-8 -*-
"""
GStools subpackage providing generators for spatial random fields.

.. currentmodule:: gstools.field.generator

The following classes are provided

.. autosummary::
   RandMeth
   IncomprRandMeth
"""
# pylint: disable=C0103

from copy import deepcopy as dcp
import numpy as np
from gstools.covmodel.base import CovModel
from gstools.random.rng import RNG
from gstools.field.summator import summate, summate_incompr


__all__ = ["RandMeth", "IncomprRandMeth"]


SAMPLING = ["auto", "inversion", "mcmc"]


class RandMeth:
    r"""Randomization method for calculating isotropic spatial random fields.

    Parameters
    ----------
    model : :any:`CovModel`
        Covariance model
    mode_no : :class:`int`, optional
        Number of Fourier modes. Default: ``1000``
    seed : :class:`int` or :any:`None`, optional
        The seed of the random number generator.
        If "None", a random seed is used. Default: :any:`None`
    verbose : :class:`bool`, optional
        Be chatty during the generation.
        Default: :any:`False`
    sampling :class:`str`, optional
        Sampling strategy. Either

            * "auto": select best strategy depending on given model
            * "inversion": use inversion method
            * "mcmc": use mcmc sampling

    **kwargs
        Placeholder for keyword-args

    Notes
    -----
    The Randomization method is used to generate isotropic
    spatial random fields characterized by a given covariance model.
    The calculation looks like:

    .. math::
       u\left(x\right)=
       \sqrt{\frac{\sigma^{2}}{N}}\cdot
       \sum_{i=1}^{N}\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:`N` : fourier mode number
        * :math:`Z_{j,i}` : random samples from a normal distribution
        * :math:`k_i` : samples from the spectral density distribution of
          the covariance model
    """

    def __init__(
        self,
        model,
        mode_no=1000,
        seed=None,
        verbose=False,
        sampling="auto",
        **kwargs
    ):
        if kwargs:
            print("gstools.RandMeth: **kwargs are ignored")
        # initialize atributes
        self._mode_no = int(mode_no)
        self._verbose = bool(verbose)
        # initialize private atributes
        self._model = None
        self._seed = None
        self._rng = None
        self._z_1 = None
        self._z_2 = None
        self._cov_sample = None
        self._value_type = "scalar"
        # set sampling strategy
        self._sampling = None
        self.sampling = sampling
        # set model and seed
        self.update(model, seed)

    def __call__(self, pos, add_nugget=True):
        """Calculate the random modes for the randomization method.

        This method  calls the `summate_*` Cython methods, which are the
        heart of the randomization method.

        Parameters
        ----------
        pos : (d, n), :class:`numpy.ndarray`
            the position tuple with d dimensions and n points.
        add_nugget : :class:`bool`
            Wheater to add nugget noise to the field.

        Returns
        -------
        :class:`numpy.ndarray`
            the random modes
        """
        pos = np.array(pos, dtype=np.double)
        summed_modes = summate(self._cov_sample, self._z_1, self._z_2, pos)
        nugget = self.get_nugget(summed_modes.shape) if add_nugget else 0.0
        return np.sqrt(self.model.var / self._mode_no) * 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):
        """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`
        """
        # 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.RandMeth: no 'model' given"
                )
        # if the user tries to trick us, we beat him!
        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._cov_sample is not None
            ):
                if self.verbose:
                    print("RandMeth.update: Nothing will be done...")
            else:
                raise ValueError(
                    "gstools.field.generator.RandMeth: "
                    + "neither 'model' nor 'seed' given!"
                )
        # wrong model type
        else:
            raise ValueError(
                "gstools.field.generator.RandMeth: 'model' is not an "
                + "instance of 'gstools.CovModel'"
            )

    def reset_seed(self, seed=np.nan):
        """
        Recalculate the random amplitudes and wave numbers 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=self._mode_no)
        self._z_2 = self._rng.random.normal(size=self._mode_no)
        # sample uniform on a sphere
        sphere_coord = self._rng.sample_sphere(self.model.dim, self._mode_no)
        # sample radii acording to radial spectral density of the model
        if self.sampling == "inversion" or (
            self.sampling == "auto" and self.model.has_ppf
        ):
            pdf, cdf, ppf = self.model.dist_func
            rad = self._rng.sample_dist(
                size=self._mode_no, pdf=pdf, cdf=cdf, ppf=ppf, a=0
            )
        else:
            rad = self._rng.sample_ln_pdf(
                ln_pdf=self.model.ln_spectral_rad_pdf,
                size=self._mode_no,
                sample_around=1.0 / self.model.len_rescaled,
            )
        # get fully spatial samples by multiplying sphere samples and radii
        self._cov_sample = rad * sphere_coord

    @property
    def sampling(self):
        """:class:`str`: Sampling strategy."""
        return self._sampling

    @sampling.setter
    def sampling(self, sampling):
        if sampling not in ["auto", "inversion", "mcmc"]:
            raise ValueError("RandMeth: sampling not in {}.".format(SAMPLING))
        self._sampling = sampling

    @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 mode_no(self):
        """:class:`int`: Number of modes in the randomization method."""
        return self._mode_no

    @mode_no.setter
    def mode_no(self, mode_no):
        if int(mode_no) != self._mode_no:
            self._mode_no = int(mode_no)
            self.reset_seed(self._seed)

    @property
    def verbose(self):
        """:class:`bool`: Verbosity of the generator."""
        return self._verbose

    @verbose.setter
    def verbose(self, verbose):
        self._verbose = bool(verbose)

    @property
    def name(self):
        """:class:`str`: Name of the generator."""
        return self.__class__.__name__

    @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 "RandMeth(model={0}, mode_no={1}, seed={2})".format(
            repr(self.model), self._mode_no, self.seed
        )


class IncomprRandMeth(RandMeth):
    r"""RandMeth for incompressible random vector fields.

    Parameters
    ----------
    model : :any:`CovModel`
        covariance model
    mean_velocity : :class:`float`, optional
        the mean velocity in x-direction
    mode_no : :class:`int`, optional
        number of Fourier modes. Default: ``1000``
    seed : :class:`int` or :any:`None`, optional
        the seed of the random number generator.
        If "None", a random seed is used. Default: :any:`None`
    verbose : :class:`bool`, optional
        State if there should be output during the generation.
        Default: :any:`False`
    sampling :class:`str`, optional
        Sampling strategy. Either

            * "auto": select best strategy depending on given model
            * "inversion": use inversion method
            * "mcmc": use mcmc sampling

    **kwargs
        Placeholder for keyword-args

    Notes
    -----
    The Randomization method is used to generate isotropic
    spatial incompressible random vector fields characterized
    by a given covariance model. The equation is:

    .. math::
       u_i\left(x\right)= \bar{u_i} \delta_{i1} +
       \bar{u_i}\sqrt{\frac{\sigma^{2}}{N}}\cdot
       \sum_{j=1}^{N}p_i(k_{j})\left(
       Z_{1,j}\cdot\cos\left(\left\langle k_{j},x\right\rangle \right)+
       Z_{2,j}\cdot\sin\left(\left\langle k_{j},x\right\rangle \right)
       \right)

    where:

        * :math:`\bar u` : mean velocity in :math:`e_1` direction
        * :math:`N` : fourier mode number
        * :math:`Z_{k,j}` : random samples from a normal distribution
        * :math:`k_j` : samples from the spectral density distribution of
          the covariance model
        * :math:`p_i(k_j) = e_1 - \frac{k_i k_1}{k^2}` : the projector
          ensuring the incompressibility
    """

    def __init__(
        self,
        model,
        mean_velocity=1.0,
        mode_no=1000,
        seed=None,
        verbose=False,
        sampling="auto",
        **kwargs
    ):
        if model.dim < 2 or model.dim > 3:
            raise ValueError(
                "Only 2D and 3D incompressible fields can be generated."
            )
        super().__init__(model, mode_no, seed, verbose, sampling, **kwargs)

        self.mean_u = mean_velocity
        self._value_type = "vector"

    def __call__(self, pos):
        """Calculate the random modes for the randomization method.

        This method  calls the `summate_incompr_*` Cython methods,
        which are the heart of the randomization method.
        In this class the method contains a projector to
        ensure the incompressibility of the vector field.

        Parameters
        ----------
        pos : (d, n), :class:`numpy.ndarray`
            the position tuple with d dimensions and n points.

        Returns
        -------
        :class:`numpy.ndarray`
            the random modes
        """
        pos = np.array(pos, dtype=np.double)
        summed_modes = summate_incompr(
            self._cov_sample, self._z_1, self._z_2, pos
        )
        nugget = self.get_nugget(summed_modes.shape)

        e1 = self._create_unit_vector(summed_modes.shape)

        return (
            self.mean_u * e1
            + self.mean_u
            * np.sqrt(self.model.var / self._mode_no)
            * summed_modes
            + nugget
        )

    def _create_unit_vector(self, broadcast_shape, axis=0):
        """Create a unit vector.

        Can be multiplied with a vector of shape broadcast_shape

        Parameters
        ----------
        broadcast_shape : :class:`tuple`
            the shape of the array with which
            the unit vector is to be multiplied
        axis : :class:`int`, optional
            the direction of the unit vector. Default: ``0``

        Returns
        -------
        :class:`numpy.ndarray`
            the unit vector
        """
        shape = np.ones(len(broadcast_shape), dtype=int)
        shape[0] = self.model.dim

        e1 = np.zeros(shape)
        e1[axis] = 1.0
        return e1

1	# -- coding: utf-8 --
2	"""	15✔
3	GStools subpackage providing generators for spatial random fields.
4
5	.. currentmodule:: gstools.field.generator
6
7	The following classes are provided
8
9	.. autosummary::
10	RandMeth
11	IncomprRandMeth
12	"""
13	# pylint: disable=C0103
14
15	from copy import deepcopy as dcp	15✔
16	import numpy as np	15✔
17	from gstools.covmodel.base import CovModel	15✔
18	from gstools.random.rng import RNG	15✔
19	from gstools.field.summator import summate, summate_incompr	15✔
20
21
22	__all__ = ["RandMeth", "IncomprRandMeth"]	15✔
23
24
25	SAMPLING = ["auto", "inversion", "mcmc"]	15✔
26
27
28	class RandMeth:	15✔
29	r"""Randomization method for calculating isotropic spatial random fields.
30
31	Parameters
32	----------
33	model : :any:`CovModel`
34	Covariance model
35	mode_no : :class:`int`, optional
36	Number of Fourier modes. Default: ``1000``
37	seed : :class:`int` or :any:`None`, optional
38	The seed of the random number generator.
39	If "None", a random seed is used. Default: :any:`None`
40	verbose : :class:`bool`, optional
41	Be chatty during the generation.
42	Default: :any:`False`
43	sampling :class:`str`, optional
44	Sampling strategy. Either
45
46	* "auto": select best strategy depending on given model
47	* "inversion": use inversion method
48	* "mcmc": use mcmc sampling
49
50	**kwargs
51	Placeholder for keyword-args
52
53	Notes
54	-----
55	The Randomization method is used to generate isotropic
56	spatial random fields characterized by a given covariance model.
57	The calculation looks like:
58
59	.. math::
60	u\left(x\right)=
61	\sqrt{\frac{\sigma^{2}}{N}}\cdot
62	\sum_{i=1}^{N}\left(
63	Z_{1,i}\cdot\cos\left(\left\langle k_{i},x\right\rangle \right)+
64	Z_{2,i}\cdot\sin\left(\left\langle k_{i},x\right\rangle \right)
65	\right)
66
67	where:
68
69	* :math:`N` : fourier mode number
70	* :math:`Z_{j,i}` : random samples from a normal distribution
71	* :math:`k_i` : samples from the spectral density distribution of
72	the covariance model
73	"""
74
75	def __init__(	15✔
76	self,
77	model,
78	mode_no=1000,
79	seed=None,
80	verbose=False,
81	sampling="auto",
82	**kwargs
83	):
84	if kwargs:	15✔
85	print("gstools.RandMeth: **kwargs are ignored")	×
86	# initialize atributes
87	self._mode_no = int(mode_no)	15✔
88	self._verbose = bool(verbose)	15✔
89	# initialize private atributes
90	self._model = None	15✔
91	self._seed = None	15✔
92	self._rng = None	15✔
93	self._z_1 = None	15✔
94	self._z_2 = None	15✔
95	self._cov_sample = None	15✔
96	self._value_type = "scalar"	15✔
97	# set sampling strategy
98	self._sampling = None	15✔
99	self.sampling = sampling	15✔
100	# set model and seed
101	self.update(model, seed)	15✔
102
103	def __call__(self, pos, add_nugget=True):	15✔
104	"""Calculate the random modes for the randomization method.
105
106	This method calls the `summate_*` Cython methods, which are the
107	heart of the randomization method.
108
109	Parameters
110	----------
111	pos : (d, n), :class:`numpy.ndarray`
112	the position tuple with d dimensions and n points.
113	add_nugget : :class:`bool`
114	Wheater to add nugget noise to the field.
115
116	Returns
117	-------
118	:class:`numpy.ndarray`
119	the random modes
120	"""
121	pos = np.array(pos, dtype=np.double)	15✔
122	summed_modes = summate(self._cov_sample, self._z_1, self._z_2, pos)	15✔
123	nugget = self.get_nugget(summed_modes.shape) if add_nugget else 0.0	15✔
124	return np.sqrt(self.model.var / self._mode_no) * summed_modes + nugget	15✔
125
126	def get_nugget(self, shape):	15✔
127	"""
128	Generate normal distributed values for the nugget simulation.
129
130	Parameters
131	----------
132	shape : :class:`tuple`
133	the shape of the summed modes
134
135	Returns
136	-------
137	nugget : :class:`numpy.ndarray`
138	the nugget in the same shape as the summed modes
139	"""
140	if self.model.nugget > 0:	15✔
141	nugget = np.sqrt(self.model.nugget) * self._rng.random.normal(	×
142	size=shape
143	)
144	else:
145	nugget = 0.0	15✔
146	return nugget	15✔
147
148	def update(self, model=None, seed=np.nan):	15✔
149	"""Update the model and the seed.
150
151	If model and seed are not different, nothing will be done.
152
153	Parameters
154	----------
155	model : :any:`CovModel` or :any:`None`, optional
156	covariance model. Default: :any:`None`
157	seed : :class:`int` or :any:`None` or :any:`numpy.nan`, optional
158	the seed of the random number generator.
159	If :any:`None`, a random seed is used. If :any:`numpy.nan`,
160	the actual seed will be kept. Default: :any:`numpy.nan`
161	"""
162	# check if a new model is given
163	if isinstance(model, CovModel):	15✔
164	if self.model != model:	15✔
165	self._model = dcp(model)	15✔
166	if seed is None or not np.isnan(seed):	15✔
167	self.reset_seed(seed)	15✔
168	else:
169	self.reset_seed(self._seed)	15✔
170	# just update the seed, if its a new one
171	elif seed is None or not np.isnan(seed):	15✔
172	self.seed = seed	15✔
173	# or just update the seed, when no model is given
174	elif model is None and (seed is None or not np.isnan(seed)):	×
175	if isinstance(self._model, CovModel):	×
176	self.seed = seed	×
177	else:
178	raise ValueError(	×
179	"gstools.field.generator.RandMeth: no 'model' given"
180	)
181	# if the user tries to trick us, we beat him!
182	elif model is None and np.isnan(seed):	×
183	if (	×
184	isinstance(self._model, CovModel)
185	and self._z_1 is not None
186	and self._z_2 is not None
187	and self._cov_sample is not None
188	):
189	if self.verbose:	×
190	print("RandMeth.update: Nothing will be done...")	×
191	else:
192	raise ValueError(	×
193	"gstools.field.generator.RandMeth: "
194	+ "neither 'model' nor 'seed' given!"
195	)
196	# wrong model type
197	else:
198	raise ValueError(	×
199	"gstools.field.generator.RandMeth: 'model' is not an "
200	+ "instance of 'gstools.CovModel'"
201	)
202
203	def reset_seed(self, seed=np.nan):	15✔
204	"""
205	Recalculate the random amplitudes and wave numbers with the given seed.
206
207	Parameters
208	----------
209	seed : :class:`int` or :any:`None` or :any:`numpy.nan`, optional
210	the seed of the random number generator.
211	If :any:`None`, a random seed is used. If :any:`numpy.nan`,
212	the actual seed will be kept. Default: :any:`numpy.nan`
213
214	Notes
215	-----
216	Even if the given seed is the present one, modes will be recalculated.
217	"""
218	if seed is None or not np.isnan(seed):	15✔
219	self._seed = seed	15✔
220	self._rng = RNG(self._seed)	15✔
221	# normal distributed samples for randmeth
222	self._z_1 = self._rng.random.normal(size=self._mode_no)	15✔
223	self._z_2 = self._rng.random.normal(size=self._mode_no)	15✔
224	# sample uniform on a sphere
225	sphere_coord = self._rng.sample_sphere(self.model.dim, self._mode_no)	15✔
226	# sample radii acording to radial spectral density of the model
227	if self.sampling == "inversion" or (	15✔
228	self.sampling == "auto" and self.model.has_ppf
229	):
230	pdf, cdf, ppf = self.model.dist_func	15✔
231	rad = self._rng.sample_dist(	15✔
232	size=self._mode_no, pdf=pdf, cdf=cdf, ppf=ppf, a=0
233	)
234	else:
235	rad = self._rng.sample_ln_pdf(	15✔
236	ln_pdf=self.model.ln_spectral_rad_pdf,
237	size=self._mode_no,
238	sample_around=1.0 / self.model.len_rescaled,
239	)
240	# get fully spatial samples by multiplying sphere samples and radii
241	self._cov_sample = rad * sphere_coord	15✔
242
243	@property	15✔
244	def sampling(self):	6✔
245	""":class:`str`: Sampling strategy."""
246	return self._sampling	15✔
247
248	@sampling.setter	15✔
249	def sampling(self, sampling):	6✔
250	if sampling not in ["auto", "inversion", "mcmc"]:	15✔
251	raise ValueError("RandMeth: sampling not in {}.".format(SAMPLING))	×
252	self._sampling = sampling	15✔
253
254	@property	15✔
255	def seed(self):	6✔
256	""":class:`int`: Seed of the master RNG.
257
258	Notes
259	-----
260	If a new seed is given, the setter property not only saves the
261	new seed, but also creates new random modes with the new seed.
262	"""
263	return self._seed	15✔
264
265	@seed.setter	15✔
266	def seed(self, new_seed):	6✔
267	if new_seed is not self._seed:	15✔
268	self.reset_seed(new_seed)	15✔
269
270	@property	15✔
271	def model(self):	6✔
272	""":any:`CovModel`: Covariance model of the spatial random field."""
273	return self._model	15✔
274
275	@model.setter	15✔
276	def model(self, model):	6✔
277	self.update(model)	15✔
278
279	@property	15✔
280	def mode_no(self):	6✔
281	""":class:`int`: Number of modes in the randomization method."""
282	return self._mode_no	×
283
284	@mode_no.setter	15✔
285	def mode_no(self, mode_no):	6✔
286	if int(mode_no) != self._mode_no:	15✔
287	self._mode_no = int(mode_no)	15✔
288	self.reset_seed(self._seed)	15✔
289
290	@property	15✔
291	def verbose(self):	6✔
292	""":class:`bool`: Verbosity of the generator."""
293	return self._verbose	×
294
295	@verbose.setter	15✔
296	def verbose(self, verbose):	6✔
297	self._verbose = bool(verbose)	×
298
299	@property	15✔
300	def name(self):	6✔
301	""":class:`str`: Name of the generator."""
302	return self.__class__.__name__	×
303
304	@property	15✔
305	def value_type(self):	6✔
306	""":class:`str`: Type of the field values (scalar, vector)."""
307	return self._value_type	15✔
308
309	def __repr__(self):
310	"""Return String representation."""
311	return "RandMeth(model={0}, mode_no={1}, seed={2})".format(
312	repr(self.model), self._mode_no, self.seed
313	)
314
315
316	class IncomprRandMeth(RandMeth):	15✔
317	r"""RandMeth for incompressible random vector fields.
318
319	Parameters
320	----------
321	model : :any:`CovModel`
322	covariance model
323	mean_velocity : :class:`float`, optional
324	the mean velocity in x-direction
325	mode_no : :class:`int`, optional
326	number of Fourier modes. Default: ``1000``
327	seed : :class:`int` or :any:`None`, optional
328	the seed of the random number generator.
329	If "None", a random seed is used. Default: :any:`None`
330	verbose : :class:`bool`, optional
331	State if there should be output during the generation.
332	Default: :any:`False`
333	sampling :class:`str`, optional
334	Sampling strategy. Either
335
336	* "auto": select best strategy depending on given model
337	* "inversion": use inversion method
338	* "mcmc": use mcmc sampling
339
340	**kwargs
341	Placeholder for keyword-args
342
343	Notes
344	-----
345	The Randomization method is used to generate isotropic
346	spatial incompressible random vector fields characterized
347	by a given covariance model. The equation is:
348
349	.. math::
350	u_i\left(x\right)= \bar{u_i} \delta_{i1} +
351	\bar{u_i}\sqrt{\frac{\sigma^{2}}{N}}\cdot
352	\sum_{j=1}^{N}p_i(k_{j})\left(
353	Z_{1,j}\cdot\cos\left(\left\langle k_{j},x\right\rangle \right)+
354	Z_{2,j}\cdot\sin\left(\left\langle k_{j},x\right\rangle \right)
355	\right)
356
357	where:
358
359	* :math:`\bar u` : mean velocity in :math:`e_1` direction
360	* :math:`N` : fourier mode number
361	* :math:`Z_{k,j}` : random samples from a normal distribution
362	* :math:`k_j` : samples from the spectral density distribution of
363	the covariance model
364	* :math:`p_i(k_j) = e_1 - \frac{k_i k_1}{k^2}` : the projector
365	ensuring the incompressibility
366	"""
367
368	def __init__(	15✔
369	self,
370	model,
371	mean_velocity=1.0,
372	mode_no=1000,
373	seed=None,
374	verbose=False,
375	sampling="auto",
376	**kwargs
377	):
378	if model.dim < 2 or model.dim > 3:	15✔
379	raise ValueError(	15✔
380	"Only 2D and 3D incompressible fields can be generated."
381	)
382	super().__init__(model, mode_no, seed, verbose, sampling, **kwargs)	15✔
383
384	self.mean_u = mean_velocity	15✔
385	self._value_type = "vector"	15✔
386
387	def __call__(self, pos):	15✔
388	"""Calculate the random modes for the randomization method.
389
390	This method calls the `summate_incompr_*` Cython methods,
391	which are the heart of the randomization method.
392	In this class the method contains a projector to
393	ensure the incompressibility of the vector field.
394
395	Parameters
396	----------
397	pos : (d, n), :class:`numpy.ndarray`
398	the position tuple with d dimensions and n points.
399
400	Returns
401	-------
402	:class:`numpy.ndarray`
403	the random modes
404	"""
405	pos = np.array(pos, dtype=np.double)	15✔
406	summed_modes = summate_incompr(	15✔
407	self._cov_sample, self._z_1, self._z_2, pos
408	)
409	nugget = self.get_nugget(summed_modes.shape)	15✔
410
411	e1 = self._create_unit_vector(summed_modes.shape)	15✔
412
413	return (	15✔
414	self.mean_u * e1
415	+ self.mean_u
416	* np.sqrt(self.model.var / self._mode_no)
417	* summed_modes
418	+ nugget
419	)
420
421	def _create_unit_vector(self, broadcast_shape, axis=0):	15✔
422	"""Create a unit vector.
423
424	Can be multiplied with a vector of shape broadcast_shape
425
426	Parameters
427	----------
428	broadcast_shape : :class:`tuple`
429	the shape of the array with which
430	the unit vector is to be multiplied
431	axis : :class:`int`, optional
432	the direction of the unit vector. Default: ``0``
433
434	Returns
435	-------
436	:class:`numpy.ndarray`
437	the unit vector
438	"""
439	shape = np.ones(len(broadcast_shape), dtype=int)	15✔
440	shape[0] = self.model.dim	15✔
441
442	e1 = np.zeros(shape)	15✔
443	e1[axis] = 1.0	15✔
444	return e1	15✔

GeoStat-Framework / GSTools / 483994410

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous