diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml index 79993388a..da78b6457 100644 --- a/.github/workflows/main.yml +++ b/.github/workflows/main.yml @@ -55,7 +55,7 @@ jobs: fail-fast: false matrix: os: - [ubuntu-latest, ubuntu-24.04-arm, windows-latest, macos-13, macos-14] + [ubuntu-latest, ubuntu-24.04-arm, windows-latest, macos-15-intel, macos-latest] # https://github.com/scipy/oldest-supported-numpy/blob/main/setup.cfg ver: - { py: "3.8", np: "==1.20.0", sp: "==1.5.4" } @@ -64,15 +64,16 @@ jobs: - { py: "3.11", np: "==1.23.2", sp: "==1.9.2" } - { py: "3.12", np: "==1.26.2", sp: "==1.11.2" } - { py: "3.13", np: "==2.1.0", sp: "==1.14.1" } - - { py: "3.13", np: ">=2.1.0", sp: ">=1.14.1" } + - { py: "3.14", np: "==2.3.2", sp: "==1.16.1" } + - { py: "3.14", np: ">=2.3.2", sp: ">=1.16.1" } exclude: - os: ubuntu-24.04-arm ver: { py: "3.8", np: "==1.20.0", sp: "==1.5.4" } - - os: macos-14 + - os: macos-latest ver: { py: "3.8", np: "==1.20.0", sp: "==1.5.4" } - - os: macos-14 + - os: macos-latest ver: { py: "3.9", np: "==1.20.0", sp: "==1.5.4" } - - os: macos-14 + - os: macos-latest ver: { py: "3.10", np: "==1.21.6", sp: "==1.7.2" } steps: - uses: actions/checkout@v4 diff --git a/examples/02_cov_model/07_roughness.py b/examples/02_cov_model/07_roughness.py new file mode 100644 index 000000000..e392a6b79 --- /dev/null +++ b/examples/02_cov_model/07_roughness.py @@ -0,0 +1,162 @@ +r""" +Roughness +========= + +The roughness describes the power-law behavior of a variogram at the origin +([Wu2016]_): + +.. math:: + \gamma(r) \sim c \cdot r^\alpha \quad (r \to 0) + +The exponent :math:`\alpha` is the roughness information, bounded by +:math:`0 \le \alpha \le 2`. +A value of 0 corresponds to a nugget effect and 2 is the smooth limit for random fields. +You can access it via :any:`CovModel.roughness`. +On a log-log plot, the slope of the variogram near the origin approaches this value. +""" + +import numpy as np +import matplotlib.pyplot as plt + +import gstools as gs + +############################################################################### +# Variogram behavior near the origin +# ---------------------------------- +# +# Compare variograms near the origin for models with different roughness. + +models = [ + gs.Exponential(len_scale=1.0), + gs.Gaussian(len_scale=1.0), + gs.Stable(len_scale=1.0, alpha=0.7), +] + +############################################################################### +# Use a small-lag grid and fit the slope on a log-log scale to estimate the +# roughness numerically. + +r = np.logspace(-4, -1, 200) +fig, ax = plt.subplots() + +for model in models: + gamma = model.variogram(r) + slope = np.polyfit(np.log(r[:20]), np.log(gamma[:20]), 1)[0] + ax.loglog( + r, gamma, label=rf"{model.name} ($\alpha={model.roughness:.2f}$)" + ) + print(f"{model.name}: roughness={model.roughness:.2f}, fit={slope:.2f}") + +ax.set_xlabel("r") +ax.set_ylabel(r"$\gamma(r)$") +ax.set_title("Variogram near the origin (log-log)") +ax.legend() + +############################################################################### +# A nugget masks the power-law behavior at the origin, so roughness is 0. + +nugget_model = gs.Gaussian(nugget=1.0) +print("Gaussian with nugget roughness:", nugget_model.roughness) + +############################################################################### +# Roughness and random fields +# --------------------------- +# +# From the theory of fractal stochastic processes (Mandelbrot and Van Ness +# 1968) ([Mandelbrot1968]_), the roughness can be interpreted as: +# +# 1. Persistent (:math:`1 < \alpha \le 2`): smooth behavior, slowly increasing +# variogram, long memory (e.g. Gaussian-like). +# 2. Antipersistent (:math:`0 < \alpha < 1`): rough behavior, steep variogram +# near the origin, short memory. +# 3. No memory (:math:`\alpha = 1`): linear slope at the origin (Exponential). +# +# The Integral model lets us control the roughness with its parameter ``nu``. +# For this model, the roughness is given by :math:`\alpha = \min(2, \nu)`. + +############################################################################### +# Set up a common grid and plotting scales so the realizations are comparable. + +sep = 100 # separation point (mid field) +ext = 10 # field extend +imext = 2 * [0, ext] +x = y = np.linspace(0, ext, 2 * sep + 1) +xm = np.linspace(0, 5, 1001) +vmin, vmax = -3, 3 + +############################################################################### +# Create Integral models with the same integral scale but different roughness. + +rough = [0.1, 1, 10] +mod = [gs.Integral(dim=2, integral_scale=1, nu=nu) for nu in rough] + +############################################################################### +# .. note:: +# +# Rough fields require a higher ``mode_no`` so the randomization method +# samples the high-frequency part of the spectrum sufficiently well. + +############################################################################### +# Generate random field realizations with a shared seed for fair comparison. + +srf = [] +for m in mod: + mode_no = int(1000 / m.roughness) # increase mode_no for rough fields + s = gs.SRF(m, seed=20110101, mode_no=mode_no) + s.structured((x, y)) + srf.append(s) + +############################################################################### +# Plot variograms, fields, and cross sections column-wise by roughness. + +fig, axes = plt.subplots(3, 3, figsize=(10, 10)) +for i in range(3): + label = rf"$\alpha(\gamma)={mod[i].roughness:.2f}$" + axes[0, i].plot(xm, mod[i].variogram(xm), label=label) + axes[0, i].legend(loc="lower right") + axes[0, i].set_ylim([-0.05, 1.05]) + axes[0, i].set_xlim([-0.25, 5.25]) + axes[0, i].grid() + axes[1, i].imshow( + srf[i].field.T, + origin="lower", + interpolation="bicubic", + vmin=vmin, + vmax=vmax, + extent=imext, + ) + axes[1, i].axhline(y=5, color="k") + axes[2, i].plot(x, srf[i].field[:, sep]) + axes[2, i].set_ylim([vmin, vmax]) + axes[2, i].set_xlim([0, ext]) + axes[2, i].grid() + +axes[0, 0].set_ylabel(r"$\gamma$") +axes[1, 0].set_ylabel(r"$y$") +axes[2, 0].set_ylabel(r"$f$") +axes[2, 0].set_xlabel(r"$x$") +axes[2, 1].set_xlabel(r"$x$") +axes[2, 2].set_xlabel(r"$x$") +fig.tight_layout() +# sphinx_gallery_thumbnail_number = 2 + +############################################################################### +# Illustration of the impact of the roughness information on spatial random +# fields. Each column shows the variogram on top, a single field realization in +# the middle and the highlighted cross section of the field on the bottom. The +# left column shows a very rough (antipersistent) field with a steep increase +# of the variogram at the origin, the middle column shows an Exponential like +# variogram with a linear slope at the origin resulting in a less rough (no +# memory) field and the right column shows a Gaussian like variogram together +# with a very smooth (persistent) field. All variograms have the same integral +# scale and x- and y-axis are given in multiples of the integral scale. + +############################################################################### +# References +# ---------- +# .. [Wu2016] Wu, W.-Y., and C. Y. Lim. 2016. "ESTIMATION OF SMOOTHNESS OF A +# STATIONARY GAUSSIAN RANDOM FIELD." Statistica Sinica 26 (4): 1729-1745. +# https://doi.org/10.5705/ss.202014.0109. +# .. [Mandelbrot1968] Mandelbrot, B. B., and J. W. Van Ness. 1968. "Fractional +# Brownian Motions, Fractional Noises and Applications." SIAM Review 10 (4): +# 422-437. https://doi.org/10.1137/1010093. diff --git a/pyproject.toml b/pyproject.toml index 885bcd774..009e1217a 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -12,7 +12,7 @@ authors = [ { name = "Sebastian Müller, Lennart Schüler", email = "info@geostat-framework.org" }, ] readme = "README.md" -license = "LGPL-3.0" +license = "LGPL-3.0-or-later" dynamic = ["version"] classifiers = [ "Development Status :: 5 - Production/Stable", @@ -34,6 +34,7 @@ classifiers = [ "Programming Language :: Python :: 3.11", "Programming Language :: Python :: 3.12", "Programming Language :: Python :: 3.13", + "Programming Language :: Python :: 3.14", "Topic :: Scientific/Engineering", "Topic :: Scientific/Engineering :: GIS", "Topic :: Scientific/Engineering :: Hydrology", diff --git a/src/gstools/covmodel/base.py b/src/gstools/covmodel/base.py index 8e3b2ae3f..4e28bcd79 100644 --- a/src/gstools/covmodel/base.py +++ b/src/gstools/covmodel/base.py @@ -52,6 +52,7 @@ pos2latlon, rotated_main_axes, ) +from gstools.tools.misc import derivative __all__ = ["CovModel"] @@ -1167,6 +1168,53 @@ def is_isotropic(self): """:any:`bool`: State if a model is isotropic.""" return np.all(np.isclose(self.anis, 1.0)) + @property + def roughness(self): + """:any:`float`: roughness information of the model. Zero for any present nugget.""" + if self.nugget > 0: + return 0.0 + if np.isclose(self.var, 0): + return np.inf + if hasattr(self, "_roughness"): + return self._roughness() + return self.calc_roughness() + + def calc_roughness(self, x=1e-4, dx=1e-4): + """Calculate the roughness of the model. + + This ignores the nugget of the model. + + Parameters + ---------- + x : :class:`float`, optional + Point at which the derivative is calculated. + Default: ``1e-4`` + dx : :class:`float`, optional + Step size for the derivative calculation. + Default: ``1e-4`` + + Returns + ------- + roughness : :class:`float` + Roughness of the model. + + Notes + ----- + The roughness is defined as the derivative of the log-log + correlation function at zero: + + * ``roughness = d( log(1 - cor(r)) ) / d( log(r) ) | r->0`` + + This is calculated numerically by evaluating the derivative + at a small distance `x`. + """ + + def f(h): + """Function for derivative calculation.""" + return np.log(1 - self.cor(np.exp(h))) + + return derivative(f, np.log(x), dx=dx, order=1) + def __eq__(self, other): """Compare CovModels.""" if not isinstance(other, CovModel): diff --git a/src/gstools/covmodel/models.py b/src/gstools/covmodel/models.py index 3c47b351c..d0bdcbf0f 100644 --- a/src/gstools/covmodel/models.py +++ b/src/gstools/covmodel/models.py @@ -188,6 +188,9 @@ def _has_ppf(self): def calc_integral_scale(self): # noqa: D102 return self.len_rescaled * np.sqrt(np.pi) / 2.0 + def _roughness(self): + return 2.0 + class Exponential(CovModel): r"""The Exponential covariance model. @@ -272,6 +275,9 @@ def _has_ppf(self): def calc_integral_scale(self): # noqa: D102 return self.len_rescaled + def _roughness(self): + return 1.0 + class Stable(CovModel): r"""The stable covariance model. @@ -347,6 +353,9 @@ def cor(self, h): def calc_integral_scale(self): # noqa: D102 return self.len_rescaled * sps.gamma(1.0 + 1.0 / self.alpha) + def _roughness(self): + return self.alpha + class Matern(CovModel): r"""The Matérn covariance model. @@ -457,6 +466,9 @@ def calc_integral_scale(self): # noqa: D102 / sps.beta(self.nu, 0.5) ) + def _roughness(self): + return min(2.0, 2 * self.nu) + class Integral(CovModel): r"""The Exponential Integral covariance model. @@ -548,6 +560,9 @@ def calc_integral_scale(self): # noqa: D102 self.len_rescaled * self.nu * np.sqrt(np.pi) / (2 * self.nu + 2.0) ) + def _roughness(self): + return min(2.0, self.nu) + class Rational(CovModel): r"""The rational quadratic covariance model. @@ -620,6 +635,9 @@ def calc_integral_scale(self): # noqa: D102 / 2.0 ) + def _roughness(self): + return 2.0 + class Cubic(CovModel): r"""The Cubic covariance model. @@ -656,6 +674,9 @@ def cor(self, h): h = np.minimum(np.abs(h, dtype=np.double), 1.0) return 1.0 - 7 * h**2 + 8.75 * h**3 - 3.5 * h**5 + 0.75 * h**7 + def _roughness(self): + return 2.0 + class Linear(CovModel): r"""The bounded linear covariance model. @@ -692,6 +713,9 @@ def check_dim(self, dim): """Linear model is only valid in 1D.""" return dim < 2 + def _roughness(self): + return 1.0 + class Circular(CovModel): r"""The circular covariance model. @@ -741,6 +765,9 @@ def check_dim(self, dim): """Circular model is only valid in 1D and 2D.""" return dim < 3 + def _roughness(self): + return 1.0 + class Spherical(CovModel): r"""The Spherical covariance model. @@ -780,6 +807,9 @@ def check_dim(self, dim): """Spherical model is only valid in 1D, 2D and 3D.""" return dim < 4 + def _roughness(self): + return 1.0 + class HyperSpherical(CovModel): r"""The Hyper-Spherical covariance model. @@ -841,6 +871,9 @@ def spectral_density(self, k): # noqa: D102 ) return res + def _roughness(self): + return 1.0 + class SuperSpherical(CovModel): r"""The Super-Spherical covariance model. @@ -917,6 +950,9 @@ def cor(self, h): ) return res + def _roughness(self): + return 1.0 + class JBessel(CovModel): r"""The J-Bessel hole model. @@ -1004,12 +1040,7 @@ def check_opt_arg(self): def cor(self, h): """J-Bessel correlation.""" h = np.asarray(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 + return sps.hyp0f1(self.nu + 1, -0.25 * h**2) def spectral_density(self, k): # noqa: D102 k = np.asarray(k, dtype=np.double) @@ -1025,3 +1056,14 @@ def spectral_density(self, k): # noqa: D102 * (1.0 - (kk * self.len_rescaled) ** 2) ** (self.nu - self.dim / 2) ) return res + + def calc_integral_scale(self): # noqa: D102 + return ( + self.len_rescaled + * np.sqrt(np.pi) + * sps.gamma(self.nu + 1) + / sps.gamma(self.nu + 0.5) + ) + + def _roughness(self): + return 2.0 diff --git a/src/gstools/covmodel/tpl_models.py b/src/gstools/covmodel/tpl_models.py index 46472d1b7..bc894c776 100644 --- a/src/gstools/covmodel/tpl_models.py +++ b/src/gstools/covmodel/tpl_models.py @@ -211,6 +211,11 @@ def spectral_density(self, k): # noqa: D102 k, self.dim, self.len_rescaled, self.hurst, self.len_low_rescaled ) + def _roughness(self): + if self.len_low > 0: + return 2.0 # roughness of gaussian model + return 2 * self.hurst + class TPLExponential(TPLCovModel): r"""Truncated-Power-Law with Exponential modes. @@ -344,6 +349,11 @@ def spectral_density(self, k): # noqa: D102 k, self.dim, self.len_rescaled, self.hurst, self.len_low_rescaled ) + def _roughness(self): + if self.len_low > 0: + return 1.0 # roughness of exponential model + return 2 * self.hurst + class TPLStable(TPLCovModel): r"""Truncated-Power-Law with Stable modes. @@ -499,6 +509,11 @@ def correlation(self, r): - self.len_low_rescaled ** (2 * self.hurst) ) + def _roughness(self): + if self.len_low > 0: + return self.alpha # roughness of stable model + return 2 * self.hurst + class TPLSimple(CovModel): r"""The simply truncated power law model. @@ -573,3 +588,6 @@ def default_opt_arg_bounds(self): def cor(self, h): """TPL Simple - normalized correlation function.""" return np.maximum(1 - np.abs(h, dtype=np.double), 0.0) ** self.nu + + def _roughness(self): + return 1.0 diff --git a/src/gstools/normalizer/base.py b/src/gstools/normalizer/base.py index 2aedd4d95..0fa6f0b46 100644 --- a/src/gstools/normalizer/base.py +++ b/src/gstools/normalizer/base.py @@ -14,10 +14,7 @@ import numpy as np import scipy.optimize as spo - -def _derivative(f, x, dx=1e-6): - """Central difference formula.""" - return (f(x + dx) - f(x - dx)) / (2 * dx) +from gstools.tools.misc import derivative class Normalizer: @@ -60,7 +57,7 @@ def _normalize(self, data): return data def _derivative(self, data): - return _derivative(self._normalize, data, dx=self._dx) + return derivative(self._normalize, data, dx=self._dx) def _loglikelihood(self, data): add = -0.5 * np.size(data) * (np.log(2 * np.pi) + 1) diff --git a/src/gstools/tools/misc.py b/src/gstools/tools/misc.py index 3fe44678e..8bc74eb75 100755 --- a/src/gstools/tools/misc.py +++ b/src/gstools/tools/misc.py @@ -15,7 +15,33 @@ from gstools.tools.geometric import format_struct_pos_dim, generate_grid -__all__ = ["get_fig_ax", "list_format", "eval_func"] +__all__ = ["derivative", "get_fig_ax", "list_format", "eval_func"] + + +def derivative(f, x, dx=1e-6, order=2): + """Central difference formula. + + Parameters + ---------- + f : :any:`callable` + Function to differentiate. + x : array_like + Point(s) where to evaluate the derivative. + dx : :class:`float`, optional + Step size for the central difference. The default is 1e-6. + order : :class:`int`, optional + Order of the derivative approximation. Either 1 (forward difference) or + 2 (central difference). The default is 2. + + Returns + ------- + :class:`numpy.ndarray` + Derivative of f at x. + """ + x = np.asarray(x, dtype=np.double) + if order == 1: + return (f(x + dx) - f(x)) / dx + return (f(x + dx) - f(x - dx)) / (2 * dx) def get_fig_ax(fig=None, ax=None, ax_name="rectilinear"): # pragma: no cover diff --git a/tests/test_covmodel.py b/tests/test_covmodel.py index 62ea8a1c2..bf0794729 100644 --- a/tests/test_covmodel.py +++ b/tests/test_covmodel.py @@ -377,6 +377,24 @@ class Gau_err(CovModel): self.assertAlmostEqual(cor, model_cov.cor(2.5)) self.assertAlmostEqual(cor, model_cor.cor(2.5)) + def test_roughness(self): + self.assertAlmostEqual(Gaussian().roughness, 2.0) + self.assertAlmostEqual(Exponential().roughness, 1.0) + self.assertAlmostEqual(Stable(alpha=1.2).roughness, 1.2) + self.assertAlmostEqual(Matern(nu=0.4).roughness, 0.8) + self.assertAlmostEqual(Integral(nu=0.6).roughness, 0.6) + self.assertAlmostEqual(TPLGaussian(hurst=0.4).roughness, 0.8) + self.assertAlmostEqual( + TPLGaussian(hurst=0.4, len_low=1.0).roughness, 2.0 + ) + self.assertAlmostEqual( + TPLStable(hurst=0.4, alpha=1.2, len_low=1.0).roughness, + 1.2, + ) + self.assertAlmostEqual(Gaussian(nugget=0.5).roughness, 0.0) + self.assertTrue(np.isinf(Gaussian(var=0.0).roughness)) + self.assertAlmostEqual(Gau_cor().roughness, 2.0, places=2) + def test_rescale(self): model1 = Exponential() model2 = Exponential(rescale=2.1)