Merge pull request #243 from GeoStat-Framework/20220620_add_expint

Add: (Exponential-)Integral model

Author: MuellerSeb

Diff for: examples/02_cov_model/README.rst

@@ -62,6 +62,7 @@ The following standard covariance models are provided by GSTools
     Gaussian
     Exponential
     Matern
+    Integral
     Stable
     Rational
     Cubic

Diff for: src/gstools/__init__.py

@@ -66,6 +66,7 @@
    Gaussian
    Exponential
    Matern
+   Integral
    Stable
    Rational
    Cubic
@@ -145,6 +146,7 @@
     Exponential,
     Gaussian,
     HyperSpherical,
+    Integral,
     JBessel,
     Linear,
     Matern,
@@ -193,6 +195,7 @@
     "Gaussian",
     "Exponential",
     "Matern",
+    "Integral",
     "Stable",
     "Rational",
     "Cubic",

Diff for: src/gstools/covmodel/__init__.py

@@ -31,6 +31,7 @@
    Gaussian
    Exponential
    Matern
+   Integral
    Stable
    Rational
    Cubic
@@ -59,6 +60,7 @@
     Exponential,
     Gaussian,
     HyperSpherical,
+    Integral,
     JBessel,
     Linear,
     Matern,
@@ -79,6 +81,7 @@
     "Gaussian",
     "Exponential",
     "Matern",
+    "Integral",
     "Stable",
     "Rational",
     "Cubic",

Diff for: src/gstools/covmodel/models.py

@@ -10,6 +10,7 @@
    Gaussian
    Exponential
    Matern
+   Integral
    Stable
    Rational
    Cubic
@@ -28,11 +29,13 @@
 
 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",
+    "Integral",
     "Stable",
     "Rational",
     "Cubic",
@@ -363,23 +366,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))
@@ -394,6 +391,97 @@ def calc_integral_scale(self):  # noqa: D102
         )
 
 
+class Integral(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):
+        """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)
+
+    def spectral_density(self, k):  # noqa: D102
+        k = np.asarray(k, dtype=np.double)
+        x = (k * self.len_rescaled / 2.0) ** 2
+        # 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)
+                * self.nu
+                / (self.nu + self.dim)
+                * (1.0 + 2 * x / (self.nu + self.dim + 2))
+            )
+        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)
+        )
+
+
 class Rational(CovModel):
     r"""The rational quadratic covariance model.
 

Diff for: tests/test_covmodel.py

@@ -13,6 +13,7 @@
     Exponential,
     Gaussian,
     HyperSpherical,
+    Integral,
     JBessel,
     Linear,
     Matern,
@@ -82,6 +83,7 @@ def setUp(self):
             SuperSpherical,
             JBessel,
             TPLSimple,
+            Integral,
         ]
         self.tpl_cov_models = [
             TPLGaussian,
@@ -396,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 Integral converge to gaussian
+        model0 = Integral(rescale=0.5)
+        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), 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)
         # stable model gets unstable for alpha < 0.3