diff --git a/CHANGELOG.md b/CHANGELOG.md index 95ad2e2c..82b9afd4 100755 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -3,6 +3,15 @@ All notable changes to **GSTools** will be documented in this file. +## [1.3.4] - Pure Pink ? + +### Enhancements +- add GStools-Core as optional dependency [#215](https://github.com/GeoStat-Framework/GSTools/pull/215) + +### Changes +- remove unnecessary `dim` argument in Cython code [#216](https://github.com/GeoStat-Framework/GSTools/issues/216) + + ## [1.3.3] - Pure Pink - 2021-08 ### Enhancements diff --git a/README.md b/README.md index 7ff21bf8..f68e49a9 100644 --- a/README.md +++ b/README.md @@ -345,6 +345,7 @@ in memory for immediate 3D plotting in Python. ### Optional +- [GSTools-Core >= 0.2.0](https://github.com/GeoStat-Framework/GSTools-Core) - [matplotlib](https://matplotlib.org) - [pyvista](https://docs.pyvista.org/) diff --git a/docs/source/index.rst b/docs/source/index.rst index f66ca738..8ead5056 100644 --- a/docs/source/index.rst +++ b/docs/source/index.rst @@ -70,6 +70,8 @@ To get the latest development version you can install it directly from GitHub: If something went wrong during installation, try the :code:`-I` `flag from pip `_. +**Speeding up GSTools by parallelization** + To enable the OpenMP support, you have to provide a C compiler and OpenMP. Parallel support is controlled by an environment variable ``GSTOOLS_BUILD_PARALLEL``, that can be ``0`` or ``1`` (interpreted as ``0`` if not present). @@ -89,6 +91,33 @@ For the development version, you can do almost the same: export GSTOOLS_BUILD_PARALLEL=1 pip install git+git://github.com/GeoStat-Framework/GSTools.git@main +**Using experimental GSTools-Core for even more speed** + +You can install the optional dependency `GSTools-Core `_, +which is a re-implementation of the main algorithms used in GSTools. The new +package uses the language Rust and it should be faster (in some cases by orders +of magnitude), safer, and it will potentially completely replace the current +standard implementation in Cython. Once the package GSTools-Core is available +on your machine, it will be used by default. In case you want to switch back to +the Cython implementation, you can set :code:`gstools.config.USE_RUST=False` in +your code. This also works at runtime. You can install the optional dependency +e.g. by + +.. code-block:: none + + pip install gstools[rust] + +or by manually installing the package + +.. code-block:: none + + pip install gstools-core + +GSTools-Core will automatically use all your cores in parallel, without having +to use OpenMP or a local C compiler. +In case you want to restrict the number of threads used, you can set the +environment variable ``RAYON_NUM_THREADS`` to the desired amount. + Citation ======== @@ -387,6 +416,7 @@ Requirements Optional -------- +- `GSTools-Core >= 0.2.0 `_ - `matplotlib `_ - `pyvista `_ diff --git a/gstools/__init__.py b/gstools/__init__.py index 77281d8d..a82839cc 100644 --- a/gstools/__init__.py +++ b/gstools/__init__.py @@ -127,6 +127,7 @@ """ # Hooray! from gstools import ( + config, covmodel, field, krige, diff --git a/gstools/config.py b/gstools/config.py new file mode 100644 index 00000000..38f3e7fa --- /dev/null +++ b/gstools/config.py @@ -0,0 +1,14 @@ +# -*- coding: utf-8 -*- +""" +GStools subpackage providing global variables. + +.. currentmodule:: gstools.config + +""" +# pylint: disable=W0611 +try: # pragma: no cover + import gstools_core + + USE_RUST = True +except ImportError: + USE_RUST = False diff --git a/gstools/field/generator.py b/gstools/field/generator.py index 66f11e99..31cc7eea 100644 --- a/gstools/field/generator.py +++ b/gstools/field/generator.py @@ -10,16 +10,22 @@ 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 import config from gstools.covmodel.base import CovModel -from gstools.field.summator import summate, summate_incompr from gstools.random.rng import RNG +if config.USE_RUST: # pragma: no cover + # pylint: disable=E0401 + from gstools_core import summate, summate_incompr +else: + from gstools.field.summator import summate, summate_incompr + __all__ = ["RandMeth", "IncomprRandMeth"] diff --git a/gstools/krige/base.py b/gstools/krige/base.py index 1e217185..8accf508 100755 --- a/gstools/krige/base.py +++ b/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 import config from gstools.field.base import Field -from gstools.krige.krigesum import ( - calc_field_krige, - calc_field_krige_and_variance, -) 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: # pragma: no cover + # 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"] diff --git a/gstools/variogram/estimator.pyx b/gstools/variogram/estimator.pyx index 8c149b41..6c5a67df 100644 --- a/gstools/variogram/estimator.pyx +++ b/gstools/variogram/estimator.pyx @@ -178,7 +178,6 @@ cdef _normalization_func_vec choose_estimator_normalization_vec(str estimator_ty def directional( - const int dim, const double[:,:] f, const double[:] bin_edges, const double[:,:] pos, @@ -203,6 +202,7 @@ def directional( choose_estimator_normalization_vec(estimator_type) ) + cdef int dim = pos.shape[0] cdef int d_max = direction.shape[0] cdef int i_max = bin_edges.shape[0] - 1 cdef int j_max = pos.shape[1] - 1 @@ -237,13 +237,13 @@ def directional( return np.asarray(variogram), np.asarray(counts) def unstructured( - const int dim, const double[:,:] f, const double[:] bin_edges, const double[:,:] pos, str estimator_type='m', str distance_type='e', ): + cdef int dim = pos.shape[0] cdef _dist_func distance if distance_type == 'e': diff --git a/gstools/variogram/variogram.py b/gstools/variogram/variogram.py index 7c107695..7acce928 100644 --- a/gstools/variogram/variogram.py +++ b/gstools/variogram/variogram.py @@ -10,8 +10,10 @@ vario_estimate vario_estimate_axis """ +# pylint: disable=C0412 import numpy as np +from gstools import config from gstools.normalizer.tools import remove_trend_norm_mean from gstools.tools.geometric import ( ang2dir, @@ -20,12 +22,20 @@ generate_grid, ) from gstools.variogram.binning import standard_bins -from gstools.variogram.estimator import ( - directional, - ma_structured, - structured, - unstructured, -) + +if config.USE_RUST: # pragma: no cover + # pylint: disable=E0401 + from gstools_core import variogram_directional as directional + from gstools_core import variogram_ma_structured as ma_structured + from gstools_core import variogram_structured as structured + from gstools_core import variogram_unstructured as unstructured +else: + from gstools.variogram.estimator import ( + directional, + ma_structured, + structured, + unstructured, + ) __all__ = [ "vario_estimate", @@ -59,7 +69,8 @@ def _separate_dirs_test(direction, angles_tol): for j in range(i + 1, direction.shape[0]): s_prod = np.minimum(np.abs(np.dot(direction[i], direction[j])), 1) separate_dirs &= np.arccos(s_prod) >= 2 * angles_tol - return separate_dirs + # gstools-core doesn't like the type `numpy.bool_` + return bool(separate_dirs) def vario_estimate( @@ -341,7 +352,6 @@ def vario_estimate( # "h"aversine or "e"uclidean distance type distance_type = "h" if latlon else "e" estimates, counts = unstructured( - dim, field, bin_edges, pos, @@ -350,7 +360,6 @@ def vario_estimate( ) else: estimates, counts = directional( - dim, field, bin_edges, pos, @@ -441,7 +450,9 @@ def vario_estimate_axis( field = np.ma.array(field, ndmin=1, dtype=np.double) if missing: field.mask = np.logical_or(field.mask, missing_mask) - mask = np.asarray(np.ma.getmaskarray(field), dtype=np.int32) + mask = np.ma.getmaskarray(field) + if not config.USE_RUST: + mask = np.asarray(mask, dtype=np.int32) else: field = np.array(field, ndmin=1, dtype=np.double, copy=False) missing_mask = None # free space diff --git a/setup.cfg b/setup.cfg index efc2fd4c..5a59f87c 100644 --- a/setup.cfg +++ b/setup.cfg @@ -66,6 +66,8 @@ doc = plotting = matplotlib>=3,<4 pyvista>=0.29,<1 +rust = + gstools_core>=0.2.0,<1 test = coverage[toml]>=5.2.1,<6 pytest>=6.0,<7