"""
GStools subpackage providing different covariance models.
.. currentmodule:: gstools.covmodel.models
The following classes are provided
.. autosummary::
Gaussian
Exponential
Matern
Integral
Stable
Rational
Cubic
Linear
Circular
Spherical
HyperSpherical
SuperSpherical
JBessel
"""
# pylint: disable=C0103, E1101, R0201
import warnings
import numpy as np
from scipy import special as sps
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",
"Linear",
"Circular",
"Spherical",
"HyperSpherical",
"SuperSpherical",
"JBessel",
]
[docs]class Gaussian(CovModel):
r"""The Gaussian covariance model.
Notes
-----
This model is given by the following variogram [Webster2007]_:
.. math::
\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}`.
References
----------
.. [Webster2007] Webster, R. and Oliver, M. A.
"Geostatistics for environmental scientists.",
John Wiley & Sons. (2007)
"""
[docs] def cor(self, h):
"""Gaussian normalized correlation function."""
return np.exp(-(h**2))
[docs] def default_rescale(self):
"""Gaussian rescaling factor to result in integral scale."""
return np.sqrt(np.pi) / 2.0
[docs] def spectral_density(self, k): # noqa: D102
k = np.asarray(k, dtype=np.double)
return (self.len_rescaled / 2.0 / np.sqrt(np.pi)) ** self.dim * np.exp(
-((k * self.len_rescaled / 2.0) ** 2)
)
[docs] def spectral_rad_cdf(self, r):
"""Gaussian radial spectral cdf."""
r = np.asarray(r, dtype=np.double)
if self.dim == 1:
return sps.erf(r * self.len_rescaled / 2.0)
if self.dim == 2:
return 1.0 - np.exp(-((r * self.len_rescaled / 2.0) ** 2))
if self.dim == 3:
return sps.erf(
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 # pragma: no cover
[docs] def spectral_rad_ppf(self, u):
"""Gaussian radial spectral ppf.
Notes
-----
Not defined for 3D.
"""
u = np.asarray(u, dtype=np.double)
if self.dim == 1:
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]
def _has_ppf(self):
return self.dim in [1, 2]
[docs] def calc_integral_scale(self): # noqa: D102
return self.len_rescaled * np.sqrt(np.pi) / 2.0
[docs]class Exponential(CovModel):
r"""The Exponential covariance model.
Notes
-----
This model is given by the following variogram [Webster2007]_:
.. math::
\gamma(r)=
\sigma^{2}
\left(1-\exp\left(-s\cdot\frac{r}{\ell}\right)\right)+n
Where the standard rescale factor is :math:`s=1`.
References
----------
.. [Webster2007] Webster, R. and Oliver, M. A.
"Geostatistics for environmental scientists.",
John Wiley & Sons. (2007)
"""
[docs] def cor(self, h):
"""Exponential normalized correlation function."""
return np.exp(-h)
[docs] def spectral_density(self, k): # noqa: D102
k = np.asarray(k, dtype=np.double)
return (
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)
)
[docs] def spectral_rad_cdf(self, r):
"""Exponential radial spectral cdf."""
r = np.asarray(r, dtype=np.double)
if self.dim == 1:
return np.arctan(r * self.len_rescaled) * 2.0 / np.pi
if self.dim == 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_rescaled)
- r
* self.len_rescaled
/ (1.0 + (r * self.len_rescaled) ** 2)
)
* 2.0
/ np.pi
)
return None # pragma: no cover
[docs] def spectral_rad_ppf(self, u):
"""Exponential radial spectral ppf.
Notes
-----
Not defined for 3D.
"""
u = np.asarray(u, dtype=np.double)
if self.dim == 1:
return np.tan(np.pi / 2 * u) / self.len_rescaled
if self.dim == 2:
u_power = np.divide(
1,
u**2,
out=np.full_like(u, np.inf),
where=np.logical_not(np.isclose(u, 0)),
)
return np.sqrt(u_power - 1.0) / self.len_rescaled
return None # pragma: no cover
def _has_cdf(self):
return self.dim in [1, 2, 3]
def _has_ppf(self):
return self.dim in [1, 2]
[docs] def calc_integral_scale(self): # noqa: D102
return self.len_rescaled
[docs]class Stable(CovModel):
r"""The stable covariance model.
Notes
-----
This model is given by the following correlation function
[Wackernagel2003]_:
.. math::
\rho(r) =
\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]`
References
----------
.. [Wackernagel2003] Wackernagel, H. "Multivariate geostatistics",
Springer, Berlin, Heidelberg (2003)
Other Parameters
----------------
alpha : :class:`float`, optional
Shape parameter. Standard range: ``(0, 2]``
Default: ``1.5``
"""
[docs] def default_opt_arg(self):
"""Defaults for the optional arguments.
* ``{"alpha": 1.5}``
Returns
-------
:class:`dict`
Defaults for optional arguments
"""
return {"alpha": 1.5}
[docs] def default_opt_arg_bounds(self):
"""Defaults for boundaries of the optional arguments.
* ``{"alpha": [0, 2, "oc"]}``
Returns
-------
:class:`dict`
Boundaries for optional arguments
"""
return {"alpha": [0, 2, "oc"]}
[docs] def check_opt_arg(self):
"""Check the optional arguments.
Warns
-----
alpha
If alpha is < 0.3, the model tends to a nugget model and gets
numerically unstable.
"""
if self.alpha < 0.3:
warnings.warn(
"Stable: parameter 'alpha' is < 0.3, "
"count with unstable results",
AttributeWarning,
)
[docs] def cor(self, h):
r"""Stable normalized correlation function."""
return np.exp(-np.power(h, self.alpha))
[docs] def calc_integral_scale(self): # noqa: D102
return self.len_rescaled * sps.gamma(1.0 + 1.0 / self.alpha)
[docs]class Matern(CovModel):
r"""The Matérn covariance model.
Notes
-----
This model is given by the following correlation function [Rasmussen2003]_:
.. math::
\rho(r) =
\frac{2^{1-\nu}}{\Gamma\left(\nu\right)} \cdot
\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 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.
If :math:`\nu > 20`, a gaussian model is used, since it represents
the limiting case:
.. math::
\rho(r) =
\exp\left(-\left(s\cdot\frac{r}{2\ell}\right)^2\right)
References
----------
.. [Rasmussen2003] Rasmussen, C. E.,
"Gaussian processes in machine learning." Summer school on
machine learning. Springer, Berlin, Heidelberg, (2003)
Other Parameters
----------------
nu : :class:`float`, optional
Shape parameter. Standard range: ``[0.2, 30]``
Default: ``1.0``
"""
[docs] def default_opt_arg(self):
"""Defaults for the optional arguments.
* ``{"nu": 1.0}``
Returns
-------
:class:`dict`
Defaults for optional arguments
"""
return {"nu": 1.0}
[docs] def default_opt_arg_bounds(self):
"""Defaults for boundaries of the optional arguments.
* ``{"nu": [0.2, 30.0, "cc"]}``
Returns
-------
:class:`dict`
Boundaries for optional arguments
"""
return {"nu": [0.2, 30.0, "cc"]}
[docs] def cor(self, h):
"""Matérn normalized correlation function."""
h = np.asarray(np.abs(h), dtype=np.double)
# for nu > 20 we just use the gaussian model
if self.nu > 20.0:
return np.exp(-((h / 2.0) ** 2))
# calculate by log-transformation to prevent numerical errors
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) * 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 positive
res = np.maximum(res, 0.0)
return res
[docs] 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(-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 + x / self.nu)
+ sps.loggamma(self.nu + self.dim / 2.0)
- sps.loggamma(self.nu)
- self.dim * np.log(np.sqrt(self.nu))
)
[docs] def calc_integral_scale(self): # noqa: D102
return (
self.len_rescaled
* np.pi
/ np.sqrt(self.nu)
/ sps.beta(self.nu, 0.5)
)
[docs]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``
"""
[docs] def default_opt_arg(self):
"""Defaults for the optional arguments.
* ``{"nu": 1.0}``
Returns
-------
:class:`dict`
Defaults for optional arguments
"""
return {"nu": 1.0}
[docs] 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"]}
[docs] 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)
[docs] 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)
)
[docs] def calc_integral_scale(self): # noqa: D102
return (
self.len_rescaled * self.nu * np.sqrt(np.pi) / (2 * self.nu + 2.0)
)
[docs]class Rational(CovModel):
r"""The rational quadratic covariance model.
Notes
-----
This model is given by the following correlation function [Rasmussen2003]_:
.. 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)
References
----------
.. [Rasmussen2003] Rasmussen, C. E.,
"Gaussian processes in machine learning." Summer school on
machine learning. Springer, Berlin, Heidelberg, (2003)
Other Parameters
----------------
alpha : :class:`float`, optional
Shape parameter. Standard range: ``[0.5, 50]``
Default: ``1.0``
"""
[docs] def default_opt_arg(self):
"""Defaults for the optional arguments.
* ``{"alpha": 1.0}``
Returns
-------
:class:`dict`
Defaults for optional arguments
"""
return {"alpha": 1.0}
[docs] def default_opt_arg_bounds(self):
"""Defaults for boundaries of the optional arguments.
* ``{"alpha": [0.5, 50.0]}``
Returns
-------
:class:`dict`
Boundaries for optional arguments
"""
return {"alpha": [0.5, 50.0]}
[docs] def cor(self, h):
"""Rational normalized correlation function."""
return np.power(1 + h**2 / self.alpha, -self.alpha)
[docs] 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
)
[docs]class Cubic(CovModel):
r"""The Cubic covariance model.
A model with reverse curvature near the origin and a finite range of
correlation.
Notes
-----
This model is given by the following correlation function [Chiles2009]_:
.. math::
\rho(r) =
\begin{cases}
1- 7 \left(s\cdot\frac{r}{\ell}\right)^{2}
+ \frac{35}{4} \left(s\cdot\frac{r}{\ell}\right)^{3}
- \frac{7}{2} \left(s\cdot\frac{r}{\ell}\right)^{5}
+ \frac{3}{4} \left(s\cdot\frac{r}{\ell}\right)^{7}
& r<\frac{\ell}{s}\\
0 & r\geq\frac{\ell}{s}
\end{cases}
Where the standard rescale factor is :math:`s=1`.
References
----------
.. [Chiles2009] Chiles, J. P., & Delfiner, P.,
"Geostatistics: modeling spatial uncertainty" (Vol. 497),
John Wiley & Sons. (2009)
"""
[docs] def cor(self, h):
"""Spherical normalized correlation function."""
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
[docs]class Linear(CovModel):
r"""The bounded linear covariance model.
This model is derived from the relative intersection area of
two lines in 1D, where the middle points have a distance of :math:`r`
and the line lengths are :math:`\ell`.
Notes
-----
This model is given by the following correlation function [Webster2007]_:
.. math::
\rho(r) =
\begin{cases}
1-s\cdot\frac{r}{\ell} & r<\frac{\ell}{s}\\
0 & r\geq\frac{\ell}{s}
\end{cases}
Where the standard rescale factor is :math:`s=1`.
References
----------
.. [Webster2007] Webster, R. and Oliver, M. A.
"Geostatistics for environmental scientists.",
John Wiley & Sons. (2007)
"""
[docs] def cor(self, h):
"""Linear normalized correlation function."""
return np.maximum(1 - np.abs(h, dtype=np.double), 0.0)
[docs] def check_dim(self, dim):
"""Linear model is only valid in 1D."""
return dim < 2
[docs]class Circular(CovModel):
r"""The circular covariance model.
This model is derived as the relative intersection area of
two discs in 2D, where the middle points have a distance of :math:`r`
and the diameters are given by :math:`\ell`.
Notes
-----
This model is given by the following correlation function [Webster2007]_:
.. math::
\rho(r) =
\begin{cases}
\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<\frac{\ell}{s}\\
0 & r\geq\frac{\ell}{s}
\end{cases}
Where the standard rescale factor is :math:`s=1`.
References
----------
.. [Webster2007] Webster, R. and Oliver, M. A.
"Geostatistics for environmental scientists.",
John Wiley & Sons. (2007)
"""
[docs] def cor(self, h):
"""Circular normalized correlation function."""
h = np.asarray(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
[docs] def check_dim(self, dim):
"""Circular model is only valid in 1D and 2D."""
return dim < 3
[docs]class Spherical(CovModel):
r"""The Spherical covariance model.
This model is derived from the relative intersection area of
two spheres in 3D, where the middle points have a distance of :math:`r`
and the diameters are given by :math:`\ell`.
Notes
-----
This model is given by the following correlation function [Webster2007]_:
.. math::
\rho(r) =
\begin{cases}
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}
Where the standard rescale factor is :math:`s=1`.
References
----------
.. [Webster2007] Webster, R. and Oliver, M. A.
"Geostatistics for environmental scientists.",
John Wiley & Sons. (2007)
"""
[docs] 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
[docs] def check_dim(self, dim):
"""Spherical model is only valid in 1D, 2D and 3D."""
return dim < 4
[docs]class HyperSpherical(CovModel):
r"""The Hyper-Spherical covariance model.
This model is derived from the relative intersection area of
two d-dimensional hyperspheres,
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 the Circular model
and in 3D the Spherical model.
Notes
-----
This model is given by the following correlation function [Matern1960]_:
.. math::
\rho(r) =
\begin{cases}
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{d-1}{2},\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:`d` is the dimension.
References
----------
.. [Matern1960] Matern B., "Spatial Variation",
Swedish National Institute for Forestry Research, (1960)
"""
[docs] def cor(self, h):
"""Hyper-Spherical normalized correlation function."""
h = np.asarray(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
[docs] def spectral_density(self, k): # noqa: D102
k = np.asarray(k, dtype=np.double)
res = np.empty_like(k)
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
[docs]class SuperSpherical(CovModel):
r"""The Super-Spherical covariance model.
This model is derived from the relative intersection area of
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.
By default it coincides with the Hyper-Spherical model.
Notes
-----
This model is given by the following correlation function [Matern1960]_:
.. 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.
References
----------
.. [Matern1960] Matern B., "Spatial Variation",
Swedish National Institute for Forestry Research, (1960)
Other Parameters
----------------
nu : :class:`float`, optional
Shape parameter. Standard range: ``[(dim-1)/2, 50]``
Default: ``(dim-1)/2``
"""
[docs] def default_opt_arg(self):
"""Defaults for the optional arguments.
* ``{"nu": (dim-1)/2}``
Returns
-------
:class:`dict`
Defaults for optional arguments
"""
return {"nu": (self.dim - 1) / 2}
[docs] def default_opt_arg_bounds(self):
"""Defaults for boundaries of the optional arguments.
* ``{"nu": [(dim-1)/2, 50.0]}``
Returns
-------
:class:`dict`
Boundaries for optional arguments
"""
return {"nu": [(self.dim - 1) / 2, 50.0]}
[docs] def cor(self, h):
"""Super-Spherical normalized correlation function."""
h = np.asarray(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
[docs]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 function [Chiles2009]_:
.. 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}}
References
----------
.. [Chiles2009] Chiles, J. P., & Delfiner, P.,
"Geostatistics: modeling spatial uncertainty" (Vol. 497),
John Wiley & Sons. (2009)
Other Parameters
----------------
nu : :class:`float`, optional
Shape parameter. Standard range: ``[dim/2 - 1, 50]``
Default: ``dim/2``
"""
[docs] 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}
[docs] def default_opt_arg_bounds(self):
"""Defaults for boundaries of the optional arguments.
* ``{"nu": [dim/2 - 1, 50.0]}``
Returns
-------
:class:`dict`
Boundaries for optional arguments
"""
return {"nu": [self.dim / 2 - 1, 50.0]}
[docs] 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",
AttributeWarning,
)
[docs] 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
[docs] def spectral_density(self, k): # noqa: D102
k = np.asarray(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