From 42c4aa58816eeb4c305adbdf3c97b889c9279578 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Mon, 20 Jun 2022 17:16:11 +0200 Subject: [PATCH 1/8] Models: add the ExpInt model --- examples/02_cov_model/README.rst | 2 + src/gstools/__init__.py | 6 ++ src/gstools/covmodel/__init__.py | 6 ++ src/gstools/covmodel/models.py | 97 ++++++++++++++++++++++++++++++++ 4 files changed, 111 insertions(+) diff --git a/examples/02_cov_model/README.rst b/examples/02_cov_model/README.rst index 6f8d3dad2..986e62406 100644 --- a/examples/02_cov_model/README.rst +++ b/examples/02_cov_model/README.rst @@ -62,6 +62,8 @@ The following standard covariance models are provided by GSTools Gaussian Exponential Matern + ExpInt + Mueller Stable Rational Cubic diff --git a/src/gstools/__init__.py b/src/gstools/__init__.py index 2d0e72327..770ac27c9 100644 --- a/src/gstools/__init__.py +++ b/src/gstools/__init__.py @@ -64,6 +64,8 @@ Gaussian Exponential Matern + ExpInt + Mueller Stable Rational Cubic @@ -141,12 +143,14 @@ Circular, CovModel, Cubic, + ExpInt, Exponential, Gaussian, HyperSpherical, JBessel, Linear, Matern, + Mueller, Rational, Spherical, Stable, @@ -192,6 +196,8 @@ "Gaussian", "Exponential", "Matern", + "ExpInt", + "Mueller", "Stable", "Rational", "Cubic", diff --git a/src/gstools/covmodel/__init__.py b/src/gstools/covmodel/__init__.py index a793ae6a0..5f95b408d 100644 --- a/src/gstools/covmodel/__init__.py +++ b/src/gstools/covmodel/__init__.py @@ -29,6 +29,8 @@ Gaussian Exponential Matern + ExpInt + Mueller Stable Rational Cubic @@ -54,12 +56,14 @@ from gstools.covmodel.models import ( Circular, Cubic, + ExpInt, Exponential, Gaussian, HyperSpherical, JBessel, Linear, Matern, + Mueller, Rational, Spherical, Stable, @@ -77,6 +81,8 @@ "Gaussian", "Exponential", "Matern", + "ExpInt", + "Mueller", "Stable", "Rational", "Cubic", diff --git a/src/gstools/covmodel/models.py b/src/gstools/covmodel/models.py index 102467e2b..149a266c4 100644 --- a/src/gstools/covmodel/models.py +++ b/src/gstools/covmodel/models.py @@ -10,6 +10,8 @@ Gaussian Exponential Matern + ExpInt + Mueller Stable Rational Cubic @@ -28,11 +30,14 @@ from gstools.covmodel.base import CovModel from gstools.covmodel.tools import AttributeWarning +from gstools.tools.special import exp_int, inc_gamma_low __all__ = [ "Gaussian", "Exponential", "Matern", + "ExpInt", + "Mueller", "Stable", "Rational", "Cubic", @@ -394,6 +399,98 @@ def calc_integral_scale(self): # noqa: D102 ) +class Mueller(CovModel): + r"""The Exponential Integral covariance model. + + Notes + ----- + This model is given by the following correlation function [Mueller2021]_: + + .. math:: + \rho(r) = + \frac{\nu}{2}\cdot + E_{1+\frac{\nu}{2}}\left( \left( s\cdot\frac{r}{\ell} \right)^2 \right) + + Where the standard rescale factor is :math:`s=1`. + :math:`E_s(x)` is the exponential integral. + + :math:`\nu` is a shape parameter (1 by default). + + For :math:`\nu \to \infty`, a gaussian model is approached, since it represents + the limiting case: + + .. math:: + \rho(r) = + \exp\left(-\left(s\cdot\frac{r}{\ell}\right)^2\right) + + References + ---------- + .. [Mueller2021] Müller, S., Heße, F., Attinger, S., and Zech, A., + "The extended generalized radial flow model and effective + conductivity for truncated power law variograms", + Adv. Water Resour., 156, 104027, (2021) + + Other Parameters + ---------------- + nu : :class:`float`, optional + Shape parameter. Standard range: ``(0.0, 50]`` + Default: ``1.0`` + """ + + def default_opt_arg(self): + """Defaults for the optional arguments. + + * ``{"nu": 1.0}`` + + Returns + ------- + :class:`dict` + Defaults for optional arguments + """ + return {"nu": 1.0} + + def default_opt_arg_bounds(self): + """Defaults for boundaries of the optional arguments. + + * ``{"nu": [0.0, 50.0, "oc"]}`` + + Returns + ------- + :class:`dict` + Boundaries for optional arguments + """ + return {"nu": [0.0, 50.0, "oc"]} + + def cor(self, h): + """ExpInt normalized correlation function.""" + h = np.asarray(h, dtype=np.double) + return 0.5 * self.nu * exp_int(1.0 + 0.5 * self.nu, h**2) + + def spectral_density(self, k): # noqa: D102 + k = np.asarray(k, dtype=np.double) + x = (k * self.len_rescaled / 2.0) ** 2 + # for nu > 20 we just use an approximation of the gaussian model + if self.nu > 20.0: + return ( + (0.5 * self.len_rescaled / np.sqrt(np.pi)) ** self.dim + * np.exp(-x) + * (1.0 + x / (self.nu * 0.5 + 1.0)) + ) + return ( + self.nu + / (x ** (self.nu * 0.5) * 2 * (k * np.sqrt(np.pi)) ** self.dim) + * inc_gamma_low((self.nu + self.dim) / 2.0, x) + ) + + def calc_integral_scale(self): # noqa: D102 + return ( + self.len_rescaled * self.nu * np.sqrt(np.pi) / (2 * self.nu + 2.0) + ) + + +ExpInt = Mueller + + class Rational(CovModel): r"""The rational quadratic covariance model. From be2041e9b73816fb6e9d71cac581288b9897de86 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Mon, 20 Jun 2022 17:17:09 +0200 Subject: [PATCH 2/8] add ExpInt to tests --- tests/test_covmodel.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/tests/test_covmodel.py b/tests/test_covmodel.py index 415bc691a..6d1df0de1 100644 --- a/tests/test_covmodel.py +++ b/tests/test_covmodel.py @@ -16,6 +16,7 @@ JBessel, Linear, Matern, + Mueller, Rational, Spherical, Stable, @@ -82,6 +83,7 @@ def setUp(self): SuperSpherical, JBessel, TPLSimple, + Mueller, ] self.tpl_cov_models = [ TPLGaussian, From 27bfe1bcfd3723a29467dad5a29275f77be3f6cb Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Mon, 20 Jun 2022 17:27:29 +0200 Subject: [PATCH 3/8] Tests: add test for covergence of expint model --- tests/test_covmodel.py | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/tests/test_covmodel.py b/tests/test_covmodel.py index 6d1df0de1..0f0a55ce4 100644 --- a/tests/test_covmodel.py +++ b/tests/test_covmodel.py @@ -398,11 +398,16 @@ def test_rescale(self): self.assertAlmostEqual(model1.integral_scale, model3.integral_scale) def test_special_models(self): - # matern converges to gaussian + # matern and mueller converge to gaussian + model0 = Mueller(rescale=0.5) + model0.set_arg_bounds(nu=[0, 101]) + model0.nu = 100 model1 = Matern() model1.set_arg_bounds(nu=[0, 101]) model1.nu = 100 model2 = Gaussian(rescale=0.5) + self.assertAlmostEqual(model0.variogram(1), model2.variogram(1)) + self.assertAlmostEqual(model0.spectrum(1), model2.spectrum(1), 2) 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 From b6aca764d64e51cd96d9ca8ed8fd1413236bf8db Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Mon, 20 Jun 2022 23:04:06 +0200 Subject: [PATCH 4/8] Test: update precision --- 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 0f0a55ce4..82d663448 100644 --- a/tests/test_covmodel.py +++ b/tests/test_covmodel.py @@ -400,13 +400,13 @@ def test_rescale(self): def test_special_models(self): # matern and mueller converge to gaussian model0 = Mueller(rescale=0.5) - model0.set_arg_bounds(nu=[0, 101]) - model0.nu = 100 + model0.set_arg_bounds(nu=[0, 1001]) + model0.nu = 1000 model1 = Matern() model1.set_arg_bounds(nu=[0, 101]) model1.nu = 100 model2 = Gaussian(rescale=0.5) - self.assertAlmostEqual(model0.variogram(1), model2.variogram(1)) + self.assertAlmostEqual(model0.variogram(1), model2.variogram(1), 2) self.assertAlmostEqual(model0.spectrum(1), model2.spectrum(1), 2) self.assertAlmostEqual(model1.variogram(1), model2.variogram(1)) self.assertAlmostEqual(model1.spectrum(1), model2.spectrum(1), 2) From eebd160120f33a10c567670c3c75b3a9e3518797 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Mon, 20 Jun 2022 23:55:05 +0200 Subject: [PATCH 5/8] ExpInt: increase threshold for nu in spectral_density --- src/gstools/covmodel/models.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/gstools/covmodel/models.py b/src/gstools/covmodel/models.py index 149a266c4..946a80c36 100644 --- a/src/gstools/covmodel/models.py +++ b/src/gstools/covmodel/models.py @@ -469,8 +469,8 @@ def cor(self, h): def spectral_density(self, k): # noqa: D102 k = np.asarray(k, dtype=np.double) x = (k * self.len_rescaled / 2.0) ** 2 - # for nu > 20 we just use an approximation of the gaussian model - if self.nu > 20.0: + # for nu > 50 we just use an approximation of the gaussian model + if self.nu > 50.0: return ( (0.5 * self.len_rescaled / np.sqrt(np.pi)) ** self.dim * np.exp(-x) From f704deeb5fdbfc3150be58a10998aacd60ac1f8c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Tue, 21 Jun 2022 11:04:32 +0200 Subject: [PATCH 6/8] ExpInt+Matern: corrections for spectrum approximation nu -> inf --- src/gstools/covmodel/models.py | 20 ++++++++------------ 1 file changed, 8 insertions(+), 12 deletions(-) diff --git a/src/gstools/covmodel/models.py b/src/gstools/covmodel/models.py index 946a80c36..f44ba3201 100644 --- a/src/gstools/covmodel/models.py +++ b/src/gstools/covmodel/models.py @@ -368,23 +368,17 @@ def cor(self, h): def spectral_density(self, k): # noqa: D102 k = np.asarray(k, dtype=np.double) + x = (k * self.len_rescaled) ** 2 # for nu > 20 we just use an approximation of the gaussian model if self.nu > 20.0: return ( (self.len_rescaled / np.sqrt(np.pi)) ** self.dim - * np.exp(-((k * self.len_rescaled) ** 2)) - * ( - 1 - + ( - ((k * self.len_rescaled) ** 2 - self.dim / 2.0) ** 2 - - self.dim / 2.0 - ) - / self.nu - ) + * np.exp(-x) + * (1 + 0.5 * x**2 / self.nu) + * np.sqrt(1 + x / self.nu) ** (-self.dim) ) 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_rescaled) ** 2 / self.nu) + -(self.nu + self.dim / 2.0) * np.log(1.0 + x / self.nu) + sps.loggamma(self.nu + self.dim / 2.0) - sps.loggamma(self.nu) - self.dim * np.log(np.sqrt(self.nu)) @@ -474,7 +468,9 @@ def spectral_density(self, k): # noqa: D102 return ( (0.5 * self.len_rescaled / np.sqrt(np.pi)) ** self.dim * np.exp(-x) - * (1.0 + x / (self.nu * 0.5 + 1.0)) + * self.nu + / (self.nu + self.dim) + * (1.0 + 2 * x / (self.nu + self.dim + 2)) ) return ( self.nu From dfdea63baba0e247bd9bb7f28ea5e6cc4bf249d6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Wed, 19 Oct 2022 13:28:06 +0200 Subject: [PATCH 7/8] ExpInt: rename to 'Integral' --- examples/02_cov_model/README.rst | 3 +-- src/gstools/__init__.py | 9 +++------ src/gstools/covmodel/models.py | 13 ++++--------- tests/test_covmodel.py | 8 ++++---- 4 files changed, 12 insertions(+), 21 deletions(-) diff --git a/examples/02_cov_model/README.rst b/examples/02_cov_model/README.rst index 986e62406..73704183c 100644 --- a/examples/02_cov_model/README.rst +++ b/examples/02_cov_model/README.rst @@ -62,8 +62,7 @@ The following standard covariance models are provided by GSTools Gaussian Exponential Matern - ExpInt - Mueller + Integral Stable Rational Cubic diff --git a/src/gstools/__init__.py b/src/gstools/__init__.py index 770ac27c9..96fc75f27 100644 --- a/src/gstools/__init__.py +++ b/src/gstools/__init__.py @@ -64,8 +64,7 @@ Gaussian Exponential Matern - ExpInt - Mueller + Integral Stable Rational Cubic @@ -143,14 +142,13 @@ Circular, CovModel, Cubic, - ExpInt, Exponential, Gaussian, HyperSpherical, + Integral, JBessel, Linear, Matern, - Mueller, Rational, Spherical, Stable, @@ -196,8 +194,7 @@ "Gaussian", "Exponential", "Matern", - "ExpInt", - "Mueller", + "Integral", "Stable", "Rational", "Cubic", diff --git a/src/gstools/covmodel/models.py b/src/gstools/covmodel/models.py index f44ba3201..61243486f 100644 --- a/src/gstools/covmodel/models.py +++ b/src/gstools/covmodel/models.py @@ -10,8 +10,7 @@ Gaussian Exponential Matern - ExpInt - Mueller + Integral Stable Rational Cubic @@ -36,8 +35,7 @@ "Gaussian", "Exponential", "Matern", - "ExpInt", - "Mueller", + "Integral", "Stable", "Rational", "Cubic", @@ -393,7 +391,7 @@ def calc_integral_scale(self): # noqa: D102 ) -class Mueller(CovModel): +class Integral(CovModel): r"""The Exponential Integral covariance model. Notes @@ -456,7 +454,7 @@ def default_opt_arg_bounds(self): return {"nu": [0.0, 50.0, "oc"]} def cor(self, h): - """ExpInt normalized correlation function.""" + """Exponential Integral normalized correlation function.""" h = np.asarray(h, dtype=np.double) return 0.5 * self.nu * exp_int(1.0 + 0.5 * self.nu, h**2) @@ -484,9 +482,6 @@ def calc_integral_scale(self): # noqa: D102 ) -ExpInt = Mueller - - class Rational(CovModel): r"""The rational quadratic covariance model. diff --git a/tests/test_covmodel.py b/tests/test_covmodel.py index 82d663448..3b9b5eaaa 100644 --- a/tests/test_covmodel.py +++ b/tests/test_covmodel.py @@ -13,10 +13,10 @@ Exponential, Gaussian, HyperSpherical, + Integral, JBessel, Linear, Matern, - Mueller, Rational, Spherical, Stable, @@ -83,7 +83,7 @@ def setUp(self): SuperSpherical, JBessel, TPLSimple, - Mueller, + Integral, ] self.tpl_cov_models = [ TPLGaussian, @@ -398,8 +398,8 @@ def test_rescale(self): self.assertAlmostEqual(model1.integral_scale, model3.integral_scale) def test_special_models(self): - # matern and mueller converge to gaussian - model0 = Mueller(rescale=0.5) + # Matern and Integral converge to gaussian + model0 = Integral(rescale=0.5) model0.set_arg_bounds(nu=[0, 1001]) model0.nu = 1000 model1 = Matern() From 3b1ac655d7a2fae8f2ecd89e14236043b6d5d8d4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Wed, 19 Oct 2022 13:31:08 +0200 Subject: [PATCH 8/8] CovModel: fix wrong imports --- src/gstools/covmodel/__init__.py | 9 +++------ 1 file changed, 3 insertions(+), 6 deletions(-) diff --git a/src/gstools/covmodel/__init__.py b/src/gstools/covmodel/__init__.py index 5f95b408d..87affce9d 100644 --- a/src/gstools/covmodel/__init__.py +++ b/src/gstools/covmodel/__init__.py @@ -29,8 +29,7 @@ Gaussian Exponential Matern - ExpInt - Mueller + Integral Stable Rational Cubic @@ -56,14 +55,13 @@ from gstools.covmodel.models import ( Circular, Cubic, - ExpInt, Exponential, Gaussian, HyperSpherical, + Integral, JBessel, Linear, Matern, - Mueller, Rational, Spherical, Stable, @@ -81,8 +79,7 @@ "Gaussian", "Exponential", "Matern", - "ExpInt", - "Mueller", + "Integral", "Stable", "Rational", "Cubic",