From f8d92eb95e850df089e643028839d3f0dce22083 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Sebastian=20M=C3=BCller?= <mueller.seb@posteo.de>
Date: Tue, 4 Feb 2025 11:38:20 +0100
Subject: [PATCH 1/9] Cython: remove pyx files and import routines from
 gstools_cython

---
 src/gstools/field/generator.py      |   6 +-
 src/gstools/field/summator.pyx      | 130 ----------
 src/gstools/krige/base.py           |   8 +-
 src/gstools/krige/krigesum.pyx      |  84 -------
 src/gstools/variogram/estimator.pyx | 376 ----------------------------
 src/gstools/variogram/variogram.py  |   8 +-
 6 files changed, 11 insertions(+), 601 deletions(-)
 delete mode 100755 src/gstools/field/summator.pyx
 delete mode 100644 src/gstools/krige/krigesum.pyx
 delete mode 100644 src/gstools/variogram/estimator.pyx

diff --git a/src/gstools/field/generator.py b/src/gstools/field/generator.py
index d4cb0dea5..eef41f47e 100755
--- a/src/gstools/field/generator.py
+++ b/src/gstools/field/generator.py
@@ -20,12 +20,12 @@
 from copy import deepcopy as dcp
 
 import numpy as np
+from gstools_cython.field import summate as summate_c
+from gstools_cython.field import summate_fourier as summate_fourier_c
+from gstools_cython.field import summate_incompr as summate_incompr_c
 
 from gstools import config
 from gstools.covmodel.base import CovModel
-from gstools.field.summator import summate as summate_c
-from gstools.field.summator import summate_fourier as summate_fourier_c
-from gstools.field.summator import summate_incompr as summate_incompr_c
 from gstools.random.rng import RNG
 from gstools.tools.geometric import generate_grid
 
diff --git a/src/gstools/field/summator.pyx b/src/gstools/field/summator.pyx
deleted file mode 100755
index ce79ab7a0..000000000
--- a/src/gstools/field/summator.pyx
+++ /dev/null
@@ -1,130 +0,0 @@
-# cython: language_level=3, boundscheck=False, wraparound=False, cdivision=True
-"""
-This is the randomization method summator, implemented in cython.
-"""
-
-import numpy as np
-from cython.parallel import prange
-
-if OPENMP:
-    cimport openmp
-
-cimport numpy as np
-from libc.math cimport cos, sin
-
-
-def set_num_threads(num_threads):
-    cdef int num_threads_c = 1
-    if num_threads is None:
-        # OPENMP set during setup
-        if OPENMP:
-            num_threads_c = openmp.omp_get_num_procs()
-        else:
-            ...
-    else:
-        num_threads_c = num_threads
-    return num_threads_c
-
-
-def summate(
-    const double[:, :] cov_samples,
-    const double[:] z_1,
-    const double[:] z_2,
-    const double[:, :] pos,
-    num_threads=None,
-):
-    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)
-
-    cdef int num_threads_c = set_num_threads(num_threads)
-
-    for i in prange(X_len, nogil=True, num_threads=num_threads_c):
-        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,
-    num_threads=None,
-):
-    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)
-
-
-def summate_fourier(
-    const double[:] spectrum_factor,
-    const double[:, :] modes,
-    const double[:] z_1,
-    const double[:] z_2,
-    const double[:, :] pos,
-    num_threads=None,
-):
-    cdef int i, j, d
-    cdef double phase
-    cdef int dim = pos.shape[0]
-
-    cdef int X_len = pos.shape[1]
-    cdef int N = modes.shape[1]
-
-    cdef double[:] summed_modes = np.zeros(X_len, dtype=float)
-
-    cdef int num_threads_c = set_num_threads(num_threads)
-
-    for i in prange(X_len, nogil=True, num_threads=num_threads_c):
-        for j in range(N):
-            phase = 0.
-            for d in range(dim):
-                phase += modes[d, j] * pos[d, i]
-            summed_modes[i] += (
-                spectrum_factor[j] * (z_1[j] * cos(phase) + z_2[j] * sin(phase))
-            )
-
-    return np.asarray(summed_modes)
diff --git a/src/gstools/krige/base.py b/src/gstools/krige/base.py
index b30190c2f..e14830fed 100755
--- a/src/gstools/krige/base.py
+++ b/src/gstools/krige/base.py
@@ -14,14 +14,14 @@
 
 import numpy as np
 import scipy.linalg as spl
+from gstools_cython.krige import calc_field_krige as calc_field_krige_c
+from gstools_cython.krige import (
+    calc_field_krige_and_variance as calc_field_krige_and_variance_c,
+)
 from scipy.spatial.distance import cdist
 
 from gstools import config
 from gstools.field.base import Field
-from gstools.krige.krigesum import calc_field_krige as calc_field_krige_c
-from gstools.krige.krigesum import (
-    calc_field_krige_and_variance as calc_field_krige_and_variance_c,
-)
 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
diff --git a/src/gstools/krige/krigesum.pyx b/src/gstools/krige/krigesum.pyx
deleted file mode 100644
index b32c70a74..000000000
--- a/src/gstools/krige/krigesum.pyx
+++ /dev/null
@@ -1,84 +0,0 @@
-# cython: language_level=3, boundscheck=False, wraparound=False, cdivision=True
-"""
-This is a summator for the kriging routines
-"""
-
-import numpy as np
-from cython.parallel import prange
-
-if OPENMP:
-    cimport openmp
-
-cimport numpy as np
-
-
-def set_num_threads(num_threads):
-    cdef int num_threads_c = 1
-    if num_threads is None:
-        # OPENMP set during setup
-        if OPENMP:
-            num_threads_c = openmp.omp_get_num_procs()
-        else:
-            ...
-    else:
-        num_threads_c = num_threads
-    return num_threads_c
-
-
-def calc_field_krige_and_variance(
-    const double[:, :] krig_mat,
-    const double[:, :] krig_vecs,
-    const double[:] cond,
-    num_threads=None,
-):
-
-    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
-
-    cdef int num_threads_c = set_num_threads(num_threads)
-
-    # error = krig_vecs * krig_mat * krig_vecs
-    # field = cond * krig_mat * krig_vecs
-    for k in prange(res_i, nogil=True, num_threads=num_threads_c):
-        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,
-    num_threads=None,
-):
-
-    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
-
-    cdef int num_threads_c = set_num_threads(num_threads)
-
-    # field = cond * krig_mat * krig_vecs
-    for k in prange(res_i, nogil=True, num_threads=num_threads_c):
-        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)
diff --git a/src/gstools/variogram/estimator.pyx b/src/gstools/variogram/estimator.pyx
deleted file mode 100644
index 9387aa82a..000000000
--- a/src/gstools/variogram/estimator.pyx
+++ /dev/null
@@ -1,376 +0,0 @@
-# cython: language_level=3, boundscheck=False, wraparound=False, cdivision=True
-# distutils: language = c++
-"""
-This is the variogram estimater, implemented in cython.
-"""
-
-import numpy as np
-from cython.parallel import parallel, prange
-
-if OPENMP:
-    cimport openmp
-
-cimport numpy as np
-from libc.math cimport M_PI, acos, atan2, cos, fabs, isnan, pow, sin, sqrt
-
-# numpy's "bool"
-ctypedef unsigned char uint8
-
-
-def set_num_threads(num_threads):
-    cdef int num_threads_c = 1
-    if num_threads is None:
-        # OPENMP set during setup
-        if OPENMP:
-            num_threads_c = openmp.omp_get_num_procs()
-        else:
-            ...
-    else:
-        num_threads_c = num_threads
-    return num_threads_c
-
-
-cdef inline double dist_euclid(
-    const int dim,
-    const double[:, :] pos,
-    const int i,
-    const int j,
-) nogil:
-    cdef int d
-    cdef double dist_squared = 0.0
-    for d in range(dim):
-        dist_squared += ((pos[d, i] - pos[d, j]) * (pos[d, i] - pos[d, j]))
-    return sqrt(dist_squared)
-
-
-cdef inline double dist_haversine(
-    const int dim,
-    const double[:, :] pos,
-    const int i,
-    const int j,
-) nogil:
-    # pos holds lat-lon in deg
-    cdef double deg_2_rad = M_PI / 180.0
-    cdef double diff_lat = (pos[0, j] - pos[0, i]) * deg_2_rad
-    cdef double diff_lon = (pos[1, j] - pos[1, i]) * deg_2_rad
-    cdef double arg = (
-        pow(sin(diff_lat/2.0), 2) +
-        cos(pos[0, i]*deg_2_rad) *
-        cos(pos[0, j]*deg_2_rad) *
-        pow(sin(diff_lon/2.0), 2)
-    )
-    return 2.0 * atan2(sqrt(arg), sqrt(1.0-arg))
-
-
-ctypedef double (*_dist_func)(
-    const int,
-    const double[:, :],
-    const int,
-    const int,
-) nogil
-
-
-cdef inline bint dir_test(
-    const int dim,
-    const double[:, :] pos,
-    const double dist,
-    const double[:, :] direction,
-    const double angles_tol,
-    const double bandwidth,
-    const int i,
-    const int j,
-    const int d,
-) nogil:
-    cdef double s_prod = 0.0  # scalar product
-    cdef double b_dist = 0.0  # band-distance
-    cdef double tmp  # temporary variable
-    cdef int k
-    cdef bint in_band = True
-    cdef bint in_angle = True
-
-    # scalar-product calculation for bandwidth projection and angle calculation
-    for k in range(dim):
-        s_prod += (pos[k, i] - pos[k, j]) * direction[d, k]
-
-    # calculate band-distance by projection of point-pair-vec to direction line
-    if bandwidth > 0.0:
-        for k in range(dim):
-            tmp = (pos[k, i] - pos[k, j]) - s_prod * direction[d, k]
-            b_dist += tmp * tmp
-        in_band = sqrt(b_dist) < bandwidth
-
-    # allow repeating points (dist = 0)
-    if dist > 0.0:
-        # use smallest angle by taking absolute value for arccos angle formula
-        tmp = fabs(s_prod) / dist
-        if tmp < 1.0:  # else same direction (prevent numerical errors)
-            in_angle = acos(tmp) < angles_tol
-
-    return in_band and in_angle
-
-
-cdef inline double estimator_matheron(const double f_diff) nogil:
-    return f_diff * f_diff
-
-cdef inline double estimator_cressie(const double f_diff) nogil:
-    return sqrt(fabs(f_diff))
-
-ctypedef double (*_estimator_func)(const double) nogil
-
-cdef inline void normalization_matheron(
-    double[:] variogram,
-    np.int64_t[:] counts,
-):
-    cdef int i
-    for i in range(variogram.shape[0]):
-        # avoid division by zero
-        variogram[i] /= (2. * max(counts[i], 1))
-
-cdef inline void normalization_cressie(
-    double[:] variogram,
-    np.int64_t[:] counts,
-):
-    cdef int i
-    cdef np.int64_t cnt
-    for i in range(variogram.shape[0]):
-        # avoid division by zero
-        cnt = max(counts[i], 1)
-        variogram[i] = (
-            0.5 * (1./cnt * variogram[i])**4 /
-            (0.457 + 0.494 / cnt + 0.045 / cnt**2)
-        )
-
-ctypedef void (*_normalization_func)(
-    double[:],
-    np.int64_t[:],
-)
-
-cdef inline void normalization_matheron_vec(
-    double[:, :] variogram,
-    np.int64_t[:, :] counts,
-):
-    cdef int d
-    for d in range(variogram.shape[0]):
-        normalization_matheron(variogram[d, :], counts[d, :])
-
-cdef inline void normalization_cressie_vec(
-    double[:, :] variogram,
-    np.int64_t[:, :] counts,
-):
-    cdef int d
-    for d in range(variogram.shape[0]):
-        normalization_cressie(variogram[d, :], counts[d, :])
-
-ctypedef void (*_normalization_func_vec)(
-    double[:, :],
-    np.int64_t[:, :],
-)
-
-cdef _estimator_func choose_estimator_func(str estimator_type):
-    cdef _estimator_func estimator_func
-    if estimator_type == 'm':
-        estimator_func = estimator_matheron
-    else:  # estimator_type == 'c'
-        estimator_func = estimator_cressie
-    return estimator_func
-
-cdef _normalization_func choose_estimator_normalization(str estimator_type):
-    cdef _normalization_func normalization_func
-    if estimator_type == 'm':
-        normalization_func = normalization_matheron
-    else:  # estimator_type == 'c'
-        normalization_func = normalization_cressie
-    return normalization_func
-
-cdef _normalization_func_vec choose_estimator_normalization_vec(str estimator_type):
-    cdef _normalization_func_vec normalization_func_vec
-    if estimator_type == 'm':
-        normalization_func_vec = normalization_matheron_vec
-    else:  # estimator_type == 'c'
-        normalization_func_vec = normalization_cressie_vec
-    return normalization_func_vec
-
-
-def directional(
-    const double[:, :] f,
-    const double[:] bin_edges,
-    const double[:, :] pos,
-    const double[:, :] direction,  # should be normed
-    const double angles_tol=M_PI/8.0,
-    const double bandwidth=-1.0,  # negative values to turn of bandwidth search
-    const bint separate_dirs=False,  # whether the direction bands don't overlap
-    str estimator_type='m',
-    num_threads=None,
-):
-    if pos.shape[1] != f.shape[1]:
-        raise ValueError(f'len(pos) = {pos.shape[1]} != len(f) = {f.shape[1])}')
-
-    if bin_edges.shape[0] < 2:
-        raise ValueError('len(bin_edges) too small')
-
-    if angles_tol <= 0:
-        raise ValueError('tolerance for angle search masks must be > 0')
-
-    cdef _estimator_func estimator_func = choose_estimator_func(estimator_type)
-    cdef _normalization_func_vec normalization_func_vec = (
-        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
-    cdef int k_max = pos.shape[1]
-    cdef int f_max = f.shape[0]
-
-    cdef double[:, :] variogram = np.zeros((d_max, len(bin_edges)-1))
-    cdef np.int64_t[:, :] counts = np.zeros((d_max, len(bin_edges)-1), dtype=np.int64)
-    cdef int i, j, k, m, d
-    cdef double dist
-
-    cdef int num_threads_c = set_num_threads(num_threads)
-
-    for i in prange(i_max, nogil=True, num_threads=num_threads_c):
-        for j in range(j_max):
-            for k in range(j+1, k_max):
-                dist = dist_euclid(dim, pos, j, k)
-                if dist < bin_edges[i] or dist >= bin_edges[i+1]:
-                    continue  # skip if not in current bin
-                for d in range(d_max):
-                    if not dir_test(
-                      dim, pos, dist, direction, angles_tol, bandwidth, k, j, d
-                    ):
-                        continue  # skip if not in current direction
-                    for m in range(f_max):
-                        # skip no data values
-                        if not (isnan(f[m, k]) or isnan(f[m, j])):
-                            counts[d, i] += 1
-                            variogram[d, i] += estimator_func(f[m, k] - f[m, j])
-                    # once we found a fitting direction
-                    # break the search if directions are separated
-                    if separate_dirs:
-                        break
-
-    normalization_func_vec(variogram, counts)
-    return np.asarray(variogram), np.asarray(counts)
-
-
-def unstructured(
-    const double[:, :] f,
-    const double[:] bin_edges,
-    const double[:, :] pos,
-    str estimator_type='m',
-    str distance_type='e',
-    num_threads=None,
-):
-    cdef int dim = pos.shape[0]
-    cdef _dist_func distance
-
-    if distance_type == 'e':
-        distance = dist_euclid
-    else:
-        distance = dist_haversine
-        if dim != 2:
-            raise ValueError(f'Haversine: dim = {dim} != 2')
-
-    if pos.shape[1] != f.shape[1]:
-        raise ValueError(f'len(pos) = {pos.shape[1]} != len(f) = {f.shape[1])}')
-
-    if bin_edges.shape[0] < 2:
-        raise ValueError('len(bin_edges) too small')
-
-    cdef _estimator_func estimator_func = choose_estimator_func(estimator_type)
-    cdef _normalization_func normalization_func = (
-        choose_estimator_normalization(estimator_type)
-    )
-
-    cdef int i_max = bin_edges.shape[0] - 1
-    cdef int j_max = pos.shape[1] - 1
-    cdef int k_max = pos.shape[1]
-    cdef int f_max = f.shape[0]
-
-    cdef double[:] variogram = np.zeros(len(bin_edges)-1)
-    cdef np.int64_t[:] counts = np.zeros(len(bin_edges)-1, dtype=np.int64)
-    cdef int i, j, k, m
-    cdef double dist
-
-    cdef int num_threads_c = set_num_threads(num_threads)
-
-    for i in prange(i_max, nogil=True, num_threads=num_threads_c):
-        for j in range(j_max):
-            for k in range(j+1, k_max):
-                dist = distance(dim, pos, j, k)
-                if dist < bin_edges[i] or dist >= bin_edges[i+1]:
-                    continue  # skip if not in current bin
-                for m in range(f_max):
-                    # skip no data values
-                    if not (isnan(f[m, k]) or isnan(f[m, j])):
-                        counts[i] += 1
-                        variogram[i] += estimator_func(f[m, k] - f[m, j])
-
-    normalization_func(variogram, counts)
-    return np.asarray(variogram), np.asarray(counts)
-
-
-def structured(
-    const double[:, :] f,
-    str estimator_type='m',
-    num_threads=None,
-):
-    cdef _estimator_func estimator_func = choose_estimator_func(estimator_type)
-    cdef _normalization_func normalization_func = (
-        choose_estimator_normalization(estimator_type)
-    )
-
-    cdef int i_max = f.shape[0] - 1
-    cdef int j_max = f.shape[1]
-    cdef int k_max = i_max + 1
-
-    cdef double[:] variogram = np.zeros(k_max)
-    cdef np.int64_t[:] counts = np.zeros(k_max, dtype=np.int64)
-    cdef int i, j, k
-
-    cdef int num_threads_c = set_num_threads(num_threads)
-
-    with nogil, parallel(num_threads=num_threads_c):
-        for i in range(i_max):
-            for j in range(j_max):
-                for k in prange(1, k_max-i):
-                    counts[k] += 1
-                    variogram[k] += estimator_func(f[i, j] - f[i+k, j])
-
-    normalization_func(variogram, counts)
-    return np.asarray(variogram)
-
-
-def ma_structured(
-    const double[:, :] f,
-    uint8[:, :] mask,  # numpy's bools are 8bit vars
-    str estimator_type='m',
-    num_threads=None,
-):
-    cdef _estimator_func estimator_func = choose_estimator_func(estimator_type)
-    cdef _normalization_func normalization_func = (
-        choose_estimator_normalization(estimator_type)
-    )
-
-    cdef int i_max = f.shape[0] - 1
-    cdef int j_max = f.shape[1]
-    cdef int k_max = i_max + 1
-
-    cdef double[:] variogram = np.zeros(k_max)
-    cdef np.int64_t[:] counts = np.zeros(k_max, dtype=np.int64)
-    cdef int i, j, k
-
-    cdef int num_threads_c = set_num_threads(num_threads)
-
-    with nogil, parallel(num_threads=num_threads_c):
-        for i in range(i_max):
-            for j in range(j_max):
-                for k in prange(1, k_max-i):
-                    if mask[i, j] == 0 and mask[i+k, j] == 0:
-                        counts[k] += 1
-                        variogram[k] += estimator_func(f[i, j] - f[i+k, j])
-
-    normalization_func(variogram, counts)
-    return np.asarray(variogram)
diff --git a/src/gstools/variogram/variogram.py b/src/gstools/variogram/variogram.py
index 6c51f1726..dfa062b11 100644
--- a/src/gstools/variogram/variogram.py
+++ b/src/gstools/variogram/variogram.py
@@ -12,6 +12,10 @@
 
 # pylint: disable=C0412
 import numpy as np
+from gstools_cython.variogram import directional as directional_c
+from gstools_cython.variogram import ma_structured as ma_structured_c
+from gstools_cython.variogram import structured as structured_c
+from gstools_cython.variogram import unstructured as unstructured_c
 
 from gstools import config
 from gstools.normalizer.tools import remove_trend_norm_mean
@@ -23,10 +27,6 @@
     generate_grid,
 )
 from gstools.variogram.binning import standard_bins
-from gstools.variogram.estimator import directional as directional_c
-from gstools.variogram.estimator import ma_structured as ma_structured_c
-from gstools.variogram.estimator import structured as structured_c
-from gstools.variogram.estimator import unstructured as unstructured_c
 
 if config._GSTOOLS_CORE_AVAIL:  # pylint: disable=W0212; # pragma: no cover
     # pylint: disable=E0401

From 47f8a59f3ba4599078b12ac5e8dd66a17401959f Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Sebastian=20M=C3=BCller?= <mueller.seb@posteo.de>
Date: Tue, 4 Feb 2025 11:39:21 +0100
Subject: [PATCH 2/9] Switch to hatchling as build-system

---
 .github/workflows/main.yml | 43 ++++++----------------------------
 MANIFEST.in                |  4 ----
 pyproject.toml             | 48 +++++++++++++++++++-------------------
 setup.py                   | 34 ---------------------------
 4 files changed, 31 insertions(+), 98 deletions(-)
 delete mode 100644 MANIFEST.in
 delete mode 100644 setup.py

diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml
index 8ff7b1f80..6efa0ad5d 100644
--- a/.github/workflows/main.yml
+++ b/.github/workflows/main.yml
@@ -49,10 +49,6 @@ jobs:
         run: |
           python -m pylint src/gstools/
 
-      - name: cython-lint check
-        run: |
-          cython-lint src/gstools/
-
       - name: coveralls check
         run: |
           python -m pytest --cov gstools --cov-report term-missing -v tests/
@@ -60,30 +56,6 @@ jobs:
         env:
           GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
 
-  build_wheels:
-    name: wheels for ${{ matrix.os }}
-    runs-on: ${{ matrix.os }}
-    strategy:
-      fail-fast: false
-      matrix:
-        # macos-13 is an intel runner, macos-14 is apple silicon
-        os: [ubuntu-latest, windows-latest, macos-13, macos-14]
-
-    steps:
-      - uses: actions/checkout@v4
-        with:
-          fetch-depth: "0"
-
-      - name: Build wheels
-        uses: pypa/cibuildwheel@v2.22.0
-        with:
-          output-dir: dist-wheel-${{ matrix.os }}
-
-      - uses: actions/upload-artifact@v4
-        with:
-          name: dist-wheel-${{ matrix.os }}
-          path: ./dist-wheel-${{ matrix.os }}/*.whl
-
   build_sdist:
     name: sdist on ${{ matrix.os }} with py ${{ matrix.ver.py }} numpy${{ matrix.ver.np }} scipy${{ matrix.ver.sp }}
     runs-on: ${{ matrix.os }}
@@ -123,8 +95,9 @@ jobs:
           pip install build
 
       - name: Install GSTools
-        env:
-          GSTOOLS_BUILD_PARALLEL: 1
+        # env:
+        #   GSTOOLS_BUILD_PARALLEL: 1
+        #   pip install --no-binary gstools-cython -v --editable .[test]
         run: |
           pip install -v --editable .[test]
 
@@ -136,13 +109,13 @@ jobs:
       - name: Build sdist
         run: |
           # PEP 517 package builder from pypa
-          python -m build --sdist --outdir dist-sdist .
+          python -m build .
 
       - uses: actions/upload-artifact@v4
         if: matrix.os == 'ubuntu-latest' && matrix.ver.py == '3.11'
         with:
-          name: dist-sdist
-          path: dist-sdist/*.tar.gz
+          name: dist
+          path: dist/*
 
   upload_to_pypi:
     needs: [build_wheels, build_sdist]
@@ -151,9 +124,7 @@ jobs:
     steps:
       - uses: actions/download-artifact@v4
         with:
-          pattern: dist-*
-          merge-multiple: true
-          path: dist
+          name: dist
 
       - name: Publish to Test PyPI
         # only if working on main
diff --git a/MANIFEST.in b/MANIFEST.in
deleted file mode 100644
index 24184482a..000000000
--- a/MANIFEST.in
+++ /dev/null
@@ -1,4 +0,0 @@
-prune **
-recursive-include tests *.py
-recursive-include src/gstools *.py *.pyx
-include AUTHORS.md LICENSE README.md pyproject.toml setup.py
diff --git a/pyproject.toml b/pyproject.toml
index 4e98c50d2..8b82c2bf1 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -1,13 +1,11 @@
+
+
 [build-system]
 requires = [
-    "setuptools>=64",
-    "setuptools_scm>=7",
-    "numpy>=2.0.0rc1,<2.3; python_version >= '3.9'",
-    "oldest-supported-numpy; python_version < '3.9'",
-    "Cython>=3.0.10,<3.1.0",
-    "extension-helpers>=1",
+    "hatchling>=1.8.0",
+    "hatch-vcs",
 ]
-build-backend = "setuptools.build_meta"
+build-backend = "hatchling.build"
 
 [project]
 requires-python = ">=3.8"
@@ -17,7 +15,7 @@ authors = [
     {name = "Sebastian Müller, Lennart Schüler", email = "info@geostat-framework.org"},
 ]
 readme = "README.md"
-license = {text = "LGPL-3.0"}
+license = "LGPL-3.0"
 dynamic = ["version"]
 classifiers = [
     "Development Status :: 5 - Production/Stable",
@@ -53,6 +51,7 @@ dependencies = [
     "numpy>=1.20.0",
     "pyevtk>=1.1.1",
     "scipy>=1.1.0",
+    "gstools-cython>=1.0.0rc1",
 ]
 
 [project.optional-dependencies]
@@ -78,7 +77,6 @@ lint = [
     "black>=24",
     "pylint",
     "isort[colors]",
-    "cython-lint",
 ]
 
 [project.urls]
@@ -89,14 +87,25 @@ Homepage = "https://geostat-framework.org/#gstools"
 Source = "https://github.com/GeoStat-Framework/GSTools"
 Tracker = "https://github.com/GeoStat-Framework/GSTools/issues"
 
-[tool.setuptools]
-license-files = ["LICENSE"]
+[tool.hatch.version]
+source = "vcs"
+fallback_version = "0.0.0.dev0"
 
-[tool.setuptools_scm]
-write_to = "src/gstools/_version.py"
-write_to_template = "__version__ = '{version}'"
+[tool.hatch.version.raw-options]
 local_scheme = "no-local-version"
-fallback_version = "0.0.0.dev0"
+
+[tool.hatch.build.hooks.vcs]
+version-file = "src/gstools/_version.py"
+template = "__version__ = '{version}'"
+
+[tool.hatch.build.targets.sdist]
+include = [
+  "/src",
+  "/tests",
+]
+
+[tool.hatch.build.targets.wheel]
+packages = ["src/gstools"]
 
 [tool.isort]
 profile = "black"
@@ -160,12 +169,3 @@ target-version = [
     max-attributes = 25
     max-public-methods = 80
     max-positional-arguments=20
-
-[tool.cibuildwheel]
-# Switch to using build
-build-frontend = "build"
-# Disable building PyPy wheels on all platforms, 32bit for py3.10/11/12/13, musllinux builds, py3.6/7
-skip = ["cp36-*", "cp37-*", "pp*", "*-win32", "*-manylinux_i686", "*-musllinux_*"]
-# Run the package tests using `pytest`
-test-extras = "test"
-test-command = "pytest -v {package}/tests"
diff --git a/setup.py b/setup.py
deleted file mode 100644
index b27548a94..000000000
--- a/setup.py
+++ /dev/null
@@ -1,34 +0,0 @@
-"""GSTools: A geostatistical toolbox."""
-
-import os
-
-import numpy as np
-from Cython.Build import cythonize
-from extension_helpers import add_openmp_flags_if_available
-from setuptools import Extension, setup
-
-# cython extensions
-CY_MODULES = [
-    Extension(
-        name=f"gstools.{ext}",
-        sources=[os.path.join("src", "gstools", *ext.split(".")) + ".pyx"],
-        include_dirs=[np.get_include()],
-        define_macros=[("NPY_NO_DEPRECATED_API", "NPY_1_7_API_VERSION")],
-    )
-    for ext in ["field.summator", "variogram.estimator", "krige.krigesum"]
-]
-# you can set GSTOOLS_BUILD_PARALLEL=0 or GSTOOLS_BUILD_PARALLEL=1
-open_mp = False
-if int(os.getenv("GSTOOLS_BUILD_PARALLEL", "0")):
-    added = [add_openmp_flags_if_available(mod) for mod in CY_MODULES]
-    if any(added):
-        open_mp = True
-    print(f"## GSTools setup: OpenMP used: {open_mp}")
-else:
-    print("## GSTools setup: OpenMP not wanted by the user.")
-
-# setup - do not include package data to ignore .pyx files in wheels
-setup(
-    ext_modules=cythonize(CY_MODULES, compile_time_env={"OPENMP": open_mp}),
-    include_package_data=False,
-)

From 13f073da60609a30de715c3f0ce735ee977a38d9 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Sebastian=20M=C3=BCller?= <mueller.seb@posteo.de>
Date: Tue, 4 Feb 2025 11:47:24 +0100
Subject: [PATCH 3/9] update docs for cython backend usage

---
 docs/source/index.rst | 9 +++++----
 1 file changed, 5 insertions(+), 4 deletions(-)

diff --git a/docs/source/index.rst b/docs/source/index.rst
index bba2309e8..37a081bdf 100644
--- a/docs/source/index.rst
+++ b/docs/source/index.rst
@@ -90,21 +90,22 @@ When using conda, the parallel version of GSTools is installed per default.
 To enable the OpenMP support in Cython when using pip, 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). GSTools then needs to be installed from source:
+if not present). The GSTools-Cython backend then needs to be installed from source:
 
 .. code-block:: none
 
     export GSTOOLS_BUILD_PARALLEL=1
-    pip install --no-binary=gstools gstools
+    pip install --no-binary=gstools-cython gstools
 
-Note, that the ``--no-binary=gstools`` option forces pip to not use a wheel
-for GSTools.
+Note, that the ``--no-binary=gstools-cython`` option forces pip to not use a wheel
+for the GSTools-Cython backend.
 
 For the development version, you can do almost the same:
 
 .. code-block:: none
 
     export GSTOOLS_BUILD_PARALLEL=1
+    pip install git+git://github.com/GeoStat-Framework/GSTools-Cython.git@main
     pip install git+git://github.com/GeoStat-Framework/GSTools.git@main
 
 

From d4db5b398c791f4a5feb156d0e3fcca149ec4d5f Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Sebastian=20M=C3=BCller?= <mueller.seb@posteo.de>
Date: Tue, 4 Feb 2025 11:50:06 +0100
Subject: [PATCH 4/9] CI: fix

---
 .github/workflows/main.yml | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml
index 6efa0ad5d..1f8c887c8 100644
--- a/.github/workflows/main.yml
+++ b/.github/workflows/main.yml
@@ -56,7 +56,7 @@ jobs:
         env:
           GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
 
-  build_sdist:
+  build:
     name: sdist on ${{ matrix.os }} with py ${{ matrix.ver.py }} numpy${{ matrix.ver.np }} scipy${{ matrix.ver.sp }}
     runs-on: ${{ matrix.os }}
     strategy:
@@ -106,7 +106,7 @@ jobs:
           pip install "numpy${{ matrix.ver.np }}" "scipy${{ matrix.ver.sp }}"
           python -m pytest -v tests/
 
-      - name: Build sdist
+      - name: Build dist
         run: |
           # PEP 517 package builder from pypa
           python -m build .
@@ -118,7 +118,7 @@ jobs:
           path: dist/*
 
   upload_to_pypi:
-    needs: [build_wheels, build_sdist]
+    needs: [build]
     runs-on: ubuntu-latest
 
     steps:

From 20967ca1d732ea44599702434d6854ff820cc5d2 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Sebastian=20M=C3=BCller?= <mueller.seb@posteo.de>
Date: Fri, 7 Feb 2025 21:51:31 +0100
Subject: [PATCH 5/9] CI: rename test job

---
 .github/workflows/main.yml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml
index 1f8c887c8..b8d573c6f 100644
--- a/.github/workflows/main.yml
+++ b/.github/workflows/main.yml
@@ -57,7 +57,7 @@ jobs:
           GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
 
   build:
-    name: sdist on ${{ matrix.os }} with py ${{ matrix.ver.py }} numpy${{ matrix.ver.np }} scipy${{ matrix.ver.sp }}
+    name: test on ${{ matrix.os }} with py ${{ matrix.ver.py }} numpy${{ matrix.ver.np }} scipy${{ matrix.ver.sp }}
     runs-on: ${{ matrix.os }}
     strategy:
       fail-fast: false

From 093d03d7de3d59c82d31f0f6a8664cac1760c989 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Sebastian=20M=C3=BCller?= <mueller.seb@posteo.de>
Date: Mon, 10 Feb 2025 11:01:12 +0100
Subject: [PATCH 6/9] Dep: gstools-cython now released, cap with v2

---
 pyproject.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/pyproject.toml b/pyproject.toml
index 8b82c2bf1..88fbffebe 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -45,13 +45,13 @@ classifiers = [
     "Topic :: Utilities",
 ]
 dependencies = [
+    "gstools-cython>=1,<2",
     "emcee>=3.0.0",
     "hankel>=1.0.0",
     "meshio>=5.1.0",
     "numpy>=1.20.0",
     "pyevtk>=1.1.1",
     "scipy>=1.1.0",
-    "gstools-cython>=1.0.0rc1",
 ]
 
 [project.optional-dependencies]

From ca3f28d32f1fec9db4c03496cb81ab1d13a44337 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Sebastian=20M=C3=BCller?= <mueller.seb@posteo.de>
Date: Mon, 10 Feb 2025 11:12:53 +0100
Subject: [PATCH 7/9] CI: test on ubuntu-24.04-arm

---
 .github/workflows/main.yml | 5 +----
 1 file changed, 1 insertion(+), 4 deletions(-)

diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml
index b8d573c6f..cf9a3e996 100644
--- a/.github/workflows/main.yml
+++ b/.github/workflows/main.yml
@@ -62,7 +62,7 @@ jobs:
     strategy:
       fail-fast: false
       matrix:
-        os: [ubuntu-latest, windows-latest, macos-13, macos-14]
+        os: [ubuntu-latest, ubuntu-24.04-arm, windows-latest, macos-13, macos-14]
         # https://github.com/scipy/oldest-supported-numpy/blob/main/setup.cfg
         ver:
           - { py: "3.8", np: "==1.20.0", sp: "==1.5.4" }
@@ -95,9 +95,6 @@ jobs:
           pip install build
 
       - name: Install GSTools
-        # env:
-        #   GSTOOLS_BUILD_PARALLEL: 1
-        #   pip install --no-binary gstools-cython -v --editable .[test]
         run: |
           pip install -v --editable .[test]
 

From 80765b253c5f4f43612da353403a7ea0db98194e Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Sebastian=20M=C3=BCller?= <mueller.seb@posteo.de>
Date: Mon, 10 Feb 2025 11:28:09 +0100
Subject: [PATCH 8/9] CI: skip py38 on linux arm since gstools-cython wheels
 will be for >py39 there

---
 .github/workflows/main.yml | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml
index cf9a3e996..4665fbf8f 100644
--- a/.github/workflows/main.yml
+++ b/.github/workflows/main.yml
@@ -73,6 +73,8 @@ jobs:
           - { py: "3.13", np: "==2.1.0", sp: "==1.14.1" }
           - { py: "3.13", np: ">=2.1.0", sp: ">=1.14.1" }
         exclude:
+          - os: ubuntu-24.04-arm
+            ver: { py: "3.8", np: "==1.20.0", sp: "==1.5.4" }
           - os: macos-14
             ver: { py: "3.8", np: "==1.20.0", sp: "==1.5.4" }
           - os: macos-14

From d5a6338d89d31dd5b9c49a2bab2619bf5e44f96b Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Sebastian=20M=C3=BCller?= <mueller.seb@posteo.de>
Date: Mon, 28 Apr 2025 09:41:21 +0200
Subject: [PATCH 9/9] docs: polish section for parallel support

---
 docs/source/index.rst | 26 ++++++++++++--------------
 1 file changed, 12 insertions(+), 14 deletions(-)

diff --git a/docs/source/index.rst b/docs/source/index.rst
index 37a081bdf..23f3b179f 100644
--- a/docs/source/index.rst
+++ b/docs/source/index.rst
@@ -87,16 +87,17 @@ When using conda, the parallel version of GSTools is installed per default.
 
 ***Parallelizing Cython***
 
-To enable the OpenMP support in Cython when using pip, 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). The GSTools-Cython backend then needs to be installed from source:
+For parallel support, the `GSTools-Cython <https://github.com/GeoStat-Framework/GSTools-Cython>`_
+backend needs to be compiled from source the following way:
 
 .. code-block:: none
 
     export GSTOOLS_BUILD_PARALLEL=1
     pip install --no-binary=gstools-cython gstools
 
+You have to provide a C compiler and OpenMP to compile GSTools-Cython with parallel support.
+The feature is controlled by the environment variable
+``GSTOOLS_BUILD_PARALLEL``, that can be ``0`` or ``1`` (interpreted as ``0`` if not present).
 Note, that the ``--no-binary=gstools-cython`` option forces pip to not use a wheel
 for the GSTools-Cython backend.
 
@@ -112,14 +113,7 @@ For the development version, you can do almost the same:
 ***Using GSTools-Core for parallelization and even more speed***
 
 You can install the optional dependency `GSTools-Core <https://github.com/GeoStat-Framework/GSTools-Core>`_,
-which is a re-implementation of the 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_GSTOOLS_CORE=False` in your code. This also works at
-runtime. You can install the optional dependency e.g. by
+which is a re-implementation of GSTools-Cython:
 
 .. code-block:: none
 
@@ -131,8 +125,12 @@ or by manually installing the package
 
     pip install gstools-core
 
-GSTools-Core will automatically run in parallel, without having to use provide
-OpenMP or a local C compiler.
+The new package uses the language Rust and it should be safer and faster (in some cases by orders of magnitude).
+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_GSTOOLS_CORE=False` in your code. This also works at runtime.
+
+GSTools-Core will automatically run in parallel, without having to provide OpenMP or a local C compiler.
 
 
 Citation