1099186911

Build Type

github

Committed by GitHub

Commit Message

<a href="https://github.com/GeoStat-Framework/GSTools/commit/<a class=hub.com/GeoStat-Framework/GSTools/commit/1eabd8b1adc36ea239da4d7151294349ba72e649">1eabd8b1a<a href="https://github.com/GeoStat-Framework/GSTools/commit/1eabd8b1adc36ea239da4d7151294349ba72e649">">Merge </a><a class="double-link" href="https://github.com/GeoStat-Framework/GSTools/commit/<a class="double-link" href="https://github.com/GeoStat-Framework/GSTools/commit/2808a1bcb3228c02316ef0f53ed45f0c07e646d8">2808a1bcb</a>">2808a1bcb</a><a href="https://github.com/GeoStat-Framework/GSTools/commit/1eabd8b1adc36ea239da4d7151294349ba72e649"> into 21c97fc74">21c97fc74</a>

Pull Request Pull Request #197: Refactor: transform submodule; Field storage control; CondSRF optimization

Run Details

443 of 443 new or added lines in 20 files covered. (100.0%)

2770 of 2957 relevant lines covered (93.68%)

11.24 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

76.54

/gstools/tools/special.py

# -*- coding: utf-8 -*-
"""
GStools subpackage providing special functions.

.. currentmodule:: gstools.tools.special

The following functions are provided

.. autosummary::
   inc_gamma
   inc_gamma_low
   exp_int
   inc_beta
   tplstable_cor
   tpl_exp_spec_dens
   tpl_gau_spec_dens
"""
# pylint: disable=C0103, E1101
import numpy as np
from scipy import special as sps

__all__ = [
    "confidence_scaling",
    "inc_gamma",
    "inc_gamma_low",
    "exp_int",
    "inc_beta",
    "tplstable_cor",
    "tpl_exp_spec_dens",
    "tpl_gau_spec_dens",
]


# special functions ###########################################################


def confidence_scaling(per=0.95):
    """
    Scaling of standard deviation to get the desired confidence interval.

    Parameters
    ----------
    per : :class:`float`, optional
        Confidence level. The default is 0.95.

    Returns
    -------
    :class:`float`
        Scale to multiply the standard deviation with.
    """
    return np.sqrt(2) * sps.erfinv(per)


def inc_gamma(s, x):
    r"""Calculate the (upper) incomplete gamma function.

    Given by: :math:`\Gamma(s,x) = \int_x^{\infty} t^{s-1}\,e^{-t}\,{\rm d}t`

    Parameters
    ----------
    s : :class:`float`
        exponent in the integral
    x : :class:`numpy.ndarray`
        input values
    """
    if np.isclose(s, 0):
        return sps.exp1(x)
    if np.isclose(s, np.around(s)) and s < -0.5:
        return x ** s * sps.expn(int(1 - np.around(s)), x)
    if s < 0:
        return (inc_gamma(s + 1, x) - x ** s * np.exp(-x)) / s
    return sps.gamma(s) * sps.gammaincc(s, x)


def inc_gamma_low(s, x):
    r"""Calculate the lower incomplete gamma function.

    Given by: :math:`\gamma(s,x) = \int_0^x t^{s-1}\,e^{-t}\,{\rm d}t`

    Parameters
    ----------
    s : :class:`float`
        exponent in the integral
    x : :class:`numpy.ndarray`
        input values
    """
    if np.isclose(s, np.around(s)) and s < 0.5:
        return np.full_like(x, np.inf, dtype=np.double)
    if s < 0:
        return (inc_gamma_low(s + 1, x) + x ** s * np.exp(-x)) / s
    return sps.gamma(s) * sps.gammainc(s, x)


def exp_int(s, x):
    r"""Calculate the exponential integral :math:`E_s(x)`.

    Given by: :math:`E_s(x) = \int_1^\infty \frac{e^{-xt}}{t^s}\,\mathrm dt`

    Parameters
    ----------
    s : :class:`float`
        exponent in the integral (should be > -100)
    x : :class:`numpy.ndarray`
        input values
    """
    if np.isclose(s, 1):
        return sps.exp1(x)
    if np.isclose(s, np.around(s)) and s > -0.5:
        return sps.expn(int(np.around(s)), x)
    x = np.asarray(x, dtype=np.double)
    x_neg = x < 0
    x = np.abs(x)
    x_compare = x ** min((10, max(((1 - s), 1))))
    res = np.empty_like(x)
    # use asymptotic behavior for zeros
    x_zero = np.isclose(x_compare, 0, atol=1e-20)
    x_inf = x > max(30, -s / 2)  # function is like exp(-x)*(1/x + s/x^2)
    x_fin = np.logical_not(np.logical_or(x_zero, x_inf))
    x_fin_pos = np.logical_and(x_fin, np.logical_not(x_neg))
    if s > 1.0:  # limit at x=+0
        res[x_zero] = 1.0 / (s - 1.0)
    else:
        res[x_zero] = np.inf
    res[x_inf] = np.exp(-x[x_inf]) * (x[x_inf] ** -1 - s * x[x_inf] ** -2)
    res[x_fin_pos] = inc_gamma(1 - s, x[x_fin_pos]) * x[x_fin_pos] ** (s - 1)
    res[x_neg] = np.nan  # nan for x < 0
    return res


def inc_beta(a, b, x):
    r"""Calculate the incomplete Beta function.

    Given by: :math:`B(a,b;\,x) = \int_0^x t^{a-1}\,(1-t)^{b-1}\,dt`

    Parameters
    ----------
    a : :class:`float`
        first exponent in the integral
    b : :class:`float`
        second exponent in the integral
    x : :class:`numpy.ndarray`
        input values
    """
    return sps.betainc(a, b, x) * sps.beta(a, b)


def tplstable_cor(r, len_scale, hurst, alpha):
    r"""Calculate the correlation function of the TPLStable model.

    Given by the following correlation function:

    .. math::
       \rho(r) =
       \frac{2H}{\alpha} \cdot
       E_{1+\frac{2H}{\alpha}}
       \left(\left(\frac{r}{\ell}\right)^{\alpha} \right)


    Parameters
    ----------
    r : :class:`numpy.ndarray`
        input values
    len_scale : :class:`float`
        length-scale of the model.
    hurst : :class:`float`
        Hurst coefficient of the power law.
    alpha : :class:`float`, optional
        Shape parameter of the stable model.
    """
    r = np.asarray(np.abs(r / len_scale), dtype=np.double)
    r[np.isclose(r, 0)] = 0  # hack to prevent numerical errors
    res = np.ones_like(r)
    res[r > 0] = (2 * hurst / alpha) * exp_int(
        1 + 2 * hurst / alpha, (r[r > 0]) ** alpha
    )
    return res


def tpl_exp_spec_dens(k, dim, len_scale, hurst, len_low=0.0):
    r"""
    Spectal density of the TPLExponential covariance model.

    Parameters
    ----------
    k : :class:`float`
        Radius of the phase: :math:`k=\left\Vert\mathbf{k}\right\Vert`
    dim : :class:`int`
        Dimension of the model.
    len_scale : :class:`float`
        Length scale of the model.
    hurst : :class:`float`
        Hurst coefficient of the power law.
    len_low : :class:`float`, optional
        The lower length scale truncation of the model.
        Default: 0.0

    Returns
    -------
    :class:`float`
        spectal density of the TPLExponential model
    """
    if np.isclose(len_low, 0.0):
        k = np.asarray(k, dtype=np.double)
        z = (k * len_scale) ** 2
        a = hurst + dim / 2.0
        b = hurst + 0.5
        c = hurst + dim / 2.0 + 1.0
        d = dim / 2.0 + 0.5
        fac = len_scale ** dim * hurst * sps.gamma(d) / (np.pi ** d * a)
        return fac / (1.0 + z) ** a * sps.hyp2f1(a, b, c, z / (1.0 + z))
    fac_up = (len_scale + len_low) ** (2 * hurst)
    spec_up = tpl_exp_spec_dens(k, dim, len_scale + len_low, hurst)
    fac_low = len_low ** (2 * hurst)
    spec_low = tpl_exp_spec_dens(k, dim, len_low, hurst)
    return (fac_up * spec_up - fac_low * spec_low) / (fac_up - fac_low)


def tpl_gau_spec_dens(k, dim, len_scale, hurst, len_low=0.0):
    r"""
    Spectal density of the TPLGaussian covariance model.

    Parameters
    ----------
    k : :class:`float`
        Radius of the phase: :math:`k=\left\Vert\mathbf{k}\right\Vert`
    dim : :class:`int`
        Dimension of the model.
    len_scale : :class:`float`
        Length scale of the model.
    hurst : :class:`float`
        Hurst coefficient of the power law.
    len_low : :class:`float`, optional
        The lower length scale truncation of the model.
        Default: 0.0

    Returns
    -------
    :class:`float`
        spectal density of the TPLExponential model
    """
    if np.isclose(len_low, 0.0):
        k = np.asarray(k, dtype=np.double)
        z = np.array((k * len_scale / 2.0) ** 2)
        res = np.empty_like(z)
        z_gz = z > 0.1  # greater zero
        z_nz = np.logical_not(z_gz)  # near zero
        a = hurst + dim / 2.0
        fac = (len_scale / 2.0) ** dim * hurst / np.pi ** (dim / 2.0)
        res[z_gz] = fac * inc_gamma_low(a, z[z_gz]) / z[z_gz] ** a
        # first order approximation for z near zero
        res[z_nz] = fac * (1.0 / a - z[z_nz] / (a + 1.0))
        return res
    fac_up = (len_scale + len_low) ** (2 * hurst)
    spec_up = tpl_gau_spec_dens(k, dim, len_scale + len_low, hurst)
    fac_low = len_low ** (2 * hurst)
    spec_low = tpl_gau_spec_dens(k, dim, len_low, hurst)
    return (fac_up * spec_up - fac_low * spec_low) / (fac_up - fac_low)

1	# -- coding: utf-8 --
2	"""	12✔
3	GStools subpackage providing special functions.
4
5	.. currentmodule:: gstools.tools.special
6
7	The following functions are provided
8
9	.. autosummary::
10	inc_gamma
11	inc_gamma_low
12	exp_int
13	inc_beta
14	tplstable_cor
15	tpl_exp_spec_dens
16	tpl_gau_spec_dens
17	"""
18	# pylint: disable=C0103, E1101
19	import numpy as np	12✔
20	from scipy import special as sps	12✔
21
22	__all__ = [	12✔
23	"confidence_scaling",
24	"inc_gamma",
25	"inc_gamma_low",
26	"exp_int",
27	"inc_beta",
28	"tplstable_cor",
29	"tpl_exp_spec_dens",
30	"tpl_gau_spec_dens",
31	]
32
33
34	# special functions ###########################################################
35
36
37	def confidence_scaling(per=0.95):	12✔
38	"""
39	Scaling of standard deviation to get the desired confidence interval.
40
41	Parameters
42	----------
43	per : :class:`float`, optional
44	Confidence level. The default is 0.95.
45
46	Returns
47	-------
48	:class:`float`
49	Scale to multiply the standard deviation with.
50	"""
51	return np.sqrt(2) * sps.erfinv(per)	×
52
53
54	def inc_gamma(s, x):	12✔
55	r"""Calculate the (upper) incomplete gamma function.
56
57	Given by: :math:`\Gamma(s,x) = \int_x^{\infty} t^{s-1}\,e^{-t}\,{\rm d}t`
58
59	Parameters
60	----------
61	s : :class:`float`
62	exponent in the integral
63	x : :class:`numpy.ndarray`
64	input values
65	"""
66	if np.isclose(s, 0):	12✔
67	return sps.exp1(x)	×
68	if np.isclose(s, np.around(s)) and s < -0.5:	12✔
69	return x ** s * sps.expn(int(1 - np.around(s)), x)	×
70	if s < 0:	12✔
71	return (inc_gamma(s + 1, x) - x ** s * np.exp(-x)) / s	12✔
72	return sps.gamma(s) * sps.gammaincc(s, x)	12✔
73
74
75	def inc_gamma_low(s, x):	12✔
76	r"""Calculate the lower incomplete gamma function.
77
78	Given by: :math:`\gamma(s,x) = \int_0^x t^{s-1}\,e^{-t}\,{\rm d}t`
79
80	Parameters
81	----------
82	s : :class:`float`
83	exponent in the integral
84	x : :class:`numpy.ndarray`
85	input values
86	"""
87	if np.isclose(s, np.around(s)) and s < 0.5:	12✔
88	return np.full_like(x, np.inf, dtype=np.double)	×
89	if s < 0:	12✔
90	return (inc_gamma_low(s + 1, x) + x ** s * np.exp(-x)) / s	×
91	return sps.gamma(s) * sps.gammainc(s, x)	12✔
92
93
94	def exp_int(s, x):	12✔
95	r"""Calculate the exponential integral :math:`E_s(x)`.
96
97	Given by: :math:`E_s(x) = \int_1^\infty \frac{e^{-xt}}{t^s}\,\mathrm dt`
98
99	Parameters
100	----------
101	s : :class:`float`
102	exponent in the integral (should be > -100)
103	x : :class:`numpy.ndarray`
104	input values
105	"""
106	if np.isclose(s, 1):	12✔
107	return sps.exp1(x)	×
108	if np.isclose(s, np.around(s)) and s > -0.5:	12✔
109	return sps.expn(int(np.around(s)), x)	×
110	x = np.asarray(x, dtype=np.double)	12✔
111	x_neg = x < 0	12✔
112	x = np.abs(x)	12✔
113	x_compare = x ** min((10, max(((1 - s), 1))))	12✔
114	res = np.empty_like(x)	12✔
115	# use asymptotic behavior for zeros
116	x_zero = np.isclose(x_compare, 0, atol=1e-20)	12✔
117	x_inf = x > max(30, -s / 2) # function is like exp(-x)*(1/x + s/x^2)	12✔
118	x_fin = np.logical_not(np.logical_or(x_zero, x_inf))	12✔
119	x_fin_pos = np.logical_and(x_fin, np.logical_not(x_neg))	12✔
120	if s > 1.0: # limit at x=+0	12✔
121	res[x_zero] = 1.0 / (s - 1.0)	12✔
122	else:
123	res[x_zero] = np.inf	×
124	res[x_inf] = np.exp(-x[x_inf]) * (x[x_inf] ** -1 - s * x[x_inf] ** -2)	12✔
125	res[x_fin_pos] = inc_gamma(1 - s, x[x_fin_pos]) * x[x_fin_pos] ** (s - 1)	12✔
126	res[x_neg] = np.nan # nan for x < 0	12✔
127	return res	12✔
128
129
130	def inc_beta(a, b, x):	12✔
131	r"""Calculate the incomplete Beta function.
132
133	Given by: :math:`B(a,b;\,x) = \int_0^x t^{a-1}\,(1-t)^{b-1}\,dt`
134
135	Parameters
136	----------
137	a : :class:`float`
138	first exponent in the integral
139	b : :class:`float`
140	second exponent in the integral
141	x : :class:`numpy.ndarray`
142	input values
143	"""
144	return sps.betainc(a, b, x) * sps.beta(a, b)	×
145
146
147	def tplstable_cor(r, len_scale, hurst, alpha):	12✔
148	r"""Calculate the correlation function of the TPLStable model.
149
150	Given by the following correlation function:
151
152	.. math::
153	\rho(r) =
154	\frac{2H}{\alpha} \cdot
155	E_{1+\frac{2H}{\alpha}}
156	\left(\left(\frac{r}{\ell}\right)^{\alpha} \right)
157
158
159	Parameters
160	----------
161	r : :class:`numpy.ndarray`
162	input values
163	len_scale : :class:`float`
164	length-scale of the model.
165	hurst : :class:`float`
166	Hurst coefficient of the power law.
167	alpha : :class:`float`, optional
168	Shape parameter of the stable model.
169	"""
170	r = np.asarray(np.abs(r / len_scale), dtype=np.double)	12✔
171	r[np.isclose(r, 0)] = 0 # hack to prevent numerical errors	12✔
172	res = np.ones_like(r)	12✔
173	res[r > 0] = (2 * hurst / alpha) * exp_int(	12✔
174	1 + 2 * hurst / alpha, (r[r > 0]) ** alpha
175	)
176	return res	12✔
177
178
179	def tpl_exp_spec_dens(k, dim, len_scale, hurst, len_low=0.0):	12✔
180	r"""
181	Spectal density of the TPLExponential covariance model.
182
183	Parameters
184	----------
185	k : :class:`float`
186	Radius of the phase: :math:`k=\left\Vert\mathbf{k}\right\Vert`
187	dim : :class:`int`
188	Dimension of the model.
189	len_scale : :class:`float`
190	Length scale of the model.
191	hurst : :class:`float`
192	Hurst coefficient of the power law.
193	len_low : :class:`float`, optional
194	The lower length scale truncation of the model.
195	Default: 0.0
196
197	Returns
198	-------
199	:class:`float`
200	spectal density of the TPLExponential model
201	"""
202	if np.isclose(len_low, 0.0):	12✔
203	k = np.asarray(k, dtype=np.double)	12✔
204	z = (k * len_scale) ** 2	12✔
205	a = hurst + dim / 2.0	12✔
206	b = hurst + 0.5	12✔
207	c = hurst + dim / 2.0 + 1.0	12✔
208	d = dim / 2.0 + 0.5	12✔
209	fac = len_scale ** dim * hurst * sps.gamma(d) / (np.pi ** d * a)	12✔
210	return fac / (1.0 + z) ** a * sps.hyp2f1(a, b, c, z / (1.0 + z))	12✔
211	fac_up = (len_scale + len_low) ** (2 * hurst)	×
212	spec_up = tpl_exp_spec_dens(k, dim, len_scale + len_low, hurst)	×
213	fac_low = len_low ** (2 * hurst)	×
214	spec_low = tpl_exp_spec_dens(k, dim, len_low, hurst)	×
215	return (fac_up * spec_up - fac_low * spec_low) / (fac_up - fac_low)	×
216
217
218	def tpl_gau_spec_dens(k, dim, len_scale, hurst, len_low=0.0):	12✔
219	r"""
220	Spectal density of the TPLGaussian covariance model.
221
222	Parameters
223	----------
224	k : :class:`float`
225	Radius of the phase: :math:`k=\left\Vert\mathbf{k}\right\Vert`
226	dim : :class:`int`
227	Dimension of the model.
228	len_scale : :class:`float`
229	Length scale of the model.
230	hurst : :class:`float`
231	Hurst coefficient of the power law.
232	len_low : :class:`float`, optional
233	The lower length scale truncation of the model.
234	Default: 0.0
235
236	Returns
237	-------
238	:class:`float`
239	spectal density of the TPLExponential model
240	"""
241	if np.isclose(len_low, 0.0):	12✔
242	k = np.asarray(k, dtype=np.double)	12✔
243	z = np.array((k * len_scale / 2.0) ** 2)	12✔
244	res = np.empty_like(z)	12✔
245	z_gz = z > 0.1 # greater zero	12✔
246	z_nz = np.logical_not(z_gz) # near zero	12✔
247	a = hurst + dim / 2.0	12✔
248	fac = (len_scale / 2.0) ** dim * hurst / np.pi ** (dim / 2.0)	12✔
249	res[z_gz] = fac * inc_gamma_low(a, z[z_gz]) / z[z_gz] ** a	12✔
250	# first order approximation for z near zero
251	res[z_nz] = fac * (1.0 / a - z[z_nz] / (a + 1.0))	12✔
252	return res	12✔
253	fac_up = (len_scale + len_low) ** (2 * hurst)	×
254	spec_up = tpl_gau_spec_dens(k, dim, len_scale + len_low, hurst)	×
255	fac_low = len_low ** (2 * hurst)	×
256	spec_low = tpl_gau_spec_dens(k, dim, len_low, hurst)	×
257	return (fac_up * spec_up - fac_low * spec_low) / (fac_up - fac_low)	×

GeoStat-Framework / GSTools / 1099186911

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous