Merge pull request #75 from GeoStat-Framework/develop

Bugfix 1.2.1

CHANGELOG.md

@@ -2,6 +2,14 @@
 
 All notable changes to **GSTools** will be documented in this file.
 
+## [1.2.1] - Volatile Violet - 2020-04-14
+
+### Bugfixes
+- `ModuleNotFoundError` is not present in py35
+- Fixing Cressie-Bug #76
+- Adding analytical formula for integral scales of rational and stable model
+- remove prange from IncomprRandMeth summators to prevent errors on Win and macOS
+
 
 ## [1.2.0] - Volatile Violet - 2020-03-20
 
@@ -114,7 +122,8 @@ All notable changes to **GSTools** will be documented in this file.
 First release of GSTools.
 
 
-[Unreleased]: https://github.com/GeoStat-Framework/gstools/compare/v1.2.0...HEAD
+[Unreleased]: https://github.com/GeoStat-Framework/gstools/compare/v1.2.1...HEAD
+[1.2.1]: https://github.com/GeoStat-Framework/gstools/compare/v1.2.0...v1.2.1
 [1.2.0]: https://github.com/GeoStat-Framework/gstools/compare/v1.1.1...v1.2.0
 [1.1.1]: https://github.com/GeoStat-Framework/gstools/compare/v1.1.0...v1.1.1
 [1.1.0]: https://github.com/GeoStat-Framework/gstools/compare/v1.0.1...v1.1.0

README.md

@@ -3,7 +3,7 @@
 [![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.1313628.svg)](https://doi.org/10.5281/zenodo.1313628)
 [![PyPI version](https://badge.fury.io/py/gstools.svg)](https://badge.fury.io/py/gstools)
 [![Conda Version](https://img.shields.io/conda/vn/conda-forge/gstools.svg)](https://anaconda.org/conda-forge/gstools)
-[![Build Status](https://travis-ci.org/GeoStat-Framework/GSTools.svg?branch=master)](https://travis-ci.org/GeoStat-Framework/GSTools)
+[![Build Status](https://travis-ci.com/GeoStat-Framework/GSTools.svg?branch=master)](https://travis-ci.com/GeoStat-Framework/GSTools)
 [![Coverage Status](https://coveralls.io/repos/github/GeoStat-Framework/GSTools/badge.svg?branch=master)](https://coveralls.io/github/GeoStat-Framework/GSTools?branch=master)
 [![Documentation Status](https://readthedocs.org/projects/gstools/badge/?version=stable)](https://geostat-framework.readthedocs.io/projects/gstools/en/stable/?badge=stable)
 [![Code style: black](https://img.shields.io/badge/code%20style-black-000000.svg)](https://github.com/ambv/black)

examples/07_transformations/02_discrete.py

@@ -1,6 +1,6 @@
 """
-discrete fields
--------------
+Discrete fields
+---------------
 
 Here we transform a field to a discrete field with values.
 If we do not give thresholds, the pairwise means of the given

gstools/__init__.py

@@ -131,7 +131,7 @@
 
 try:
     from gstools._version import __version__
-except ModuleNotFoundError:  # pragma: nocover
+except ImportError:  # pragma: nocover
     # package is not installed
     __version__ = "0.0.0.dev0"
 

gstools/covmodel/models.py

@@ -66,7 +66,7 @@ def correlation(self, r):
     def spectral_density(self, k):  # noqa: D102
         k = np.array(k, dtype=np.double)
         return (self.len_scale / np.pi) ** self.dim * np.exp(
-            -(k * self.len_scale) ** 2 / np.pi
+            -((k * self.len_scale) ** 2) / np.pi
         )
 
     def spectral_rad_cdf(self, r):
@@ -75,12 +75,12 @@ def spectral_rad_cdf(self, r):
         if self.dim == 1:
             return sps.erf(self.len_scale * r / np.sqrt(np.pi))
         if self.dim == 2:
-            return 1.0 - np.exp(-(r * self.len_scale) ** 2 / np.pi)
+            return 1.0 - np.exp(-((r * self.len_scale) ** 2) / np.pi)
         if self.dim == 3:
             return sps.erf(
                 self.len_scale * r / np.sqrt(np.pi)
             ) - 2 * r * self.len_scale / np.pi * np.exp(
-                -(r * self.len_scale) ** 2 / np.pi
+                -((r * self.len_scale) ** 2) / np.pi
             )
         return None
 
@@ -250,6 +250,14 @@ def correlation(self, r):
             1 + 0.5 / self.alpha * (r / self.len_scale) ** 2, -self.alpha
         )
 
+    def calc_integral_scale(self):  # noqa: D102
+        return (
+            self.len_scale
+            * np.sqrt(np.pi * self.alpha * 0.5)
+            * sps.gamma(self.alpha - 0.5)
+            / sps.gamma(self.alpha)
+        )
+
 
 # Stable Model ################################################################
 
@@ -325,6 +333,9 @@ def correlation(self, r):
         r = np.array(np.abs(r), dtype=np.double)
         return np.exp(-np.power(r / self.len_scale, self.alpha))
 
+    def calc_integral_scale(self):  # noqa: D102
+        return self.len_scale * sps.gamma(1.0 + 1.0 / self.alpha)
+
 
 # Matérn Model ################################################################
 
@@ -399,7 +410,7 @@ def correlation(self, r):
         r = np.array(np.abs(r), dtype=np.double)
         # for nu > 20 we just use the gaussian model
         if self.nu > 20.0:
-            return np.exp(-(r / self.len_scale) ** 2 / 4)
+            return np.exp(-((r / self.len_scale) ** 2) / 4)
         # calculate by log-transformation to prevent numerical errors
         r_gz = r[r > 0.0]
         res = np.ones_like(r)
@@ -421,7 +432,7 @@ def spectral_density(self, k):  # noqa: D102
         if self.nu > 20.0:
             return (
                 (self.len_scale / np.sqrt(np.pi)) ** self.dim
-                * np.exp(-(k * self.len_scale) ** 2)
+                * np.exp(-((k * self.len_scale) ** 2))
                 * (
                     1
                     + (

gstools/field/summator.pyx

@@ -163,7 +163,7 @@ def summate_incompr_unstruct(
 
     cdef double[:,:] summed_modes = np.zeros((dim, X_len), dtype=DTYPE)
 
-    for i in prange(X_len, nogil=True):
+    for i in range(X_len):
         for j in range(N):
             k_2 = abs_square(cov_samples[:,j])
             phase = 0.
@@ -210,7 +210,7 @@ def summate_incompr_struct_2d(
 
     cdef double[:,:,:] summed_modes = np.zeros((dim, X_len, Y_len), dtype=DTYPE)
 
-    for i in prange(X_len, nogil=True):
+    for i in range(X_len):
         for j in range(Y_len):
             for k in range(N):
                 k_2 = abs_square(cov_samples[:,k])
@@ -245,7 +245,7 @@ def summate_incompr_struct_3d(
 
     cdef double[:,:,:,:] summed_modes = np.zeros((dim, X_len, Y_len, Z_len), dtype=DTYPE)
 
-    for i in prange(X_len, nogil=True):
+    for i in range(X_len):
         for j in range(Y_len):
             for k in range(Z_len):
                 for l in range(N):

gstools/transform/field.py

@@ -80,9 +80,9 @@ def discrete(fld, values, thresholds="arithmetic"):
     fld : :any:`Field`
         Spatial Random Field class containing a generated field.
         Field will be transformed inplace.
-    values : :any:`np.ndarray`
+    values : :any:`numpy.ndarray`
         The discrete values the field will take
-    thresholds : :class:`str` or :any:`np.ndarray`, optional
+    thresholds : :class:`str` or :any:`numpy.ndarray`, optional
         the thresholds, where the value classes are separated
         possible values are:
         * "arithmetic": the mean of the 2 neighbouring values

gstools/variogram/estimator.pyx

@@ -77,7 +77,7 @@ cdef inline void normalization_cressie(
         if counts[i] == 0:
             counts[i] = 1
         variogram[i] = (
-            (1./counts[i] * variogram[i])**4 /
+            0.5 * (1./counts[i] * variogram[i])**4 /
             (0.457 + 0.494 / counts[i] + 0.045 / counts[i]**2)
         )
 

gstools/variogram/variogram.py

@@ -21,7 +21,7 @@
 
 
 def _set_estimator(estimator):
-    """Translate the verbose Python estimator identifier to single char"""
+    """Translate the verbose Python estimator identifier to single char."""
     if estimator.lower() == "matheron":
         cython_estimator = "m"
     elif estimator.lower() == "cressie":
@@ -50,16 +50,18 @@ def vario_estimate_unstructured(
        \gamma(r_k) = \frac{1}{2 N(r_k)} \sum_{i=1}^{N(r_k)} (z(\mathbf x_i) -
        z(\mathbf x_i'))^2 \; ,
 
-    with :math:`r_k \leq \| \mathbf x_i - \mathbf x_i' \| < r_{k+1}` being the bins.
+    with :math:`r_k \leq \| \mathbf x_i - \mathbf x_i' \| < r_{k+1}`
+    being the bins.
 
     Or if the estimator "cressie" was chosen:
 
     .. math::
-       \gamma(r_k) = \frac{\left(\frac{1}{N(r_k)} \sum_{i=1}^{N(r_k)}
+       \gamma(r_k) = \frac{\frac{1}{2}\left(\frac{1}{N(r_k)}\sum_{i=1}^{N(r_k)}
        \left|z(\mathbf x_i) - z(\mathbf x_i')\right|^{0.5}\right)^4}
        {0.457 + 0.494 / N(r_k) + 0.045 / N^2(r_k)} \; ,
 
-    with :math:`r_k \leq \| \mathbf x_i - \mathbf x_i' \| < r_{k+1}` being the bins.
+    with :math:`r_k \leq \| \mathbf x_i - \mathbf x_i' \| < r_{k+1}`
+    being the bins.
     The Cressie estimator is more robust to outliers.
 
     Notes
@@ -133,16 +135,18 @@ def vario_estimate_structured(field, direction="x", estimator="matheron"):
        \gamma(r_k) = \frac{1}{2 N(r_k)} \sum_{i=1}^{N(r_k)} (z(\mathbf x_i) -
        z(\mathbf x_i'))^2 \; ,
 
-    with :math:`r_k \leq \| \mathbf x_i - \mathbf x_i' \| < r_{k+1}` being the bins.
+    with :math:`r_k \leq \| \mathbf x_i - \mathbf x_i' \| < r_{k+1}`
+    being the bins.
 
     Or if the estimator "cressie" was chosen:
 
     .. math::
-       \gamma(r_k) = \frac{\left(\frac{1}{N(r_k)} \sum_{i=1}^{N(r_k)}
+       \gamma(r_k) = \frac{\frac{1}{2}\left(\frac{1}{N(r_k)}\sum_{i=1}^{N(r_k)}
        \left|z(\mathbf x_i) - z(\mathbf x_i')\right|^{0.5}\right)^4}
        {0.457 + 0.494 / N(r_k) + 0.045 / N^2(r_k)} \; ,
 
-    with :math:`r_k \leq \| \mathbf x_i - \mathbf x_i' \| < r_{k+1}` being the bins.
+    with :math:`r_k \leq \| \mathbf x_i - \mathbf x_i' \| < r_{k+1}`
+    being the bins.
     The Cressie estimator is more robust to outliers.
 
     Warnings

tests/test_variogram_structured.py

@@ -51,7 +51,7 @@ def test_list(self):
     def test_cressie_1d(self):
         z = [41.2, 40.2, 39.7, 39.2, 40.1, 38.3, 39.1, 40.0, 41.1, 40.3]
         gamma = variogram.vario_estimate_structured(z, estimator="cressie")
-        self.assertAlmostEqual(gamma[1], 1.546, places=3)
+        self.assertAlmostEqual(gamma[1], 1.546 / 2.0, places=3)
 
     def test_1d(self):
         # literature values
@@ -169,8 +169,7 @@ def test_uncorrelated_cressie_2d(self):
             field, direction="y", estimator="cressie"
         )
 
-        # TODO figure out what is going on here
-        var = 0.177
+        var = 1.0 / 12.0
         self.assertAlmostEqual(gamma_x[0], 0.0, places=1)
         self.assertAlmostEqual(gamma_x[len(gamma_x) // 2], var, places=1)
         self.assertAlmostEqual(gamma_y[0], 0.0, places=1)