API: reduce size of wheels:

- merge cython source files (reduce disk overhead from generated boilerplate code)
- exclude unused Python implementation of core algorithm from distributions

.gitignore

@@ -10,6 +10,6 @@ __pycache__
 MANIFEST
 build
 dist
-src/gpgi/clib/*.html
-src/gpgi/clib/*.so
-src/gpgi/clib/*.c
+src/gpgi/_lib.html
+src/gpgi/_lib.*.so
+src/gpgi/_lib.c

MANIFEST.in

@@ -1,7 +1,8 @@
 include LICENCE
 include src/gpgi/py.typed
-include src/gpgi/clib/*.pyx
-exclude src/gpgi/clib/*.c
+include src/gpgi/_lib.pyx
+exclude src/gpgi/_lib.c
+exclude src/gpgi/_lib.py
 
 recursive-exclude tests *
 exclude .pre-commit-config.yaml

setup.py

@@ -1,33 +1,19 @@
-import os
 from distutils.extension import Extension
 
 import numpy
 from Cython.Build import cythonize
 from setuptools import setup
 
-
-def make_ext(path: str) -> Extension:
-    name, _ = os.path.splitext(path)
-    name = name.replace("/", ".")
-
-    return Extension(
-        name,
-        sources=[f"src/{path}"],
-        include_dirs=[numpy.get_include()],
-        define_macros=[
-            # keep in sync with runtime requirements (pyproject.toml)
-            ("NPY_TARGET_VERSION", "NPY_1_23_API_VERSION"),
-            ("NPY_NO_DEPRECATED_API", "NPY_1_7_API_VERSION"),
-        ],
-    )
+lib_extension = Extension(
+    "gpgi._lib",
+    sources=["src/gpgi/_lib.pyx"],
+    include_dirs=[numpy.get_include()],
+    define_macros=[
+        # keep in sync with runtime requirements (pyproject.toml)
+        ("NPY_TARGET_VERSION", "NPY_1_23_API_VERSION"),
+        ("NPY_NO_DEPRECATED_API", "NPY_1_7_API_VERSION"),
+    ],
+)
 
 
-setup(
-    ext_modules=cythonize(
-        [
-            make_ext("gpgi/clib/_indexing.pyx"),
-            make_ext("gpgi/clib/_deposition_methods.pyx"),
-        ],
-        compiler_directives={"language_level": 3},
-    ),
-)
+setup(ext_modules=cythonize([lib_extension], compiler_directives={"language_level": 3}))

src/gpgi/lib/_deposition_methods.py → src/gpgi/_lib.py

@@ -1,6 +1,80 @@
+"""Python implementations of deposition algorithms (dev only)."""
+
 import numpy as np
 
 
+def _index_particles(
+    cell_edges_x1,
+    cell_edges_x2,
+    cell_edges_x3,
+    particles_x1,
+    particles_x2,
+    particles_x3,
+    dx,
+    out,
+):
+    if cell_edges_x3.shape[0] > 1:
+        ndim = 3
+    elif cell_edges_x2.shape[0] > 1:
+        ndim = 2
+    else:
+        ndim = 1
+
+    particle_count = particles_x1.shape[0]
+
+    for ipart in range(particle_count):
+        x = particles_x1[ipart]
+        if dx[0] > 0:
+            out[ipart, 0] = int((x - cell_edges_x1[0]) // dx[0])
+        else:
+            iL = 1
+            iR = cell_edges_x1.shape[0] - 2
+            idx = (iL + iR) // 2
+            while idx != iL:
+                if cell_edges_x1[idx] > x:
+                    iR = idx
+                else:
+                    iL = idx
+                idx = (iL + iR) // 2
+            out[ipart, 0] = idx
+
+        if ndim < 2:
+            continue
+
+        x = particles_x2[ipart]
+        if dx[1] > 0:
+            out[ipart, 1] = int((x - cell_edges_x2[0]) // dx[1])
+        else:
+            iL = 1
+            iR = cell_edges_x2.shape[00] - 2
+            idx = (iL + iR) // 2
+            while idx != iL:
+                if cell_edges_x2[idx] > x:
+                    iR = idx
+                else:
+                    iL = idx
+                idx = (iL + iR) // 2
+            out[ipart, 1] = idx
+
+        if ndim < 3:
+            continue
+
+        x = particles_x3[ipart]
+        if dx[2] > 0:
+            out[ipart, 2] = int((x - cell_edges_x3[0]) // dx[2])
+        else:
+            iL = 1
+            iR = cell_edges_x3.shape[0] - 2
+            idx = (iL + iR) // 2
+            while idx != iL:
+                if cell_edges_x3[idx] > x:
+                    iR = idx
+                else:
+                    iL = idx
+                idx = (iL + iR) // 2
+            out[ipart, 2] = idx
+
+
 def _deposit_ngp(
     cell_edges_x1,
     cell_edges_x2,

src/gpgi/clib/_deposition_methods.pyx → src/gpgi/_lib.pyx

@@ -7,6 +7,104 @@ cdef fused real:
     np.float32_t
 
 
+@cython.cdivision(True)
+@cython.boundscheck(False)
+@cython.wraparound(False)
+def _index_particles(
+    np.ndarray[real, ndim=1] cell_edges_x1,
+    np.ndarray[real, ndim=1] cell_edges_x2,
+    np.ndarray[real, ndim=1] cell_edges_x3,
+    np.ndarray[real, ndim=1] particles_x1,
+    np.ndarray[real, ndim=1] particles_x2,
+    np.ndarray[real, ndim=1] particles_x3,
+    np.ndarray[real, ndim=1] dx,
+    np.ndarray[np.uint16_t, ndim=2] out,
+):
+    r"""
+    compute the host cell index (*out*)
+    The result is a mapping from particle index to containing cell index.
+    Note that cell edges are to be treated as ordained from left to right
+    and are padded with ghost layers (one for each side and direction).
+    Particles on the other hand are not assumed to be sorted, but they *are*
+    assumed to be all contained in the active domain.
+
+    dx is a 3 element array that contains the constant step used in each direction,
+    if any, or -1 as a filler value otherwise.
+
+    In directions where a constant step is used, the result is immediate (O(1))
+    Otherwise we use a bisection search (O(log(N)))
+    """
+    cdef Py_ssize_t ipart, particle_count
+    cdef np.uint16_t iL, iR, idx
+
+    # reduce the chance of floating point arithmetic errors
+    # by forcing positions be read as double precision
+    # (instead of reading as the input type)
+    cdef np.float64_t x
+
+    cdef int ndim
+    if cell_edges_x3.shape[0] > 1:
+        ndim = 3
+    elif cell_edges_x2.shape[0] > 1:
+        ndim = 2
+    else:
+        ndim = 1
+
+    particle_count = particles_x1.shape[0]
+
+    for ipart in range(particle_count):
+        x = particles_x1[ipart]
+        if dx[0] > 0:
+            out[ipart, 0] = int((x - cell_edges_x1[0]) // dx[0])
+        else:
+            iL = 1
+            iR = cell_edges_x1.shape[0] - 2
+            idx = (iL + iR) // 2
+            while idx != iL:
+                if cell_edges_x1[idx] > x:
+                    iR = idx
+                else:
+                    iL = idx
+                idx = (iL + iR) // 2
+            out[ipart, 0] = idx
+
+        if ndim < 2:
+            continue
+
+        x = particles_x2[ipart]
+        if dx[1] > 0:
+            out[ipart, 1] = int((x - cell_edges_x2[0]) // dx[1])
+        else:
+            iL = 1
+            iR = cell_edges_x2.shape[0] - 2
+            idx = (iL + iR) // 2
+            while idx != iL:
+                if cell_edges_x2[idx] > x:
+                    iR = idx
+                else:
+                    iL = idx
+                idx = (iL + iR) // 2
+            out[ipart, 1] = idx
+
+        if ndim < 3:
+            continue
+
+        x = particles_x3[ipart]
+        if dx[2] > 0:
+            out[ipart, 2] = int((x - cell_edges_x3[0]) // dx[2])
+        else:
+            iL = 1
+            iR = cell_edges_x3.shape[0] - 2
+            idx = (iL + iR) // 2
+            while idx != iL:
+                if cell_edges_x3[idx] > x:
+                    iR = idx
+                else:
+                    iL = idx
+                idx = (iL + iR) // 2
+            out[ipart, 2] = idx
+
+
 @cython.boundscheck(False)
 @cython.wraparound(False)
 def _deposit_ngp_1D(