Restore Cython code and make Rust impl. optional

Now, the GSTools-Core is used, if the package can be imported, but it
can also be switched off, by setting the global var.
`gstools.config.USE_RUST=False` during the runtime.

MANIFEST.in

@@ -1,6 +1,6 @@
 prune *
 graft tests
-recursive-include gstools *.py
+recursive-include gstools *.py *.pyx
 recursive-exclude gstools *.c *.cpp
 include LICENSE README.md pyproject.toml setup.py setup.cfg
 exclude CHANGELOG.md CONTRIBUTING.md AUTHORS.md

README.md

@@ -338,14 +338,14 @@ in memory for immediate 3D plotting in Python.
 
 - [NumPy >= 1.14.5](https://www.numpy.org)
 - [SciPy >= 1.1.0](https://www.scipy.org/scipylib)
-- [GSTools-Core >= 0.1.0](https://github.com/GeoStat-Framework/GSTools-Core)
 - [hankel >= 1.0.2](https://github.com/steven-murray/hankel)
 - [emcee >= 3.0.0](https://github.com/dfm/emcee)
 - [pyevtk >= 1.1.1](https://github.com/pyscience-projects/pyevtk)
 - [meshio>=4.0.3, <5.0](https://github.com/nschloe/meshio)
 
 ### Optional
 
+- [GSTools-Core >= 0.1.0](https://github.com/GeoStat-Framework/GSTools-Core)
 - [matplotlib](https://matplotlib.org)
 - [pyvista](https://docs.pyvista.org/)
 

docs/source/index.rst

@@ -378,7 +378,6 @@ Requirements
 
 - `Numpy >= 1.14.5 <http://www.numpy.org>`_
 - `SciPy >= 1.1.0 <http://www.scipy.org>`_
-- `GSTools-Core >= 0.1.0 <https://github.com/GeoStat-Framework/GSTools-Core>`_
 - `hankel >= 1.0.2 <https://github.com/steven-murray/hankel>`_
 - `emcee >= 3.0.0 <https://github.com/dfm/emcee>`_
 - `pyevtk >= 1.1.1 <https://github.com/pyscience-projects/pyevtk>`_
@@ -388,6 +387,7 @@ Requirements
 Optional
 --------
 
+- `GSTools-Core >= 0.1.0 <https://github.com/GeoStat-Framework/GSTools-Core>`_
 - `matplotlib <https://matplotlib.org>`_
 - `pyvista <https://docs.pyvista.org>`_
 

gstools/__init__.py

@@ -127,6 +127,7 @@
 """
 # Hooray!
 from gstools import (
+    config,
     covmodel,
     field,
     krige,

gstools/config.py

@@ -0,0 +1,14 @@
+# -*- coding: utf-8 -*-
+"""
+GStools subpackage providing global variables.
+
+.. currentmodule:: gstools.config
+
+"""
+# pylint: disable=W0611
+try:
+    import gstools_core
+
+    USE_RUST = True
+except ImportError:
+    USE_RUST = False

gstools/field/generator.py

@@ -10,16 +10,23 @@
    RandMeth
    IncomprRandMeth
 """
-# pylint: disable=C0103, W0222
+# pylint: disable=C0103, W0222, C0412
 import warnings
 from copy import deepcopy as dcp
 
 import numpy as np
 
-from gstools_core import summate, summate_incompr
+from gstools import config
 from gstools.covmodel.base import CovModel
 from gstools.random.rng import RNG
 
+if config.USE_RUST:
+    # pylint: disable=E0401
+    from gstools_core import summate, summate_incompr
+else:
+    # pylint: disable=C0412
+    from gstools.field.summator import summate, summate_incompr
+
 __all__ = ["RandMeth", "IncomprRandMeth"]
 
 

gstools/field/summator.pyx

@@ -0,0 +1,82 @@
+#cython: language_level=3, boundscheck=False, wraparound=False, cdivision=True
+# -*- coding: utf-8 -*-
+"""
+This is the randomization method summator, implemented in cython.
+"""
+
+import numpy as np
+
+cimport cython
+
+from cython.parallel import prange
+
+cimport numpy as np
+from libc.math cimport cos, sin
+
+
+def summate(
+    const double[:, :] cov_samples,
+    const double[:] z_1,
+    const double[:] z_2,
+    const double[:, :] pos
+    ):
+    cdef int i, j, d
+    cdef double phase
+    cdef int dim = pos.shape[0]
+
+    cdef int X_len = pos.shape[1]
+    cdef int N = cov_samples.shape[1]
+
+    cdef double[:] summed_modes = np.zeros(X_len, dtype=float)
+
+    for i in prange(X_len, nogil=True):
+        for j in range(N):
+            phase = 0.
+            for d in range(dim):
+                phase += cov_samples[d, j] * pos[d, i]
+            summed_modes[i] += z_1[j] * cos(phase) + z_2[j] * sin(phase)
+
+    return np.asarray(summed_modes)
+
+
+cdef (double) abs_square(const double[:] vec) nogil:
+    cdef int i
+    cdef double r = 0.
+
+    for i in range(vec.shape[0]):
+        r += vec[i]**2
+
+    return r
+
+
+def summate_incompr(
+    const double[:, :] cov_samples,
+    const double[:] z_1,
+    const double[:] z_2,
+    const double[:, :] pos
+    ):
+    cdef int i, j, d
+    cdef double phase
+    cdef double k_2
+    cdef int dim = pos.shape[0]
+
+    cdef double[:] e1 = np.zeros(dim, dtype=float)
+    e1[0] = 1.
+    cdef double[:] proj = np.empty(dim)
+
+    cdef int X_len = pos.shape[1]
+    cdef int N = cov_samples.shape[1]
+
+    cdef double[:, :] summed_modes = np.zeros((dim, X_len), dtype=float)
+
+    for i in range(X_len):
+        for j in range(N):
+            k_2 = abs_square(cov_samples[:, j])
+            phase = 0.
+            for d in range(dim):
+                phase += cov_samples[d, j] * pos[d, i]
+            for d in range(dim):
+                proj[d] = e1[d] - cov_samples[d, j] * cov_samples[0, j] / k_2
+                summed_modes[d, i] += proj[d] * (z_1[j] * cos(phase) + z_2[j] * sin(phase))
+
+    return np.asarray(summed_modes)

gstools/krige/base.py

@@ -9,23 +9,29 @@
 .. autosummary::
    Krige
 """
-# pylint: disable=C0103, W0221, E1102, R0201
+# pylint: disable=C0103, W0221, E1102, R0201, C0412
 import collections
 
 import numpy as np
 import scipy.linalg as spl
 from scipy.spatial.distance import cdist
 
-from gstools_core import (
-    calc_field_krige,
-    calc_field_krige_and_variance,
-)
+from gstools import config
 from gstools.field.base import Field
 from gstools.krige.tools import get_drift_functions, set_condition
 from gstools.tools.geometric import rotated_main_axes
 from gstools.tools.misc import eval_func
 from gstools.variogram import vario_estimate
 
+if config.USE_RUST:
+    # pylint: disable=E0401
+    from gstools_core import calc_field_krige, calc_field_krige_and_variance
+else:
+    from gstools.krige.krigesum import (
+        calc_field_krige,
+        calc_field_krige_and_variance,
+    )
+
 __all__ = ["Krige"]
 
 

gstools/krige/krigesum.pyx

@@ -0,0 +1,64 @@
+#cython: language_level=3, boundscheck=False, wraparound=False, cdivision=True
+# -*- coding: utf-8 -*-
+"""
+This is a summator for the kriging routines
+"""
+
+import numpy as np
+
+cimport cython
+from cython.parallel import prange
+cimport numpy as np
+
+
+def calc_field_krige_and_variance(
+    const double[:, :] krig_mat,
+    const double[:, :] krig_vecs,
+    const double[:] cond
+):
+
+    cdef int mat_i = krig_mat.shape[0]
+    cdef int res_i = krig_vecs.shape[1]
+
+    cdef double[:] field = np.zeros(res_i)
+    cdef double[:] error = np.zeros(res_i)
+    cdef double krig_fac
+
+    cdef int i, j, k
+
+    # error = krig_vecs * krig_mat * krig_vecs
+    # field = cond * krig_mat * krig_vecs
+    for k in prange(res_i, nogil=True):
+        for i in range(mat_i):
+            krig_fac = 0.0
+            for j in range(mat_i):
+                krig_fac += krig_mat[i, j] * krig_vecs[j, k]
+            error[k] += krig_vecs[i, k] * krig_fac
+            field[k] += cond[i] * krig_fac
+
+    return np.asarray(field), np.asarray(error)
+
+
+def calc_field_krige(
+    const double[:, :] krig_mat,
+    const double[:, :] krig_vecs,
+    const double[:] cond
+):
+
+    cdef int mat_i = krig_mat.shape[0]
+    cdef int res_i = krig_vecs.shape[1]
+
+    cdef double[:] field = np.zeros(res_i)
+    cdef double krig_fac
+
+    cdef int i, j, k
+
+    # field = cond * krig_mat * krig_vecs
+    for k in prange(res_i, nogil=True):
+        for i in range(mat_i):
+            krig_fac = 0.0
+            for j in range(mat_i):
+                krig_fac += krig_mat[i, j] * krig_vecs[j, k]
+            field[k] += cond[i] * krig_fac
+
+    return np.asarray(field)