From 50e51cd2a32c21f1b29d89f4e1fae5f2fe9f6709 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Tue, 10 Nov 2020 21:26:53 +0100 Subject: [PATCH 01/69] CovModel: add 'rescale' argument to rescale the len_scale --- gstools/covmodel/base.py | 56 +++++++++++++++++++++------------------- 1 file changed, 29 insertions(+), 27 deletions(-) diff --git a/gstools/covmodel/base.py b/gstools/covmodel/base.py index b8dacd51..551484ff 100644 --- a/gstools/covmodel/base.py +++ b/gstools/covmodel/base.py @@ -82,7 +82,13 @@ class CovModel(metaclass=InitSubclassMeta): integral_scale : :class:`float` or :class:`list` or :any:`None`, optional If given, ``len_scale`` will be ignored and recalculated, so that the integral scale of the model matches the given one. - Default: ``None`` + Default: :any:`None` + rescale : :class:`float` or :any:`None`, optional + Optional rescaling factor to divide the length scale with. + This could be used for unit convertion or rescaling the length scale + to coincide with e.g. the integral scale. + Will be set by each model individually. + Default: :any:`None` var_raw : :class:`float` or :any:`None`, optional raw variance of the model which will be multiplied with :any:`CovModel.var_factor` to result in the actual variance. @@ -95,18 +101,6 @@ class CovModel(metaclass=InitSubclassMeta): used for the spectrum calculation. Use with caution (Better: Don't!). ``None`` is equivalent to ``{"a": -1, "b": 1, "N": 1000, "h": 0.001}``. Default: :any:`None` - - Examples - -------- - >>> from gstools import CovModel - >>> import numpy as np - >>> class Gau(CovModel): - ... def cor(self, h): - ... return np.exp(-h**2) - ... - >>> model = Gau() - >>> model.spectrum(2) - 0.00825830126008459 """ def __init__( @@ -118,6 +112,7 @@ def __init__( anis=1.0, angles=0.0, integral_scale=None, + rescale=None, var_raw=None, hankel_kw=None, **opt_arg @@ -182,6 +177,9 @@ def __init__( self.set_arg_bounds(check_args=False, **bounds) # set parameters + self._rescale = ( + self.default_rescale() if rescale is None else abs(float(rescale)) + ) self._nugget = nugget self._angles = set_angles(self.dim, angles) self._len_scale, self._anis = set_len_anis(self.dim, len_scale, anis) @@ -265,7 +263,7 @@ def correlation_from_cor(self, r): :math:`\rho(0)=1` and :math:`\rho(\infty)=0`. """ r = np.array(np.abs(r), dtype=np.double) - return self.cor(r / self.len_scale) + return self.cor(r * self.rescale / self.len_scale) def cor_from_correlation(self, h): r"""Normalziled correlation function taking a normalized range. @@ -273,11 +271,12 @@ def cor_from_correlation(self, h): Given by: :math:`\mathrm{cor}\left(r/\ell\right) = \rho(r)` """ h = np.array(np.abs(h), dtype=np.double) - return self.correlation(h * self.len_scale) + return self.correlation(h * self.len_scale / self.rescale) abstract = True if hasattr(cls, "cor"): - cls.correlation = correlation_from_cor + if not hasattr(cls, "correlation"): + cls.correlation = correlation_from_cor abstract = False else: cls.cor = cor_from_correlation @@ -298,8 +297,8 @@ def cor_from_correlation(self, h): "Can't instantiate class '" + cls.__name__ + "', " - + "without overriding at least on of the methods " - + "'variogram', 'covariance' or 'correlation'." + + "without providing at least one of the methods " + + "'cor', 'variogram', 'covariance' or 'correlation'." ) # modify the docstrings @@ -307,11 +306,11 @@ def cor_from_correlation(self, h): # class docstring gets attributes added if cls.__doc__ is None: cls.__doc__ = ( - "User defined GSTools Covariance-Model " - + CovModel.__doc__[44:-296] + "User defined GSTools Covariance-Model." + + CovModel.__doc__[45:] ) else: - cls.__doc__ += CovModel.__doc__[44:-296] + cls.__doc__ += CovModel.__doc__[45:] # overridden functions get standard doc if no new doc was created ignore = ["__", "variogram", "covariance", "correlation"] for attr in cls.__dict__: @@ -533,6 +532,10 @@ def var_factor(self): """Factor for the variance.""" return 1.0 + def default_rescale(self): + """Provide default rescaling factor.""" + return 1.0 + # calculation of different scales def calc_integral_scale(self): @@ -1054,6 +1057,11 @@ def len_scale(self, len_scale): ) self.check_arg_bounds() + @property + def rescale(self): + """:class:`float`: Rescale factor for the length scale of the model.""" + return self._rescale + @property def anis(self): """:class:`numpy.ndarray`: The anisotropy factors of the model.""" @@ -1263,9 +1271,3 @@ def __repr__(self): + opt_str + ")" ) - - -if __name__ == "__main__": - import doctest - - doctest.testmod() From f9ec092d052cd864baa250809bb4dc968d2b6e48 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Tue, 10 Nov 2020 21:27:26 +0100 Subject: [PATCH 02/69] Cleanup: remove doctest calls --- gstools/field/base.py | 6 ------ gstools/field/generator.py | 6 ------ gstools/field/srf.py | 6 ------ gstools/krige/base.py | 6 ------ gstools/krige/methods.py | 6 ------ gstools/random/rng.py | 6 ------ gstools/random/tools.py | 6 ------ 7 files changed, 42 deletions(-) diff --git a/gstools/field/base.py b/gstools/field/base.py index 2606a9c7..b560e947 100755 --- a/gstools/field/base.py +++ b/gstools/field/base.py @@ -394,9 +394,3 @@ def __str__(self): def __repr__(self): """Return String representation.""" return "Field(model={0}, mean={1})".format(self.model, self.mean) - - -if __name__ == "__main__": - import doctest - - doctest.testmod() diff --git a/gstools/field/generator.py b/gstools/field/generator.py index 67ed7ae9..a5a9c2ab 100644 --- a/gstools/field/generator.py +++ b/gstools/field/generator.py @@ -500,9 +500,3 @@ def _set_dtype(x, y=None, z=None, dtype=np.double): if z is not None: z = z.astype(dtype, copy=False) return x, y, z - - -if __name__ == "__main__": - import doctest - - doctest.testmod() diff --git a/gstools/field/srf.py b/gstools/field/srf.py index bc6052fa..6fee323b 100644 --- a/gstools/field/srf.py +++ b/gstools/field/srf.py @@ -286,9 +286,3 @@ def __repr__(self): return "SRF(model={0}, mean={1}, generator={2}".format( self.model, self.mean, self.generator ) - - -if __name__ == "__main__": - import doctest - - doctest.testmod() diff --git a/gstools/krige/base.py b/gstools/krige/base.py index 9c3560bd..57cd4121 100755 --- a/gstools/krige/base.py +++ b/gstools/krige/base.py @@ -640,9 +640,3 @@ def __repr__(self): ) + ")" ) - - -if __name__ == "__main__": - import doctest - - doctest.testmod() diff --git a/gstools/krige/methods.py b/gstools/krige/methods.py index c139ef89..616f3d90 100644 --- a/gstools/krige/methods.py +++ b/gstools/krige/methods.py @@ -412,9 +412,3 @@ def __init__( pseudo_inv=pseudo_inv, pseudo_inv_type=pseudo_inv_type, ) - - -if __name__ == "__main__": - import doctest - - doctest.testmod() diff --git a/gstools/random/rng.py b/gstools/random/rng.py index 8f6fd578..f49fed62 100644 --- a/gstools/random/rng.py +++ b/gstools/random/rng.py @@ -221,9 +221,3 @@ def __str__(self): def __repr__(self): """Return String representation.""" return "RNG(seed={})".format(self.seed) - - -if __name__ == "__main__": - import doctest - - doctest.testmod() diff --git a/gstools/random/tools.py b/gstools/random/tools.py index 5a4b0d9d..2b577096 100644 --- a/gstools/random/tools.py +++ b/gstools/random/tools.py @@ -186,9 +186,3 @@ def _cdf(self, x, *args): def _ppf(self, q, *args): return self.ppf_in(q) - - -if __name__ == "__main__": - import doctest - - doctest.testmod() From 333321b3f406b0e6c7810b5cc77faa33c1758706 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Wed, 11 Nov 2020 00:01:46 +0100 Subject: [PATCH 03/69] covmodel: add attribute 'len_rescaled' which is the rescaled len_scale --- gstools/covmodel/base.py | 72 +++++++++------------------------------- 1 file changed, 16 insertions(+), 56 deletions(-) diff --git a/gstools/covmodel/base.py b/gstools/covmodel/base.py index 551484ff..ceb65627 100644 --- a/gstools/covmodel/base.py +++ b/gstools/covmodel/base.py @@ -204,74 +204,29 @@ def __init__( # one of these functions needs to be overridden def __init_subclass__(cls): - r"""Initialize gstools covariance model. - - Warnings - -------- - Don't instantiate ``CovModel`` directly. You need to inherit a - child class which overrides one of the following methods: - - * ``model.variogram(r)`` - :math:`\gamma\left(r\right)= - \sigma^2\cdot\left(1-\rho\left(r\right)\right)+n` - * ``model.covariance(r)`` - :math:`C\left(r\right)= - \sigma^2\cdot\rho\left(r\right)` - * ``model.correlation(r)`` - :math:`\rho\left(r\right)` - - Best practice is to use the ``correlation`` function, or the ``cor`` - function. The latter one takes the dimensionles distance h=r/l. - """ + """Initialize gstools covariance model.""" def variogram(self, r): - r"""Isotropic variogram of the model. - - Given by: :math:`\gamma\left(r\right)= - \sigma^2\cdot\left(1-\rho\left(r\right)\right)+n` - - Where :math:`\rho(r)` is the correlation function. - """ + """Isotropic variogram of the model.""" return self.var - self.covariance(r) + self.nugget def covariance(self, r): - r"""Covariance of the model. - - Given by: :math:`C\left(r\right)= - \sigma^2\cdot\rho\left(r\right)` - - Where :math:`\rho(r)` is the correlation function. - """ + """Covariance of the model.""" return self.var * self.correlation(r) def correlation(self, r): - r"""Correlation function (or normalized covariance) of the model. - - Given by: :math:`\rho\left(r\right)` - - It has to be a monotonic decreasing function with - :math:`\rho(0)=1` and :math:`\rho(\infty)=0`. - """ + """Correlation function of the model.""" return 1.0 - (self.variogram(r) - self.nugget) / self.var def correlation_from_cor(self, r): - r"""Correlation function (or normalized covariance) of the model. - - Given by: :math:`\rho\left(r\right)` - - It has to be a monotonic decreasing function with - :math:`\rho(0)=1` and :math:`\rho(\infty)=0`. - """ + """Correlation function of the model.""" r = np.array(np.abs(r), dtype=np.double) - return self.cor(r * self.rescale / self.len_scale) + return self.cor(r / self.len_rescaled) def cor_from_correlation(self, h): - r"""Normalziled correlation function taking a normalized range. - - Given by: :math:`\mathrm{cor}\left(r/\ell\right) = \rho(r)` - """ + """Correlation taking a non-dimensional range.""" h = np.array(np.abs(h), dtype=np.double) - return self.correlation(h * self.len_scale / self.rescale) + return self.correlation(h * self.len_rescaled) abstract = True if hasattr(cls, "cor"): @@ -312,7 +267,7 @@ def cor_from_correlation(self, h): else: cls.__doc__ += CovModel.__doc__[45:] # overridden functions get standard doc if no new doc was created - ignore = ["__", "variogram", "covariance", "correlation"] + ignore = ["__", "variogram", "covariance", "cor"] for attr in cls.__dict__: if any( [attr.startswith(ign) for ign in ignore] @@ -570,8 +525,8 @@ def spectrum(self, k): This is given by: - .. math:: S(k) = \left(\frac{1}{2\pi}\right)^n - \int C(r) e^{i b\mathbf{k}\cdot\mathbf{r}} d^n\mathbf{r} + .. math:: S(\mathbf{k}) = \left(\frac{1}{2\pi}\right)^n + \int C(r) e^{i \mathbf{k}\cdot\mathbf{r}} d^n\mathbf{r} Internally, this is calculated by the hankel transformation: @@ -1062,6 +1017,11 @@ def rescale(self): """:class:`float`: Rescale factor for the length scale of the model.""" return self._rescale + @property + def len_rescaled(self): + """:class:`float`: The rescaled main length scale of the model.""" + return self._len_scale / self._rescale + @property def anis(self): """:class:`numpy.ndarray`: The anisotropy factors of the model.""" From 474e71ad0546d3ee8ed0c3d3798fbac74bc90a20 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Wed, 11 Nov 2020 00:32:14 +0100 Subject: [PATCH 04/69] CovModel: rewrite Gaussian covmodel to use 'cor' and 'rescale' --- gstools/covmodel/models.py | 53 +++++++++++++++++--------------------- 1 file changed, 23 insertions(+), 30 deletions(-) diff --git a/gstools/covmodel/models.py b/gstools/covmodel/models.py index 8cce6919..b665a1ec 100644 --- a/gstools/covmodel/models.py +++ b/gstools/covmodel/models.py @@ -45,47 +45,46 @@ class Gaussian(CovModel): Notes ----- - This model is given by the following correlation function: + This model is given by the following variogram: .. math:: - \rho(r) = - \exp\left(- \frac{\pi}{4} \cdot \left(\frac{r}{\ell}\right)^2\right) + \gamma(r)= + \sigma^{2} + \left(1-\exp\left(-\left(s\cdot\frac{r}{\ell}\right)^{2}\right)\right)+n + Where the standard rescale factor is :math:`s=\frac{\sqrt{\pi}}{2}`. """ - def correlation(self, r): - r"""Gaussian correlation function. + def cor(self, h): + """Gaussian normalized correlation function.""" + return np.exp(-(h ** 2)) - .. math:: - \rho(r) = - \exp\left(- \frac{\pi}{4}\cdot \left(\frac{r}{\ell}\right)^2\right) - """ - r = np.array(np.abs(r), dtype=np.double) - return np.exp(-np.pi / 4 * (r / self.len_scale) ** 2) + def default_rescale(self): + """Gaussian rescaling factor to result in integral scale.""" + return np.sqrt(np.pi) / 2.0 def spectral_density(self, k): # noqa: D102 k = np.array(k, dtype=np.double) - return (self.len_scale / np.pi) ** self.dim * np.exp( - -((k * self.len_scale) ** 2) / np.pi + return (self.len_rescaled / 2.0 / np.sqrt(np.pi)) ** self.dim * np.exp( + -((k * self.len_rescaled / 2.0) ** 2) ) def spectral_rad_cdf(self, r): - """Radial spectral cdf.""" + """Gaussian radial spectral cdf.""" r = np.array(r, dtype=np.double) if self.dim == 1: - return sps.erf(self.len_scale * r / np.sqrt(np.pi)) + return sps.erf(r * self.len_rescaled / 2.0) if self.dim == 2: - return 1.0 - np.exp(-((r * self.len_scale) ** 2) / np.pi) + return 1.0 - np.exp(-((r * self.len_rescaled / 2.0) ** 2)) if self.dim == 3: return sps.erf( - self.len_scale * r / np.sqrt(np.pi) - ) - 2 * r * self.len_scale / np.pi * np.exp( - -((r * self.len_scale) ** 2) / np.pi + r * self.len_rescaled / 2.0 + ) - r * self.len_rescaled / np.sqrt(np.pi) * np.exp( + -((r * self.len_rescaled / 2.0) ** 2) ) - return None def spectral_rad_ppf(self, u): - """Radial spectral ppf. + """Gaussian radial spectral ppf. Notes ----- @@ -93,10 +92,9 @@ def spectral_rad_ppf(self, u): """ u = np.array(u, dtype=np.double) if self.dim == 1: - return sps.erfinv(u) * np.sqrt(np.pi) / self.len_scale + return 2.0 / self.len_rescaled * sps.erfinv(u) if self.dim == 2: - return np.sqrt(np.pi) / self.len_scale * np.sqrt(-np.log(1.0 - u)) - return None + return 2.0 / self.len_rescaled * np.sqrt(-np.log(1.0 - u)) def _has_cdf(self): return self.dim in [1, 2, 3] @@ -105,7 +103,7 @@ def _has_ppf(self): return self.dim in [1, 2] def calc_integral_scale(self): # noqa: D102 - return self.len_scale + return self.len_rescaled * np.sqrt(np.pi) / 2.0 # Exponential Model ########################################################### @@ -121,7 +119,6 @@ class Exponential(CovModel): .. math:: \rho(r) = \exp\left(- \frac{r}{\ell} \right) - """ def correlation(self, r): @@ -481,7 +478,6 @@ class Linear(CovModel): & r<\ell\\ 0 & r\geq\ell \end{cases} - """ def correlation(self, r): @@ -527,7 +523,6 @@ class Circular(CovModel): & r<\ell\\ 0 & r\geq\ell \end{cases} - """ def correlation(self, r): @@ -584,7 +579,6 @@ class Spherical(CovModel): & r<\ell\\ 0 & r\geq\ell \end{cases} - """ def correlation(self, r): @@ -659,7 +653,6 @@ class Intersection(CovModel): & r<\ell\\ 0 & r\geq\ell \end{cases} - """ def correlation(self, r): # noqa: D102 From 60c92aa1b73b42fe2ae37699a5febff9b0b4b238 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Wed, 11 Nov 2020 00:33:02 +0100 Subject: [PATCH 05/69] CovModel: move **opt_arg doc to main class --- gstools/covmodel/base.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/gstools/covmodel/base.py b/gstools/covmodel/base.py index ceb65627..a61b0b96 100644 --- a/gstools/covmodel/base.py +++ b/gstools/covmodel/base.py @@ -101,6 +101,9 @@ class CovModel(metaclass=InitSubclassMeta): used for the spectrum calculation. Use with caution (Better: Don't!). ``None`` is equivalent to ``{"a": -1, "b": 1, "N": 1000, "h": 0.001}``. Default: :any:`None` + **opt_arg + Optional arguments are covered by these keyword arguments. + If present, they are described in the section `Other Parameters`. """ def __init__( From 795b3d27a12fccf99415d02077f0bcb0e0fc964b Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Wed, 11 Nov 2020 00:42:50 +0100 Subject: [PATCH 06/69] CovModel: rewrite Exponential covmodel to use 'cor' and 'rescale' --- gstools/covmodel/models.py | 51 +++++++++++++++++++------------------- 1 file changed, 25 insertions(+), 26 deletions(-) diff --git a/gstools/covmodel/models.py b/gstools/covmodel/models.py index b665a1ec..57afe4bd 100644 --- a/gstools/covmodel/models.py +++ b/gstools/covmodel/models.py @@ -114,52 +114,51 @@ class Exponential(CovModel): Notes ----- - This model is given by the following correlation function: + This model is given by the following variogram: .. math:: - \rho(r) = - \exp\left(- \frac{r}{\ell} \right) - """ + \gamma(r)= + \sigma^{2} + \left(1-\exp\left(-s\cdot\frac{r}{\ell}\right)\right)+n - def correlation(self, r): - r"""Exponential correlation function. + Where the standard rescale factor is :math:`s=1`. + """ - .. math:: - \rho(r) = - \exp\left(- \frac{r}{\ell} \right) - """ - r = np.array(np.abs(r), dtype=np.double) - return np.exp(-1 * r / self.len_scale) + def cor(self, h): + """Exponential normalized correlation function.""" + return np.exp(-h) def spectral_density(self, k): # noqa: D102 k = np.array(k, dtype=np.double) return ( - self.len_scale ** self.dim - * sps.gamma((self.dim + 1) / 2) - / (np.pi * (1.0 + (k * self.len_scale) ** 2)) - ** ((self.dim + 1) / 2) + self.len_rescaled ** self.dim + * sps.gamma((self.dim + 1) / 2.0) + / (np.pi * (1.0 + (k * self.len_rescaled) ** 2)) + ** ((self.dim + 1) / 2.0) ) def spectral_rad_cdf(self, r): - """Radial spectral cdf.""" + """Exponential radial spectral cdf.""" r = np.array(r, dtype=np.double) if self.dim == 1: - return np.arctan(r * self.len_scale) * 2 / np.pi + return np.arctan(r * self.len_rescaled) * 2.0 / np.pi if self.dim == 2: - return 1.0 - 1 / np.sqrt(1 + (r * self.len_scale) ** 2) + return 1.0 - 1.0 / np.sqrt(1.0 + (r * self.len_rescaled) ** 2) if self.dim == 3: return ( ( - np.arctan(r * self.len_scale) - - r * self.len_scale / (1 + (r * self.len_scale) ** 2) + np.arctan(r * self.len_rescaled) + - r + * self.len_rescaled + / (1.0 + (r * self.len_rescaled) ** 2) ) - * 2 + * 2.0 / np.pi ) return None def spectral_rad_ppf(self, u): - """Radial spectral ppf. + """Exponential radial spectral ppf. Notes ----- @@ -167,7 +166,7 @@ def spectral_rad_ppf(self, u): """ u = np.array(u, dtype=np.double) if self.dim == 1: - return np.tan(np.pi / 2 * u) / self.len_scale + return np.tan(np.pi / 2 * u) / self.len_rescaled if self.dim == 2: u_power = np.divide( 1, @@ -175,7 +174,7 @@ def spectral_rad_ppf(self, u): out=np.full_like(u, np.inf), where=np.logical_not(np.isclose(u, 0)), ) - return np.sqrt(u_power - 1.0) / self.len_scale + return np.sqrt(u_power - 1.0) / self.len_rescaled return None def _has_cdf(self): @@ -185,7 +184,7 @@ def _has_ppf(self): return self.dim in [1, 2] def calc_integral_scale(self): # noqa: D102 - return self.len_scale + return self.len_rescaled # Rational Model ############################################################## From d9daac64aceae0014422c4ea62c3d25ff940b0e6 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Wed, 11 Nov 2020 00:51:55 +0100 Subject: [PATCH 07/69] CovModel: rewrite Rational and Stable covmodel to use 'cor' and 'rescale' --- gstools/covmodel/models.py | 41 +++++++++++--------------------------- 1 file changed, 12 insertions(+), 29 deletions(-) diff --git a/gstools/covmodel/models.py b/gstools/covmodel/models.py index 57afe4bd..ca693e21 100644 --- a/gstools/covmodel/models.py +++ b/gstools/covmodel/models.py @@ -200,14 +200,13 @@ class Rational(CovModel): .. math:: \rho(r) = \left(1 + \frac{1}{2\alpha} \cdot - \left(\frac{r}{\ell}\right)^2\right)^{-\alpha} + \left(s\cdot\frac{r}{\ell}\right)^2\right)^{-\alpha} + Where the standard rescale factor is :math:`s=1`. :math:`\alpha` is a shape parameter and should be > 0.5. Other Parameters ---------------- - **opt_arg - The following parameters are covered by these keyword arguments alpha : :class:`float`, optional Shape parameter. Standard range: ``(0, inf)`` Default: ``1.0`` @@ -237,22 +236,13 @@ def default_opt_arg_bounds(self): """ return {"alpha": [0.5, np.inf]} - def correlation(self, r): - r"""Rational correlation function. - - .. math:: - \rho(r) = - \left(1 + \frac{1}{2\alpha} \cdot - \left(\frac{r}{\ell}\right)^2\right)^{-\alpha} - """ - r = np.array(np.abs(r), dtype=np.double) - return np.power( - 1 + 0.5 / self.alpha * (r / self.len_scale) ** 2, -self.alpha - ) + def cor(self, h): + """Rational normalized correlation function.""" + return np.power(1 + 0.5 / self.alpha * h ** 2, -self.alpha) def calc_integral_scale(self): # noqa: D102 return ( - self.len_scale + self.len_rescaled * np.sqrt(np.pi * self.alpha * 0.5) * sps.gamma(self.alpha - 0.5) / sps.gamma(self.alpha) @@ -271,14 +261,13 @@ class Stable(CovModel): .. math:: \rho(r) = - \exp\left(- \left(\frac{r}{\ell}\right)^{\alpha}\right) + \exp\left(- \left(s\cdot\frac{r}{\ell}\right)^{\alpha}\right) + Where the standard rescale factor is :math:`s=1`. :math:`\alpha` is a shape parameter with :math:`\alpha\in(0,2]` Other Parameters ---------------- - **opt_arg - The following parameters are covered by these keyword arguments alpha : :class:`float`, optional Shape parameter. Standard range: ``(0, 2]`` Default: ``1.5`` @@ -323,18 +312,12 @@ def check_opt_arg(self): + "count with unstable results" ) - def correlation(self, r): - r"""Stable correlation function. - - .. math:: - \rho(r) = - \exp\left(- \left(\frac{r}{\ell}\right)^{\alpha}\right) - """ - r = np.array(np.abs(r), dtype=np.double) - return np.exp(-np.power(r / self.len_scale, self.alpha)) + def cor(self, h): + r"""Stable normalized correlation function.""" + return np.exp(-np.power(h, self.alpha)) def calc_integral_scale(self): # noqa: D102 - return self.len_scale * sps.gamma(1.0 + 1.0 / self.alpha) + return self.len_rescaled * sps.gamma(1.0 + 1.0 / self.alpha) # Matérn Model ################################################################ From b34c37228ada6dec6abbb6b06df4f51b5bfb9292 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Wed, 11 Nov 2020 01:02:53 +0100 Subject: [PATCH 08/69] CovModel: rewrite Matern covmodel to use 'cor' and 'rescale' --- gstools/covmodel/models.py | 52 +++++++++++++++++--------------------- 1 file changed, 23 insertions(+), 29 deletions(-) diff --git a/gstools/covmodel/models.py b/gstools/covmodel/models.py index ca693e21..6e6386ba 100644 --- a/gstools/covmodel/models.py +++ b/gstools/covmodel/models.py @@ -333,10 +333,11 @@ class Matern(CovModel): .. math:: \rho(r) = \frac{2^{1-\nu}}{\Gamma\left(\nu\right)} \cdot - \left(\sqrt{\nu}\cdot\frac{r}{\ell}\right)^{\nu} \cdot - \mathrm{K}_{\nu}\left(\sqrt{\nu}\cdot\frac{r}{\ell}\right) + \left(\sqrt{\nu}\cdot s\cdot\frac{r}{\ell}\right)^{\nu} \cdot + \mathrm{K}_{\nu}\left(\sqrt{\nu}\cdot s\cdot\frac{r}{\ell}\right) - Where :math:`\Gamma` is the gamma function and :math:`\mathrm{K}_{\nu}` + Where the standard rescale factor is :math:`s=1`. + :math:`\Gamma` is the gamma function and :math:`\mathrm{K}_{\nu}` is the modified Bessel function of the second kind. :math:`\nu` is a shape parameter and should be >= 0.2. @@ -346,12 +347,10 @@ class Matern(CovModel): .. math:: \rho(r) = - \exp\left(- \frac{1}{4} \cdot \left(\frac{r}{\ell}\right)^2\right) + \exp\left(-\left(\frac{r}{2\ell}\right)^2\right) Other Parameters ---------------- - **opt_arg - The following parameters are covered by these keyword arguments nu : :class:`float`, optional Shape parameter. Standard range: ``[0.2, 30]`` Default: ``1.0`` @@ -381,28 +380,20 @@ def default_opt_arg_bounds(self): """ return {"nu": [0.2, 30.0, "cc"]} - def correlation(self, r): - r"""Matérn correlation function. - - .. math:: - \rho(r) = - \frac{2^{1-\nu}}{\Gamma\left(\nu\right)} \cdot - \left(\sqrt{\nu}\cdot\frac{r}{\ell}\right)^{\nu} \cdot - \mathrm{K}_{\nu}\left(\sqrt{\nu}\cdot\frac{r}{\ell}\right) - """ - r = np.array(np.abs(r), dtype=np.double) + def cor(self, h): + """Matérn normalized correlation function.""" + h = np.array(np.abs(h), dtype=np.double) # for nu > 20 we just use the gaussian model if self.nu > 20.0: - return np.exp(-((r / self.len_scale) ** 2) / 4) + return np.exp(-((h / 2.0) ** 2)) # calculate by log-transformation to prevent numerical errors - r_gz = r[r > 0.0] - res = np.ones_like(r) - # with np.errstate(over="ignore", invalid="ignore"): - res[r > 0.0] = np.exp( + h_gz = h[h > 0.0] + res = np.ones_like(h) + res[h > 0.0] = np.exp( (1.0 - self.nu) * np.log(2) - sps.loggamma(self.nu) - + self.nu * np.log(np.sqrt(self.nu) * r_gz / self.len_scale) - ) * sps.kv(self.nu, np.sqrt(self.nu) * r_gz / self.len_scale) + + self.nu * np.log(np.sqrt(self.nu) * h_gz) + ) * sps.kv(self.nu, np.sqrt(self.nu) * h_gz) # if nu >> 1 we get errors for the farfield, there 0 is approached res[np.logical_not(np.isfinite(res))] = 0.0 # covariance is positiv @@ -414,20 +405,20 @@ def spectral_density(self, k): # noqa: D102 # for nu > 20 we just use an approximation of the gaussian model if self.nu > 20.0: return ( - (self.len_scale / np.sqrt(np.pi)) ** self.dim - * np.exp(-((k * self.len_scale) ** 2)) + (self.len_rescaled / np.sqrt(np.pi)) ** self.dim + * np.exp(-((k * self.len_rescaled) ** 2)) * ( 1 + ( - ((k * self.len_scale) ** 2 - self.dim / 2.0) ** 2 + ((k * self.len_rescaled) ** 2 - self.dim / 2.0) ** 2 - self.dim / 2.0 ) / self.nu ) ) - return (self.len_scale / np.sqrt(np.pi)) ** self.dim * np.exp( + return (self.len_rescaled / np.sqrt(np.pi)) ** self.dim * np.exp( -(self.nu + self.dim / 2.0) - * np.log(1.0 + (k * self.len_scale) ** 2 / self.nu) + * np.log(1.0 + (k * self.len_rescaled) ** 2 / self.nu) + sps.loggamma(self.nu + self.dim / 2.0) - sps.loggamma(self.nu) - self.dim * np.log(np.sqrt(self.nu)) @@ -435,7 +426,10 @@ def spectral_density(self, k): # noqa: D102 def calc_integral_scale(self): # noqa: D102 return ( - self.len_scale * np.pi / np.sqrt(self.nu) / sps.beta(self.nu, 0.5) + self.len_rescaled + * np.pi + / np.sqrt(self.nu) + / sps.beta(self.nu, 0.5) ) From f324d2a7c565418f13a828d078423d97b15c18cd Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Wed, 11 Nov 2020 01:27:16 +0100 Subject: [PATCH 09/69] CovModel: rewrite Linear, Circular and Spherical covmodel to use 'cor' and 'rescale' --- gstools/covmodel/models.py | 113 +++++++++++-------------------------- 1 file changed, 34 insertions(+), 79 deletions(-) diff --git a/gstools/covmodel/models.py b/gstools/covmodel/models.py index 6e6386ba..d350d032 100644 --- a/gstools/covmodel/models.py +++ b/gstools/covmodel/models.py @@ -450,29 +450,16 @@ class Linear(CovModel): .. math:: \rho(r) = \begin{cases} - 1-\frac{r}{\ell} - & r<\ell\\ - 0 & r\geq\ell + 1-s\cdot\frac{r}{\ell} & r<\frac{\ell}{s}\\ + 0 & r\geq\frac{\ell}{s} \end{cases} - """ - def correlation(self, r): - r"""Linear correlation function. + Where the standard rescale factor is :math:`s=1`. + """ - .. math:: - \rho(r) = - \begin{cases} - 1-\frac{r}{\ell} - & r<\ell\\ - 0 & r\geq\ell - \end{cases} - """ - r = np.array(np.abs(r), dtype=np.double) - res = np.zeros_like(r) - r_ll = r < self.len_scale - r_low = r[r_ll] - res[r_ll] = 1.0 - r_low / self.len_scale - return res + def cor(self, h): + """Linear normalized correlation function.""" + return np.maximum(1 - np.abs(h, dtype=np.double), 0.0) # Circular Model ############################################################## @@ -492,43 +479,27 @@ class Circular(CovModel): .. math:: \rho(r) = \begin{cases} - \frac{2}{\pi}\cdot\left( - \cos^{-1}\left(\frac{r}{\ell}\right) - - \frac{r}{\ell}\cdot\sqrt{1-\left(\frac{r}{\ell}\right)^{2}} + \frac{2}{\pi}\cdot + \left( + \cos^{-1}\left(s\cdot\frac{r}{\ell}\right) - + s\cdot\frac{r}{\ell}\cdot\sqrt{1-\left(s\cdot\frac{r}{\ell}\right)^{2}} \right) - & r<\ell\\ - 0 & r\geq\ell + & r<\frac{\ell}{s}\\ + 0 & r\geq\frac{\ell}{s} \end{cases} - """ - - def correlation(self, r): - r"""Circular correlation function. - .. math:: - \rho(r) = - \begin{cases} - \frac{2}{\pi}\cdot\left( - \cos^{-1}\left(\frac{r}{\ell}\right) - - \frac{r}{\ell}\cdot\sqrt{1-\left(\frac{r}{\ell}\right)^{2}} - \right) - & r<\ell\\ - 0 & r\geq\ell - \end{cases} + Where the standard rescale factor is :math:`s=1`. + """ - """ - r = np.array(np.abs(r), dtype=np.double) - res = np.zeros_like(r) - r_ll = r < self.len_scale - r_low = r[r_ll] - res[r_ll] = ( - 2 - / np.pi - * ( - np.arccos(r_low / self.len_scale) - - r_low - / self.len_scale - * np.sqrt(1 - (r_low / self.len_scale) ** 2) - ) + def cor(self, h): + """Circular normalized correlation function.""" + h = np.array(np.abs(h), dtype=np.double) + res = np.zeros_like(h) + # arccos is instable around h=1 + h_l1 = h < 1.0 + h_low = h[h_l1] + res[h_l1] = ( + 2 / np.pi * (np.arccos(h_low) - h_low * np.sqrt(1 - h_low ** 2)) ) return res @@ -550,35 +521,19 @@ class Spherical(CovModel): .. math:: \rho(r) = \begin{cases} - 1-\frac{3}{2}\cdot\frac{r}{\ell} + - \frac{1}{2}\cdot\left(\frac{r}{\ell}\right)^{3} - & r<\ell\\ - 0 & r\geq\ell + 1-\frac{3}{2}\cdot s\cdot\frac{r}{\ell} + + \frac{1}{2}\cdot\left(s\cdot\frac{r}{\ell}\right)^{3} + & r<\frac{\ell}{s}\\ + 0 & r\geq\frac{\ell}{s} \end{cases} - """ - def correlation(self, r): - r"""Spherical correlation function. + Where the standard rescale factor is :math:`s=1`. + """ - .. math:: - \rho(r) = - \begin{cases} - 1-\frac{3}{2}\cdot\frac{r}{\ell} + - \frac{1}{2}\cdot\left(\frac{r}{\ell}\right)^{3} - & r<\ell\\ - 0 & r\geq\ell - \end{cases} - """ - r = np.array(np.abs(r), dtype=np.double) - res = np.zeros_like(r) - r_ll = r < self.len_scale - r_low = r[r_ll] - res[r_ll] = ( - 1.0 - - 3.0 / 2.0 * r_low / self.len_scale - + 1.0 / 2.0 * (r_low / self.len_scale) ** 3 - ) - return res + def cor(self, h): + """Spherical normalized correlation function.""" + h = np.minimum(np.abs(h, dtype=np.double), 1.0) + return 1.0 - 1.5 * h + 0.5 * h ** 3 class Intersection(CovModel): From 6722f638cf1574a7d5abcc7690515ee14b66bdac Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Wed, 11 Nov 2020 01:40:47 +0100 Subject: [PATCH 10/69] CovModel: replace the 'Intersection' model with the 'HyperSpherical' model, since it's a better implementation --- gstools/covmodel/models.py | 146 ++++++++----------------------------- 1 file changed, 31 insertions(+), 115 deletions(-) diff --git a/gstools/covmodel/models.py b/gstools/covmodel/models.py index d350d032..406c26b3 100644 --- a/gstools/covmodel/models.py +++ b/gstools/covmodel/models.py @@ -536,141 +536,57 @@ def cor(self, h): return 1.0 - 1.5 * h + 0.5 * h ** 3 -class Intersection(CovModel): - r"""The Intersection covariance model. +class HyperSpherical(CovModel): + r"""The Hyper-Spherical covariance model. This model is derived from the relative intersection area of - two d-dimensional spheres, + two d-dimensional hyper spheres, where the middle points have a distance of :math:`r` and the diameters are given by :math:`\ell`. - In 1D this is the Linear model, in 2D this is the Circular model - and in 3D this is the Spherical model. + In 1D this is the Linear model, in 2D the Circular model + and in 3D the Spherical model. Notes ----- This model is given by the following correlation functions. - In 1D: - - .. math:: - \rho(r) = - \begin{cases} - 1-\frac{r}{\ell} - & r<\ell\\ - 0 & r\geq\ell - \end{cases} - - In 2D: - .. math:: \rho(r) = \begin{cases} - \frac{2}{\pi}\cdot\left( - \cos^{-1}\left(\frac{r}{\ell}\right) - - \frac{r}{\ell}\cdot\sqrt{1-\left(\frac{r}{\ell}\right)^{2}} - \right) - & r<\ell\\ - 0 & r\geq\ell + 1-\frac{ + _{2}F_{1}\left(\frac{1}{2},\frac{1-d}{2},\frac{3}{2}, + \left(s\cdot\frac{r}{\ell}\right)^{2}\right)} + {_{2}F_{1}\left(\frac{1}{2},\frac{1-d}{2},\frac{3}{2},1\right)} + & r<\frac{\ell}{s}\\ + 0 & r\geq\frac{\ell}{s} \end{cases} - In >=3D: - - .. math:: - \rho(r) = - \begin{cases} - 1-\frac{3}{2}\cdot\frac{r}{\ell} + - \frac{1}{2}\cdot\left(\frac{r}{\ell}\right)^{3} - & r<\ell\\ - 0 & r\geq\ell - \end{cases} + Where the standard rescale factor is :math:`s=1`. + :math:`d` is the dimension. """ - def correlation(self, r): # noqa: D102 - r""" - Intersection correlation function. - - In 1D: - - .. math:: - \rho(r) = - \begin{cases} - 1-\frac{r}{\ell} - & r<\ell\\ - 0 & r\geq\ell - \end{cases} - - In 2D: - - .. math:: - \rho(r) = - \begin{cases} - \frac{2}{\pi}\cdot\left( - \cos^{-1}\left(\frac{r}{\ell}\right) - - \frac{r}{\ell}\cdot\sqrt{1-\left(\frac{r}{\ell}\right)^{2}} - \right) - & r<\ell\\ - 0 & r\geq\ell - \end{cases} - - In >=3D: - - .. math:: - \rho(r) = - \begin{cases} - 1-\frac{3}{2}\cdot\frac{r}{\ell} + - \frac{1}{2}\cdot\left(\frac{r}{\ell}\right)^{3} - & r<\ell\\ - 0 & r\geq\ell - \end{cases} - """ - r = np.array(np.abs(r), dtype=np.double) - res = np.zeros_like(r) - r_ll = r < self.len_scale - r_low = r[r_ll] - if self.dim == 1: - res[r_ll] = 1.0 - r_low / self.len_scale - elif self.dim == 2: - res[r_ll] = ( - 2 - / np.pi - * ( - np.arccos(r_low / self.len_scale) - - r_low - / self.len_scale - * np.sqrt(1 - (r_low / self.len_scale) ** 2) - ) - ) - else: - res[r_ll] = ( - 1.0 - - 3.0 / 2.0 * r_low / self.len_scale - + 1.0 / 2.0 * (r_low / self.len_scale) ** 3 - ) + def cor(self, h): + """Hyper-Spherical normalized correlation function.""" + h = np.array(h, dtype=np.double) + res = np.zeros_like(h) + h_l1 = h < 1 + nu = (self.dim - 1.0) / 2.0 + fac = 1.0 / sps.hyp2f1(0.5, -nu, 1.5, 1) + res[h_l1] = 1 - h[h_l1] * fac * sps.hyp2f1(0.5, -nu, 1.5, h[h_l1] ** 2) return res def spectral_density(self, k): # noqa: D102 k = np.array(k, dtype=np.double) res = np.empty_like(k) - kl = k * self.len_scale - kl_gz = kl > 0 - # for k=0 we calculate the limit by hand - if self.dim == 1: - res[kl_gz] = (1.0 - np.cos(kl[kl_gz])) / ( - np.pi * k[kl_gz] * kl[kl_gz] - ) - res[np.logical_not(kl_gz)] = self.len_scale / 2.0 / np.pi - elif self.dim == 2: - res[kl_gz] = sps.j1(kl[kl_gz] / 2.0) ** 2 / np.pi / k[kl_gz] ** 2 - res[np.logical_not(kl_gz)] = self.len_scale ** 2 / 16.0 / np.pi - else: - res[kl_gz] = -( - 12 * kl[kl_gz] * np.sin(kl[kl_gz]) - + (12 - 3 * kl[kl_gz] ** 2) * np.cos(kl[kl_gz]) - - 3 * kl[kl_gz] ** 2 - - 12 - ) / (2 * np.pi ** 2 * kl[kl_gz] ** 3 * k[kl_gz] ** 3) - res[np.logical_not(kl_gz)] = ( - self.len_scale ** 3 / 48.0 / np.pi ** 2 - ) + kl = k * self.len_rescaled + kl_gz = np.logical_not(np.isclose(k, 0)) + res[kl_gz] = sps.gamma(self.dim / 2 + 1) / np.sqrt(np.pi) ** self.dim + res[kl_gz] *= sps.jv(self.dim / 2, kl[kl_gz] / 2) ** 2 + res[kl_gz] /= k[kl_gz] ** self.dim + res[np.logical_not(kl_gz)] = ( + (self.len_rescaled / 4) ** self.dim + / sps.gamma(self.dim / 2 + 1) + / np.sqrt(np.pi) ** self.dim + ) return res From a1b0b2a13a471c30148268a7ba36df63817d8e73 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Wed, 11 Nov 2020 01:43:51 +0100 Subject: [PATCH 11/69] CovModel: replace the 'Intersection' model with the 'HyperSpherical' model in the package imports and documentation --- examples/02_cov_model/README.rst | 2 +- gstools/__init__.py | 6 +++--- gstools/covmodel/__init__.py | 6 +++--- gstools/covmodel/models.py | 4 ++-- tests/test_covmodel.py | 6 +++--- 5 files changed, 12 insertions(+), 12 deletions(-) diff --git a/examples/02_cov_model/README.rst b/examples/02_cov_model/README.rst index 6ea0032a..a8d34023 100644 --- a/examples/02_cov_model/README.rst +++ b/examples/02_cov_model/README.rst @@ -54,7 +54,7 @@ The following standard covariance models are provided by GSTools Linear Circular Spherical - Intersection + HyperSpherical As a special feature, we also provide truncated power law (TPL) covariance models diff --git a/gstools/__init__.py b/gstools/__init__.py index 1bafdfe4..7eadf470 100644 --- a/gstools/__init__.py +++ b/gstools/__init__.py @@ -60,7 +60,7 @@ Linear Circular Spherical - Intersection + HyperSpherical Truncated Power Law Covariance Models ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -123,7 +123,7 @@ Linear, Circular, Spherical, - Intersection, + HyperSpherical, TPLGaussian, TPLExponential, TPLStable, @@ -148,7 +148,7 @@ "Linear", "Circular", "Spherical", - "Intersection", + "HyperSpherical", "TPLGaussian", "TPLExponential", "TPLStable", diff --git a/gstools/covmodel/__init__.py b/gstools/covmodel/__init__.py index f638f2cf..4d569f53 100644 --- a/gstools/covmodel/__init__.py +++ b/gstools/covmodel/__init__.py @@ -34,7 +34,7 @@ Linear Circular Spherical - Intersection + HyperSpherical Truncated Power Law Covariance Models @@ -56,7 +56,7 @@ Linear, Circular, Spherical, - Intersection, + HyperSpherical, ) from gstools.covmodel.tpl_models import TPLGaussian, TPLExponential, TPLStable @@ -70,7 +70,7 @@ "Linear", "Circular", "Spherical", - "Intersection", + "HyperSpherical", "TPLGaussian", "TPLExponential", "TPLStable", diff --git a/gstools/covmodel/models.py b/gstools/covmodel/models.py index 406c26b3..aab98123 100644 --- a/gstools/covmodel/models.py +++ b/gstools/covmodel/models.py @@ -15,7 +15,7 @@ Linear Circular Spherical - Intersection + HyperSpherical """ # pylint: disable=C0103, E1101, E1137 @@ -33,7 +33,7 @@ "Linear", "Circular", "Spherical", - "Intersection", + "HyperSpherical", ] diff --git a/tests/test_covmodel.py b/tests/test_covmodel.py index 138d6b39..372d77b7 100644 --- a/tests/test_covmodel.py +++ b/tests/test_covmodel.py @@ -14,7 +14,7 @@ Linear, Circular, Spherical, - Intersection, + HyperSpherical, TPLGaussian, TPLExponential, TPLStable, @@ -32,7 +32,7 @@ def setUp(self): Linear, Circular, Spherical, - Intersection, + HyperSpherical, TPLGaussian, TPLExponential, TPLStable, @@ -46,7 +46,7 @@ def setUp(self): Linear, Circular, Spherical, - Intersection, + HyperSpherical, ] self.dims = range(1, 4) self.lens = [[10, 5, 2]] From a03a2e5c47610d69a99e758f5b80d7091d5c52fe Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Wed, 11 Nov 2020 01:48:59 +0100 Subject: [PATCH 12/69] covmodel: comment update --- gstools/covmodel/models.py | 26 +------------------------- 1 file changed, 1 insertion(+), 25 deletions(-) diff --git a/gstools/covmodel/models.py b/gstools/covmodel/models.py index aab98123..df9af9ae 100644 --- a/gstools/covmodel/models.py +++ b/gstools/covmodel/models.py @@ -4,7 +4,7 @@ .. currentmodule:: gstools.covmodel.models -The following classes and functions are provided +The following classes are provided .. autosummary:: Gaussian @@ -37,9 +37,6 @@ ] -# Gaussian Model ############################################################## - - class Gaussian(CovModel): r"""The Gaussian covariance model. @@ -106,9 +103,6 @@ def calc_integral_scale(self): # noqa: D102 return self.len_rescaled * np.sqrt(np.pi) / 2.0 -# Exponential Model ########################################################### - - class Exponential(CovModel): r"""The Exponential covariance model. @@ -187,9 +181,6 @@ def calc_integral_scale(self): # noqa: D102 return self.len_rescaled -# Rational Model ############################################################## - - class Rational(CovModel): r"""The rational quadratic covariance model. @@ -249,9 +240,6 @@ def calc_integral_scale(self): # noqa: D102 ) -# Stable Model ################################################################ - - class Stable(CovModel): r"""The stable covariance model. @@ -320,9 +308,6 @@ def calc_integral_scale(self): # noqa: D102 return self.len_rescaled * sps.gamma(1.0 + 1.0 / self.alpha) -# Matérn Model ################################################################ - - class Matern(CovModel): r"""The Matérn covariance model. @@ -433,9 +418,6 @@ def calc_integral_scale(self): # noqa: D102 ) -# Bounded linear Model ######################################################## - - class Linear(CovModel): r"""The bounded linear covariance model. @@ -462,9 +444,6 @@ def cor(self, h): return np.maximum(1 - np.abs(h, dtype=np.double), 0.0) -# Circular Model ############################################################## - - class Circular(CovModel): r"""The circular covariance model. @@ -504,9 +483,6 @@ def cor(self, h): return res -# Spherical Model ############################################################# - - class Spherical(CovModel): r"""The Spherical covariance model. From 7392ea04075745cbfc506fe268e6d9d1404854be Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Wed, 11 Nov 2020 02:04:28 +0100 Subject: [PATCH 13/69] CovModel: add the Super-Spherical model --- examples/02_cov_model/README.rst | 1 + gstools/__init__.py | 3 ++ gstools/covmodel/__init__.py | 3 ++ gstools/covmodel/models.py | 80 ++++++++++++++++++++++++++++++-- tests/test_covmodel.py | 3 ++ 5 files changed, 87 insertions(+), 3 deletions(-) diff --git a/examples/02_cov_model/README.rst b/examples/02_cov_model/README.rst index a8d34023..68155956 100644 --- a/examples/02_cov_model/README.rst +++ b/examples/02_cov_model/README.rst @@ -55,6 +55,7 @@ The following standard covariance models are provided by GSTools Circular Spherical HyperSpherical + SuperSpherical As a special feature, we also provide truncated power law (TPL) covariance models diff --git a/gstools/__init__.py b/gstools/__init__.py index 7eadf470..3d9f12bb 100644 --- a/gstools/__init__.py +++ b/gstools/__init__.py @@ -61,6 +61,7 @@ Circular Spherical HyperSpherical + SuperSpherical Truncated Power Law Covariance Models ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -124,6 +125,7 @@ Circular, Spherical, HyperSpherical, + SuperSpherical, TPLGaussian, TPLExponential, TPLStable, @@ -149,6 +151,7 @@ "Circular", "Spherical", "HyperSpherical", + "SuperSpherical", "TPLGaussian", "TPLExponential", "TPLStable", diff --git a/gstools/covmodel/__init__.py b/gstools/covmodel/__init__.py index 4d569f53..9b36401c 100644 --- a/gstools/covmodel/__init__.py +++ b/gstools/covmodel/__init__.py @@ -35,6 +35,7 @@ Circular Spherical HyperSpherical + SuperSpherical Truncated Power Law Covariance Models @@ -57,6 +58,7 @@ Circular, Spherical, HyperSpherical, + SuperSpherical, ) from gstools.covmodel.tpl_models import TPLGaussian, TPLExponential, TPLStable @@ -71,6 +73,7 @@ "Circular", "Spherical", "HyperSpherical", + "SuperSpherical", "TPLGaussian", "TPLExponential", "TPLStable", diff --git a/gstools/covmodel/models.py b/gstools/covmodel/models.py index df9af9ae..af9dc8a5 100644 --- a/gstools/covmodel/models.py +++ b/gstools/covmodel/models.py @@ -16,6 +16,7 @@ Circular Spherical HyperSpherical + SuperSpherical """ # pylint: disable=C0103, E1101, E1137 @@ -34,6 +35,7 @@ "Circular", "Spherical", "HyperSpherical", + "SuperSpherical", ] @@ -530,10 +532,10 @@ class HyperSpherical(CovModel): .. math:: \rho(r) = \begin{cases} - 1-\frac{ - _{2}F_{1}\left(\frac{1}{2},\frac{1-d}{2},\frac{3}{2}, + 1-s\cdot\frac{r}{\ell}\cdot\frac{ + _{2}F_{1}\left(\frac{1}{2},-\frac{d-1}{2},\frac{3}{2}, \left(s\cdot\frac{r}{\ell}\right)^{2}\right)} - {_{2}F_{1}\left(\frac{1}{2},\frac{1-d}{2},\frac{3}{2},1\right)} + {_{2}F_{1}\left(\frac{1}{2},-\frac{d-1}{2},\frac{3}{2},1\right)} & r<\frac{\ell}{s}\\ 0 & r\geq\frac{\ell}{s} \end{cases} @@ -566,3 +568,75 @@ def spectral_density(self, k): # noqa: D102 / np.sqrt(np.pi) ** self.dim ) return res + + +class SuperSpherical(CovModel): + r"""The Super-Spherical covariance model. + + This model is derived from the relative intersection area of + two d-dimensional hyper spheres, + where the middle points have a distance of :math:`r` + and the diameters are given by :math:`\ell`. + It is than valid in all lower dimensions. + By default it coincides with the Hyper-Spherical model. + + Notes + ----- + This model is given by the following correlation functions. + + .. math:: + \rho(r) = + \begin{cases} + 1-s\cdot\frac{r}{\ell}\cdot\frac{ + _{2}F_{1}\left(\frac{1}{2},-\nu,\frac{3}{2}, + \left(s\cdot\frac{r}{\ell}\right)^{2}\right)} + {_{2}F_{1}\left(\frac{1}{2},-\nu,\frac{3}{2},1\right)} + & r<\frac{\ell}{s}\\ + 0 & r\geq\frac{\ell}{s} + \end{cases} + + Where the standard rescale factor is :math:`s=1`. + :math:`\nu\geq\frac{d-1}{2}` is a shape parameter. + + Other Parameters + ---------------- + nu : :class:`float`, optional + Shape parameter. Standard range: ``[(dim-1)/2, inf]`` + Default: ``(dim-1)/2`` + """ + + def cor(self, h): + """Super-Spherical normalized correlation function.""" + h = np.array(h, dtype=np.double) + res = np.zeros_like(h) + h_l1 = h < 1 + nu = self.nu + fac = 1.0 / sps.hyp2f1(0.5, -nu, 1.5, 1.0) + res[h_l1] = 1.0 - h[h_l1] * fac * sps.hyp2f1( + 0.5, -nu, 1.5, h[h_l1] ** 2 + ) + return res + + def default_opt_arg(self): + """Defaults for the optional arguments. + + * ``{"nu": 1.0}`` + + Returns + ------- + :class:`dict` + Defaults for optional arguments + """ + return {"nu": (self.dim - 1) / 2} + + def default_opt_arg_bounds(self): + """Defaults for boundaries of the optional arguments. + + * ``{"nu": [0.5, 30.0, "cc"]}`` + + Returns + ------- + :class:`dict` + Boundaries for optional arguments + """ + return {"nu": [(self.dim - 1) / 2, np.inf, "co"]} diff --git a/tests/test_covmodel.py b/tests/test_covmodel.py index 372d77b7..a9bb1e00 100644 --- a/tests/test_covmodel.py +++ b/tests/test_covmodel.py @@ -15,6 +15,7 @@ Circular, Spherical, HyperSpherical, + SuperSpherical, TPLGaussian, TPLExponential, TPLStable, @@ -33,6 +34,7 @@ def setUp(self): Circular, Spherical, HyperSpherical, + SuperSpherical, TPLGaussian, TPLExponential, TPLStable, @@ -47,6 +49,7 @@ def setUp(self): Circular, Spherical, HyperSpherical, + SuperSpherical, ] self.dims = range(1, 4) self.lens = [[10, 5, 2]] From fb28c337426b3d14dc04ef98db479151e4bf9418 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Wed, 11 Nov 2020 02:05:48 +0100 Subject: [PATCH 14/69] SuperSpherical: doc fix --- gstools/covmodel/models.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/gstools/covmodel/models.py b/gstools/covmodel/models.py index af9dc8a5..b719d8ff 100644 --- a/gstools/covmodel/models.py +++ b/gstools/covmodel/models.py @@ -620,7 +620,7 @@ def cor(self, h): def default_opt_arg(self): """Defaults for the optional arguments. - * ``{"nu": 1.0}`` + * ``{"nu": (dim-1)/2}`` Returns ------- @@ -632,7 +632,7 @@ def default_opt_arg(self): def default_opt_arg_bounds(self): """Defaults for boundaries of the optional arguments. - * ``{"nu": [0.5, 30.0, "cc"]}`` + * ``{"nu": [(dim-1)/2, inf, "co"]}`` Returns ------- From 61bdee237bb791a02cdf0c124f017816720a27b0 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Wed, 11 Nov 2020 13:11:12 +0100 Subject: [PATCH 15/69] CovModel: add the JBessel hole model --- examples/02_cov_model/README.rst | 1 + gstools/__init__.py | 3 + gstools/covmodel/__init__.py | 3 + gstools/covmodel/models.py | 112 +++++++++++++++++++++++++++++-- tests/test_covmodel.py | 3 + 5 files changed, 118 insertions(+), 4 deletions(-) diff --git a/examples/02_cov_model/README.rst b/examples/02_cov_model/README.rst index 68155956..c5477449 100644 --- a/examples/02_cov_model/README.rst +++ b/examples/02_cov_model/README.rst @@ -56,6 +56,7 @@ The following standard covariance models are provided by GSTools Spherical HyperSpherical SuperSpherical + JBessel As a special feature, we also provide truncated power law (TPL) covariance models diff --git a/gstools/__init__.py b/gstools/__init__.py index 3d9f12bb..b7cc17db 100644 --- a/gstools/__init__.py +++ b/gstools/__init__.py @@ -62,6 +62,7 @@ Spherical HyperSpherical SuperSpherical + JBessel Truncated Power Law Covariance Models ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -126,6 +127,7 @@ Spherical, HyperSpherical, SuperSpherical, + JBessel, TPLGaussian, TPLExponential, TPLStable, @@ -152,6 +154,7 @@ "Spherical", "HyperSpherical", "SuperSpherical", + "JBessel", "TPLGaussian", "TPLExponential", "TPLStable", diff --git a/gstools/covmodel/__init__.py b/gstools/covmodel/__init__.py index 9b36401c..5a478364 100644 --- a/gstools/covmodel/__init__.py +++ b/gstools/covmodel/__init__.py @@ -36,6 +36,7 @@ Spherical HyperSpherical SuperSpherical + JBessel Truncated Power Law Covariance Models @@ -59,6 +60,7 @@ Spherical, HyperSpherical, SuperSpherical, + JBessel, ) from gstools.covmodel.tpl_models import TPLGaussian, TPLExponential, TPLStable @@ -74,6 +76,7 @@ "Spherical", "HyperSpherical", "SuperSpherical", + "JBessel", "TPLGaussian", "TPLExponential", "TPLStable", diff --git a/gstools/covmodel/models.py b/gstools/covmodel/models.py index b719d8ff..26ba038e 100644 --- a/gstools/covmodel/models.py +++ b/gstools/covmodel/models.py @@ -17,6 +17,7 @@ Spherical HyperSpherical SuperSpherical + JBessel """ # pylint: disable=C0103, E1101, E1137 @@ -36,6 +37,7 @@ "Spherical", "HyperSpherical", "SuperSpherical", + "JBessel", ] @@ -329,12 +331,12 @@ class Matern(CovModel): :math:`\nu` is a shape parameter and should be >= 0.2. - If :math:`\nu > 20`, a gaussian model is used, since it is the limit - case: + If :math:`\nu > 20`, a gaussian model is used, since it represents + the limiting case: .. math:: \rho(r) = - \exp\left(-\left(\frac{r}{2\ell}\right)^2\right) + \exp\left(-\left(s\cdot\frac{r}{2\ell}\right)^2\right) Other Parameters ---------------- @@ -601,7 +603,7 @@ class SuperSpherical(CovModel): Other Parameters ---------------- nu : :class:`float`, optional - Shape parameter. Standard range: ``[(dim-1)/2, inf]`` + Shape parameter. Standard range: ``[(dim-1)/2, inf)`` Default: ``(dim-1)/2`` """ @@ -640,3 +642,105 @@ def default_opt_arg_bounds(self): Boundaries for optional arguments """ return {"nu": [(self.dim - 1) / 2, np.inf, "co"]} + + +class JBessel(CovModel): + r"""The J-Bessel hole model. + + This covariance model is a valid hole model, meaning it has areas + of negative correlation but a valid spectral density. + + Notes + ----- + This model is given by the following correlation functions. + + .. math:: + \rho(r) = + \Gamma(\nu+1) \cdot + \frac{\mathrm{J}_{\nu}\left(s\cdot\frac{r}{\ell}\right)} + {\left(s\cdot\frac{r}{2\ell}\right)^{\nu}} + + Where the standard rescale factor is :math:`s=1`. + :math:`\Gamma` is the gamma function and :math:`\mathrm{J}_{\nu}` + is the Bessel functions of the first kind. + :math:`\nu\geq\frac{d}{2}-1` is a shape parameter, + which defaults to :math:`\nu=\frac{d}{2}`, + since the spectrum of the model gets instable for + :math:`\nu\to\frac{d}{2}-1`. + + For :math:`\nu=\frac{1}{2}` (valid in d=1,2,3) + we get the so-called 'Wave' model: + + .. math:: + \rho(r) = + \frac{\sin\left(s\cdot\frac{r}{\ell}\right)}{s\cdot\frac{r}{\ell}} + + Other Parameters + ---------------- + nu : :class:`float`, optional + Shape parameter. Standard range: ``[dim/2 - 1, inf)`` + Default: ``dim/2`` + """ + + def default_opt_arg(self): + """Defaults for the optional arguments. + + * ``{"nu": dim/2}`` + + Returns + ------- + :class:`dict` + Defaults for optional arguments + """ + return {"nu": self.dim / 2} + + def default_opt_arg_bounds(self): + """Defaults for boundaries of the optional arguments. + + * ``{"nu": [dim/2 - 1, inf, "co"]}`` + + Returns + ------- + :class:`dict` + Boundaries for optional arguments + """ + return {"nu": [self.dim / 2 - 1, np.inf, "co"]} + + def check_opt_arg(self): + """Check the optional arguments. + + Warns + ----- + nu + If nu is close to dim/2 - 1, the model tends to get unstable. + """ + if abs(self.nu - self.dim / 2 + 1) < 0.01: + warnings.warn( + "JBessel: parameter 'nu' is close to d/2-1, " + + "count with unstable results" + ) + + def cor(self, h): + """J-Bessel correlation.""" + h = np.array(h, dtype=np.double) + h_gz = np.logical_not(np.isclose(h, 0)) + hh = h[h_gz] + res = np.ones_like(h) + nu = self.nu + res[h_gz] = sps.gamma(nu + 1) * sps.jv(nu, hh) / (hh / 2.0) ** nu + return res + + def spectral_density(self, k): # noqa: D102 + k = np.array(k, dtype=np.double) + k_ll = k < 1.0 / self.len_rescaled + kk = k[k_ll] + res = np.zeros_like(k) + # the model is degenerated for nu=d/2-1, so we tweak the spectral pdf + # and cut of the divisor at nu-(d/2-1)=0.01 (gamma(0.01) about 100) + res[k_ll] = ( + (self.len_rescaled / np.sqrt(np.pi)) ** self.dim + * sps.gamma(self.nu + 1.0) + / np.minimum(sps.gamma(self.nu - self.dim / 2 + 1), 100.0) + * (1.0 - (kk * self.len_rescaled) ** 2) ** (self.nu - self.dim / 2) + ) + return res diff --git a/tests/test_covmodel.py b/tests/test_covmodel.py index a9bb1e00..05eb77c6 100644 --- a/tests/test_covmodel.py +++ b/tests/test_covmodel.py @@ -16,6 +16,7 @@ Spherical, HyperSpherical, SuperSpherical, + JBessel, TPLGaussian, TPLExponential, TPLStable, @@ -35,6 +36,7 @@ def setUp(self): Spherical, HyperSpherical, SuperSpherical, + JBessel, TPLGaussian, TPLExponential, TPLStable, @@ -50,6 +52,7 @@ def setUp(self): Spherical, HyperSpherical, SuperSpherical, + JBessel, ] self.dims = range(1, 4) self.lens = [[10, 5, 2]] From 97c3487ec0e7bc7998a40a0984b70d6b534a060b Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Wed, 11 Nov 2020 14:06:29 +0100 Subject: [PATCH 16/69] TPLmodels: refactor to use 'cor' and 'rescale'; use a new base-class for TPL models derived from super-positioning --- gstools/covmodel/tpl_models.py | 275 +++++++++++---------------------- 1 file changed, 87 insertions(+), 188 deletions(-) diff --git a/gstools/covmodel/tpl_models.py b/gstools/covmodel/tpl_models.py index 27a9e0dc..d88e25e2 100644 --- a/gstools/covmodel/tpl_models.py +++ b/gstools/covmodel/tpl_models.py @@ -26,10 +26,45 @@ __all__ = ["TPLGaussian", "TPLExponential", "TPLStable"] +class TPLCovModel(CovModel): + """Truncated-Power-Law Covariance Model base class for super-position.""" + + @property + def len_up(self): + """:class:`float`: Upper length scale truncation of the model. + + * ``len_up = len_low + len_scale`` + """ + return self.len_low + self.len_scale + + @property + def len_up_rescaled(self): + """:class:`float`: Upper length scale truncation rescaled. + + * ``len_up_rescaled = (len_low + len_scale) / rescale`` + """ + return (self.len_low + self.len_scale) / self.rescale + + @property + def len_low_rescaled(self): + """:class:`float`: Lower length scale truncation rescaled. + + * ``len_low_rescaled = len_low / rescale`` + """ + return self.len_low / self.rescale + + def var_factor(self): + """Factor for C (intensity of variation) to result in variance.""" + return ( + self.len_up_rescaled ** (2 * self.hurst) + - self.len_low_rescaled ** (2 * self.hurst) + ) / (2 * self.hurst) + + # Truncated power law ######################################################### -class TPLGaussian(CovModel): +class TPLGaussian(TPLCovModel): r"""Truncated-Power-Law with Gaussian modes. Notes @@ -77,9 +112,10 @@ class TPLGaussian(CovModel): If you want to define an upper scale truncation, you should set ``len_low`` and ``len_scale`` accordingly. - The following Parameters occure: + The following Parameters occur: - * :math:`C>0` : scaling factor from the Power-Law + * :math:`C>0` : + scaling factor from the Power-Law (intensity of variation) This parameter will be calculated internally by the given variance. You can access C directly by ``model.var_raw`` * :math:`00` : scaling factor from the Power-Law + * :math:`C>0` : + scaling factor from the Power-Law (intensity of variation) This parameter will be calculated internally by the given variance. You can access C directly by ``model.var_raw`` * :math:`00` : scaling factor from the Power-Law + * :math:`C>0` : + scaling factor from the Power-Law (intensity of variation) This parameter will be calculated internally by the given variance. You can access C directly by ``model.var_raw`` * :math:`0 Date: Wed, 11 Nov 2020 14:06:48 +0100 Subject: [PATCH 17/69] Doc: minor fixes --- examples/05_kriging/README.rst | 2 +- gstools/covmodel/base.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/05_kriging/README.rst b/examples/05_kriging/README.rst index 834efa95..d2bbe098 100644 --- a/examples/05_kriging/README.rst +++ b/examples/05_kriging/README.rst @@ -42,7 +42,7 @@ the features you want: one can set `exact` to `False` and provide either individual measurement errors for each point or set the nugget as a constant measurement error everywhere. * `pseudo_inv`: Sometimes the inversion of the kriging matrix can be numerically unstable. - This occures for examples in cases of redundant input values. In this case we provide a switch to + This occurs for examples in cases of redundant input values. In this case we provide a switch to use the pseudo-inverse of the matrix. Then redundant conditional values will automatically be averaged. diff --git a/gstools/covmodel/base.py b/gstools/covmodel/base.py index a61b0b96..0f1a7d01 100644 --- a/gstools/covmodel/base.py +++ b/gstools/covmodel/base.py @@ -450,7 +450,7 @@ def pykrige_kwargs(self): kwargs.update(add_kwargs) return kwargs - # methods for optional arguments (can be overridden) + # methods for optional/default arguments (can be overridden) def default_opt_arg(self): """Provide default optional arguments by the user. From 550e02e5c36f92c19cdd9f42b6b991288290a0e9 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Wed, 11 Nov 2020 14:17:25 +0100 Subject: [PATCH 18/69] TPLCovModel: tweak base class to provide cor and correlation like derived classes --- gstools/covmodel/tpl_models.py | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/gstools/covmodel/tpl_models.py b/gstools/covmodel/tpl_models.py index d88e25e2..d138c7e2 100644 --- a/gstools/covmodel/tpl_models.py +++ b/gstools/covmodel/tpl_models.py @@ -60,6 +60,14 @@ def var_factor(self): - self.len_low_rescaled ** (2 * self.hurst) ) / (2 * self.hurst) + def cor(self, h): + """TPL - normalized correlation function.""" + return 1.0 + + def correlation(self, r): + """TPL - correlation function.""" + return 1.0 + # Truncated power law ######################################################### @@ -320,7 +328,7 @@ def spectral_density(self, k): # noqa: D102 ) -class TPLStable(CovModel): +class TPLStable(TPLCovModel): r"""Truncated-Power-Law with Stable modes. Notes From 64017b3edae0b0952c55db23b35e4f2b00d95036 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Wed, 11 Nov 2020 17:05:40 +0100 Subject: [PATCH 19/69] CovModel: add TPLSimple model --- examples/02_cov_model/README.rst | 1 + gstools/__init__.py | 3 ++ gstools/covmodel/__init__.py | 9 +++- gstools/covmodel/tpl_models.py | 71 +++++++++++++++++++++++++++++++- tests/test_covmodel.py | 2 + 5 files changed, 84 insertions(+), 2 deletions(-) diff --git a/examples/02_cov_model/README.rst b/examples/02_cov_model/README.rst index c5477449..5825b8cb 100644 --- a/examples/02_cov_model/README.rst +++ b/examples/02_cov_model/README.rst @@ -64,6 +64,7 @@ As a special feature, we also provide truncated power law (TPL) covariance model TPLGaussian TPLExponential TPLStable + TPLSimple Gallery ------- diff --git a/gstools/__init__.py b/gstools/__init__.py index 61612148..5873a84d 100644 --- a/gstools/__init__.py +++ b/gstools/__init__.py @@ -72,6 +72,7 @@ TPLGaussian TPLExponential TPLStable + TPLSimple Functions ========= @@ -133,6 +134,7 @@ TPLGaussian, TPLExponential, TPLStable, + TPLSimple, ) try: @@ -160,6 +162,7 @@ "TPLGaussian", "TPLExponential", "TPLStable", + "TPLSimple", ] __all__ += [ diff --git a/gstools/covmodel/__init__.py b/gstools/covmodel/__init__.py index 5a478364..fef7ff4b 100644 --- a/gstools/covmodel/__init__.py +++ b/gstools/covmodel/__init__.py @@ -46,6 +46,7 @@ TPLGaussian TPLExponential TPLStable + TPLSimple """ from gstools.covmodel.base import CovModel @@ -62,7 +63,12 @@ SuperSpherical, JBessel, ) -from gstools.covmodel.tpl_models import TPLGaussian, TPLExponential, TPLStable +from gstools.covmodel.tpl_models import ( + TPLGaussian, + TPLExponential, + TPLStable, + TPLSimple, +) __all__ = [ "CovModel", @@ -80,4 +86,5 @@ "TPLGaussian", "TPLExponential", "TPLStable", + "TPLSimple", ] diff --git a/gstools/covmodel/tpl_models.py b/gstools/covmodel/tpl_models.py index d138c7e2..8deae0b6 100644 --- a/gstools/covmodel/tpl_models.py +++ b/gstools/covmodel/tpl_models.py @@ -10,6 +10,7 @@ TPLGaussian TPLExponential TPLStable + TPLSimple """ # pylint: disable=C0103, E1101 @@ -23,7 +24,7 @@ ) -__all__ = ["TPLGaussian", "TPLExponential", "TPLStable"] +__all__ = ["TPLGaussian", "TPLExponential", "TPLStable", "TPLSimple"] class TPLCovModel(CovModel): @@ -480,3 +481,71 @@ def correlation(self, r): self.len_up_rescaled ** (2 * self.hurst) - self.len_low_rescaled ** (2 * self.hurst) ) + + +class TPLSimple(CovModel): + r"""The simply truncated power law model. + + This model describes a simple truncated power law + with a finite length scale. In contrast to other models, + this one is not derived from super-postioning modes. + + Notes + ----- + This model is given by the following correlation function: + + .. math:: + \rho(r) = + \begin{cases} + \left(1-s\cdot\frac{r}{\ell}\right)^{\nu} & r<\frac{\ell}{s}\\ + 0 & r\geq\frac{\ell}{s} + \end{cases} + + Where the standard rescale factor is :math:`s=1`. + :math:`\nu\geq\frac{d+1}{2}` is a shape parameter, + which defaults to :math:`\nu=\frac{d+1}{2}`, + + For :math:`\nu=1` (valid only in d=1) + this coincides with the truncated linear model: + + .. math:: + \rho(r) = + \begin{cases} + 1-s\cdot\frac{r}{\ell} & r<\frac{\ell}{s}\\ + 0 & r\geq\frac{\ell}{s} + \end{cases} + + Other Parameters + ---------------- + nu : :class:`float`, optional + Shape parameter. Standard range: ``[(dim+1)/2, inf)`` + Default: ``dim/2`` + """ + + def default_opt_arg(self): + """Defaults for the optional arguments. + + * ``{"nu": dim/2}`` + + Returns + ------- + :class:`dict` + Defaults for optional arguments + """ + return {"nu": (self.dim + 1) / 2} + + def default_opt_arg_bounds(self): + """Defaults for boundaries of the optional arguments. + + * ``{"nu": [dim/2 - 1, inf, "co"]}`` + + Returns + ------- + :class:`dict` + Boundaries for optional arguments + """ + return {"nu": [(self.dim + 1) / 2, np.inf, "co"]} + + def cor(self, h): + """TPL Simple - normalized correlation function.""" + return np.maximum(1 - np.abs(h, dtype=np.double), 0.0) ** self.nu diff --git a/tests/test_covmodel.py b/tests/test_covmodel.py index 05eb77c6..c9f7794f 100644 --- a/tests/test_covmodel.py +++ b/tests/test_covmodel.py @@ -20,6 +20,7 @@ TPLGaussian, TPLExponential, TPLStable, + TPLSimple, ) @@ -40,6 +41,7 @@ def setUp(self): TPLGaussian, TPLExponential, TPLStable, + TPLSimple, ] self.std_cov_models = [ Gaussian, From 82ab29c8bfec7481af5cc1f1df30a4d33bdbaccb Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Fri, 13 Nov 2020 01:33:30 +0100 Subject: [PATCH 20/69] field: unify rotation and anisotropy; simplify code --- gstools/field/tools.py | 20 +-- gstools/tools/__init__.py | 48 +++++-- gstools/tools/geometric.py | 276 +++++++++++++++++++++++++++++++++---- 3 files changed, 292 insertions(+), 52 deletions(-) diff --git a/gstools/field/tools.py b/gstools/field/tools.py index e9f289a2..00beb8fd 100644 --- a/gstools/field/tools.py +++ b/gstools/field/tools.py @@ -21,7 +21,7 @@ # pylint: disable=C0103 import numpy as np -from gstools.tools.geometric import r3d_x, r3d_y, r3d_z +from gstools.tools.geometric import matrix_rotate, matrix_derotate __all__ = [ "reshape_input", @@ -137,24 +137,21 @@ def unrotate_mesh(dim, angles, x, y, z): """Rotate axes in order to implement rotation. for 3d: yaw, pitch, and roll angles are alpha, beta, and gamma, - of intrinsic rotation rotation whose Tait-Bryan angles are + of intrinsic rotation whose Tait-Bryan angles are alpha, beta, gamma about axes x, y, z. """ if dim == 1: return x, y, z if dim == 2: # extract 2d rotation matrix - rot_mat = r3d_z(-angles[0])[0:2, 0:2] + rot_mat = matrix_derotate(dim, angles) pos_tuple = np.vstack((x, y)) pos_tuple = np.vsplit(np.dot(rot_mat, pos_tuple), 2) x = pos_tuple[0].reshape(np.shape(x)) y = pos_tuple[1].reshape(np.shape(y)) return x, y, z if dim == 3: - alpha = -angles[0] - beta = -angles[1] - gamma = -angles[2] - rot_mat = np.dot(np.dot(r3d_z(alpha), r3d_y(beta)), r3d_x(gamma)) + rot_mat = matrix_derotate(dim, angles) pos_tuple = np.vstack((x, y, z)) pos_tuple = np.vsplit(np.dot(rot_mat, pos_tuple), 3) x = pos_tuple[0].reshape(np.shape(x)) @@ -168,24 +165,21 @@ def rotate_mesh(dim, angles, x, y, z): """Rotate axes. for 3d: yaw, pitch, and roll angles are alpha, beta, and gamma, - of intrinsic rotation rotation whose Tait-Bryan angles are + of intrinsic rotation whose Tait-Bryan angles are alpha, beta, gamma about axes x, y, z. """ if dim == 1: return x, y, z if dim == 2: # extract 2d rotation matrix - rot_mat = r3d_z(angles[0])[0:2, 0:2] + rot_mat = matrix_rotate(dim, angles) pos_tuple = np.vstack((x, y)) pos_tuple = np.vsplit(np.dot(rot_mat, pos_tuple), 2) x = pos_tuple[0].reshape(np.shape(x)) y = pos_tuple[1].reshape(np.shape(y)) return x, y, z if dim == 3: - alpha = angles[0] - beta = angles[1] - gamma = angles[2] - rot_mat = np.dot(np.dot(r3d_x(gamma), r3d_y(beta)), r3d_z(alpha)) + rot_mat = matrix_rotate(dim, angles) pos_tuple = np.vstack((x, y, z)) pos_tuple = np.vsplit(np.dot(rot_mat, pos_tuple), 3) x = pos_tuple[0].reshape(np.shape(x)) diff --git a/gstools/tools/__init__.py b/gstools/tools/__init__.py index 368cc900..d431b6b2 100644 --- a/gstools/tools/__init__.py +++ b/gstools/tools/__init__.py @@ -31,12 +31,20 @@ ^^^^^^^^^ .. autosummary:: - xyz2pos + set_angles + set_anis + no_of_angles + rotation_planes + givens_rotation + matrix_rotate + matrix_derotate + matrix_isotropify + matrix_anisotropify + matrix_isometrize + matrix_anisometrize pos2xyz + xyz2pos ang2dir - r3d_x - r3d_y - r3d_z ---- """ @@ -61,11 +69,19 @@ ) from gstools.tools.geometric import ( - r3d_x, - r3d_y, - r3d_z, - xyz2pos, + set_angles, + set_anis, + no_of_angles, + rotation_planes, + givens_rotation, + matrix_rotate, + matrix_derotate, + matrix_isotropify, + matrix_anisotropify, + matrix_isometrize, + matrix_anisometrize, pos2xyz, + xyz2pos, ang2dir, ) @@ -83,10 +99,18 @@ "tplstable_cor", "tpl_exp_spec_dens", "tpl_gau_spec_dens", - "xyz2pos", + "set_angles", + "set_anis", + "no_of_angles", + "rotation_planes", + "givens_rotation", + "matrix_rotate", + "matrix_derotate", + "matrix_isotropify", + "matrix_anisotropify", + "matrix_isometrize", + "matrix_anisometrize", "pos2xyz", + "xyz2pos", "ang2dir", - "r3d_x", - "r3d_y", - "r3d_z", ] diff --git a/gstools/tools/geometric.py b/gstools/tools/geometric.py index dc34edba..788b621b 100644 --- a/gstools/tools/geometric.py +++ b/gstools/tools/geometric.py @@ -7,9 +7,17 @@ The following functions are provided .. autosummary:: - r3d_x - r3d_y - r3d_z + set_angles + set_anis + no_of_angles + rotation_planes + givens_rotation + matrix_rotate + matrix_derotate + matrix_isotropify + matrix_anisotropify + matrix_isometrize + matrix_anisometrize pos2xyz xyz2pos ang2dir @@ -18,64 +26,278 @@ import numpy as np -__all__ = ["r3d_x", "r3d_y", "r3d_z", "pos2xyz", "xyz2pos", "ang2dir"] +__all__ = [ + "set_angles", + "set_anis", + "no_of_angles", + "rotation_planes", + "givens_rotation", + "matrix_rotate", + "matrix_derotate", + "matrix_isotropify", + "matrix_anisotropify", + "matrix_isometrize", + "matrix_anisometrize", + "pos2xyz", + "xyz2pos", + "ang2dir", +] # Geometric functions ######################################################### -def r3d_x(theta): - """Rotation matrix about x axis. +def set_angles(dim, angles): + """Set the angles for the given dimension. Parameters ---------- - theta : :class:`float` - Rotation angle + dim : :class:`int` + spatial dimension + angles : :class:`float` or :class:`list` + the angles of the SRF + + Returns + ------- + angles : :class:`float` + the angles fitting to the dimension + + Notes + ----- + If too few angles are given, they are filled up with `0`. + """ + out_angles = np.atleast_1d(angles)[: no_of_angles(dim)] + # fill up the rotation angle array with zeros + out_angles = np.pad( + out_angles, + (0, no_of_angles(dim) - len(out_angles)), + "constant", + constant_values=0.0, + ) + return out_angles + + +def set_anis(dim, anis): + """Set the anisotropy ratios for the given dimension. + + Parameters + ---------- + dim : :class:`int` + spatial dimension + anis : :class:`list` of :class:`float` + the anisotropy of length scales along the transversal directions + + Returns + ------- + anis : :class:`list` of :class:`float` + the anisotropy of length scales fitting the dimensions + + Notes + ----- + If too few anisotrpy ratios are given, they are filled up with `1`. + """ + out_anis = np.atleast_1d(anis)[: dim - 1] + if len(out_anis) < dim - 1: + # fill up the anisotropies with ones, such that len()==dim-1 + out_anis = np.pad( + out_anis, + (dim - len(out_anis) - 1, 0), + "constant", + constant_values=1.0, + ) + return out_anis + + +def no_of_angles(dim): + """Calculate number of rotation angles depending on the dimension. + + Parameters + ---------- + dim : :class:`int` + spatial dimension + + Returns + ------- + :class:`int` + Number of angles. + """ + return (dim * (dim - 1)) // 2 + + +def rotation_planes(dim): + """Get all 2D sub-planes for rotation. + + Parameters + ---------- + dim : :class:`int` + spatial dimension + + Returns + ------- + :class:`list` of :class:`tuple` of :class:`int` + All 2D sub-planes for rotation. + """ + return [(i, j) for j in range(1, dim) for i in range(j)] + + +def givens_rotation(dim, plane, angle): + """Givens rotation matrix in arbitrary dimensions. + + Parameters + ---------- + dim : :class:`int` + spatial dimension + plane : :class:`list` of :class:`int` + the plane to rotate in, given by the indices of the two defining axes. + For example the xy plane is defined by `(0,1)` + angle : :class:`float` or :class:`list` + the rotation angle in the given plane + + Returns + ------- + :class:`numpy.ndarray` + Rotation matrix. + """ + result = np.eye(dim, dtype=np.double) + result[plane[0], plane[0]] = np.cos(angle) + result[plane[1], plane[1]] = np.cos(angle) + result[plane[0], plane[1]] = -np.sin(angle) + result[plane[1], plane[0]] = np.sin(angle) + return result + + +def matrix_rotate(dim, angles): + """Create a matrix to rotate points to the target coordinate-system. + + Parameters + ---------- + dim : :class:`int` + spatial dimension + angles : :class:`float` or :class:`list` + the rotation angles of the target coordinate-system + + Returns + ------- + :class:`numpy.ndarray` + Rotation matrix. + """ + angles = set_angles(dim, angles) + planes = rotation_planes(dim) + result = np.eye(dim, dtype=np.double) + for i, (angle, plane) in enumerate(zip(angles, planes)): + # angles have alternating signs to match tait bryan + result = np.matmul( + givens_rotation(dim, plane, (-1) ** i * angle), result + ) + return result + + +def matrix_derotate(dim, angles): + """Create a matrix to derotate points to the initial coordinate-system. + + Parameters + ---------- + dim : :class:`int` + spatial dimension + angles : :class:`float` or :class:`list` + the rotation angles of the target coordinate-system Returns ------- :class:`numpy.ndarray` Rotation matrix. """ - sin = np.sin(theta) - cos = np.cos(theta) - return np.array(((1.0, +0.0, +0.0), (0.0, cos, -sin), (0.0, sin, cos))) + angles = -set_angles(dim, angles) + planes = rotation_planes(dim) + result = np.eye(dim, dtype=np.double) + for i, (angle, plane) in enumerate(zip(angles, planes)): + # angles have alternating signs to match tait bryan + result = np.matmul( + result, givens_rotation(dim, plane, (-1) ** i * angle) + ) + return result -def r3d_y(theta): - """Rotation matrix about y axis. +def matrix_isotropify(dim, anis): + """Create a stretching matrix to make things isotrope. Parameters ---------- - theta : :class:`float` - Rotation angle + dim : :class:`int` + spatial dimension + anis : :class:`list` of :class:`float` + the anisotropy of length scales along the transversal directions Returns ------- :class:`numpy.ndarray` - Rotation matrix. + Stretching matrix. """ - sin = np.sin(theta) - cos = np.cos(theta) - return np.array(((+cos, 0.0, sin), (+0.0, 1.0, +0.0), (-sin, 0.0, cos))) + anis = set_anis(dim, anis) + return np.diag(np.concatenate(([1.0], 1.0 / anis))) -def r3d_z(theta): - """Rotation matrix about z axis. +def matrix_anisotropify(dim, anis): + """Create a stretching matrix to make things anisotrope. Parameters ---------- - theta : :class:`float` - Rotation angle + dim : :class:`int` + spatial dimension + anis : :class:`list` of :class:`float` + the anisotropy of length scales along the transversal directions Returns ------- :class:`numpy.ndarray` - Rotation matrix. + Stretching matrix. + """ + anis = set_anis(dim, anis) + return np.diag(np.concatenate(([1.0], anis))) + + +def matrix_isometrize(dim, angles, anis): + """Create a matrix to derotate points and make them isotrope. + + Parameters + ---------- + dim : :class:`int` + spatial dimension + angles : :class:`float` or :class:`list` + the rotation angles of the target coordinate-system + anis : :class:`list` of :class:`float` + the anisotropy of length scales along the transversal directions + + Returns + ------- + :class:`numpy.ndarray` + Transformation matrix. + """ + return np.matmul( + matrix_isotropify(dim, anis), matrix_derotate(dim, angles) + ) + + +def matrix_anisometrize(dim, angles, anis): + """Create a matrix to rotate points and make them anisotrope. + + Parameters + ---------- + dim : :class:`int` + spatial dimension + angles : :class:`float` or :class:`list` + the rotation angles of the target coordinate-system + anis : :class:`list` of :class:`float` + the anisotropy of length scales along the transversal directions + + Returns + ------- + :class:`numpy.ndarray` + Transformation matrix. """ - sin = np.sin(theta) - cos = np.cos(theta) - return np.array(((cos, -sin, 0.0), (sin, +cos, 0.0), (+0.0, +0.0, 1.0))) + return np.matmul( + matrix_rotate(dim, angles), matrix_anisotropify(dim, anis) + ) # conversion ################################################################## From 5b29c37282d157d3c56a5eff14ea3586d0ede167 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Fri, 13 Nov 2020 01:35:16 +0100 Subject: [PATCH 21/69] CovModel: add [an]isometrize method; simplify code --- gstools/covmodel/base.py | 33 +++++++++----- gstools/covmodel/tools.py | 93 ++++----------------------------------- 2 files changed, 30 insertions(+), 96 deletions(-) diff --git a/gstools/covmodel/base.py b/gstools/covmodel/base.py index 0f1a7d01..e512cdb5 100644 --- a/gstools/covmodel/base.py +++ b/gstools/covmodel/base.py @@ -17,13 +17,15 @@ from scipy.integrate import quad as integral from scipy.optimize import root from hankel import SymmetricFourierTransform as SFT -from gstools.field.tools import make_isotropic, unrotate_mesh -from gstools.tools.geometric import pos2xyz +from gstools.tools.geometric import ( + set_angles, + matrix_anisometrize, + matrix_isometrize, +) from gstools.covmodel.tools import ( InitSubclassMeta, rad_fac, set_len_anis, - set_angles, check_bounds, check_arg_in_bounds, default_arg_from_bounds, @@ -283,14 +285,6 @@ def cor_from_correlation(self, h): # special variogram functions - def _get_iso_rad(self, pos): - x, y, z = pos2xyz(pos, max_dim=self.dim) - if self.do_rotation: - x, y, z = unrotate_mesh(self.dim, self.angles, x, y, z) - if not self.is_isotropic: - y, z = make_isotropic(self.dim, self.anis, y, z) - return np.linalg.norm((x, y, z)[: self.dim], axis=0) - def vario_spatial(self, pos): r"""Spatial variogram respecting anisotropy and rotation.""" return self.variogram(self._get_iso_rad(pos)) @@ -595,6 +589,23 @@ def _has_ppf(self): """State if a ppf is defined with 'spectral_rad_ppf'.""" return hasattr(self, "spectral_rad_ppf") + # spatial routines + + def isometrize(self, pos): + """Make a position tuple ready for isotropic operations.""" + pos = np.array(pos, dtype=np.double).reshape((self.dim, -1)) + return np.dot(matrix_isometrize(self.dim, self.angles, self.anis), pos) + + def anisometrize(self, pos): + """Bring a position tuple into the anisotropic coordinate-system.""" + pos = np.array(pos, dtype=np.double).reshape((self.dim, -1)) + return np.dot( + matrix_anisometrize(self.dim, self.angles, self.anis), pos + ) + + def _get_iso_rad(self, pos): + return np.linalg.norm(self.isometrize(pos), axis=0) + # fitting routine def fit_variogram( diff --git a/gstools/covmodel/tools.py b/gstools/covmodel/tools.py index 0065dece..ed1ebfff 100644 --- a/gstools/covmodel/tools.py +++ b/gstools/covmodel/tools.py @@ -10,9 +10,6 @@ InitSubclassMeta rad_fac set_len_anis - set_angles - no_of_angles - rot_planes check_bounds check_arg_in_bounds default_arg_from_bounds @@ -21,21 +18,19 @@ # pylint: disable=C0103 import numpy as np from scipy import special as sps +from gstools.tools.geometric import set_anis __all__ = [ "InitSubclassMeta", "rad_fac", "set_len_anis", - "set_angles", - "no_of_angles", - "rot_planes", "check_bounds", "check_arg_in_bounds", "default_arg_from_bounds", ] -# __init_subclass__ hack ###################################################### +# __init_subclass__ hack for Python 3.5 if hasattr(object, "__init_subclass__"): InitSubclassMeta = type @@ -109,16 +104,16 @@ def set_len_anis(dim, len_scale, anis): dim : :class:`int` spatial dimension len_scale : :class:`float` or :class:`list` - the length scale of the SRF in x direction or in x- (y-, z-) direction - anis : :class:`float`/list - the anisotropy of length scales along the y- and z-directions + the length scale of the SRF in x direction or in x- (y-, ...) direction + anis : :class:`float` or :class:`list` + the anisotropy of length scales along the transversal axes Returns ------- len_scale : :class:`float` the main length scale of the SRF in x direction - anis : :class:`float`/list, optional - the anisotropy of length scales along the y- and z-directions + anis : :class:`list`, optional + the anisotropy of length scales along the transversal axes Notes ----- @@ -137,15 +132,7 @@ def set_len_anis(dim, len_scale, anis): out_len_scale = ls_tmp[0] # set the anisotropies in y- and z-direction according to the input if len(ls_tmp) == 1: - out_anis = np.atleast_1d(anis)[: dim - 1] - if len(out_anis) < dim - 1: - # fill up the anisotropies with ones, such that len()==dim-1 - out_anis = np.pad( - out_anis, - (dim - len(out_anis) - 1, 0), - "constant", - constant_values=1.0, - ) + out_anis = set_anis(dim, anis) else: # fill up length-scales with the latter len_scale, such that len()==dim if len(ls_tmp) < dim: @@ -163,70 +150,6 @@ def set_len_anis(dim, len_scale, anis): return out_len_scale, out_anis -def set_angles(dim, angles): - """Set the angles for the given dimension. - - Parameters - ---------- - dim : :class:`int` - spatial dimension - angles : :class:`float` or :class:`list` - the angles of the SRF - - Returns - ------- - angles : :class:`float` - the angles fitting to the dimension - - Notes - ----- - If too few angles are given, they are filled up with `0`. - """ - out_angles = np.atleast_1d(angles)[: no_of_angles(dim)] - # fill up the rotation angle array with zeros - out_angles = np.pad( - out_angles, - (0, no_of_angles(dim) - len(out_angles)), - "constant", - constant_values=0.0, - ) - return out_angles - - -def no_of_angles(dim): - """ - Calculate number of rotation angles depending on the dimension. - - Parameters - ---------- - dim : :class:`int` - spatial dimension - - Returns - ------- - :class:`int` - Number of angles. - """ - return (dim * (dim - 1)) // 2 - - -def rot_planes(dim): - """ - Get all 2D sub-planes for rotation. - - Parameters - ---------- - dim : :class:`int` - spatial dimension - - Returns - ------- - :class:`list` - All 2D sub-planes for rotation. - """ - return [(i, j) for i in range(dim - 1) for j in range(i + 1, dim)] - - def check_bounds(bounds): """ Check if given bounds are valid. From b31bc8a3d1cb7ce902b1a0d8347e84d304d78111 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Fri, 13 Nov 2020 01:35:54 +0100 Subject: [PATCH 22/69] tests: adopt tests: don't allow high-dim input in lower dimensions in spatial routines in CovModel --- tests/test_covmodel.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/tests/test_covmodel.py b/tests/test_covmodel.py index c9f7794f..faa974d7 100644 --- a/tests/test_covmodel.py +++ b/tests/test_covmodel.py @@ -102,10 +102,10 @@ def cor(self, h): model.var * model.correlation(1), ) self.assertAlmostEqual( - model.vario_spatial(([1], [2], [3]))[0], + model.vario_spatial(([1], [2], [3])[:dim])[0], model.var + model.nugget - - model.cov_spatial(([1], [2], [3]))[0], + - model.cov_spatial(([1], [2], [3])[:dim])[0], ) self.assertAlmostEqual( model.cov_nugget(0), model.sill @@ -119,7 +119,7 @@ def cor(self, h): model.vario_nugget(1), model.variogram(1) ) # check if callable - model.vario_spatial((1, 2, 3)) + model.vario_spatial((1, 2, 3)[:dim]) model.spectral_density([0, 1]) model.spectrum([0, 1]) model.spectral_rad_pdf([0, 1]) From 5481b17db3411053d9d5994db68cb818ac7a0642 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Fri, 13 Nov 2020 11:14:11 +0100 Subject: [PATCH 23/69] CovModel: more details in documentation --- examples/02_cov_model/README.rst | 29 +++++++++++++++++++++-------- gstools/covmodel/models.py | 29 +++++++++++++++++------------ 2 files changed, 38 insertions(+), 20 deletions(-) diff --git a/examples/02_cov_model/README.rst b/examples/02_cov_model/README.rst index 5825b8cb..15cee254 100644 --- a/examples/02_cov_model/README.rst +++ b/examples/02_cov_model/README.rst @@ -8,7 +8,6 @@ class, which allows you to easily define arbitrary covariance models by yourself. The resulting models provide a bunch of nice features to explore the covariance models. - A covariance model is used to characterize the `semi-variogram `_, denoted by :math:`\gamma`, of a spatial random field. @@ -16,24 +15,38 @@ In GSTools, we use the following form for an isotropic and stationary field: .. math:: \gamma\left(r\right)= - \sigma^2\cdot\left(1-\rho\left(r\right)\right)+n + \sigma^2\cdot\left(1-\mathrm{cor}\left(s\cdot\frac{r}{\ell}\right)\right)+n Where: - - :math:`\rho(r)` is the so called - `correlation `_ - function depending on the distance :math:`r` + - :math:`r` is the lag distance + - :math:`\ell` is the main correlation length + - :math:`s` is a scaling factor for unit conversion or normalization - :math:`\sigma^2` is the variance - :math:`n` is the nugget (subscale variance) + - :math:`\mathrm{cor}(h)` is normalized correlation function depending on + the non-dimensional distance :math:`h=s\cdot\frac{r}{\ell}` + +Depending on the normalized correlation function, all covariance models in +GSTools are providing the following functions: + + - :math:`\rho(r)=\mathrm{cor}\left(s\cdot\frac{r}{\ell}\right)` + is the so called + `correlation `_ + function + - :math:`C(r)=\sigma^2\cdot\rho(r)` is the so called + `covariance `_ + function, which gives the name for our GSTools class .. note:: We are not limited to isotropic models. GSTools supports anisotropy ratios for length scales in orthogonal transversal directions like: - - :math:`x` (main direction) - - :math:`y` (1. transversal direction) - - :math:`z` (2. transversal direction) + - :math:`x_0` (main direction) + - :math:`x_1` (1. transversal direction) + - :math:`x_2` (2. transversal direction) + - ... These main directions can also be rotated. Just have a look at the corresponding examples. diff --git a/gstools/covmodel/models.py b/gstools/covmodel/models.py index 26ba038e..27ee03d3 100644 --- a/gstools/covmodel/models.py +++ b/gstools/covmodel/models.py @@ -200,6 +200,12 @@ class Rational(CovModel): Where the standard rescale factor is :math:`s=1`. :math:`\alpha` is a shape parameter and should be > 0.5. + For :math:`\alpha\to\infty` this model converges to the Gaussian model: + + .. math:: + \rho(r)= + \exp\left(-\frac{\pi}{4}\left(s\cdot\frac{r}{\ell}\right)^{2}\right) + Other Parameters ---------------- alpha : :class:`float`, optional @@ -607,18 +613,6 @@ class SuperSpherical(CovModel): Default: ``(dim-1)/2`` """ - def cor(self, h): - """Super-Spherical normalized correlation function.""" - h = np.array(h, dtype=np.double) - res = np.zeros_like(h) - h_l1 = h < 1 - nu = self.nu - fac = 1.0 / sps.hyp2f1(0.5, -nu, 1.5, 1.0) - res[h_l1] = 1.0 - h[h_l1] * fac * sps.hyp2f1( - 0.5, -nu, 1.5, h[h_l1] ** 2 - ) - return res - def default_opt_arg(self): """Defaults for the optional arguments. @@ -643,6 +637,17 @@ def default_opt_arg_bounds(self): """ return {"nu": [(self.dim - 1) / 2, np.inf, "co"]} + def cor(self, h): + """Super-Spherical normalized correlation function.""" + h = np.array(h, dtype=np.double) + res = np.zeros_like(h) + h_l1 = h < 1 + fac = 1.0 / sps.hyp2f1(0.5, -self.nu, 1.5, 1.0) + res[h_l1] = 1.0 - h[h_l1] * fac * sps.hyp2f1( + 0.5, -self.nu, 1.5, h[h_l1] ** 2 + ) + return res + class JBessel(CovModel): r"""The J-Bessel hole model. From ec7ba7035b0c71b162cb0ed60c07f1ae02efba15 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Fri, 13 Nov 2020 11:40:36 +0100 Subject: [PATCH 24/69] Tools: add rotated_main_axes routine to determine axes from rotation angles --- gstools/__init__.py | 15 ++++++++++----- gstools/tools/__init__.py | 3 +++ gstools/tools/geometric.py | 23 ++++++++++++++++++++++- 3 files changed, 35 insertions(+), 6 deletions(-) diff --git a/gstools/__init__.py b/gstools/__init__.py index 5873a84d..bbf0ff14 100644 --- a/gstools/__init__.py +++ b/gstools/__init__.py @@ -85,13 +85,16 @@ .. autosummary:: vtk_export - vtk_export_structured - vtk_export_unstructured to_vtk - to_vtk_structured - to_vtk_unstructured -variogram estimation +Geometric +^^^^^^^^^ +Some convenient functions for geometric operations + +.. autosummary:: + rotated_main_axes + +Variogram Estimation ^^^^^^^^^^^^^^^^^^^^ Estimate the variogram of a given field @@ -105,6 +108,7 @@ from gstools import field, variogram, random, covmodel, tools, krige, transform from gstools.field import SRF from gstools.tools import ( + rotated_main_axes, vtk_export, vtk_export_structured, vtk_export_unstructured, @@ -174,6 +178,7 @@ __all__ += [ "SRF", + "rotated_main_axes", "vtk_export", "vtk_export_structured", "vtk_export_unstructured", diff --git a/gstools/tools/__init__.py b/gstools/tools/__init__.py index d431b6b2..bd80462c 100644 --- a/gstools/tools/__init__.py +++ b/gstools/tools/__init__.py @@ -31,6 +31,7 @@ ^^^^^^^^^ .. autosummary:: + rotated_main_axes set_angles set_anis no_of_angles @@ -80,6 +81,7 @@ matrix_anisotropify, matrix_isometrize, matrix_anisometrize, + rotated_main_axes, pos2xyz, xyz2pos, ang2dir, @@ -110,6 +112,7 @@ "matrix_anisotropify", "matrix_isometrize", "matrix_anisometrize", + "rotated_main_axes", "pos2xyz", "xyz2pos", "ang2dir", diff --git a/gstools/tools/geometric.py b/gstools/tools/geometric.py index 788b621b..79f81d9d 100644 --- a/gstools/tools/geometric.py +++ b/gstools/tools/geometric.py @@ -18,6 +18,7 @@ matrix_anisotropify matrix_isometrize matrix_anisometrize + rotated_main_axes pos2xyz xyz2pos ang2dir @@ -38,6 +39,7 @@ "matrix_anisotropify", "matrix_isometrize", "matrix_anisometrize", + "rotated_main_axes", "pos2xyz", "xyz2pos", "ang2dir", @@ -185,7 +187,7 @@ def matrix_rotate(dim, angles): planes = rotation_planes(dim) result = np.eye(dim, dtype=np.double) for i, (angle, plane) in enumerate(zip(angles, planes)): - # angles have alternating signs to match tait bryan + # angles have alternating signs to match tait-bryan result = np.matmul( givens_rotation(dim, plane, (-1) ** i * angle), result ) @@ -207,6 +209,7 @@ def matrix_derotate(dim, angles): :class:`numpy.ndarray` Rotation matrix. """ + # derotating by taking negative angles angles = -set_angles(dim, angles) planes = rotation_planes(dim) result = np.eye(dim, dtype=np.double) @@ -300,6 +303,24 @@ def matrix_anisometrize(dim, angles, anis): ) +def rotated_main_axes(dim, angles): + """Create list of the main axis defined by the given system rotations. + + Parameters + ---------- + dim : :class:`int` + spatial dimension + angles : :class:`float` or :class:`list` + the rotation angles of the target coordinate-system + + Returns + ------- + :class:`numpy.ndarray` + Main axes of the target coordinate-system. + """ + return matrix_rotate(dim, angles).T + + # conversion ################################################################## From b552859f391ecb39fc25174baaf2e0f310c998ed Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Fri, 13 Nov 2020 12:39:15 +0100 Subject: [PATCH 25/69] CovModel: simplify Rational model --- gstools/covmodel/models.py | 131 +++++++++++++++++++------------------ 1 file changed, 66 insertions(+), 65 deletions(-) diff --git a/gstools/covmodel/models.py b/gstools/covmodel/models.py index 27ee03d3..c8066544 100644 --- a/gstools/covmodel/models.py +++ b/gstools/covmodel/models.py @@ -185,71 +185,6 @@ def calc_integral_scale(self): # noqa: D102 return self.len_rescaled -class Rational(CovModel): - r"""The rational quadratic covariance model. - - Notes - ----- - This model is given by the following correlation function: - - .. math:: - \rho(r) = - \left(1 + \frac{1}{2\alpha} \cdot - \left(s\cdot\frac{r}{\ell}\right)^2\right)^{-\alpha} - - Where the standard rescale factor is :math:`s=1`. - :math:`\alpha` is a shape parameter and should be > 0.5. - - For :math:`\alpha\to\infty` this model converges to the Gaussian model: - - .. math:: - \rho(r)= - \exp\left(-\frac{\pi}{4}\left(s\cdot\frac{r}{\ell}\right)^{2}\right) - - Other Parameters - ---------------- - alpha : :class:`float`, optional - Shape parameter. Standard range: ``(0, inf)`` - Default: ``1.0`` - """ - - def default_opt_arg(self): - """Defaults for the optional arguments. - - * ``{"alpha": 1.0}`` - - Returns - ------- - :class:`dict` - Defaults for optional arguments - """ - return {"alpha": 1.0} - - def default_opt_arg_bounds(self): - """Defaults for boundaries of the optional arguments. - - * ``{"alpha": [0.5, inf]}`` - - Returns - ------- - :class:`dict` - Boundaries for optional arguments - """ - return {"alpha": [0.5, np.inf]} - - def cor(self, h): - """Rational normalized correlation function.""" - return np.power(1 + 0.5 / self.alpha * h ** 2, -self.alpha) - - def calc_integral_scale(self): # noqa: D102 - return ( - self.len_rescaled - * np.sqrt(np.pi * self.alpha * 0.5) - * sps.gamma(self.alpha - 0.5) - / sps.gamma(self.alpha) - ) - - class Stable(CovModel): r"""The stable covariance model. @@ -428,6 +363,72 @@ def calc_integral_scale(self): # noqa: D102 ) +class Rational(CovModel): + r"""The rational quadratic covariance model. + + Notes + ----- + This model is given by the following correlation function: + + .. math:: + \rho(r) = + \left(1 + \frac{1}{\alpha} \cdot + \left(s\cdot\frac{r}{\ell}\right)^2\right)^{-\alpha} + + Where the standard rescale factor is :math:`s=1`. + :math:`\alpha` is a shape parameter and should be > 0.5. + + For :math:`\alpha\to\infty` this model converges to the Gaussian model: + + .. math:: + \rho(r)= + \exp\left(-\left(s\cdot\frac{r}{\ell}\right)^{2}\right) + + Other Parameters + ---------------- + alpha : :class:`float`, optional + Shape parameter. Standard range: ``(0, inf)`` + Default: ``1.0`` + """ + + def default_opt_arg(self): + """Defaults for the optional arguments. + + * ``{"alpha": 1.0}`` + + Returns + ------- + :class:`dict` + Defaults for optional arguments + """ + return {"alpha": 1.0} + + def default_opt_arg_bounds(self): + """Defaults for boundaries of the optional arguments. + + * ``{"alpha": [0.5, inf]}`` + + Returns + ------- + :class:`dict` + Boundaries for optional arguments + """ + return {"alpha": [0.5, np.inf]} + + def cor(self, h): + """Rational normalized correlation function.""" + return np.power(1 + h ** 2 / self.alpha, -self.alpha) + + def calc_integral_scale(self): # noqa: D102 + return ( + self.len_rescaled + * np.sqrt(np.pi * self.alpha) + * sps.gamma(self.alpha - 0.5) + / sps.gamma(self.alpha) + / 2.0 + ) + + class Linear(CovModel): r"""The bounded linear covariance model. From 305a0d4ce225ba2e8f87175908b3a06bb990d9bf Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Fri, 13 Nov 2020 12:40:45 +0100 Subject: [PATCH 26/69] Tools: make use of rotated_main_axes in CovModel and Examples --- examples/02_cov_model/README.rst | 2 +- examples/03_variogram/04_directional_2d.py | 2 +- examples/03_variogram/05_directional_3d.py | 22 ++++++++++------------ gstools/covmodel/base.py | 6 ++++++ 4 files changed, 18 insertions(+), 14 deletions(-) diff --git a/examples/02_cov_model/README.rst b/examples/02_cov_model/README.rst index 15cee254..8e69698b 100644 --- a/examples/02_cov_model/README.rst +++ b/examples/02_cov_model/README.rst @@ -24,7 +24,7 @@ Where: - :math:`s` is a scaling factor for unit conversion or normalization - :math:`\sigma^2` is the variance - :math:`n` is the nugget (subscale variance) - - :math:`\mathrm{cor}(h)` is normalized correlation function depending on + - :math:`\mathrm{cor}(h)` is the normalized correlation function depending on the non-dimensional distance :math:`h=s\cdot\frac{r}{\ell}` Depending on the normalized correlation function, all covariance models in diff --git a/examples/03_variogram/04_directional_2d.py b/examples/03_variogram/04_directional_2d.py index b1a42920..fe65264b 100755 --- a/examples/03_variogram/04_directional_2d.py +++ b/examples/03_variogram/04_directional_2d.py @@ -27,7 +27,7 @@ bins = range(0, 40, 3) bin_c, vario, cnt = gs.vario_estimate( *((x, y), field, bins), - angles=(angle, angle + np.pi / 2), # main dir. and transversal dir. + direction=gs.rotated_main_axes(dim=2, angles=angle), angles_tol=np.pi / 16, bandwidth=1.0, mesh_type="structured", diff --git a/examples/03_variogram/05_directional_3d.py b/examples/03_variogram/05_directional_3d.py index 521cacf9..bd70b81d 100755 --- a/examples/03_variogram/05_directional_3d.py +++ b/examples/03_variogram/05_directional_3d.py @@ -22,17 +22,15 @@ field = srf.structured((x, y, z)) ############################################################################### -# Here we generate the rotated coordinate system to get an impression what -# the rotation angles do. +# Here we generate the axes of the rotated coordinate system +# to get an impression what the rotation angles do. -x1, x2, x3 = (1, 0, 0), (0, 1, 0), (0, 0, 1) -ret = np.array(gs.field.tools.rotate_mesh(dim, angles, x1, x2, x3)) -dir0 = ret[:, 0] # main direction -dir1 = ret[:, 1] # first lateral direction -dir2 = ret[:, 2] # second lateral direction +# All 3 axes of the rotated coordinate-system +main_axes = gs.rotated_main_axes(dim, angles) +axis1, axis2, axis3 = main_axes ############################################################################### -# Now we estimate the variogram along the main axis. When the main axis is +# Now we estimate the variogram along the main axes. When the main axes are # unknown, one would need to sample multiple directions and look for the one # with the longest correlation length (flattest gradient). # Then check the transversal directions and so on. @@ -40,7 +38,7 @@ bins = range(0, 40, 3) bin_c, vario = gs.vario_estimate( *([x, y, z], field, bins), - direction=(dir0, dir1, dir2), + direction=main_axes, bandwidth=10, sampling_size=2000, sampling_seed=1001, @@ -58,9 +56,9 @@ srf.plot(ax=ax1) ax1.set_aspect("equal") -ax2.plot([0, dir0[0]], [0, dir0[1]], [0, dir0[2]], label="1.") -ax2.plot([0, dir1[0]], [0, dir1[1]], [0, dir1[2]], label="2.") -ax2.plot([0, dir2[0]], [0, dir2[1]], [0, dir2[2]], label="3.") +ax2.plot([0, axis1[0]], [0, axis1[1]], [0, axis1[2]], label="1.") +ax2.plot([0, axis2[0]], [0, axis2[1]], [0, axis2[2]], label="2.") +ax2.plot([0, axis3[0]], [0, axis3[1]], [0, axis3[2]], label="3.") ax2.set_xlim(-1, 1) ax2.set_ylim(-1, 1) ax2.set_zlim(-1, 1) diff --git a/gstools/covmodel/base.py b/gstools/covmodel/base.py index e512cdb5..a7f1e68a 100644 --- a/gstools/covmodel/base.py +++ b/gstools/covmodel/base.py @@ -21,6 +21,7 @@ set_angles, matrix_anisometrize, matrix_isometrize, + rotated_main_axes, ) from gstools.covmodel.tools import ( InitSubclassMeta, @@ -603,7 +604,12 @@ def anisometrize(self, pos): matrix_anisometrize(self.dim, self.angles, self.anis), pos ) + def main_axes(self): + """Axes of the rotated coordinate-system.""" + return rotated_main_axes(self.dim, self.angles) + def _get_iso_rad(self, pos): + """Isometrized radians.""" return np.linalg.norm(self.isometrize(pos), axis=0) # fitting routine From 43278f0d5dd7eafcafb1109a45c296d4713317e1 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Sat, 14 Nov 2020 02:42:55 +0100 Subject: [PATCH 27/69] CovModel: Fitting anisotropy to directional empirical variograms; variogram routines along axis added --- gstools/covmodel/base.py | 40 ++++++++++++++--- gstools/covmodel/fit.py | 96 +++++++++++++++++++++++++++++++++------- 2 files changed, 114 insertions(+), 22 deletions(-) diff --git a/gstools/covmodel/base.py b/gstools/covmodel/base.py index a7f1e68a..d5bb568f 100644 --- a/gstools/covmodel/base.py +++ b/gstools/covmodel/base.py @@ -286,6 +286,24 @@ def cor_from_correlation(self, h): # special variogram functions + def vario_axis(self, r, axis=0): + r"""Variogram along axis of anisotropy.""" + if axis == 0: + return self.variogram(r) + return self.variogram(np.abs(r) / self.anis[axis - 1]) + + def cov_axis(self, r, axis=0): + r"""Covariance along axis of anisotropy.""" + if axis == 0: + return self.covariance(r) + return self.covariance(np.abs(r) / self.anis[axis - 1]) + + def cor_axis(self, r, axis=0): + r"""Correlation along axis of anisotropy.""" + if axis == 0: + return self.correlation(r) + return self.correlation(np.abs(r) / self.anis[axis - 1]) + def vario_spatial(self, pos): r"""Spatial variogram respecting anisotropy and rotation.""" return self.variogram(self._get_iso_rad(pos)) @@ -618,6 +636,7 @@ def fit_variogram( self, x_data, y_data, + anis=True, sill=None, init_guess="default", weights=None, @@ -629,14 +648,25 @@ def fit_variogram( **para_select ): """ - Fiting the isotropic variogram-model to given data. + Fiting the variogram-model to an empirical variogram. Parameters ---------- x_data : :class:`numpy.ndarray` - The radii of the meassured variogram. + The bin-centers of the empirical variogram. y_data : :class:`numpy.ndarray` The messured variogram + If multiple are given, they are interpreted as the directional + variograms along the main axis of the associated rotated + coordinate system. + Anisotropy ratios will be estimated in that case. + anis : :class:`bool`, optional + In case of a directional variogram, you can control anisotropy + by this argument. Deselect the parameter from fitting, by setting + it "False". + You could also pass a fixed value to be set in the model. + Then the anisotropy ratios wont be altered during fitting. + Default: True sill : :class:`float` or :class:`bool`, optional Here you can provide a fixed sill for the variogram. It needs to be in a fitting range for the var and nugget bounds. @@ -742,6 +772,7 @@ def fit_variogram( model=self, x_data=x_data, y_data=y_data, + anis=anis, sill=sill, init_guess=init_guess, weights=weights, @@ -773,9 +804,8 @@ def set_arg_bounds(self, check_args=True, **kwargs): Parameters ---------- check_args : bool, optional - Whether to check if the arguments need to resetted to be in the - given bounds. In case, a propper default value will be determined - from the given bounds. + Whether to check if the arguments are in their valid bounds. + In case not, a propper default value will be determined. Default: True **kwargs Parameter name as keyword ("var", "len_scale", "nugget", ) diff --git a/gstools/covmodel/fit.py b/gstools/covmodel/fit.py index cbf1aecf..9896cd7b 100755 --- a/gstools/covmodel/fit.py +++ b/gstools/covmodel/fit.py @@ -26,6 +26,7 @@ def fit_variogram( model, x_data, y_data, + anis=True, sill=None, init_guess="default", weights=None, @@ -37,17 +38,28 @@ def fit_variogram( **para_select ): """ - Fiting the isotropic variogram-model to given data. + Fiting a variogram-model to an empirical variogram. Parameters ---------- model : :any:`CovModel` Covariance Model to fit. x_data : :class:`numpy.ndarray` - The radii of the meassured variogram. + The bin-centers of the empirical variogram. y_data : :class:`numpy.ndarray` The messured variogram - sill : :class:`float` or :class:`bool`, optional + If multiple are given, they are interpreted as the directional + variograms along the main axis of the associated rotated + coordinate system. + Anisotropy ratios will be estimated in that case. + anis : :class:`bool`, optional + In case of a directional variogram, you can control anisotropy + by this argument. Deselect the parameter from fitting, by setting + it "False". + You could also pass a fixed value to be set in the model. + Then the anisotropy ratios wont be altered during fitting. + Default: True + sill : :class:`float` or :class:`bool` or :any:`None`, optional Here you can provide a fixed sill for the variogram. It needs to be in a fitting range for the var and nugget bounds. If variance or nugget are not selected for estimation, @@ -72,7 +84,7 @@ def fit_variogram( * dict(name: val): specified value for each parameter by name Default: "default" - weights : :class:`str`, :class:`numpy.ndarray`, :class:`callable`, optional + weights : :class:`str`, :class:`numpy.ndarray`, :class:`callable`optional Weights applied to each point in the estimation. Either: * 'inv': inverse distance ``1 / (x_data + 1)`` @@ -153,29 +165,53 @@ def fit_variogram( # check method if method not in ["trf", "dogbox"]: raise ValueError("fit: method needs to be either 'trf' or 'dogbox'") + # prepare variogram data + x_data = np.array(x_data).reshape(-1) + y_data = np.array(y_data) + # check if anisotropy should be fitted or set + if not isinstance(anis, bool): + model.anis = anis + anis = False + # if multiple variograms are given, they will be interpreted + # as directional variograms along the main rotated axes of the model + is_dir_vario = False + if model.dim > 1 and x_data.size * model.dim == y_data.size: + is_dir_vario = True + # concatenate multiple variograms + x_data = np.tile(x_data, model.dim) + y_data = y_data.reshape(-1) + elif x_data.size != y_data.size: + raise ValueError( + "CovModel.fit_variogram: Wrong number of empirical variograms! " + + "Either provide only one variogram to fit an isotropic model, " + + "or directional ones for all main axes to fit anisotropy." + ) + anis &= is_dir_vario # set weights _set_weights(weights, x_data, curve_fit_kwargs) # set the lower/upper boundaries for the variogram-parameters bounds, init_guess_list = _init_curve_fit_para( - model, para, init_guess, constrain_sill, sill + model, para, init_guess, constrain_sill, sill, anis ) # create the fitting curve - curve_fit_kwargs["f"] = _get_curve(model, para, constrain_sill, sill) + curve_fit_kwargs["f"] = _get_curve( + model, para, constrain_sill, sill, anis, is_dir_vario + ) # set the remaining kwargs for curve_fit curve_fit_kwargs["bounds"] = bounds curve_fit_kwargs["p0"] = init_guess_list - curve_fit_kwargs["xdata"] = np.array(x_data) - curve_fit_kwargs["ydata"] = np.array(y_data) + curve_fit_kwargs["xdata"] = x_data + curve_fit_kwargs["ydata"] = y_data curve_fit_kwargs["loss"] = loss curve_fit_kwargs["max_nfev"] = max_eval curve_fit_kwargs["method"] = method # fit the variogram popt, pcov = curve_fit(**curve_fit_kwargs) # convert the results - fit_para = _post_fitting(model, para, popt) + fit_para = _post_fitting(model, para, popt, anis, is_dir_vario) # calculate the r2 score if wanted if return_r2: - return fit_para, pcov, _r2_score(model, x_data, y_data) + return fit_para, pcov, _r2_score(model, x_data, y_data, is_dir_vario) return fit_para, pcov @@ -247,16 +283,16 @@ def _pre_para(model, para_select, sill): def _set_weights(weights, x_data, curve_fit_kwargs): if weights is not None: if callable(weights): - weights = 1.0 / weights(np.array(x_data)) + weights = 1.0 / weights(x_data) elif weights == "inv": - weights = 1.0 + np.array(x_data) + weights = 1.0 + x_data else: weights = 1.0 / np.array(weights) curve_fit_kwargs["sigma"] = weights curve_fit_kwargs["absolute_sigma"] = True -def _get_curve(model, para, constrain_sill, sill): +def _get_curve(model, para, constrain_sill, sill, anis, is_dir_vario): """Create the curve for scipys curve_fit.""" # we need arg1, otherwise curve_fit throws an error (bug?!) def curve(x, arg1, *args): @@ -287,12 +323,20 @@ def curve(x, arg1, *args): # set var at last because of var_factor (other parameter needed) if para["var"]: model.var = var_tmp + if is_dir_vario: + if anis: + model.anis = args[1 - model.dim :] + xs = x[: x.size // model.dim] + out = np.array([], dtype=np.double) + for i in range(model.dim): + out = np.concatenate((out, model.vario_axis(xs, axis=i))) + return out return model.variogram(x) return curve -def _init_curve_fit_para(model, para, init_guess, constrain_sill, sill): +def _init_curve_fit_para(model, para, init_guess, constrain_sill, sill, anis): """Create initial guess and bounds for fitting.""" low_bounds = [] top_bounds = [] @@ -326,6 +370,13 @@ def _init_curve_fit_para(model, para, init_guess, constrain_sill, sill): para_name=opt, ) ) + if anis: + low_bounds += [0.0] * (model.dim - 1) + top_bounds += [np.inf] * (model.dim - 1) + if init_guess == "default": + init_guess_list += [1.0] * (model.dim - 1) + else: + init_guess_list += list(model.anis) return (low_bounds, top_bounds), init_guess_list @@ -347,7 +398,7 @@ def _init_guess(bounds, current, default, typ, para_name): raise ValueError("CovModel.fit: unkwon init_guess: '{}'".format(typ)) -def _post_fitting(model, para, popt): +def _post_fitting(model, para, popt, anis, is_dir_vario): """Postprocess fitting results and application to model.""" fit_para = {} para_skip = 0 @@ -369,15 +420,26 @@ def _post_fitting(model, para, popt): opt_skip += 1 else: fit_para[opt] = getattr(model, opt) + if is_dir_vario: + if anis: + model.anis = popt[1 - model.dim :] + fit_para["anis"] = model.anis # set var at last because of var_factor (other parameter needed) if para["var"]: model.var = var_tmp return fit_para -def _r2_score(model, x_data, y_data): +def _r2_score(model, x_data, y_data, is_dir_vario): """Calculate the R2 score.""" - residuals = y_data - model.variogram(x_data) + if is_dir_vario: + xs = x_data[: x_data.size // model.dim] + vario = np.array([], dtype=np.double) + for i in range(model.dim): + vario = np.concatenate((vario, model.vario_axis(xs, axis=i))) + else: + vario = model.variogram(x_data) + residuals = y_data - vario ss_res = np.sum(residuals ** 2) ss_tot = np.sum((y_data - np.mean(y_data)) ** 2) return 1.0 - (ss_res / ss_tot) From 92287fed9322366fa73dace0835b5993defc1c1f Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Sat, 14 Nov 2020 02:43:56 +0100 Subject: [PATCH 28/69] CovModel: update plotting routines: add axis routines; forward kwargs --- gstools/covmodel/plot.py | 101 ++++++++++++++++++++++++++++++--------- 1 file changed, 79 insertions(+), 22 deletions(-) diff --git a/gstools/covmodel/plot.py b/gstools/covmodel/plot.py index cc9535f7..66084aa5 100644 --- a/gstools/covmodel/plot.py +++ b/gstools/covmodel/plot.py @@ -10,6 +10,9 @@ plot_variogram plot_covariance plot_correlation + plot_vario_axis + plot_cov_axis + plot_cor_axis plot_vario_spatial plot_cov_spatial plot_cor_spatial @@ -17,7 +20,7 @@ plot_spectral_density plot_spectral_rad_pdf """ -# pylint: disable=C0103 +# pylint: disable=C0103, C0415 import numpy as np import gstools @@ -27,6 +30,9 @@ "plot_variogram", "plot_covariance", "plot_correlation", + "plot_vario_axis", + "plot_cov_axis", + "plot_cor_axis", "plot_vario_spatial", "plot_cov_spatial", "plot_cor_spatial", @@ -111,96 +117,147 @@ def plot_cor_spatial( def plot_variogram( - model, x_min=0.0, x_max=None, fig=None, ax=None + model, x_min=0.0, x_max=None, fig=None, ax=None, **kwargs ): # pragma: no cover """Plot variogram of a given CovModel.""" fig, ax = _get_fig_ax(fig, ax) if x_max is None: x_max = 3 * model.len_scale x_s = np.linspace(x_min, x_max) - ax.plot(x_s, model.variogram(x_s), label=model.name + " variogram") + kwargs.setdefault("label", model.name + " variogram") + ax.plot(x_s, model.variogram(x_s), **kwargs) ax.legend() fig.show() return ax def plot_covariance( - model, x_min=0.0, x_max=None, fig=None, ax=None + model, x_min=0.0, x_max=None, fig=None, ax=None, **kwargs ): # pragma: no cover """Plot covariance of a given CovModel.""" fig, ax = _get_fig_ax(fig, ax) if x_max is None: x_max = 3 * model.len_scale x_s = np.linspace(x_min, x_max) - ax.plot(x_s, model.covariance(x_s), label=model.name + " covariance") + kwargs.setdefault("label", model.name + " covariance") + ax.plot(x_s, model.covariance(x_s), **kwargs) ax.legend() fig.show() return ax def plot_correlation( - model, x_min=0.0, x_max=None, fig=None, ax=None + model, x_min=0.0, x_max=None, fig=None, ax=None, **kwargs ): # pragma: no cover """Plot correlation function of a given CovModel.""" fig, ax = _get_fig_ax(fig, ax) if x_max is None: x_max = 3 * model.len_scale x_s = np.linspace(x_min, x_max) - ax.plot(x_s, model.correlation(x_s), label=model.name + " correlation") + kwargs.setdefault("label", model.name + " correlation") + ax.plot(x_s, model.correlation(x_s), **kwargs) + ax.legend() + fig.show() + return ax + + +def plot_vario_axis( + model, axis=0, x_min=0.0, x_max=None, fig=None, ax=None, **kwargs +): # pragma: no cover + """Plot variogram of a given CovModel.""" + fig, ax = _get_fig_ax(fig, ax) + if x_max is None: + x_max = 3 * model.len_scale + x_s = np.linspace(x_min, x_max) + kwargs.setdefault( + "label", model.name + " variogram on axis {}".format(axis) + ) + ax.plot(x_s, model.vario_axis(x_s, axis), **kwargs) + ax.legend() + fig.show() + return ax + + +def plot_cov_axis( + model, axis=0, x_min=0.0, x_max=None, fig=None, ax=None, **kwargs +): # pragma: no cover + """Plot variogram of a given CovModel.""" + fig, ax = _get_fig_ax(fig, ax) + if x_max is None: + x_max = 3 * model.len_scale + x_s = np.linspace(x_min, x_max) + kwargs.setdefault( + "label", model.name + " covariance on axis {}".format(axis) + ) + ax.plot(x_s, model.cov_axis(x_s, axis), **kwargs) + ax.legend() + fig.show() + return ax + + +def plot_cor_axis( + model, axis=0, x_min=0.0, x_max=None, fig=None, ax=None, **kwargs +): # pragma: no cover + """Plot variogram of a given CovModel.""" + fig, ax = _get_fig_ax(fig, ax) + if x_max is None: + x_max = 3 * model.len_scale + x_s = np.linspace(x_min, x_max) + kwargs.setdefault( + "label", model.name + " correlation on axis {}".format(axis) + ) + ax.plot(x_s, model.cor_axis(x_s, axis), **kwargs) ax.legend() fig.show() return ax def plot_spectrum( - model, x_min=0.0, x_max=None, fig=None, ax=None + model, x_min=0.0, x_max=None, fig=None, ax=None, **kwargs ): # pragma: no cover """Plot specturm of a given CovModel.""" fig, ax = _get_fig_ax(fig, ax) if x_max is None: x_max = 3 / model.len_scale x_s = np.linspace(x_min, x_max) - ax.plot( - x_s, - model.spectrum(x_s), - label=model.name + " " + str(model.dim) + "D spectrum", + kwargs.setdefault( + "label", model.name + " " + str(model.dim) + "D spectrum" ) + ax.plot(x_s, model.spectrum(x_s), **kwargs) ax.legend() fig.show() return ax def plot_spectral_density( - model, x_min=0.0, x_max=None, fig=None, ax=None + model, x_min=0.0, x_max=None, fig=None, ax=None, **kwargs ): # pragma: no cover """Plot spectral density of a given CovModel.""" fig, ax = _get_fig_ax(fig, ax) if x_max is None: x_max = 3 / model.len_scale x_s = np.linspace(x_min, x_max) - ax.plot( - x_s, - model.spectral_density(x_s), - label=model.name + " " + str(model.dim) + "D spectral-density", + kwargs.setdefault( + "label", model.name + " " + str(model.dim) + "D spectral-density" ) + ax.plot(x_s, model.spectral_density(x_s), **kwargs) ax.legend() fig.show() return ax def plot_spectral_rad_pdf( - model, x_min=0.0, x_max=None, fig=None, ax=None + model, x_min=0.0, x_max=None, fig=None, ax=None, **kwargs ): # pragma: no cover """Plot radial spectral pdf of a given CovModel.""" fig, ax = _get_fig_ax(fig, ax) if x_max is None: x_max = 3 / model.len_scale x_s = np.linspace(x_min, x_max) - ax.plot( - x_s, - model.spectral_rad_pdf(x_s), - label=model.name + " " + str(model.dim) + "D spectral-rad-pdf", + kwargs.setdefault( + "label", model.name + " " + str(model.dim) + "D spectral-rad-pdf" ) + ax.plot(x_s, model.spectral_rad_pdf(x_s), **kwargs) ax.legend() fig.show() return ax From bcd984e6d1866af9c3e767daa268803f8a9ad2db Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Sat, 14 Nov 2020 02:54:09 +0100 Subject: [PATCH 29/69] Examples: demonstrate fitting of anisotropy --- examples/03_variogram/04_directional_2d.py | 33 +++++++++++++------- examples/03_variogram/05_directional_3d.py | 35 ++++++++++++++++------ 2 files changed, 49 insertions(+), 19 deletions(-) diff --git a/examples/03_variogram/04_directional_2d.py b/examples/03_variogram/04_directional_2d.py index fe65264b..7d5dd55d 100755 --- a/examples/03_variogram/04_directional_2d.py +++ b/examples/03_variogram/04_directional_2d.py @@ -1,9 +1,11 @@ """ -Directional variogram estimation in 2D --------------------------------------- +Directional variogram estimation and fitting in 2D +-------------------------------------------------- In this example, we demonstrate how to estimate a directional variogram by setting the direction angles in 2D. + +Afterwards we will fit a model to this estimated variogram and show the result. """ import numpy as np import gstools as gs @@ -20,29 +22,40 @@ ############################################################################### # Now we are going to estimate a directional variogram with an angular -# tolerance of 11.25 degree and a bandwith of 1. -# We provide the rotation angle of the covariance model and the orthogonal -# direction by adding 90 degree. +# tolerance of 11.25 degree and a bandwith of 8. -bins = range(0, 40, 3) -bin_c, vario, cnt = gs.vario_estimate( +bins = range(0, 40, 2) +bin_center, dir_vario, counts = gs.vario_estimate( *((x, y), field, bins), direction=gs.rotated_main_axes(dim=2, angles=angle), angles_tol=np.pi / 16, - bandwidth=1.0, + bandwidth=8, mesh_type="structured", return_counts=True, ) +############################################################################### +# Afterwards we can use the estimated variogram to fit a model to it: + +print("Original:") +print(model) +model.fit_variogram(bin_center, dir_vario) +print("Fitted:") +print(model) + ############################################################################### # Plotting. fig, (ax1, ax2) = plt.subplots(1, 2, figsize=[10, 5]) -ax1.plot(bin_c, vario[0], label="emp. vario: pi/8") -ax1.plot(bin_c, vario[1], label="emp. vario: pi*5/8") +ax1.plot(bin_center, dir_vario[0], label="emp. vario: pi/8") +ax1.plot(bin_center, dir_vario[1], label="emp. vario: pi*5/8") ax1.legend(loc="lower right") +model.plot("vario_axis", axis=0, ax=ax1, x_max=40, label="fit on axis 0") +model.plot("vario_axis", axis=1, ax=ax1, x_max=40, label="fit on axis 1") +ax1.set_title("Fitting an anisotropic model") + srf.plot(ax=ax2) ax2.set_aspect("equal") diff --git a/examples/03_variogram/05_directional_3d.py b/examples/03_variogram/05_directional_3d.py index bd70b81d..0bd92775 100755 --- a/examples/03_variogram/05_directional_3d.py +++ b/examples/03_variogram/05_directional_3d.py @@ -1,9 +1,11 @@ """ -Directional variogram estimation in 3D --------------------------------------- +Directional variogram estimation and fitting in 3D +-------------------------------------------------- In this example, we demonstrate how to estimate a directional variogram by setting the estimation directions in 3D. + +Afterwards we will fit a model to this estimated variogram and show the result. """ import numpy as np import gstools as gs @@ -36,7 +38,7 @@ # Then check the transversal directions and so on. bins = range(0, 40, 3) -bin_c, vario = gs.vario_estimate( +bin_center, dir_vario = gs.vario_estimate( *([x, y, z], field, bins), direction=main_axes, bandwidth=10, @@ -45,6 +47,16 @@ mesh_type="structured" ) +############################################################################### +# Afterwards we can use the estimated variogram to fit a model to it. +# Note, that the rotation angles need to be set beforehand. + +print("Original:") +print(model) +model.fit_variogram(bin_center, dir_vario) +print("Fitted:") +print(model) + ############################################################################### # Plotting. @@ -56,9 +68,9 @@ srf.plot(ax=ax1) ax1.set_aspect("equal") -ax2.plot([0, axis1[0]], [0, axis1[1]], [0, axis1[2]], label="1.") -ax2.plot([0, axis2[0]], [0, axis2[1]], [0, axis2[2]], label="2.") -ax2.plot([0, axis3[0]], [0, axis3[1]], [0, axis3[2]], label="3.") +ax2.plot([0, axis1[0]], [0, axis1[1]], [0, axis1[2]], label="0.") +ax2.plot([0, axis2[0]], [0, axis2[1]], [0, axis2[2]], label="1.") +ax2.plot([0, axis3[0]], [0, axis3[1]], [0, axis3[2]], label="2.") ax2.set_xlim(-1, 1) ax2.set_ylim(-1, 1) ax2.set_zlim(-1, 1) @@ -68,8 +80,13 @@ ax2.set_title("Tait-Bryan main axis") ax2.legend(loc="lower left") -ax3.plot(bin_c, vario[0], label="1. axis") -ax3.plot(bin_c, vario[1], label="2. axis") -ax3.plot(bin_c, vario[2], label="3. axis") +ax3.plot(bin_center, dir_vario[0], label="0. axis") +ax3.plot(bin_center, dir_vario[1], label="1. axis") +ax3.plot(bin_center, dir_vario[2], label="2. axis") +model.plot("vario_axis", axis=0, ax=ax3, label="fit on axis 0") +model.plot("vario_axis", axis=1, ax=ax3, label="fit on axis 1") +model.plot("vario_axis", axis=2, ax=ax3, label="fit on axis 2") +ax3.set_title("Fitting an anisotropic model") ax3.legend() + plt.show() From 4bdcecaf107514ccfa789dc6db0ebfe9672ab6e7 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Sat, 14 Nov 2020 13:22:40 +0100 Subject: [PATCH 30/69] CovModel: fitting with weights respects directional variogram now; restructure+cleanup fit routine --- gstools/covmodel/fit.py | 153 ++++++++++++++++++++++------------------ 1 file changed, 83 insertions(+), 70 deletions(-) diff --git a/gstools/covmodel/fit.py b/gstools/covmodel/fit.py index 9896cd7b..eb92bda4 100755 --- a/gstools/covmodel/fit.py +++ b/gstools/covmodel/fit.py @@ -159,36 +159,21 @@ def fit_variogram( The fitted parameters will be instantly set in the model. """ # preprocess selected parameters - para, sill, constrain_sill = _pre_para(model, para_select, sill) + para, sill, constrain_sill, anis = _pre_para( + model, para_select, sill, anis + ) # check curve_fit kwargs curve_fit_kwargs = {} if curve_fit_kwargs is None else curve_fit_kwargs # check method if method not in ["trf", "dogbox"]: raise ValueError("fit: method needs to be either 'trf' or 'dogbox'") # prepare variogram data - x_data = np.array(x_data).reshape(-1) - y_data = np.array(y_data) - # check if anisotropy should be fitted or set - if not isinstance(anis, bool): - model.anis = anis - anis = False - # if multiple variograms are given, they will be interpreted - # as directional variograms along the main rotated axes of the model - is_dir_vario = False - if model.dim > 1 and x_data.size * model.dim == y_data.size: - is_dir_vario = True - # concatenate multiple variograms - x_data = np.tile(x_data, model.dim) - y_data = y_data.reshape(-1) - elif x_data.size != y_data.size: - raise ValueError( - "CovModel.fit_variogram: Wrong number of empirical variograms! " - + "Either provide only one variogram to fit an isotropic model, " - + "or directional ones for all main axes to fit anisotropy." - ) + # => concatenate directional variograms to have a 1D array for x and y + x_data, y_data, is_dir_vario = _check_vario(model, x_data, y_data) + # only fit anisotropy if a directional variogram was given anis &= is_dir_vario # set weights - _set_weights(weights, x_data, curve_fit_kwargs) + _set_weights(model, weights, x_data, curve_fit_kwargs, is_dir_vario) # set the lower/upper boundaries for the variogram-parameters bounds, init_guess_list = _init_curve_fit_para( model, para, init_guess, constrain_sill, sill, anis @@ -215,7 +200,7 @@ def fit_variogram( return fit_para, pcov -def _pre_para(model, para_select, sill): +def _pre_para(model, para_select, sill, anis): """Preprocess selected parameters.""" var_last = False for par in para_select: @@ -277,65 +262,49 @@ def _pre_para(model, para_select, sill): para.update({opt: True for opt in model.opt_arg}) # now deselect unwanted parameters para.update(para_select) - return para, sill, constrain_sill + # check if anisotropy should be fitted or set + if not isinstance(anis, bool): + model.anis = anis + anis = False + return para, sill, constrain_sill, anis + + +def _check_vario(model, x_data, y_data): + # prepare variogram data + x_data = np.array(x_data).reshape(-1) + y_data = np.array(y_data).reshape(-1) + # if multiple variograms are given, they will be interpreted + # as directional variograms along the main rotated axes of the model + is_dir_vario = False + if model.dim > 1 and x_data.size * model.dim == y_data.size: + is_dir_vario = True + # concatenate multiple variograms + x_data = np.tile(x_data, model.dim) + elif x_data.size != y_data.size: + raise ValueError( + "CovModel.fit_variogram: Wrong number of empirical variograms! " + + "Either provide only one variogram to fit an isotropic model, " + + "or directional ones for all main axes to fit anisotropy." + ) + return x_data, y_data, is_dir_vario -def _set_weights(weights, x_data, curve_fit_kwargs): +def _set_weights(model, weights, x_data, curve_fit_kwargs, is_dir_vario): if weights is not None: if callable(weights): weights = 1.0 / weights(x_data) elif weights == "inv": weights = 1.0 + x_data + elif is_dir_vario: + if weights.size * model.dim == x_data.size: + weights = np.tile(weights, model.dim) + weights = 1.0 / np.array(weights).reshape(-1) else: - weights = 1.0 / np.array(weights) + weights = 1.0 / np.array(weights).reshape(-1) curve_fit_kwargs["sigma"] = weights curve_fit_kwargs["absolute_sigma"] = True -def _get_curve(model, para, constrain_sill, sill, anis, is_dir_vario): - """Create the curve for scipys curve_fit.""" - # we need arg1, otherwise curve_fit throws an error (bug?!) - def curve(x, arg1, *args): - """Adapted Variogram function.""" - args = (arg1,) + args - para_skip = 0 - opt_skip = 0 - if para["var"]: - var_tmp = args[para_skip] - if constrain_sill: - nugget_tmp = sill - var_tmp - # punishment, if resulting nugget out of range for fixed sill - if check_arg_in_bounds(model, "nugget", nugget_tmp) > 0: - return np.full_like(x, np.inf) - # nugget estimation deselected in this case - model.nugget = nugget_tmp - para_skip += 1 - if para["len_scale"]: - model.len_scale = args[para_skip] - para_skip += 1 - if para["nugget"]: - model.nugget = args[para_skip] - para_skip += 1 - for opt in model.opt_arg: - if para[opt]: - setattr(model, opt, args[para_skip + opt_skip]) - opt_skip += 1 - # set var at last because of var_factor (other parameter needed) - if para["var"]: - model.var = var_tmp - if is_dir_vario: - if anis: - model.anis = args[1 - model.dim :] - xs = x[: x.size // model.dim] - out = np.array([], dtype=np.double) - for i in range(model.dim): - out = np.concatenate((out, model.vario_axis(xs, axis=i))) - return out - return model.variogram(x) - - return curve - - def _init_curve_fit_para(model, para, init_guess, constrain_sill, sill, anis): """Create initial guess and bounds for fitting.""" low_bounds = [] @@ -398,6 +367,50 @@ def _init_guess(bounds, current, default, typ, para_name): raise ValueError("CovModel.fit: unkwon init_guess: '{}'".format(typ)) +def _get_curve(model, para, constrain_sill, sill, anis, is_dir_vario): + """Create the curve for scipys curve_fit.""" + # we need arg1, otherwise curve_fit throws an error (bug?!) + def curve(x, arg1, *args): + """Adapted Variogram function.""" + args = (arg1,) + args + para_skip = 0 + opt_skip = 0 + if para["var"]: + var_tmp = args[para_skip] + if constrain_sill: + nugget_tmp = sill - var_tmp + # punishment, if resulting nugget out of range for fixed sill + if check_arg_in_bounds(model, "nugget", nugget_tmp) > 0: + return np.full_like(x, np.inf) + # nugget estimation deselected in this case + model.nugget = nugget_tmp + para_skip += 1 + if para["len_scale"]: + model.len_scale = args[para_skip] + para_skip += 1 + if para["nugget"]: + model.nugget = args[para_skip] + para_skip += 1 + for opt in model.opt_arg: + if para[opt]: + setattr(model, opt, args[para_skip + opt_skip]) + opt_skip += 1 + # set var at last because of var_factor (other parameter needed) + if para["var"]: + model.var = var_tmp + if is_dir_vario: + if anis: + model.anis = args[1 - model.dim :] + xs = x[: x.size // model.dim] + out = np.array([], dtype=np.double) + for i in range(model.dim): + out = np.concatenate((out, model.vario_axis(xs, axis=i))) + return out + return model.variogram(x) + + return curve + + def _post_fitting(model, para, popt, anis, is_dir_vario): """Postprocess fitting results and application to model.""" fit_para = {} From 316ca186a109b801f7264c805b9b404bff004943 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Sat, 14 Nov 2020 13:23:28 +0100 Subject: [PATCH 31/69] CovModel: always use floats for 'angles' and 'anis' --- gstools/tools/geometric.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/gstools/tools/geometric.py b/gstools/tools/geometric.py index 79f81d9d..168df1a6 100644 --- a/gstools/tools/geometric.py +++ b/gstools/tools/geometric.py @@ -68,7 +68,8 @@ def set_angles(dim, angles): ----- If too few angles are given, they are filled up with `0`. """ - out_angles = np.atleast_1d(angles)[: no_of_angles(dim)] + out_angles = np.array(angles, dtype=np.double) + out_angles = np.atleast_1d(out_angles)[: no_of_angles(dim)] # fill up the rotation angle array with zeros out_angles = np.pad( out_angles, @@ -98,7 +99,8 @@ def set_anis(dim, anis): ----- If too few anisotrpy ratios are given, they are filled up with `1`. """ - out_anis = np.atleast_1d(anis)[: dim - 1] + out_anis = np.array(anis, dtype=np.double) + out_anis = np.atleast_1d(out_anis)[: dim - 1] if len(out_anis) < dim - 1: # fill up the anisotropies with ones, such that len()==dim-1 out_anis = np.pad( From 2bea823eb232ca00dce025955da565949efc9727 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Sat, 14 Nov 2020 13:24:13 +0100 Subject: [PATCH 32/69] CovModel: always use floats for 'len_scale' --- gstools/covmodel/tools.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/gstools/covmodel/tools.py b/gstools/covmodel/tools.py index ed1ebfff..93a7d381 100644 --- a/gstools/covmodel/tools.py +++ b/gstools/covmodel/tools.py @@ -127,7 +127,8 @@ def set_len_anis(dim, len_scale, anis): If to few ``anis`` values are given, the first dimensions will be filled up with 1. (eg. anis=[e] in 3D is equal to anis=[1, e]) """ - ls_tmp = np.atleast_1d(len_scale)[:dim] + ls_tmp = np.array(len_scale, dtype=np.double) + ls_tmp = np.atleast_1d(ls_tmp)[:dim] # use just one length scale (x-direction) out_len_scale = ls_tmp[0] # set the anisotropies in y- and z-direction according to the input From 2fd07e30f695a62eb7f3125934cbf9c901e011a7 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Sat, 14 Nov 2020 13:25:26 +0100 Subject: [PATCH 33/69] Tools: add a string formatter to pretty-print lists of floats --- gstools/tools/misc.py | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) create mode 100755 gstools/tools/misc.py diff --git a/gstools/tools/misc.py b/gstools/tools/misc.py new file mode 100755 index 00000000..5c908b8a --- /dev/null +++ b/gstools/tools/misc.py @@ -0,0 +1,16 @@ +# -*- coding: utf-8 -*- +""" +GStools subpackage providing miscellaneous tools. + +.. currentmodule:: gstools.tools.misc + +The following functions are provided + +.. autosummary:: + list_format +""" + + +def list_format(lst, prec): + """Format a list of floats.""" + return "[{}]".format(", ".join(f"{x:.{prec}}" for x in lst)) From b84439793c4ca9ca388647f498ab4e7c0b8a4591 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Sat, 14 Nov 2020 13:26:26 +0100 Subject: [PATCH 34/69] CovModel: always use floats for opt_args; better formatting of CovModel string representation --- gstools/covmodel/base.py | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/gstools/covmodel/base.py b/gstools/covmodel/base.py index d5bb568f..355cd0f6 100644 --- a/gstools/covmodel/base.py +++ b/gstools/covmodel/base.py @@ -23,6 +23,7 @@ matrix_isometrize, rotated_main_axes, ) +from gstools.tools.misc import list_format from gstools.covmodel.tools import ( InitSubclassMeta, rad_fac, @@ -175,7 +176,7 @@ def __init__( + "the class. It could not be added to the model" ) # Magic happens here - setattr(self, opt_name, opt_arg[opt_name]) + setattr(self, opt_name, float(opt_arg[opt_name])) # set standard boundaries for variance, len_scale, nugget and opt_arg bounds = self.default_arg_bounds() @@ -1266,17 +1267,17 @@ def __repr__(self): """Return String representation.""" opt_str = "" for opt in self.opt_arg: - opt_str += ", " + opt + "={}".format(getattr(self, opt)) + opt_str += ", " + opt + "={.3}".format(getattr(self, opt)) return ( - "{0}(dim={1}, var={2}, len_scale={3}, " - "nugget={4}, anis={5}, angles={6}".format( + "{0}(dim={1}, var={2:.3}, len_scale={3:.3}, " + "nugget={4:.3}, anis={5}, angles={6}".format( self.name, self.dim, self.var, self.len_scale, self.nugget, - self.anis, - self.angles, + list_format(self.anis, 3), + list_format(self.angles, 3), ) + opt_str + ")" From 621db8a3fcd57d7c52665ddd977a8be2348dc887 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Sat, 14 Nov 2020 13:27:12 +0100 Subject: [PATCH 35/69] Examples: demonstrate using bin-counts as weights during variogram fitting --- examples/03_variogram/05_directional_3d.py | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) diff --git a/examples/03_variogram/05_directional_3d.py b/examples/03_variogram/05_directional_3d.py index 0bd92775..d2358a8f 100755 --- a/examples/03_variogram/05_directional_3d.py +++ b/examples/03_variogram/05_directional_3d.py @@ -38,22 +38,28 @@ # Then check the transversal directions and so on. bins = range(0, 40, 3) -bin_center, dir_vario = gs.vario_estimate( +bin_center, dir_vario, counts = gs.vario_estimate( *([x, y, z], field, bins), direction=main_axes, bandwidth=10, sampling_size=2000, sampling_seed=1001, - mesh_type="structured" + mesh_type="structured", + return_counts=True, ) ############################################################################### # Afterwards we can use the estimated variogram to fit a model to it. # Note, that the rotation angles need to be set beforehand. +# +# We can use the `counts` of data pairs per bin as weights for the fitting +# routines to give more attention to areas where more data was available. +# In order to not introduce to much offset at the origin, we disable +# fitting the nugget. print("Original:") print(model) -model.fit_variogram(bin_center, dir_vario) +model.fit_variogram(bin_center, dir_vario, weights=counts, nugget=False) print("Fitted:") print(model) From 933aea0257d4b46b08de1cb976046fa38bc3e871 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Sat, 14 Nov 2020 13:39:02 +0100 Subject: [PATCH 36/69] tools/misc: f-string not working in py35; use format --- gstools/tools/misc.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/gstools/tools/misc.py b/gstools/tools/misc.py index 5c908b8a..93ef4cc0 100755 --- a/gstools/tools/misc.py +++ b/gstools/tools/misc.py @@ -13,4 +13,6 @@ def list_format(lst, prec): """Format a list of floats.""" - return "[{}]".format(", ".join(f"{x:.{prec}}" for x in lst)) + return "[{}]".format( + ", ".join("{x:.{p}}".format(x=x, p=prec) for x in lst) + ) From 916867e7d177886b3c94fe9709c44f0a0ba6faed Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Sat, 14 Nov 2020 13:45:38 +0100 Subject: [PATCH 37/69] CovModel: use precision as privat attribute to set printing format --- gstools/covmodel/base.py | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/gstools/covmodel/base.py b/gstools/covmodel/base.py index 355cd0f6..9175b2f9 100644 --- a/gstools/covmodel/base.py +++ b/gstools/covmodel/base.py @@ -208,6 +208,8 @@ def __init__( self.check_arg_bounds() # additional checks for the optional arguments (provided by user) self.check_opt_arg() + # precision for printing + self._prec = 3 # one of these functions needs to be overridden def __init_subclass__(cls): @@ -1267,10 +1269,12 @@ def __repr__(self): """Return String representation.""" opt_str = "" for opt in self.opt_arg: - opt_str += ", " + opt + "={.3}".format(getattr(self, opt)) + opt_str += ( + ", " + opt + "={0:.{1}}".format(getattr(self, opt), self._prec) + ) return ( - "{0}(dim={1}, var={2:.3}, len_scale={3:.3}, " - "nugget={4:.3}, anis={5}, angles={6}".format( + "{0}(dim={1}, var={2:.{p}}, len_scale={3:.{p}}, " + "nugget={4:.{p}}, anis={5}, angles={6}".format( self.name, self.dim, self.var, @@ -1278,6 +1282,7 @@ def __repr__(self): self.nugget, list_format(self.anis, 3), list_format(self.angles, 3), + p=self._prec, ) + opt_str + ")" From 373251dbb7c01fafa57c834f92fd8b8e2805e431 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Sat, 14 Nov 2020 13:45:57 +0100 Subject: [PATCH 38/69] Examples: pimp variogram plot --- examples/03_variogram/05_directional_3d.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/examples/03_variogram/05_directional_3d.py b/examples/03_variogram/05_directional_3d.py index d2358a8f..9b981a05 100755 --- a/examples/03_variogram/05_directional_3d.py +++ b/examples/03_variogram/05_directional_3d.py @@ -89,9 +89,9 @@ ax3.plot(bin_center, dir_vario[0], label="0. axis") ax3.plot(bin_center, dir_vario[1], label="1. axis") ax3.plot(bin_center, dir_vario[2], label="2. axis") -model.plot("vario_axis", axis=0, ax=ax3, label="fit on axis 0") -model.plot("vario_axis", axis=1, ax=ax3, label="fit on axis 1") -model.plot("vario_axis", axis=2, ax=ax3, label="fit on axis 2") +model.plot("vario_axis", axis=0, ax=ax3, label="fit on axis 0", linestyle=":") +model.plot("vario_axis", axis=1, ax=ax3, label="fit on axis 1", linestyle=":") +model.plot("vario_axis", axis=2, ax=ax3, label="fit on axis 2", linestyle=":") ax3.set_title("Fitting an anisotropic model") ax3.legend() From 24b4f8304308570cd2a1d88f0e600e693e18b3ac Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Sat, 14 Nov 2020 14:01:04 +0100 Subject: [PATCH 39/69] Field: better printing format --- gstools/field/base.py | 4 +++- gstools/field/srf.py | 4 ++-- 2 files changed, 5 insertions(+), 3 deletions(-) diff --git a/gstools/field/base.py b/gstools/field/base.py index b560e947..c3349e0a 100755 --- a/gstools/field/base.py +++ b/gstools/field/base.py @@ -393,4 +393,6 @@ def __str__(self): def __repr__(self): """Return String representation.""" - return "Field(model={0}, mean={1})".format(self.model, self.mean) + return "Field(model={0}, mean={1:.{p}})".format( + self.model, self.mean, p=self.model._prec + ) diff --git a/gstools/field/srf.py b/gstools/field/srf.py index 6fee323b..befba22c 100644 --- a/gstools/field/srf.py +++ b/gstools/field/srf.py @@ -283,6 +283,6 @@ def upscaling(self, upscaling): def __repr__(self): """Return String representation.""" - return "SRF(model={0}, mean={1}, generator={2}".format( - self.model, self.mean, self.generator + return "SRF(model={0}, mean={1:.{p}}, generator={2}".format( + self.model, self.mean, self.generator, p=self.model._prec ) From 3f1c0bad13a2a2631965feefe8ee1e2909d96c64 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Sat, 14 Nov 2020 14:23:13 +0100 Subject: [PATCH 40/69] CovModel.check_arg_in_bounds: allow checking lists --- gstools/covmodel/tools.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/gstools/covmodel/tools.py b/gstools/covmodel/tools.py index 93a7d381..1c4fa70a 100644 --- a/gstools/covmodel/tools.py +++ b/gstools/covmodel/tools.py @@ -184,20 +184,21 @@ def check_arg_in_bounds(model, arg, val=None): raise ValueError("check bounds: unknown argument: {}".format(arg)) bnd = list(model.arg_bounds[arg]) val = getattr(model, arg) if val is None else val + val = np.array(val) error_case = 0 if len(bnd) == 2: bnd.append("cc") # use closed intervals by default if bnd[2][0] == "c": - if val < bnd[0]: + if np.any(val < bnd[0]): error_case = 1 else: - if val <= bnd[0]: + if np.any(val <= bnd[0]): error_case = 2 if bnd[2][1] == "c": - if val > bnd[1]: + if np.any(val > bnd[1]): error_case = 3 else: - if val >= bnd[1]: + if np.any(val >= bnd[1]): error_case = 4 return error_case From 602a90a3a7220cc3317e58ae685bf69980cc364a Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Sun, 15 Nov 2020 17:14:00 +0100 Subject: [PATCH 41/69] CovModel: add bounds for 'anis' since we can fit it now --- gstools/covmodel/base.py | 37 ++++++++++++++++++++++++++++++++++++- gstools/covmodel/tools.py | 2 +- 2 files changed, 37 insertions(+), 2 deletions(-) diff --git a/gstools/covmodel/base.py b/gstools/covmodel/base.py index 9175b2f9..54f608b1 100644 --- a/gstools/covmodel/base.py +++ b/gstools/covmodel/base.py @@ -142,6 +142,7 @@ def __init__( self._var_bounds = None self._len_scale_bounds = None self._nugget_bounds = None + self._anis_bounds = None self._opt_arg_bounds = {} # SFT class will be created within dim.setter but needs hankel_kw self.hankel_kw = hankel_kw @@ -798,6 +799,7 @@ def default_arg_bounds(self): "var": (0.0, np.inf, "oo"), "len_scale": (0.0, np.inf, "oo"), "nugget": (0.0, np.inf, "co"), + "anis": (0.0, np.inf, "oo"), } return res @@ -838,12 +840,18 @@ def set_arg_bounds(self, check_args=True, **kwargs): self.len_scale_bounds = kwargs[arg] elif arg == "nugget": self.nugget_bounds = kwargs[arg] + elif arg == "anis": + self.anis_bounds = kwargs[arg] else: raise ValueError( "set_arg_bounds: unknown argument '{}'".format(arg) ) if check_args and check_arg_in_bounds(self, arg) > 0: - setattr(self, arg, default_arg_from_bounds(kwargs[arg])) + def_arg = default_arg_from_bounds(kwargs[arg]) + if arg == "anis": + setattr(self, arg, [def_arg] * (self.dim - 1)) + else: + setattr(self, arg, def_arg) # set var last like allways if var_bnds: self.var_bounds = var_bnds @@ -954,6 +962,32 @@ def nugget_bounds(self, bounds): ) self._nugget_bounds = bounds + @property + def anis_bounds(self): + """:class:`list`: Bounds for the nugget. + + Notes + ----- + Is a list of 2 or 3 values: + + * ``[a, b]`` or + * ``[a, b, ]`` + + is one of ``"oo"``, ``"cc"``, ``"oc"`` or ``"co"`` + to define if the bounds are open ("o") or closed ("c"). + """ + return self._anis_bounds + + @anis_bounds.setter + def anis_bounds(self, bounds): + if not check_bounds(bounds): + raise ValueError( + "Given bounds for 'anis' are not " + + "valid, got: " + + str(bounds) + ) + self._anis_bounds = bounds + @property def opt_arg_bounds(self): """:class:`dict`: Bounds for the optional arguments. @@ -988,6 +1022,7 @@ def arg_bounds(self): "var": self.var_bounds, "len_scale": self.len_scale_bounds, "nugget": self.nugget_bounds, + "anis": self.anis_bounds, } res.update(self.opt_arg_bounds) return res diff --git a/gstools/covmodel/tools.py b/gstools/covmodel/tools.py index 1c4fa70a..5e7752c8 100644 --- a/gstools/covmodel/tools.py +++ b/gstools/covmodel/tools.py @@ -142,7 +142,7 @@ def set_len_anis(dim, len_scale, anis): out_anis = np.zeros(dim - 1, dtype=np.double) for i in range(1, dim): out_anis[i - 1] = ls_tmp[i] / ls_tmp[0] - + # sanity check for ani in out_anis: if not ani > 0.0: raise ValueError( From 3e2680f0fc9a80243993a6709f25727b07649ac9 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Sun, 15 Nov 2020 17:32:11 +0100 Subject: [PATCH 42/69] CovModel: add arg_list attribute to get list of values; also for isotropic arguments only --- gstools/covmodel/base.py | 67 +++++++++++++++++++++------------------- 1 file changed, 36 insertions(+), 31 deletions(-) diff --git a/gstools/covmodel/base.py b/gstools/covmodel/base.py index 54f608b1..60c40d6f 100644 --- a/gstools/covmodel/base.py +++ b/gstools/covmodel/base.py @@ -44,8 +44,6 @@ class AttributeWarning(UserWarning): """Attribute warning for CovModel class.""" - pass - class CovModel(metaclass=InitSubclassMeta): r"""Base class for the GSTools covariance models. @@ -320,34 +318,22 @@ def cor_spatial(self, pos): r"""Spatial correlation respecting anisotropy and rotation.""" return self.correlation(self._get_iso_rad(pos)) - def cov_nugget(self, r): - r"""Covariance of the model respecting the nugget at r=0. - - Given by: :math:`C\left(r\right)= - \sigma^2\cdot\rho\left(r\right)` - - Where :math:`\rho(r)` is the correlation function. - """ + def vario_nugget(self, r): + """Isotropic variogram of the model respecting the nugget at r=0.""" r = np.array(np.abs(r), dtype=np.double) r_gz = np.logical_not(np.isclose(r, 0)) res = np.empty_like(r, dtype=np.double) - res[r_gz] = self.covariance(r[r_gz]) - res[np.logical_not(r_gz)] = self.sill + res[r_gz] = self.variogram(r[r_gz]) + res[np.logical_not(r_gz)] = 0.0 return res - def vario_nugget(self, r): - r"""Isotropic variogram of the model respecting the nugget at r=0. - - Given by: :math:`\gamma\left(r\right)= - \sigma^2\cdot\left(1-\rho\left(r\right)\right)+n` - - Where :math:`\rho(r)` is the correlation function. - """ + def cov_nugget(self, r): + """Isotropic covariance of the model respecting the nugget at r=0.""" r = np.array(np.abs(r), dtype=np.double) r_gz = np.logical_not(np.isclose(r, 0)) res = np.empty_like(r, dtype=np.double) - res[r_gz] = self.variogram(r[r_gz]) - res[np.logical_not(r_gz)] = 0.0 + res[r_gz] = self.covariance(r[r_gz]) + res[np.logical_not(r_gz)] = self.sill return res def plot(self, func="variogram", **kwargs): # pragma: no cover @@ -365,6 +351,9 @@ def plot(self, func="variogram", **kwargs): # pragma: no cover * "vario_spatial" * "cov_spatial" * "cor_spatial" + * "vario_axis" + * "cov_axis" + * "cor_axis" * "spectrum" * "spectral_density" * "spectral_rad_pdf" @@ -383,13 +372,7 @@ def plot(self, func="variogram", **kwargs): # pragma: no cover # pykrige functions def pykrige_vario(self, args=None, r=0): - r"""Isotropic variogram of the model for pykrige. - - Given by: :math:`\gamma\left(r\right)= - \sigma^2\cdot\left(1-\rho\left(r\right)\right)+n` - - Where :math:`\rho(r)` is the correlation function. - """ + """Isotropic variogram of the model for pykrige.""" return self.variogram(r) @property @@ -497,7 +480,6 @@ def check_opt_arg(self): * Any return value will be ignored * This method will only be run once, when the class is initialized """ - pass def fix_dim(self): """Set a fix dimension for the model.""" @@ -1038,7 +1020,9 @@ def dim(self): def dim(self, dim): # check if a fixed dimension should be used if self.fix_dim() is not None: - print(self.name + ": using fixed dimension " + str(self.fix_dim())) + warnings.warn( + self.name + ": using fixed dimension " + str(self.fix_dim()) + ) dim = self.fix_dim() # set the dimension # TODO: add a routine to check if dimension is valid @@ -1216,6 +1200,27 @@ def arg(self): """:class:`list` of :class:`str`: Names of all arguments.""" return ["var", "len_scale", "nugget", "anis", "angles"] + self._opt_arg + @property + def arg_list(self): + """:class:`list` of :class:`float`: Values of all arguments.""" + alist = [self.var, self.len_scale, self.nugget, self.anis, self.angles] + for opt in self.opt_arg: + alist.append(getattr(self, opt)) + return alist + + @property + def iso_arg(self): + """:class:`list` of :class:`str`: Names of isotropic arguments.""" + return ["var", "len_scale", "nugget"] + self._opt_arg + + @property + def iso_arg_list(self): + """:class:`list` of :class:`float`: Values of isotropic arguments.""" + alist = [self.var, self.len_scale, self.nugget] + for opt in self.opt_arg: + alist.append(getattr(self, opt)) + return alist + @property def opt_arg(self): """:class:`list` of :class:`str`: Names of the optional arguments.""" From d9ad2e7bf824e492e8d587fd426e3e77181c7c54 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Sun, 15 Nov 2020 18:26:31 +0100 Subject: [PATCH 43/69] CovModel.fit: better handling of bounds for anis when fitting directional variogram --- gstools/covmodel/fit.py | 25 ++++++++++++++++++------- 1 file changed, 18 insertions(+), 7 deletions(-) diff --git a/gstools/covmodel/fit.py b/gstools/covmodel/fit.py index eb92bda4..a6f73801 100755 --- a/gstools/covmodel/fit.py +++ b/gstools/covmodel/fit.py @@ -14,6 +14,7 @@ import numpy as np from scipy.optimize import curve_fit from gstools.covmodel.tools import check_arg_in_bounds, default_arg_from_bounds +from gstools.tools.geometric import set_anis __all__ = ["fit_variogram"] @@ -340,12 +341,22 @@ def _init_curve_fit_para(model, para, init_guess, constrain_sill, sill, anis): ) ) if anis: - low_bounds += [0.0] * (model.dim - 1) - top_bounds += [np.inf] * (model.dim - 1) - if init_guess == "default": - init_guess_list += [1.0] * (model.dim - 1) - else: + low_bounds += [model.anis_bounds[0]] * (model.dim - 1) + top_bounds += [model.anis_bounds[1]] * (model.dim - 1) + if isinstance(init_guess, dict): + if "anis" not in init_guess: + raise ValueError("CovModel.fit: missing init guess for 'anis'") + init_guess_list += list(set_anis(model.dim, init_guess["anis"])) + elif init_guess == "default": + def_arg = default_arg_from_bounds(model.anis_bounds) + init_guess_list += [def_arg] * (model.dim - 1) + elif init_guess == "current": init_guess_list += list(model.anis) + else: + raise ValueError( + "CovModel.fit: unknown init_guess: '{}'".format(init_guess) + ) + return (low_bounds, top_bounds), init_guess_list @@ -356,7 +367,7 @@ def _init_guess(bounds, current, default, typ, para_name): return typ[para_name] # if we have a dict, all parameters need a given init_guess raise ValueError( - "CovModel.fit: missing init guess for: '{}'".format(para_name) + "CovModel.fit: missing init guess for '{}'".format(para_name) ) if typ == "default": if bounds[0] < default < bounds[1]: @@ -364,7 +375,7 @@ def _init_guess(bounds, current, default, typ, para_name): return default_arg_from_bounds(bounds) if typ == "current": return current - raise ValueError("CovModel.fit: unkwon init_guess: '{}'".format(typ)) + raise ValueError("CovModel.fit: unknown init_guess: '{}'".format(typ)) def _get_curve(model, para, constrain_sill, sill, anis, is_dir_vario): From c4e629c4d0161e4a11a466dbd313ffeb8ae30dae Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Mon, 16 Nov 2020 01:57:59 +0100 Subject: [PATCH 44/69] Examples: add overview gallery to tutorials page in Docs --- docs/source/tutorials.rst | 45 +++++++++++++++++++++++ examples/00_misc/README.rst | 2 + examples/01_random_field/README.rst | 2 + examples/02_cov_model/README.rst | 2 + examples/03_variogram/README.rst | 2 + examples/04_vector_field/README.rst | 2 + examples/05_kriging/README.rst | 2 + examples/06_conditioned_fields/README.rst | 2 + examples/07_transformations/README.rst | 2 + 9 files changed, 61 insertions(+) diff --git a/docs/source/tutorials.rst b/docs/source/tutorials.rst index 5c2b6787..9c0b11d1 100644 --- a/docs/source/tutorials.rst +++ b/docs/source/tutorials.rst @@ -20,3 +20,48 @@ explore its whole beauty and power. examples/06_conditioned_fields/index examples/07_transformations/index examples/00_misc/index + +Gallery of all Tutorial Examples +-------------------------------- + +.. include:: examples/01_random_field/index.rst + :start-after: Below is a gallery of examples + :end-line: -11 + +.. include:: examples/02_cov_model/index.rst + :start-after: Below is a gallery of examples + :end-line: -11 + +.. include:: examples/03_variogram/index.rst + :start-after: Below is a gallery of examples + :end-line: -11 + +.. include:: examples/04_vector_field/index.rst + :start-after: Below is a gallery of examples + :end-line: -11 + +.. include:: examples/05_kriging/index.rst + :start-after: Below is a gallery of examples + :end-line: -11 + +.. include:: examples/06_conditioned_fields/index.rst + :start-after: Below is a gallery of examples + :end-line: -11 + +.. include:: examples/07_transformations/index.rst + :start-after: Below is a gallery of examples + :end-line: -11 + +.. include:: examples/00_misc/index.rst + :start-after: Below is a gallery of examples + :end-line: -11 + +.. raw:: html + +
+ +.. only:: html + + .. rst-class:: sphx-glr-signature + + `Gallery generated by Sphinx-Gallery `_ \ No newline at end of file diff --git a/examples/00_misc/README.rst b/examples/00_misc/README.rst index cf021250..8b407d4f 100644 --- a/examples/00_misc/README.rst +++ b/examples/00_misc/README.rst @@ -5,3 +5,5 @@ A few miscellaneous examples Gallery ------- + +Below is a gallery of examples diff --git a/examples/01_random_field/README.rst b/examples/01_random_field/README.rst index a0452568..f3408bfd 100644 --- a/examples/01_random_field/README.rst +++ b/examples/01_random_field/README.rst @@ -15,3 +15,5 @@ GSTools supports arbitrary and non-isotropic covariance models. Gallery ------- + +Below is a gallery of examples diff --git a/examples/02_cov_model/README.rst b/examples/02_cov_model/README.rst index 8e69698b..775d2e56 100644 --- a/examples/02_cov_model/README.rst +++ b/examples/02_cov_model/README.rst @@ -81,3 +81,5 @@ As a special feature, we also provide truncated power law (TPL) covariance model Gallery ------- + +Below is a gallery of examples diff --git a/examples/03_variogram/README.rst b/examples/03_variogram/README.rst index d16c16c2..34313f64 100644 --- a/examples/03_variogram/README.rst +++ b/examples/03_variogram/README.rst @@ -12,3 +12,5 @@ by this subpackage. Gallery ------- + +Below is a gallery of examples diff --git a/examples/04_vector_field/README.rst b/examples/04_vector_field/README.rst index 1875531e..7c913f44 100644 --- a/examples/04_vector_field/README.rst +++ b/examples/04_vector_field/README.rst @@ -34,3 +34,5 @@ the resulting field is indeed incompressible. Gallery ------- + +Below is a gallery of examples diff --git a/examples/05_kriging/README.rst b/examples/05_kriging/README.rst index d2bbe098..6a9fb8ea 100644 --- a/examples/05_kriging/README.rst +++ b/examples/05_kriging/README.rst @@ -89,3 +89,5 @@ submodule :any:`gstools.krige`. Gallery ------- + +Below is a gallery of examples diff --git a/examples/06_conditioned_fields/README.rst b/examples/06_conditioned_fields/README.rst index bd5a6d9d..94596257 100644 --- a/examples/06_conditioned_fields/README.rst +++ b/examples/06_conditioned_fields/README.rst @@ -32,3 +32,5 @@ or: Gallery ------- + +Below is a gallery of examples diff --git a/examples/07_transformations/README.rst b/examples/07_transformations/README.rst index d3a2ba06..47c5d184 100644 --- a/examples/07_transformations/README.rst +++ b/examples/07_transformations/README.rst @@ -35,3 +35,5 @@ Simply import the transform submodule and apply a transformation to the srf clas Gallery ------- + +Below is a gallery of examples From 49589e784c78ad2a3af1efaba21fbc0287a42c56 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Mon, 16 Nov 2020 14:09:18 +0100 Subject: [PATCH 45/69] Docs: remove Overview gallery again: TOC-tree struggles --- docs/source/conf.py | 2 +- docs/source/tutorials.rst | 63 ++++++++++++++++++++------------------- 2 files changed, 33 insertions(+), 32 deletions(-) diff --git a/docs/source/conf.py b/docs/source/conf.py index b9214669..25b2dc41 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -201,8 +201,8 @@ def setup(app): "pointsize": "10pt", "papersize": "a4paper", "fncychap": "\\usepackage[Glenn]{fncychap}", + # 'inputenc': r'\usepackage[utf8]{inputenc}', } - # Grouping the document tree into LaTeX files. List of tuples # (source start file, target name, title, # author, documentclass [howto, manual, or own class]). diff --git a/docs/source/tutorials.rst b/docs/source/tutorials.rst index 9c0b11d1..363b3c65 100644 --- a/docs/source/tutorials.rst +++ b/docs/source/tutorials.rst @@ -21,47 +21,48 @@ explore its whole beauty and power. examples/07_transformations/index examples/00_misc/index -Gallery of all Tutorial Examples --------------------------------- +.. comment + Gallery of all Tutorial Examples + -------------------------------- -.. include:: examples/01_random_field/index.rst - :start-after: Below is a gallery of examples - :end-line: -11 + .. include:: examples/01_random_field/index.rst + :start-after: Below is a gallery of examples + :end-line: -11 -.. include:: examples/02_cov_model/index.rst - :start-after: Below is a gallery of examples - :end-line: -11 + .. include:: examples/02_cov_model/index.rst + :start-after: Below is a gallery of examples + :end-line: -11 -.. include:: examples/03_variogram/index.rst - :start-after: Below is a gallery of examples - :end-line: -11 + .. include:: examples/03_variogram/index.rst + :start-after: Below is a gallery of examples + :end-line: -11 -.. include:: examples/04_vector_field/index.rst - :start-after: Below is a gallery of examples - :end-line: -11 + .. include:: examples/04_vector_field/index.rst + :start-after: Below is a gallery of examples + :end-line: -11 -.. include:: examples/05_kriging/index.rst - :start-after: Below is a gallery of examples - :end-line: -11 + .. include:: examples/05_kriging/index.rst + :start-after: Below is a gallery of examples + :end-line: -11 -.. include:: examples/06_conditioned_fields/index.rst - :start-after: Below is a gallery of examples - :end-line: -11 + .. include:: examples/06_conditioned_fields/index.rst + :start-after: Below is a gallery of examples + :end-line: -11 -.. include:: examples/07_transformations/index.rst - :start-after: Below is a gallery of examples - :end-line: -11 + .. include:: examples/07_transformations/index.rst + :start-after: Below is a gallery of examples + :end-line: -11 -.. include:: examples/00_misc/index.rst - :start-after: Below is a gallery of examples - :end-line: -11 + .. include:: examples/00_misc/index.rst + :start-after: Below is a gallery of examples + :end-line: -11 -.. raw:: html + .. raw:: html -
+
-.. only:: html + .. only:: html - .. rst-class:: sphx-glr-signature + .. rst-class:: sphx-glr-signature - `Gallery generated by Sphinx-Gallery `_ \ No newline at end of file + `Gallery generated by Sphinx-Gallery `_ From c605cd25f240a025876705d7aeb813fab93fdbbc Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Mon, 16 Nov 2020 14:09:50 +0100 Subject: [PATCH 46/69] Doc: resolve minor problems --- examples/00_misc/README.rst | 8 +++++--- examples/01_random_field/README.rst | 8 +++++--- examples/02_cov_model/README.rst | 8 +++++--- examples/03_variogram/README.rst | 8 +++++--- examples/04_vector_field/README.rst | 8 +++++--- examples/05_kriging/README.rst | 8 +++++--- examples/06_conditioned_fields/README.rst | 8 +++++--- examples/07_transformations/README.rst | 8 +++++--- gstools/random/rng.py | 4 ++-- gstools/variogram/variogram.py | 8 +++++--- 10 files changed, 47 insertions(+), 29 deletions(-) diff --git a/examples/00_misc/README.rst b/examples/00_misc/README.rst index 8b407d4f..15492cbd 100644 --- a/examples/00_misc/README.rst +++ b/examples/00_misc/README.rst @@ -3,7 +3,9 @@ Miscellaneous A few miscellaneous examples -Gallery -------- +.. only:: html -Below is a gallery of examples + Gallery + ------- + + Below is a gallery of examples diff --git a/examples/01_random_field/README.rst b/examples/01_random_field/README.rst index f3408bfd..0682d23e 100644 --- a/examples/01_random_field/README.rst +++ b/examples/01_random_field/README.rst @@ -13,7 +13,9 @@ and its discretised modes are evaluated at random frequencies. GSTools supports arbitrary and non-isotropic covariance models. -Gallery -------- +.. only:: html -Below is a gallery of examples + Gallery + ------- + + Below is a gallery of examples diff --git a/examples/02_cov_model/README.rst b/examples/02_cov_model/README.rst index 775d2e56..7c9c413c 100644 --- a/examples/02_cov_model/README.rst +++ b/examples/02_cov_model/README.rst @@ -79,7 +79,9 @@ As a special feature, we also provide truncated power law (TPL) covariance model TPLStable TPLSimple -Gallery -------- +.. only:: html -Below is a gallery of examples + Gallery + ------- + + Below is a gallery of examples diff --git a/examples/03_variogram/README.rst b/examples/03_variogram/README.rst index 34313f64..5ba19bd7 100644 --- a/examples/03_variogram/README.rst +++ b/examples/03_variogram/README.rst @@ -10,7 +10,9 @@ The same `(semi-)variogram Date: Mon, 16 Nov 2020 14:12:04 +0100 Subject: [PATCH 47/69] Fix some typos --- gstools/covmodel/fit.py | 4 ++-- gstools/covmodel/models.py | 4 ++-- gstools/covmodel/tpl_models.py | 2 +- gstools/tools/geometric.py | 2 +- 4 files changed, 6 insertions(+), 6 deletions(-) diff --git a/gstools/covmodel/fit.py b/gstools/covmodel/fit.py index a6f73801..693a1af8 100755 --- a/gstools/covmodel/fit.py +++ b/gstools/covmodel/fit.py @@ -39,7 +39,7 @@ def fit_variogram( **para_select ): """ - Fiting a variogram-model to an empirical variogram. + Fitting a variogram-model to an empirical variogram. Parameters ---------- @@ -58,7 +58,7 @@ def fit_variogram( by this argument. Deselect the parameter from fitting, by setting it "False". You could also pass a fixed value to be set in the model. - Then the anisotropy ratios wont be altered during fitting. + Then the anisotropy ratios won't be altered during fitting. Default: True sill : :class:`float` or :class:`bool` or :any:`None`, optional Here you can provide a fixed sill for the variogram. diff --git a/gstools/covmodel/models.py b/gstools/covmodel/models.py index c8066544..9f93992e 100644 --- a/gstools/covmodel/models.py +++ b/gstools/covmodel/models.py @@ -527,7 +527,7 @@ class HyperSpherical(CovModel): r"""The Hyper-Spherical covariance model. This model is derived from the relative intersection area of - two d-dimensional hyper spheres, + two d-dimensional hyperspheres, where the middle points have a distance of :math:`r` and the diameters are given by :math:`\ell`. @@ -583,7 +583,7 @@ class SuperSpherical(CovModel): r"""The Super-Spherical covariance model. This model is derived from the relative intersection area of - two d-dimensional hyper spheres, + two d-dimensional hyperspheres, where the middle points have a distance of :math:`r` and the diameters are given by :math:`\ell`. It is than valid in all lower dimensions. diff --git a/gstools/covmodel/tpl_models.py b/gstools/covmodel/tpl_models.py index 8deae0b6..d881e441 100644 --- a/gstools/covmodel/tpl_models.py +++ b/gstools/covmodel/tpl_models.py @@ -488,7 +488,7 @@ class TPLSimple(CovModel): This model describes a simple truncated power law with a finite length scale. In contrast to other models, - this one is not derived from super-postioning modes. + this one is not derived from super-positioning modes. Notes ----- diff --git a/gstools/tools/geometric.py b/gstools/tools/geometric.py index 168df1a6..c11ad3f4 100644 --- a/gstools/tools/geometric.py +++ b/gstools/tools/geometric.py @@ -97,7 +97,7 @@ def set_anis(dim, anis): Notes ----- - If too few anisotrpy ratios are given, they are filled up with `1`. + If too few anisotropy ratios are given, they are filled up with `1`. """ out_anis = np.array(anis, dtype=np.double) out_anis = np.atleast_1d(out_anis)[: dim - 1] From 5d6815b681cd3ecdb7dc0feb9848d3c293b8f4a8 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Tue, 17 Nov 2020 10:58:36 +0100 Subject: [PATCH 48/69] CovModel: add routine 'check_dim' --- gstools/covmodel/base.py | 15 ++++++++++++--- 1 file changed, 12 insertions(+), 3 deletions(-) diff --git a/gstools/covmodel/base.py b/gstools/covmodel/base.py index 60c40d6f..cd106593 100644 --- a/gstools/covmodel/base.py +++ b/gstools/covmodel/base.py @@ -481,6 +481,10 @@ def check_opt_arg(self): * This method will only be run once, when the class is initialized """ + def check_dim(self, dim): + """Check the given dimension.""" + return True + def fix_dim(self): """Set a fix dimension for the model.""" return None @@ -1025,9 +1029,13 @@ def dim(self, dim): ) dim = self.fix_dim() # set the dimension - # TODO: add a routine to check if dimension is valid - if dim < 1 or dim > 3: - raise ValueError("Only dimensions of 1 <= d <= 3 are supported.") + if dim < 1: + raise ValueError("Only dimensions of d >= 1 are supported.") + check = self.check_dim(dim) + if not check: + warnings.warn( + "Dimension {} is not appropriate for this model.".format(dim) + ) self._dim = int(dim) # create fourier transform just once (recreate for dim change) self._sft = SFT(ndim=self.dim, **self.hankel_kw) @@ -1038,6 +1046,7 @@ def dim(self, dim): ) if self._angles is not None: self._angles = set_angles(self.dim, self._angles) + self.check_arg_bounds() @property def var(self): From 0045cb36986af3f4f144c77877e07f2649fa9572 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Tue, 17 Nov 2020 11:19:31 +0100 Subject: [PATCH 49/69] CovModel: skip bound checking, when no bounds present (e.g. during init); add AttributeWarning class to warnings --- gstools/covmodel/base.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/gstools/covmodel/base.py b/gstools/covmodel/base.py index cd106593..72fe6c2d 100644 --- a/gstools/covmodel/base.py +++ b/gstools/covmodel/base.py @@ -848,6 +848,8 @@ def check_arg_bounds(self): """Check arguments to be within the given bounds.""" # check var, len_scale, nugget and optional-arguments for arg in self.arg_bounds: + if not self.arg_bounds[arg]: + continue bnd = list(self.arg_bounds[arg]) val = getattr(self, arg) error_case = check_arg_in_bounds(self, arg) @@ -1025,7 +1027,8 @@ def dim(self, dim): # check if a fixed dimension should be used if self.fix_dim() is not None: warnings.warn( - self.name + ": using fixed dimension " + str(self.fix_dim()) + self.name + ": using fixed dimension " + str(self.fix_dim()), + AttributeWarning, ) dim = self.fix_dim() # set the dimension @@ -1034,7 +1037,8 @@ def dim(self, dim): check = self.check_dim(dim) if not check: warnings.warn( - "Dimension {} is not appropriate for this model.".format(dim) + "Dimension {} is not appropriate for this model.".format(dim), + AttributeWarning, ) self._dim = int(dim) # create fourier transform just once (recreate for dim change) From cd93f0c34f5916ec3e77f9e8cdd3948bf50c6226 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Tue, 17 Nov 2020 11:35:23 +0100 Subject: [PATCH 50/69] CovModels: add dim checks to the linear, circular and spherical model --- gstools/covmodel/base.py | 2 +- gstools/covmodel/models.py | 12 ++++++++++++ 2 files changed, 13 insertions(+), 1 deletion(-) diff --git a/gstools/covmodel/base.py b/gstools/covmodel/base.py index 72fe6c2d..e3cdbbda 100644 --- a/gstools/covmodel/base.py +++ b/gstools/covmodel/base.py @@ -849,7 +849,7 @@ def check_arg_bounds(self): # check var, len_scale, nugget and optional-arguments for arg in self.arg_bounds: if not self.arg_bounds[arg]: - continue + continue # no bounds given during init (called from self.dim) bnd = list(self.arg_bounds[arg]) val = getattr(self, arg) error_case = check_arg_in_bounds(self, arg) diff --git a/gstools/covmodel/models.py b/gstools/covmodel/models.py index 9f93992e..be04ea89 100644 --- a/gstools/covmodel/models.py +++ b/gstools/covmodel/models.py @@ -454,6 +454,10 @@ def cor(self, h): """Linear normalized correlation function.""" return np.maximum(1 - np.abs(h, dtype=np.double), 0.0) + def check_dim(self, dim): + """Linear model is only valid in 1D.""" + return dim < 2 + class Circular(CovModel): r"""The circular covariance model. @@ -493,6 +497,10 @@ def cor(self, h): ) return res + def check_dim(self, dim): + """Circular model is only valid in 1D and 2D.""" + return dim < 3 + class Spherical(CovModel): r"""The Spherical covariance model. @@ -522,6 +530,10 @@ def cor(self, h): h = np.minimum(np.abs(h, dtype=np.double), 1.0) return 1.0 - 1.5 * h + 0.5 * h ** 3 + def check_dim(self, dim): + """Spherical model is only valid in 1D, 2D and 3D.""" + return dim < 4 + class HyperSpherical(CovModel): r"""The Hyper-Spherical covariance model. From 88c6496176f7e9a2c0f430ec4dbaeeb787b3ec83 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Wed, 18 Nov 2020 13:22:14 +0100 Subject: [PATCH 51/69] CovModel: move spectral_rad_pdf to tools --- gstools/covmodel/base.py | 18 ++---------------- gstools/covmodel/tools.py | 34 ++++++++++++++++++++++++++++++++++ 2 files changed, 36 insertions(+), 16 deletions(-) diff --git a/gstools/covmodel/base.py b/gstools/covmodel/base.py index e3cdbbda..8f2a3375 100644 --- a/gstools/covmodel/base.py +++ b/gstools/covmodel/base.py @@ -26,7 +26,7 @@ from gstools.tools.misc import list_format from gstools.covmodel.tools import ( InitSubclassMeta, - rad_fac, + spectral_rad_pdf, set_len_anis, check_bounds, check_arg_in_bounds, @@ -569,21 +569,7 @@ def spectral_density(self, k): def spectral_rad_pdf(self, r): """Radial spectral density of the model.""" - r = np.array(np.abs(r), dtype=np.double) - if self.dim > 1: - r_gz = np.logical_not(np.isclose(r, 0)) - # to prevent numerical errors, we just calculate where r>0 - res = np.zeros_like(r, dtype=np.double) - res[r_gz] = rad_fac(self.dim, r[r_gz]) * np.abs( - self.spectral_density(r[r_gz]) - ) - else: - res = rad_fac(self.dim, r) * np.abs(self.spectral_density(r)) - # prevent numerical errors in hankel for small r values (set 0) - res[np.logical_not(np.isfinite(res))] = 0.0 - # prevent numerical errors in hankel for big r (set non-negative) - res = np.maximum(res, 0.0) - return res + return spectral_rad_pdf(self, r) def ln_spectral_rad_pdf(self, r): """Log radial spectral density of the model.""" diff --git a/gstools/covmodel/tools.py b/gstools/covmodel/tools.py index 5e7752c8..a6b62599 100644 --- a/gstools/covmodel/tools.py +++ b/gstools/covmodel/tools.py @@ -96,6 +96,40 @@ def rad_fac(dim, r): return fac +def spectral_rad_pdf(model, r): + """ + Spectral radians PDF of a model. + + Parameters + ---------- + model : :any:`CovModel` + The covariance model in use. + r : :class:`numpy.ndarray` + Given radii. + + Returns + ------- + :class:`numpy.ndarray` + PDF values. + + """ + r = np.array(np.abs(r), dtype=np.double) + if model.dim > 1: + r_gz = np.logical_not(np.isclose(r, 0)) + # to prevent numerical errors, we just calculate where r>0 + res = np.zeros_like(r, dtype=np.double) + res[r_gz] = rad_fac(model.dim, r[r_gz]) * np.abs( + model.spectral_density(r[r_gz]) + ) + else: + res = rad_fac(model.dim, r) * np.abs(model.spectral_density(r)) + # prevent numerical errors in hankel for small r values (set 0) + res[np.logical_not(np.isfinite(res))] = 0.0 + # prevent numerical errors in hankel for big r (set non-negative) + res = np.maximum(res, 0.0) + return res + + def set_len_anis(dim, len_scale, anis): """Set the length scale and anisotropy factors for the given dimension. From cf14a5f1c5df81e12e5dcb557a73da04dcc8e4ad Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Wed, 18 Nov 2020 13:30:30 +0100 Subject: [PATCH 52/69] CovModel: outsource percentile_scale to tools --- gstools/covmodel/base.py | 17 +----- gstools/covmodel/tools.py | 116 +++++++++++++++++++++++++++----------- 2 files changed, 85 insertions(+), 48 deletions(-) diff --git a/gstools/covmodel/base.py b/gstools/covmodel/base.py index 8f2a3375..dd3d1140 100644 --- a/gstools/covmodel/base.py +++ b/gstools/covmodel/base.py @@ -15,7 +15,6 @@ import copy import numpy as np from scipy.integrate import quad as integral -from scipy.optimize import root from hankel import SymmetricFourierTransform as SFT from gstools.tools.geometric import ( set_angles, @@ -26,11 +25,12 @@ from gstools.tools.misc import list_format from gstools.covmodel.tools import ( InitSubclassMeta, - spectral_rad_pdf, set_len_anis, check_bounds, check_arg_in_bounds, default_arg_from_bounds, + spectral_rad_pdf, + percentile_scale, ) from gstools.covmodel import plot from gstools.covmodel.fit import fit_variogram @@ -510,18 +510,7 @@ def percentile_scale(self, per=0.9): This is the distance, where the given percentile of the variance is reached by the variogram """ - # check the given percentile - if not 0.0 < per < 1.0: - raise ValueError( - "percentile needs to be within (0, 1), got: " + str(per) - ) - - # define a curve, that has its root at the wanted point - def curve(x): - return 1.0 - self.correlation(x) - per - - # take 'per * len_scale' as initial guess - return root(curve, per * self.len_scale)["x"][0] + return percentile_scale(self, per) # spectrum methods (can be overridden for speedup) diff --git a/gstools/covmodel/tools.py b/gstools/covmodel/tools.py index a6b62599..31ee54d2 100644 --- a/gstools/covmodel/tools.py +++ b/gstools/covmodel/tools.py @@ -13,10 +13,13 @@ check_bounds check_arg_in_bounds default_arg_from_bounds + spectral_rad_pdf + percentile_scale """ # pylint: disable=C0103 import numpy as np +from scipy.optimize import root from scipy import special as sps from gstools.tools.geometric import set_anis @@ -27,6 +30,8 @@ "check_bounds", "check_arg_in_bounds", "default_arg_from_bounds", + "spectral_rad_pdf", + "percentile_scale", ] @@ -96,40 +101,6 @@ def rad_fac(dim, r): return fac -def spectral_rad_pdf(model, r): - """ - Spectral radians PDF of a model. - - Parameters - ---------- - model : :any:`CovModel` - The covariance model in use. - r : :class:`numpy.ndarray` - Given radii. - - Returns - ------- - :class:`numpy.ndarray` - PDF values. - - """ - r = np.array(np.abs(r), dtype=np.double) - if model.dim > 1: - r_gz = np.logical_not(np.isclose(r, 0)) - # to prevent numerical errors, we just calculate where r>0 - res = np.zeros_like(r, dtype=np.double) - res[r_gz] = rad_fac(model.dim, r[r_gz]) * np.abs( - model.spectral_density(r[r_gz]) - ) - else: - res = rad_fac(model.dim, r) * np.abs(model.spectral_density(r)) - # prevent numerical errors in hankel for small r values (set 0) - res[np.logical_not(np.isfinite(res))] = 0.0 - # prevent numerical errors in hankel for big r (set non-negative) - res = np.maximum(res, 0.0) - return res - - def set_len_anis(dim, len_scale, anis): """Set the length scale and anisotropy factors for the given dimension. @@ -258,3 +229,80 @@ def default_arg_from_bounds(bounds): if bounds[1] < np.inf: return bounds[1] - 1.0 return 0.0 + + +# outsourced routines + + +def spectral_rad_pdf(model, r): + """ + Spectral radians PDF of a model. + + Parameters + ---------- + model : :any:`CovModel` + The covariance model in use. + r : :class:`numpy.ndarray` + Given radii. + + Returns + ------- + :class:`numpy.ndarray` + PDF values. + + """ + r = np.array(np.abs(r), dtype=np.double) + if model.dim > 1: + r_gz = np.logical_not(np.isclose(r, 0)) + # to prevent numerical errors, we just calculate where r>0 + res = np.zeros_like(r, dtype=np.double) + res[r_gz] = rad_fac(model.dim, r[r_gz]) * np.abs( + model.spectral_density(r[r_gz]) + ) + else: + res = rad_fac(model.dim, r) * np.abs(model.spectral_density(r)) + # prevent numerical errors in hankel for small r values (set 0) + res[np.logical_not(np.isfinite(res))] = 0.0 + # prevent numerical errors in hankel for big r (set non-negative) + res = np.maximum(res, 0.0) + return res + + +def percentile_scale(model, per=0.9): + """ + Calculate the percentile scale of the isotrope model. + + This is the distance, where the given percentile of the variance + is reached by the variogram + + + Parameters + ---------- + model : :any:`CovModel` + The covariance model in use. + per : float, optional + Percentile to use. The default is 0.9. + + Raises + ------ + ValueError + When percentile is not in (0, 1). + + Returns + ------- + float + Percentile scale. + + """ + # check the given percentile + if not 0.0 < per < 1.0: + raise ValueError( + "percentile needs to be within (0, 1), got: " + str(per) + ) + + # define a curve, that has its root at the wanted point + def curve(x): + return 1.0 - model.correlation(x) - per + + # take 'per * len_rescaled' as initial guess + return root(curve, per * model.len_rescaled)["x"][0] From b7462f68d497e6a432a2f2008d25510ff663413b Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Wed, 18 Nov 2020 13:41:26 +0100 Subject: [PATCH 53/69] CovModel: add 'set_opt_arg' routine to tools to shrink __init__ --- gstools/covmodel/base.py | 36 +++------------------------ gstools/covmodel/tools.py | 52 +++++++++++++++++++++++++++++++++++++++ 2 files changed, 55 insertions(+), 33 deletions(-) diff --git a/gstools/covmodel/base.py b/gstools/covmodel/base.py index dd3d1140..2f2018fc 100644 --- a/gstools/covmodel/base.py +++ b/gstools/covmodel/base.py @@ -24,7 +24,9 @@ ) from gstools.tools.misc import list_format from gstools.covmodel.tools import ( + AttributeWarning, InitSubclassMeta, + set_opt_args, set_len_anis, check_bounds, check_arg_in_bounds, @@ -41,10 +43,6 @@ HANKEL_DEFAULT = {"a": -1, "b": 1, "N": 200, "h": 0.001, "alt": True} -class AttributeWarning(UserWarning): - """Attribute warning for CovModel class.""" - - class CovModel(metaclass=InitSubclassMeta): r"""Base class for the GSTools covariance models. @@ -147,35 +145,7 @@ def __init__( self.dim = dim # optional arguments for the variogram-model - self._opt_arg = [] - # look up the defaults for the optional arguments (defined by the user) - default = self.default_opt_arg() - for opt_name in opt_arg: - if opt_name not in default: - warnings.warn( - "The given optional argument '{}' ".format(opt_name) - + "is unknown or has at least no defined standard value. " - + "Or you made a Typo... hehe.", - AttributeWarning, - ) - # add the default vaules if not specified - for def_arg in default: - if def_arg not in opt_arg: - opt_arg[def_arg] = default[def_arg] - # save names of the optional arguments (sort them by name) - self._opt_arg = list(opt_arg.keys()) - self._opt_arg.sort() - # add the optional arguments as attributes to the class - for opt_name in opt_arg: - if opt_name in dir(self): # "dir" also respects properties - raise ValueError( - "parameter '" - + opt_name - + "' has a 'bad' name, since it is already present in " - + "the class. It could not be added to the model" - ) - # Magic happens here - setattr(self, opt_name, float(opt_arg[opt_name])) + set_opt_args(self, opt_arg) # set standard boundaries for variance, len_scale, nugget and opt_arg bounds = self.default_arg_bounds() diff --git a/gstools/covmodel/tools.py b/gstools/covmodel/tools.py index 31ee54d2..a369e604 100644 --- a/gstools/covmodel/tools.py +++ b/gstools/covmodel/tools.py @@ -18,6 +18,7 @@ """ # pylint: disable=C0103 +import warnings import numpy as np from scipy.optimize import root from scipy import special as sps @@ -70,6 +71,10 @@ def __init__(cls, name, bases, ns, **kwargs): super_class.__init_subclass__.__func__(cls, **kwargs) +class AttributeWarning(UserWarning): + """Attribute warning for CovModel class.""" + + # Helping functions ########################################################### @@ -101,6 +106,53 @@ def rad_fac(dim, r): return fac +def set_opt_args(model, opt_arg): + """ + Setting optional arguments in the model class. + + Parameters + ---------- + model : :any:`CovModel` + The covariance model in use. + opt_arg : :class:`dict` + Dictionary with optional arguments. + + Raises + ------ + ValueError + When an optional argument has an already taken name. + """ + model._opt_arg = [] + # look up the defaults for the optional arguments (defined by the user) + default = model.default_opt_arg() + for opt_name in opt_arg: + if opt_name not in default: + warnings.warn( + "The given optional argument '{}' ".format(opt_name) + + "is unknown or has at least no defined standard value. " + + "Or you made a Typo... hehe.", + AttributeWarning, + ) + # add the default vaules if not specified + for def_arg in default: + if def_arg not in opt_arg: + opt_arg[def_arg] = default[def_arg] + # save names of the optional arguments (sort them by name) + model._opt_arg = list(opt_arg.keys()) + model._opt_arg.sort() + # add the optional arguments as attributes to the class + for opt_name in opt_arg: + if opt_name in dir(model): # "dir" also respects properties + raise ValueError( + "parameter '" + + opt_name + + "' has a 'bad' name, since it is already present in " + + "the class. It could not be added to the model" + ) + # Magic happens here + setattr(model, opt_name, float(opt_arg[opt_name])) + + def set_len_anis(dim, len_scale, anis): """Set the length scale and anisotropy factors for the given dimension. From e805f375c6faf45b6c81257c04a8294d62e81a75 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Wed, 18 Nov 2020 13:54:11 +0100 Subject: [PATCH 54/69] CovModel: move '[set/check]_arg_bounds' routines to tools --- gstools/covmodel/base.py | 64 ++--------------------- gstools/covmodel/tools.py | 107 +++++++++++++++++++++++++++++++++++++- 2 files changed, 111 insertions(+), 60 deletions(-) diff --git a/gstools/covmodel/base.py b/gstools/covmodel/base.py index 2f2018fc..90756407 100644 --- a/gstools/covmodel/base.py +++ b/gstools/covmodel/base.py @@ -33,6 +33,8 @@ default_arg_from_bounds, spectral_rad_pdf, percentile_scale, + set_arg_bounds, + check_arg_bounds, ) from gstools.covmodel import plot from gstools.covmodel.fit import fit_variogram @@ -753,67 +755,11 @@ def set_arg_bounds(self, check_args=True, **kwargs): is one of ``"oo"``, ``"cc"``, ``"oc"`` or ``"co"`` to define if the bounds are open ("o") or closed ("c"). """ - # if variance needs to be resetted, do this at last - var_bnds = [] - for arg in kwargs: - if not check_bounds(kwargs[arg]): - raise ValueError( - "Given bounds for '{0}' are not valid, got: {1}".format( - arg, kwargs[arg] - ) - ) - if arg in self.opt_arg: - self._opt_arg_bounds[arg] = kwargs[arg] - elif arg == "var": - var_bnds = kwargs[arg] - continue - elif arg == "len_scale": - self.len_scale_bounds = kwargs[arg] - elif arg == "nugget": - self.nugget_bounds = kwargs[arg] - elif arg == "anis": - self.anis_bounds = kwargs[arg] - else: - raise ValueError( - "set_arg_bounds: unknown argument '{}'".format(arg) - ) - if check_args and check_arg_in_bounds(self, arg) > 0: - def_arg = default_arg_from_bounds(kwargs[arg]) - if arg == "anis": - setattr(self, arg, [def_arg] * (self.dim - 1)) - else: - setattr(self, arg, def_arg) - # set var last like allways - if var_bnds: - self.var_bounds = var_bnds - if check_args and check_arg_in_bounds(self, "var") > 0: - self.var = default_arg_from_bounds(var_bnds) + return set_arg_bounds(self, check_args, **kwargs) def check_arg_bounds(self): - """Check arguments to be within the given bounds.""" - # check var, len_scale, nugget and optional-arguments - for arg in self.arg_bounds: - if not self.arg_bounds[arg]: - continue # no bounds given during init (called from self.dim) - bnd = list(self.arg_bounds[arg]) - val = getattr(self, arg) - error_case = check_arg_in_bounds(self, arg) - if error_case == 1: - raise ValueError( - "{0} needs to be >= {1}, got: {2}".format(arg, bnd[0], val) - ) - if error_case == 2: - raise ValueError( - "{0} needs to be > {1}, got: {2}".format(arg, bnd[0], val) - ) - if error_case == 3: - raise ValueError( - "{0} needs to be <= {1}, got: {2}".format(arg, bnd[1], val) - ) - if error_case == 4: - raise ValueError( - "{0} needs to be < {1}, got: {2}".format(arg, bnd[1], val) - ) + """Check arguments to be within their given bounds.""" + return check_arg_bounds(self) # bounds properties diff --git a/gstools/covmodel/tools.py b/gstools/covmodel/tools.py index a369e604..d9d39a17 100644 --- a/gstools/covmodel/tools.py +++ b/gstools/covmodel/tools.py @@ -8,13 +8,17 @@ .. autosummary:: InitSubclassMeta + AttributeWarning rad_fac + set_opt_args set_len_anis check_bounds check_arg_in_bounds default_arg_from_bounds spectral_rad_pdf percentile_scale + set_arg_bounds + check_arg_bounds """ # pylint: disable=C0103 @@ -26,13 +30,17 @@ __all__ = [ "InitSubclassMeta", + "AttributeWarning", "rad_fac", + "set_opt_args", "set_len_anis", "check_bounds", "check_arg_in_bounds", "default_arg_from_bounds", "spectral_rad_pdf", "percentile_scale", + "set_arg_bounds", + "check_arg_bounds", ] @@ -108,7 +116,7 @@ def rad_fac(dim, r): def set_opt_args(model, opt_arg): """ - Setting optional arguments in the model class. + Set optional arguments in the model class. Parameters ---------- @@ -358,3 +366,100 @@ def curve(x): # take 'per * len_rescaled' as initial guess return root(curve, per * model.len_rescaled)["x"][0] + + +def set_arg_bounds(model, check_args=True, **kwargs): + r"""Set bounds for the parameters of the model. + + Parameters + ---------- + model : :any:`CovModel` + The covariance model in use. + check_args : bool, optional + Whether to check if the arguments are in their valid bounds. + In case not, a propper default value will be determined. + Default: True + **kwargs + Parameter name as keyword ("var", "len_scale", "nugget", ) + and a list of 2 or 3 values as value: + + * ``[a, b]`` or + * ``[a, b, ]`` + + is one of ``"oo"``, ``"cc"``, ``"oc"`` or ``"co"`` + to define if the bounds are open ("o") or closed ("c"). + """ + # if variance needs to be resetted, do this at last + var_bnds = [] + for arg in kwargs: + if not check_bounds(kwargs[arg]): + raise ValueError( + "Given bounds for '{0}' are not valid, got: {1}".format( + arg, kwargs[arg] + ) + ) + if arg in model.opt_arg: + model._opt_arg_bounds[arg] = kwargs[arg] + elif arg == "var": + var_bnds = kwargs[arg] + continue + elif arg == "len_scale": + model.len_scale_bounds = kwargs[arg] + elif arg == "nugget": + model.nugget_bounds = kwargs[arg] + elif arg == "anis": + model.anis_bounds = kwargs[arg] + else: + raise ValueError( + "set_arg_bounds: unknown argument '{}'".format(arg) + ) + if check_args and check_arg_in_bounds(model, arg) > 0: + def_arg = default_arg_from_bounds(kwargs[arg]) + if arg == "anis": + setattr(model, arg, [def_arg] * (model.dim - 1)) + else: + setattr(model, arg, def_arg) + # set var last like allways + if var_bnds: + model.var_bounds = var_bnds + if check_args and check_arg_in_bounds(model, "var") > 0: + model.var = default_arg_from_bounds(var_bnds) + + +def check_arg_bounds(model): + """ + Check arguments to be within their given bounds. + + Parameters + ---------- + model : :any:`CovModel` + The covariance model in use. + + Raises + ------ + ValueError + When an argument is not in its valid bounds. + """ + # check var, len_scale, nugget and optional-arguments + for arg in model.arg_bounds: + if not model.arg_bounds[arg]: + continue # no bounds given during init (called from self.dim) + bnd = list(model.arg_bounds[arg]) + val = getattr(model, arg) + error_case = check_arg_in_bounds(model, arg) + if error_case == 1: + raise ValueError( + "{0} needs to be >= {1}, got: {2}".format(arg, bnd[0], val) + ) + if error_case == 2: + raise ValueError( + "{0} needs to be > {1}, got: {2}".format(arg, bnd[0], val) + ) + if error_case == 3: + raise ValueError( + "{0} needs to be <= {1}, got: {2}".format(arg, bnd[1], val) + ) + if error_case == 4: + raise ValueError( + "{0} needs to be < {1}, got: {2}".format(arg, bnd[1], val) + ) From ce6b8f41dd49e52fb3b7798d523be500604fbab0 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Wed, 18 Nov 2020 14:05:32 +0100 Subject: [PATCH 55/69] CovModel: add 'set_dim' to tools --- gstools/covmodel/base.py | 33 ++----------------------- gstools/covmodel/tools.py | 52 +++++++++++++++++++++++++++++++++++++-- 2 files changed, 52 insertions(+), 33 deletions(-) diff --git a/gstools/covmodel/base.py b/gstools/covmodel/base.py index 90756407..37febef5 100644 --- a/gstools/covmodel/base.py +++ b/gstools/covmodel/base.py @@ -11,7 +11,6 @@ """ # pylint: disable=C0103, R0201 -import warnings import copy import numpy as np from scipy.integrate import quad as integral @@ -24,17 +23,15 @@ ) from gstools.tools.misc import list_format from gstools.covmodel.tools import ( - AttributeWarning, InitSubclassMeta, set_opt_args, set_len_anis, check_bounds, - check_arg_in_bounds, - default_arg_from_bounds, spectral_rad_pdf, percentile_scale, set_arg_bounds, check_arg_bounds, + set_dim, ) from gstools.covmodel import plot from gstools.covmodel.fit import fit_variogram @@ -915,33 +912,7 @@ def dim(self): @dim.setter def dim(self, dim): - # check if a fixed dimension should be used - if self.fix_dim() is not None: - warnings.warn( - self.name + ": using fixed dimension " + str(self.fix_dim()), - AttributeWarning, - ) - dim = self.fix_dim() - # set the dimension - if dim < 1: - raise ValueError("Only dimensions of d >= 1 are supported.") - check = self.check_dim(dim) - if not check: - warnings.warn( - "Dimension {} is not appropriate for this model.".format(dim), - AttributeWarning, - ) - self._dim = int(dim) - # create fourier transform just once (recreate for dim change) - self._sft = SFT(ndim=self.dim, **self.hankel_kw) - # recalculate dimension related parameters - if self._anis is not None: - self._len_scale, self._anis = set_len_anis( - self.dim, self._len_scale, self._anis - ) - if self._angles is not None: - self._angles = set_angles(self.dim, self._angles) - self.check_arg_bounds() + set_dim(self, dim) @property def var(self): diff --git a/gstools/covmodel/tools.py b/gstools/covmodel/tools.py index d9d39a17..96565cab 100644 --- a/gstools/covmodel/tools.py +++ b/gstools/covmodel/tools.py @@ -19,14 +19,16 @@ percentile_scale set_arg_bounds check_arg_bounds + set_dim """ -# pylint: disable=C0103 +# pylint: disable=C0103, W0212 import warnings import numpy as np from scipy.optimize import root from scipy import special as sps -from gstools.tools.geometric import set_anis +from hankel import SymmetricFourierTransform as SFT +from gstools.tools.geometric import set_anis, set_angles __all__ = [ "InitSubclassMeta", @@ -41,6 +43,7 @@ "percentile_scale", "set_arg_bounds", "check_arg_bounds", + "set_dim", ] @@ -463,3 +466,48 @@ def check_arg_bounds(model): raise ValueError( "{0} needs to be < {1}, got: {2}".format(arg, bnd[1], val) ) + + +def set_dim(model, dim): + """ + Set the dimension in the given model. + + Parameters + ---------- + model : :any:`CovModel` + The covariance model in use. + dim : :class:`int` + dimension of the model. + + Raises + ------ + ValueError + When dimension is < 1. + """ + # check if a fixed dimension should be used + if model.fix_dim() is not None: + warnings.warn( + model.name + ": using fixed dimension " + str(model.fix_dim()), + AttributeWarning, + ) + dim = model.fix_dim() + # set the dimension + if dim < 1: + raise ValueError("Only dimensions of d >= 1 are supported.") + check = model.check_dim(dim) + if not check: + warnings.warn( + "Dimension {} is not appropriate for this model.".format(dim), + AttributeWarning, + ) + model._dim = int(dim) + # create fourier transform just once (recreate for dim change) + model._sft = SFT(ndim=model.dim, **model.hankel_kw) + # recalculate dimension related parameters + if model._anis is not None: + model._len_scale, model._anis = set_len_anis( + model.dim, model._len_scale, model._anis + ) + if model._angles is not None: + model._angles = set_angles(model.dim, model._angles) + model.check_arg_bounds() From f40769bda47e859ba82eb6df1072c9a67e8ea93c Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Wed, 18 Nov 2020 14:16:26 +0100 Subject: [PATCH 56/69] CovModel: add '_init_subclass' to tools --- gstools/covmodel/base.py | 56 ++----------------------------------- gstools/covmodel/tools.py | 58 +++++++++++++++++++++++++++++++++++++-- 2 files changed, 59 insertions(+), 55 deletions(-) diff --git a/gstools/covmodel/base.py b/gstools/covmodel/base.py index 37febef5..c4ee62dd 100644 --- a/gstools/covmodel/base.py +++ b/gstools/covmodel/base.py @@ -24,6 +24,7 @@ from gstools.tools.misc import list_format from gstools.covmodel.tools import ( InitSubclassMeta, + _init_subclass, set_opt_args, set_len_anis, check_bounds, @@ -182,60 +183,9 @@ def __init__( # one of these functions needs to be overridden def __init_subclass__(cls): """Initialize gstools covariance model.""" + _init_subclass(cls) - def variogram(self, r): - """Isotropic variogram of the model.""" - return self.var - self.covariance(r) + self.nugget - - def covariance(self, r): - """Covariance of the model.""" - return self.var * self.correlation(r) - - def correlation(self, r): - """Correlation function of the model.""" - return 1.0 - (self.variogram(r) - self.nugget) / self.var - - def correlation_from_cor(self, r): - """Correlation function of the model.""" - r = np.array(np.abs(r), dtype=np.double) - return self.cor(r / self.len_rescaled) - - def cor_from_correlation(self, h): - """Correlation taking a non-dimensional range.""" - h = np.array(np.abs(h), dtype=np.double) - return self.correlation(h * self.len_rescaled) - - abstract = True - if hasattr(cls, "cor"): - if not hasattr(cls, "correlation"): - cls.correlation = correlation_from_cor - abstract = False - else: - cls.cor = cor_from_correlation - if not hasattr(cls, "variogram"): - cls.variogram = variogram - else: - abstract = False - if not hasattr(cls, "covariance"): - cls.covariance = covariance - else: - abstract = False - if not hasattr(cls, "correlation"): - cls.correlation = correlation - else: - abstract = False - if abstract: - raise TypeError( - "Can't instantiate class '" - + cls.__name__ - + "', " - + "without providing at least one of the methods " - + "'cor', 'variogram', 'covariance' or 'correlation'." - ) - - # modify the docstrings - - # class docstring gets attributes added + # modify the docstrings: class docstring gets attributes added if cls.__doc__ is None: cls.__doc__ = ( "User defined GSTools Covariance-Model." diff --git a/gstools/covmodel/tools.py b/gstools/covmodel/tools.py index 96565cab..12ec0f39 100644 --- a/gstools/covmodel/tools.py +++ b/gstools/covmodel/tools.py @@ -47,6 +47,10 @@ ] +class AttributeWarning(UserWarning): + """Attribute warning for CovModel class.""" + + # __init_subclass__ hack for Python 3.5 if hasattr(object, "__init_subclass__"): @@ -82,8 +86,58 @@ def __init__(cls, name, bases, ns, **kwargs): super_class.__init_subclass__.__func__(cls, **kwargs) -class AttributeWarning(UserWarning): - """Attribute warning for CovModel class.""" +def _init_subclass(cls): + """Initialize gstools covariance model.""" + + def variogram(self, r): + """Isotropic variogram of the model.""" + return self.var - self.covariance(r) + self.nugget + + def covariance(self, r): + """Covariance of the model.""" + return self.var * self.correlation(r) + + def correlation(self, r): + """Correlation function of the model.""" + return 1.0 - (self.variogram(r) - self.nugget) / self.var + + def correlation_from_cor(self, r): + """Correlation function of the model.""" + r = np.array(np.abs(r), dtype=np.double) + return self.cor(r / self.len_rescaled) + + def cor_from_correlation(self, h): + """Correlation taking a non-dimensional range.""" + h = np.array(np.abs(h), dtype=np.double) + return self.correlation(h * self.len_rescaled) + + abstract = True + if hasattr(cls, "cor"): + if not hasattr(cls, "correlation"): + cls.correlation = correlation_from_cor + abstract = False + else: + cls.cor = cor_from_correlation + if not hasattr(cls, "variogram"): + cls.variogram = variogram + else: + abstract = False + if not hasattr(cls, "covariance"): + cls.covariance = covariance + else: + abstract = False + if not hasattr(cls, "correlation"): + cls.correlation = correlation + else: + abstract = False + if abstract: + raise TypeError( + "Can't instantiate class '" + + cls.__name__ + + "', " + + "without providing at least one of the methods " + + "'cor', 'variogram', 'covariance' or 'correlation'." + ) # Helping functions ########################################################### From 791147a044a6bf7d7cd51b8d1a06ed6c2db1a05d Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Wed, 18 Nov 2020 17:39:40 +0100 Subject: [PATCH 57/69] CovModel: add more tests for class functionality --- gstools/covmodel/tools.py | 2 +- tests/test_covmodel.py | 113 ++++++++++++++++++++++++++++++++++++++ 2 files changed, 114 insertions(+), 1 deletion(-) diff --git a/gstools/covmodel/tools.py b/gstools/covmodel/tools.py index 12ec0f39..276d23e6 100644 --- a/gstools/covmodel/tools.py +++ b/gstools/covmodel/tools.py @@ -161,7 +161,7 @@ def rad_fac(dim, r): fac = 2 * np.pi * r elif dim == 3: fac = 4 * np.pi * r ** 2 - else: # general solution ( for the record :D ) + else: # pragma: no cover fac = ( dim * r ** (dim - 1) diff --git a/tests/test_covmodel.py b/tests/test_covmodel.py index faa974d7..05607fc9 100644 --- a/tests/test_covmodel.py +++ b/tests/test_covmodel.py @@ -4,6 +4,7 @@ """ import numpy as np import unittest +from gstools.covmodel.tools import AttributeWarning, check_bounds from gstools import ( CovModel, Gaussian, @@ -24,6 +25,40 @@ ) +class Gau_var(CovModel): + def variogram(self, r): + h = np.abs(r) / self.len_rescaled + return self.var * (1.0 - np.exp(-(h ** 2))) + self.nugget + + +class Gau_cov(CovModel): + def covariance(self, r): + h = np.abs(r) / self.len_rescaled + return self.var * np.exp(-(h ** 2)) + + +class Gau_cor(CovModel): + def correlation(self, r): + h = np.abs(r) / self.len_rescaled + return np.exp(-(h ** 2)) + + +class Gau_fix(CovModel): + def cor(self, h): + return np.exp(-(h ** 2)) + + def fix_dim(self): + return 2 + + +class Mod_add(CovModel): + def cor(self, h): + return 1.0 + + def default_opt_arg(self): + return {"alpha": 1} + + class TestCovModel(unittest.TestCase): def setUp(self): self.cov_models = [ @@ -169,6 +204,84 @@ def test_fitting(self): ) self.assertEqual(len(ret), 3) + def test_covmodel_class(self): + model_std = Gaussian(rescale=3, var=1.1, nugget=1.2, len_scale=1.3) + model_var = Gau_var(rescale=3, var=1.1, nugget=1.2, len_scale=1.3) + model_cov = Gau_cov(rescale=3, var=1.1, nugget=1.2, len_scale=1.3) + model_cor = Gau_cor(rescale=3, var=1.1, nugget=1.2, len_scale=1.3) + var = model_std.variogram(2.5) + cov = model_std.covariance(2.5) + corr = model_std.correlation(2.5) + cor = model_std.cor(2.5) + + self.assertFalse(check_bounds(bounds=[0])) + self.assertFalse(check_bounds(bounds=[0, 1, 2, 3])) + self.assertFalse(check_bounds(bounds=[0, 1, "kk"])) + self.assertRaises(ValueError, model_std.set_arg_bounds, wrong_arg=[1]) + self.assertRaises( + ValueError, model_std.set_arg_bounds, wrong_arg=[-1, 1] + ) + + # arg in bounds check + model_std.set_arg_bounds(var=[0.5, 1.5]) + with self.assertRaises(ValueError): + model_std.var = 0.4 + with self.assertRaises(ValueError): + model_std.var = 1.6 + model_std.set_arg_bounds(var=[0.5, 1.5, "oo"]) + with self.assertRaises(ValueError): + model_std.var = 0.5 + with self.assertRaises(ValueError): + model_std.var = 1.5 + model_std.var = 1 + # std value from bounds with neg. inf and finit bound + model_add = Mod_add() + model_add.set_arg_bounds(alpha=[-np.inf, 0]) + self.assertAlmostEqual(model_add.alpha, -1) + # special treatment of anis check + model_std.set_arg_bounds(anis=[2, 4, "oo"]) + self.assertTrue(np.all(np.isclose(model_std.anis, 3))) + # dim specific checks + with self.assertWarns(AttributeWarning): + Gau_fix(dim=1) + self.assertRaises(ValueError, Gaussian, dim=0) + # check inputs + self.assertRaises(ValueError, model_std.percentile_scale, per=-1.0) + self.assertRaises(ValueError, Gaussian, anis=-1.0) + self.assertRaises(ValueError, Gaussian, len_scale=[1, -1]) + self.assertWarns(AttributeWarning, Gaussian, wrong_arg=1.0) + with self.assertWarns(AttributeWarning): + self.assertRaises(ValueError, Gaussian, len_rescaled=1.0) + + # check correct subclassing + with self.assertRaises(TypeError): + + class Gau_err(CovModel): + pass + + self.assertAlmostEqual(var, model_var.variogram(2.5)) + self.assertAlmostEqual(var, model_cov.variogram(2.5)) + self.assertAlmostEqual(var, model_cor.variogram(2.5)) + self.assertAlmostEqual(cov, model_var.covariance(2.5)) + self.assertAlmostEqual(cov, model_cov.covariance(2.5)) + self.assertAlmostEqual(cov, model_cor.covariance(2.5)) + self.assertAlmostEqual(corr, model_var.correlation(2.5)) + self.assertAlmostEqual(corr, model_cov.correlation(2.5)) + self.assertAlmostEqual(corr, model_cor.correlation(2.5)) + self.assertAlmostEqual(cor, model_var.cor(2.5)) + self.assertAlmostEqual(cor, model_cov.cor(2.5)) + self.assertAlmostEqual(cor, model_cor.cor(2.5)) + + def test_rescale(self): + model1 = Exponential() + model2 = Exponential(rescale=2.1) + model3 = Exponential(rescale=2.1, len_scale=2.1) + + self.assertAlmostEqual( + model1.integral_scale, 2.1 * model2.integral_scale + ) + self.assertAlmostEqual(model1.integral_scale, model3.integral_scale) + if __name__ == "__main__": unittest.main() From 38426da8f361622238e1a8e3cfa42a41e7d23e8b Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Wed, 18 Nov 2020 17:47:21 +0100 Subject: [PATCH 58/69] CovModel: minor simplifications for TPL models --- gstools/covmodel/tpl_models.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/gstools/covmodel/tpl_models.py b/gstools/covmodel/tpl_models.py index d881e441..e3c820e8 100644 --- a/gstools/covmodel/tpl_models.py +++ b/gstools/covmodel/tpl_models.py @@ -44,7 +44,7 @@ def len_up_rescaled(self): * ``len_up_rescaled = (len_low + len_scale) / rescale`` """ - return (self.len_low + self.len_scale) / self.rescale + return self.len_up / self.rescale @property def len_low_rescaled(self): @@ -63,11 +63,9 @@ def var_factor(self): def cor(self, h): """TPL - normalized correlation function.""" - return 1.0 def correlation(self, r): """TPL - correlation function.""" - return 1.0 # Truncated power law ######################################################### From 97f9b443271bfa7a85783cae4e2a856fbb83c97c Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Wed, 18 Nov 2020 17:57:24 +0100 Subject: [PATCH 59/69] Tests: covmodel.tools to 100% coverage --- gstools/covmodel/tools.py | 2 +- tests/test_covmodel.py | 8 +++++++- 2 files changed, 8 insertions(+), 2 deletions(-) diff --git a/gstools/covmodel/tools.py b/gstools/covmodel/tools.py index 276d23e6..96a5b0ae 100644 --- a/gstools/covmodel/tools.py +++ b/gstools/covmodel/tools.py @@ -345,7 +345,7 @@ def default_arg_from_bounds(bounds): return bounds[0] + 1.0 if bounds[1] < np.inf: return bounds[1] - 1.0 - return 0.0 + return 0.0 # pragma: no cover # outsourced routines diff --git a/tests/test_covmodel.py b/tests/test_covmodel.py index 05607fc9..a95e1dae 100644 --- a/tests/test_covmodel.py +++ b/tests/test_covmodel.py @@ -4,7 +4,11 @@ """ import numpy as np import unittest -from gstools.covmodel.tools import AttributeWarning, check_bounds +from gstools.covmodel.tools import ( + AttributeWarning, + check_bounds, + check_arg_in_bounds, +) from gstools import ( CovModel, Gaussian, @@ -215,6 +219,7 @@ def test_covmodel_class(self): cor = model_std.cor(2.5) self.assertFalse(check_bounds(bounds=[0])) + self.assertFalse(check_bounds(bounds=[1, -1])) self.assertFalse(check_bounds(bounds=[0, 1, 2, 3])) self.assertFalse(check_bounds(bounds=[0, 1, "kk"])) self.assertRaises(ValueError, model_std.set_arg_bounds, wrong_arg=[1]) @@ -249,6 +254,7 @@ def test_covmodel_class(self): self.assertRaises(ValueError, model_std.percentile_scale, per=-1.0) self.assertRaises(ValueError, Gaussian, anis=-1.0) self.assertRaises(ValueError, Gaussian, len_scale=[1, -1]) + self.assertRaises(ValueError, check_arg_in_bounds, model_std, "wrong") self.assertWarns(AttributeWarning, Gaussian, wrong_arg=1.0) with self.assertWarns(AttributeWarning): self.assertRaises(ValueError, Gaussian, len_rescaled=1.0) From 6e1b0fad2985cc1604f0acc7acb057ec0ea4f32e Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Wed, 18 Nov 2020 18:10:55 +0100 Subject: [PATCH 60/69] CovModel: test *_axis methods and cor_spatial --- tests/test_covmodel.py | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/tests/test_covmodel.py b/tests/test_covmodel.py index a95e1dae..6cd806ff 100644 --- a/tests/test_covmodel.py +++ b/tests/test_covmodel.py @@ -146,6 +146,19 @@ def cor(self, h): + model.nugget - model.cov_spatial(([1], [2], [3])[:dim])[0], ) + self.assertAlmostEqual( + model.cor_spatial(([1], [2], [3])[:dim])[0], + model.cov_spatial(([1], [2], [3])[:dim])[0] + / model.var, + ) + self.assertAlmostEqual( + model.vario_axis(1), + model.var - model.cov_axis(1) + model.nugget, + ) + self.assertAlmostEqual( + model.cor_axis(1), + model.cov_axis(1) / model.var, + ) self.assertAlmostEqual( model.cov_nugget(0), model.sill ) From def555e72517444aec29ece187c79b1290ec0787 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Wed, 18 Nov 2020 18:34:44 +0100 Subject: [PATCH 61/69] CovModel: better testing TPL models --- gstools/covmodel/fit.py | 4 ++++ tests/test_covmodel.py | 33 ++++++++++++++++++--------------- 2 files changed, 22 insertions(+), 15 deletions(-) diff --git a/gstools/covmodel/fit.py b/gstools/covmodel/fit.py index 693a1af8..b09ce0db 100755 --- a/gstools/covmodel/fit.py +++ b/gstools/covmodel/fit.py @@ -380,6 +380,7 @@ def _init_guess(bounds, current, default, typ, para_name): def _get_curve(model, para, constrain_sill, sill, anis, is_dir_vario): """Create the curve for scipys curve_fit.""" + var_save = model.var # we need arg1, otherwise curve_fit throws an error (bug?!) def curve(x, arg1, *args): """Adapted Variogram function.""" @@ -409,6 +410,9 @@ def curve(x, arg1, *args): # set var at last because of var_factor (other parameter needed) if para["var"]: model.var = var_tmp + # needs to be reset for TPL models when len_scale was changed + else: + model.var = var_save if is_dir_vario: if anis: model.anis = args[1 - model.dim :] diff --git a/tests/test_covmodel.py b/tests/test_covmodel.py index 6cd806ff..005b53f1 100644 --- a/tests/test_covmodel.py +++ b/tests/test_covmodel.py @@ -65,7 +65,7 @@ def default_opt_arg(self): class TestCovModel(unittest.TestCase): def setUp(self): - self.cov_models = [ + self.std_cov_models = [ Gaussian, Exponential, Rational, @@ -77,24 +77,14 @@ def setUp(self): HyperSpherical, SuperSpherical, JBessel, + TPLSimple, + ] + self.tpl_cov_models = [ TPLGaussian, TPLExponential, TPLStable, - TPLSimple, - ] - self.std_cov_models = [ - Gaussian, - Exponential, - Rational, - Stable, - Matern, - Linear, - Circular, - Spherical, - HyperSpherical, - SuperSpherical, - JBessel, ] + self.cov_models = self.std_cov_models + self.tpl_cov_models self.dims = range(1, 4) self.lens = [[10, 5, 2]] self.anis = [[0.5, 0.2]] @@ -190,6 +180,19 @@ def cor(self, h): self.assertAlmostEqual(model.var, 3) self.assertAlmostEqual(model.nugget, 1.5) + def test_tpl_models(self): + for Model in self.tpl_cov_models: + for dim in self.dims: + model = Model(dim=dim, len_scale=9, len_low=1, rescale=2) + self.assertAlmostEqual(model.len_up_rescaled, 5) + model.len_low = 0.0 + self.assertAlmostEqual(model.cor(2), model.correlation(9)) + # also check resetting of var when sill is given lower + model.fit_variogram( + self.gamma_x, self.gamma_y, sill=2, nugget=False + ) + self.assertAlmostEqual(model.var, 2.0, delta=1e-5) + def test_fitting(self): for Model in self.std_cov_models: for dim in self.dims: From ae2ff57f848ad6da1aed7aeb0d4b7a7ebd96007b Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Wed, 18 Nov 2020 20:06:48 +0100 Subject: [PATCH 62/69] CovModel: check special models; better check-case for TPL vario fitting --- gstools/covmodel/models.py | 15 ++++++++++----- gstools/covmodel/tpl_models.py | 6 ++++-- tests/test_covmodel.py | 22 ++++++++++++++++++++-- 3 files changed, 34 insertions(+), 9 deletions(-) diff --git a/gstools/covmodel/models.py b/gstools/covmodel/models.py index be04ea89..f7ac5a1f 100644 --- a/gstools/covmodel/models.py +++ b/gstools/covmodel/models.py @@ -19,11 +19,12 @@ SuperSpherical JBessel """ -# pylint: disable=C0103, E1101, E1137 +# pylint: disable=C0103, E1101, E1137, R0201 import warnings import numpy as np from scipy import special as sps +from gstools.covmodel.tools import AttributeWarning from gstools.covmodel import CovModel __all__ = [ @@ -83,6 +84,7 @@ def spectral_rad_cdf(self, r): ) - r * self.len_rescaled / np.sqrt(np.pi) * np.exp( -((r * self.len_rescaled / 2.0) ** 2) ) + return None # pragma: no cover def spectral_rad_ppf(self, u): """Gaussian radial spectral ppf. @@ -96,6 +98,7 @@ def spectral_rad_ppf(self, u): return 2.0 / self.len_rescaled * sps.erfinv(u) if self.dim == 2: return 2.0 / self.len_rescaled * np.sqrt(-np.log(1.0 - u)) + return None # pragma: no cover def _has_cdf(self): return self.dim in [1, 2, 3] @@ -153,7 +156,7 @@ def spectral_rad_cdf(self, r): * 2.0 / np.pi ) - return None + return None # pragma: no cover def spectral_rad_ppf(self, u): """Exponential radial spectral ppf. @@ -173,7 +176,7 @@ def spectral_rad_ppf(self, u): where=np.logical_not(np.isclose(u, 0)), ) return np.sqrt(u_power - 1.0) / self.len_rescaled - return None + return None # pragma: no cover def _has_cdf(self): return self.dim in [1, 2, 3] @@ -242,7 +245,8 @@ def check_opt_arg(self): if self.alpha < 0.3: warnings.warn( "Stable: parameter 'alpha' is < 0.3, " - + "count with unstable results" + + "count with unstable results", + AttributeWarning, ) def cor(self, h): @@ -735,7 +739,8 @@ def check_opt_arg(self): if abs(self.nu - self.dim / 2 + 1) < 0.01: warnings.warn( "JBessel: parameter 'nu' is close to d/2-1, " - + "count with unstable results" + + "count with unstable results", + AttributeWarning, ) def cor(self, h): diff --git a/gstools/covmodel/tpl_models.py b/gstools/covmodel/tpl_models.py index e3c820e8..a097e888 100644 --- a/gstools/covmodel/tpl_models.py +++ b/gstools/covmodel/tpl_models.py @@ -16,7 +16,8 @@ import warnings import numpy as np -from gstools.covmodel.base import CovModel +from gstools.covmodel import CovModel +from gstools.covmodel.tools import AttributeWarning from gstools.tools.special import ( tplstable_cor, tpl_gau_spec_dens, @@ -458,7 +459,8 @@ def check_opt_arg(self): if self.alpha < 0.3: warnings.warn( "TPLStable: parameter 'alpha' is < 0.3, " - + "count with unstable results" + + "count with unstable results", + AttributeWarning, ) def cor(self, h): diff --git a/tests/test_covmodel.py b/tests/test_covmodel.py index 005b53f1..a0c96fa0 100644 --- a/tests/test_covmodel.py +++ b/tests/test_covmodel.py @@ -189,9 +189,10 @@ def test_tpl_models(self): self.assertAlmostEqual(model.cor(2), model.correlation(9)) # also check resetting of var when sill is given lower model.fit_variogram( - self.gamma_x, self.gamma_y, sill=2, nugget=False + self.gamma_x, self.gamma_y, sill=1.1, nugget=False ) - self.assertAlmostEqual(model.var, 2.0, delta=1e-5) + print(model) + self.assertAlmostEqual(model.var, 1.1, delta=1e-5) def test_fitting(self): for Model in self.std_cov_models: @@ -304,6 +305,23 @@ def test_rescale(self): ) self.assertAlmostEqual(model1.integral_scale, model3.integral_scale) + def test_special_models(self): + # matern converges to gaussian + model1 = Matern() + model1.set_arg_bounds(nu=[0, 101]) + model1.nu = 100 + model2 = Gaussian(rescale=0.5) + self.assertAlmostEqual(model1.variogram(1), model2.variogram(1)) + self.assertAlmostEqual(model1.spectrum(1), model2.spectrum(1), 2) + # stable model gets unstable for alpha < 0.3 + with self.assertWarns(AttributeWarning): + Stable(alpha=0.2) + with self.assertWarns(AttributeWarning): + TPLStable(alpha=0.2) + # corner case for JBessel model + with self.assertWarns(AttributeWarning): + JBessel(dim=3, nu=0.5) + if __name__ == "__main__": unittest.main() From 55d41fd2088df83f50ab1a82987b0a960abb88ef Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Wed, 18 Nov 2020 20:38:36 +0100 Subject: [PATCH 63/69] CovModel: further testing of properties; default_opt_args always empty now; check var_raw with TPL models; skip pykrige routine default vals --- gstools/covmodel/base.py | 23 ++++++++++------------- tests/test_covmodel.py | 27 +++++++++++++++++++++++++-- 2 files changed, 35 insertions(+), 15 deletions(-) diff --git a/gstools/covmodel/base.py b/gstools/covmodel/base.py index c4ee62dd..0973389d 100644 --- a/gstools/covmodel/base.py +++ b/gstools/covmodel/base.py @@ -292,56 +292,56 @@ def plot(self, func="variogram", **kwargs): # pragma: no cover def pykrige_vario(self, args=None, r=0): """Isotropic variogram of the model for pykrige.""" - return self.variogram(r) + return self.variogram(r) # pragma: no cover @property def pykrige_anis(self): """2D anisotropy ratio for pykrige.""" if self.dim == 2: return 1 / self.anis[0] - return 1.0 + return 1.0 # pragma: no cover @property def pykrige_anis_y(self): """3D anisotropy ratio in y direction for pykrige.""" if self.dim >= 2: return 1 / self.anis[0] - return 1.0 + return 1.0 # pragma: no cover @property def pykrige_anis_z(self): """3D anisotropy ratio in z direction for pykrige.""" if self.dim == 3: return 1 / self.anis[1] - return 1.0 + return 1.0 # pragma: no cover @property def pykrige_angle(self): """2D rotation angle for pykrige.""" if self.dim == 2: return self.angles[0] / np.pi * 180 - return 0.0 + return 0.0 # pragma: no cover @property def pykrige_angle_z(self): """3D rotation angle around z for pykrige.""" if self.dim >= 2: return self.angles[0] / np.pi * 180 - return 0.0 + return 0.0 # pragma: no cover @property def pykrige_angle_y(self): """3D rotation angle around y for pykrige.""" if self.dim == 3: return self.angles[1] / np.pi * 180 - return 0.0 + return 0.0 # pragma: no cover @property def pykrige_angle_x(self): """3D rotation angle around x for pykrige.""" if self.dim == 3: return self.angles[2] / np.pi * 180 - return 0.0 + return 0.0 # pragma: no cover @property def pykrige_kwargs(self): @@ -374,12 +374,9 @@ def pykrige_kwargs(self): def default_opt_arg(self): """Provide default optional arguments by the user. - Should be given as a dictionary. + Should be given as a dictionary when overridden. """ - res = {} - for opt in self.opt_arg: - res[opt] = 0.0 - return res + return {} def default_opt_arg_bounds(self): """Provide default boundaries for optional arguments.""" diff --git a/tests/test_covmodel.py b/tests/test_covmodel.py index a0c96fa0..159c0154 100644 --- a/tests/test_covmodel.py +++ b/tests/test_covmodel.py @@ -191,8 +191,13 @@ def test_tpl_models(self): model.fit_variogram( self.gamma_x, self.gamma_y, sill=1.1, nugget=False ) - print(model) self.assertAlmostEqual(model.var, 1.1, delta=1e-5) + # check var_raw handling + model = Model(var_raw=1, len_low=0, integral_scale=10) + var_save = model.var + model.var_raw = 1.1 + self.assertAlmostEqual(model.var, var_save * 1.1) + self.assertAlmostEqual(model.integral_scale, 10) def test_fitting(self): for Model in self.std_cov_models: @@ -244,6 +249,15 @@ def test_covmodel_class(self): ValueError, model_std.set_arg_bounds, wrong_arg=[-1, 1] ) + # checking some properties + self.assertEqual(len(model_std.arg), len(model_std.arg_list)) + self.assertEqual(len(model_std.iso_arg), len(model_std.iso_arg_list)) + self.assertEqual(len(model_std.arg), len(model_std.iso_arg) + 2) + self.assertEqual(len(model_std.len_scale_vec), model_std.dim) + self.assertFalse(Gaussian() == Stable()) + model_std.hankel_kw = {"N": 300} + self.assertEqual(model_std.hankel_kw["N"], 300) + # arg in bounds check model_std.set_arg_bounds(var=[0.5, 1.5]) with self.assertRaises(ValueError): @@ -255,7 +269,16 @@ def test_covmodel_class(self): model_std.var = 0.5 with self.assertRaises(ValueError): model_std.var = 1.5 - model_std.var = 1 + with self.assertRaises(ValueError): + model_std.var_bounds = [1, -1] + with self.assertRaises(ValueError): + model_std.len_scale_bounds = [1, -1] + with self.assertRaises(ValueError): + model_std.nugget_bounds = [1, -1] + with self.assertRaises(ValueError): + model_std.anis_bounds = [1, -1] + # reset the standard model + model_std = Gaussian(rescale=3, var=1.1, nugget=1.2, len_scale=1.3) # std value from bounds with neg. inf and finit bound model_add = Mod_add() model_add.set_arg_bounds(alpha=[-np.inf, 0]) From d47d74f510ca5ea8a82847772da25f4a5ce8e164 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Wed, 18 Nov 2020 20:43:27 +0100 Subject: [PATCH 64/69] CovModel: test all axis specific routines --- tests/test_covmodel.py | 19 +++++++++++-------- 1 file changed, 11 insertions(+), 8 deletions(-) diff --git a/tests/test_covmodel.py b/tests/test_covmodel.py index 159c0154..fac2e76a 100644 --- a/tests/test_covmodel.py +++ b/tests/test_covmodel.py @@ -141,14 +141,17 @@ def cor(self, h): model.cov_spatial(([1], [2], [3])[:dim])[0] / model.var, ) - self.assertAlmostEqual( - model.vario_axis(1), - model.var - model.cov_axis(1) + model.nugget, - ) - self.assertAlmostEqual( - model.cor_axis(1), - model.cov_axis(1) / model.var, - ) + for d in range(dim): + self.assertAlmostEqual( + model.vario_axis(1, axis=d), + model.var + + model.nugget + - model.cov_axis(1, axis=d), + ) + self.assertAlmostEqual( + model.cor_axis(1, axis=d), + model.cov_axis(1, axis=d) / model.var, + ) self.assertAlmostEqual( model.cov_nugget(0), model.sill ) From bc2fb16994fa12f4d99d69470766ba11995616ec Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Wed, 18 Nov 2020 21:58:35 +0100 Subject: [PATCH 65/69] CovModel.fit: remove redundant 'dict' option for init_guess --- gstools/covmodel/base.py | 3 +-- gstools/covmodel/fit.py | 17 +++-------------- 2 files changed, 4 insertions(+), 16 deletions(-) diff --git a/gstools/covmodel/base.py b/gstools/covmodel/base.py index 0973389d..56b5b924 100644 --- a/gstools/covmodel/base.py +++ b/gstools/covmodel/base.py @@ -565,12 +565,11 @@ def fit_variogram( and set to the current sill of the model. Then, the procedure above is applied. Default: None - init_guess : :class:`str` or :class:`dict`, optional + init_guess : :class:`str`, optional Initial guess for the estimation. Either: * "default": using the default values of the covariance model * "current": using the current values of the covariance model - * dict(name: val): specified value for each parameter by name Default: "default" weights : :class:`str`, :class:`numpy.ndarray`, :class:`callable`, optional diff --git a/gstools/covmodel/fit.py b/gstools/covmodel/fit.py index b09ce0db..e911a27e 100755 --- a/gstools/covmodel/fit.py +++ b/gstools/covmodel/fit.py @@ -77,12 +77,11 @@ def fit_variogram( and set to the current sill of the model. Then, the procedure above is applied. Default: None - init_guess : :class:`str` or :class:`dict`, optional + init_guess : :class:`str`, optional Initial guess for the estimation. Either: * "default": using the default values of the covariance model * "current": using the current values of the covariance model - * dict(name: val): specified value for each parameter by name Default: "default" weights : :class:`str`, :class:`numpy.ndarray`, :class:`callable`optional @@ -343,11 +342,7 @@ def _init_curve_fit_para(model, para, init_guess, constrain_sill, sill, anis): if anis: low_bounds += [model.anis_bounds[0]] * (model.dim - 1) top_bounds += [model.anis_bounds[1]] * (model.dim - 1) - if isinstance(init_guess, dict): - if "anis" not in init_guess: - raise ValueError("CovModel.fit: missing init guess for 'anis'") - init_guess_list += list(set_anis(model.dim, init_guess["anis"])) - elif init_guess == "default": + if init_guess == "default": def_arg = default_arg_from_bounds(model.anis_bounds) init_guess_list += [def_arg] * (model.dim - 1) elif init_guess == "current": @@ -362,13 +357,6 @@ def _init_curve_fit_para(model, para, init_guess, constrain_sill, sill, anis): def _init_guess(bounds, current, default, typ, para_name): """Proper determination of initial guess.""" - if isinstance(typ, dict): - if para_name in typ: - return typ[para_name] - # if we have a dict, all parameters need a given init_guess - raise ValueError( - "CovModel.fit: missing init guess for '{}'".format(para_name) - ) if typ == "default": if bounds[0] < default < bounds[1]: return default @@ -381,6 +369,7 @@ def _init_guess(bounds, current, default, typ, para_name): def _get_curve(model, para, constrain_sill, sill, anis, is_dir_vario): """Create the curve for scipys curve_fit.""" var_save = model.var + # we need arg1, otherwise curve_fit throws an error (bug?!) def curve(x, arg1, *args): """Adapted Variogram function.""" From c4001bfe444dd9ddfcc431f23be3264a8b1cd63b Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Wed, 18 Nov 2020 22:03:13 +0100 Subject: [PATCH 66/69] CovModel.fit: more testing --- tests/test_covmodel.py | 57 +++++++++++++++++++++++++--- tests/test_variogram_unstructured.py | 18 +++++++++ 2 files changed, 69 insertions(+), 6 deletions(-) diff --git a/tests/test_covmodel.py b/tests/test_covmodel.py index fac2e76a..f95c4801 100644 --- a/tests/test_covmodel.py +++ b/tests/test_covmodel.py @@ -201,6 +201,9 @@ def test_tpl_models(self): model.var_raw = 1.1 self.assertAlmostEqual(model.var, var_save * 1.1) self.assertAlmostEqual(model.integral_scale, 10) + # integral scale is not setable when len_low is not 0 + with self.assertRaises(ValueError): + Model(var_raw=1, len_low=5, integral_scale=10) def test_fitting(self): for Model in self.std_cov_models: @@ -218,6 +221,8 @@ def test_fitting(self): self.gamma_x, self.gamma_y, sill=2, nugget=False ) self.assertAlmostEqual(model.var, 2.0) + model.fit_variogram(self.gamma_x, self.gamma_y, var=1) + self.assertAlmostEqual(model.var, 1) model = Model(dim=dim) model.fit_variogram( self.gamma_x, self.gamma_y, sill=2, nugget=1 @@ -230,9 +235,48 @@ def test_fitting(self): loss="linear", return_r2=True, weights="inv", + init_guess="current", ) self.assertEqual(len(ret), 3) + # treatment of sill/var/nugget by fitting + model = Stable() + model.fit_variogram( + self.gamma_x, self.gamma_y, nugget=False, var=False, sill=2 + ) + self.assertAlmostEqual(model.var, 1) + self.assertAlmostEqual(model.nugget, 1) + model.fit_variogram(self.gamma_x, self.gamma_y, var=2, sill=3) + self.assertAlmostEqual(model.var, 2) + self.assertAlmostEqual(model.nugget, 1) + model.var = 3 + model.fit_variogram( + self.gamma_x, self.gamma_y, nugget=False, var=False, sill=2 + ) + self.assertAlmostEqual(model.var, 2) + self.assertAlmostEqual(model.nugget, 0) + # check ValueErrors + with self.assertRaises(ValueError): + model.fit_variogram(self.gamma_x, self.gamma_y, sill=2, var=3) + with self.assertRaises(ValueError): + model.fit_variogram(self.gamma_x, self.gamma_y, sill=2, nugget=3) + with self.assertRaises(ValueError): + model.fit_variogram(self.gamma_x, self.gamma_y, method="wrong") + with self.assertRaises(ValueError): + model.fit_variogram(self.gamma_x, self.gamma_y, wrong=False) + model.var_bounds = [0, 1] + model.nugget_bounds = [0, 1] + with self.assertRaises(ValueError): + model.fit_variogram(self.gamma_x, self.gamma_y, sill=3) + # init guess + with self.assertRaises(ValueError): + model.fit_variogram(self.gamma_x, self.gamma_y, init_guess="wrong") + model.var_bounds = [0, np.inf] + model.fit_variogram( + self.gamma_x, np.array(self.gamma_y) + 1, sill=2, alpha=False + ) + self.assertAlmostEqual(model.var + model.nugget, 2) + def test_covmodel_class(self): model_std = Gaussian(rescale=3, var=1.1, nugget=1.2, len_scale=1.3) model_var = Gau_var(rescale=3, var=1.1, nugget=1.2, len_scale=1.3) @@ -253,13 +297,14 @@ def test_covmodel_class(self): ) # checking some properties - self.assertEqual(len(model_std.arg), len(model_std.arg_list)) - self.assertEqual(len(model_std.iso_arg), len(model_std.iso_arg_list)) - self.assertEqual(len(model_std.arg), len(model_std.iso_arg) + 2) - self.assertEqual(len(model_std.len_scale_vec), model_std.dim) + model_par = Stable() + self.assertEqual(len(model_par.arg), len(model_par.arg_list)) + self.assertEqual(len(model_par.iso_arg), len(model_par.iso_arg_list)) + self.assertEqual(len(model_par.arg), len(model_par.iso_arg) + 2) + self.assertEqual(len(model_par.len_scale_vec), model_par.dim) self.assertFalse(Gaussian() == Stable()) - model_std.hankel_kw = {"N": 300} - self.assertEqual(model_std.hankel_kw["N"], 300) + model_par.hankel_kw = {"N": 300} + self.assertEqual(model_par.hankel_kw["N"], 300) # arg in bounds check model_std.set_arg_bounds(var=[0.5, 1.5]) diff --git a/tests/test_variogram_unstructured.py b/tests/test_variogram_unstructured.py index 8e4e7ee0..4d2230a2 100644 --- a/tests/test_variogram_unstructured.py +++ b/tests/test_variogram_unstructured.py @@ -334,6 +334,24 @@ def test_mask_no_data(self): self.assertEqual(c4[0], 0) self.assertAlmostEqual(v5[0], 0.0) + def test_fit_directional(self): + model = gs.Stable(dim=3) + bins = [0, 3, 6, 9, 12] + model.len_scale_bounds = [0, 20] + bin_center, emp_vario = gs.vario_estimate( + *(self.pos, self.field, bins), + direction=gs.rotated_main_axes(3, model.angles), + mesh_type="structured", + ) + # check if this succeeds + model.fit_variogram(bin_center, emp_vario, sill=1, return_r2=True) + self.assertTrue(1 > model.anis[0] > model.anis[1]) + model.fit_variogram(bin_center, emp_vario, sill=1, anis=[0.5, 0.25]) + self.assertTrue(15 > model.len_scale) + # catch wrong dim for dir.-vario + with self.assertRaises(ValueError): + model.fit_variogram(bin_center, emp_vario[:2]) + if __name__ == "__main__": unittest.main() From 66943a575298db5cd5292ca501c82b898ec85ceb Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Wed, 18 Nov 2020 22:27:28 +0100 Subject: [PATCH 67/69] CovModel.fit: weights testing --- tests/test_covmodel.py | 8 ++++++-- tests/test_variogram_unstructured.py | 7 ++++++- 2 files changed, 12 insertions(+), 3 deletions(-) diff --git a/tests/test_covmodel.py b/tests/test_covmodel.py index f95c4801..d6730167 100644 --- a/tests/test_covmodel.py +++ b/tests/test_covmodel.py @@ -221,8 +221,6 @@ def test_fitting(self): self.gamma_x, self.gamma_y, sill=2, nugget=False ) self.assertAlmostEqual(model.var, 2.0) - model.fit_variogram(self.gamma_x, self.gamma_y, var=1) - self.assertAlmostEqual(model.var, 1) model = Model(dim=dim) model.fit_variogram( self.gamma_x, self.gamma_y, sill=2, nugget=1 @@ -255,6 +253,12 @@ def test_fitting(self): ) self.assertAlmostEqual(model.var, 2) self.assertAlmostEqual(model.nugget, 0) + model.fit_variogram(self.gamma_x, self.gamma_y, weights="inv") + len_save = model.len_scale + model.fit_variogram( + self.gamma_x, self.gamma_y, weights=lambda x: 1 / (1 + x) + ) + self.assertAlmostEqual(model.len_scale, len_save) # check ValueErrors with self.assertRaises(ValueError): model.fit_variogram(self.gamma_x, self.gamma_y, sill=2, var=3) diff --git a/tests/test_variogram_unstructured.py b/tests/test_variogram_unstructured.py index 4d2230a2..9c3ed977 100644 --- a/tests/test_variogram_unstructured.py +++ b/tests/test_variogram_unstructured.py @@ -338,16 +338,21 @@ def test_fit_directional(self): model = gs.Stable(dim=3) bins = [0, 3, 6, 9, 12] model.len_scale_bounds = [0, 20] - bin_center, emp_vario = gs.vario_estimate( + bin_center, emp_vario, counts = gs.vario_estimate( *(self.pos, self.field, bins), direction=gs.rotated_main_axes(3, model.angles), mesh_type="structured", + return_counts=True, ) # check if this succeeds model.fit_variogram(bin_center, emp_vario, sill=1, return_r2=True) self.assertTrue(1 > model.anis[0] > model.anis[1]) model.fit_variogram(bin_center, emp_vario, sill=1, anis=[0.5, 0.25]) self.assertTrue(15 > model.len_scale) + model.fit_variogram(bin_center, emp_vario, sill=1, weights=counts) + len_save = model.len_scale + model.fit_variogram(bin_center, emp_vario, sill=1, weights=counts[0]) + self.assertAlmostEqual(len_save, model.len_scale) # catch wrong dim for dir.-vario with self.assertRaises(ValueError): model.fit_variogram(bin_center, emp_vario[:2]) From 44090595a103c87fade0acd66391db29427339fd Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Wed, 18 Nov 2020 22:33:17 +0100 Subject: [PATCH 68/69] Tests: used model method for main_axes --- tests/test_variogram_unstructured.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_variogram_unstructured.py b/tests/test_variogram_unstructured.py index 9c3ed977..42eac2ca 100644 --- a/tests/test_variogram_unstructured.py +++ b/tests/test_variogram_unstructured.py @@ -340,7 +340,7 @@ def test_fit_directional(self): model.len_scale_bounds = [0, 20] bin_center, emp_vario, counts = gs.vario_estimate( *(self.pos, self.field, bins), - direction=gs.rotated_main_axes(3, model.angles), + direction=model.main_axes(), mesh_type="structured", return_counts=True, ) From 9e833843a6ffb2a722241f361337faf0fa735298 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Wed, 18 Nov 2020 23:03:09 +0100 Subject: [PATCH 69/69] Docs: remove overview gallery hack again --- docs/source/tutorials.rst | 46 --------------------------------------- 1 file changed, 46 deletions(-) diff --git a/docs/source/tutorials.rst b/docs/source/tutorials.rst index 363b3c65..5c2b6787 100644 --- a/docs/source/tutorials.rst +++ b/docs/source/tutorials.rst @@ -20,49 +20,3 @@ explore its whole beauty and power. examples/06_conditioned_fields/index examples/07_transformations/index examples/00_misc/index - -.. comment - Gallery of all Tutorial Examples - -------------------------------- - - .. include:: examples/01_random_field/index.rst - :start-after: Below is a gallery of examples - :end-line: -11 - - .. include:: examples/02_cov_model/index.rst - :start-after: Below is a gallery of examples - :end-line: -11 - - .. include:: examples/03_variogram/index.rst - :start-after: Below is a gallery of examples - :end-line: -11 - - .. include:: examples/04_vector_field/index.rst - :start-after: Below is a gallery of examples - :end-line: -11 - - .. include:: examples/05_kriging/index.rst - :start-after: Below is a gallery of examples - :end-line: -11 - - .. include:: examples/06_conditioned_fields/index.rst - :start-after: Below is a gallery of examples - :end-line: -11 - - .. include:: examples/07_transformations/index.rst - :start-after: Below is a gallery of examples - :end-line: -11 - - .. include:: examples/00_misc/index.rst - :start-after: Below is a gallery of examples - :end-line: -11 - - .. raw:: html - -
- - .. only:: html - - .. rst-class:: sphx-glr-signature - - `Gallery generated by Sphinx-Gallery `_