cytools (#382)

.gitignore

@@ -122,5 +122,6 @@ rever/
 
 # Cython byproducts
 versioned_hdf5/*.c
+versioned_hdf5/*.cpp
 versioned_hdf5/*.html
 versioned_hdf5/*.so

meson.build

@@ -9,6 +9,4 @@ project(
 py = import('python').find_installation(pure: false)
 cy = meson.get_compiler('cython')
 
-add_project_arguments('-std=c++11', language: 'cpp')
-
 subdir('versioned_hdf5')

versioned_hdf5/cytools.pxd

@@ -0,0 +1,38 @@
+# This file allows cimport'ing the functions and types declared below from other
+# Cython modules
+
+from cython cimport Py_ssize_t
+from libc.stdint cimport uint64_t
+
+# Centralized definition of hsize_t, in accordance with libhd5:
+# https://github.com/HDFGroup/hdf5/blob/6b43197b0817596f47670c6b55d26ff7f86d6bd9/src/H5public.h#L301
+#
+# versioned_hdf5 uses the same datatype for indexing as libhdf5. Notably, this differs
+# from numpy's npy_intp (same as Py_ssize_t, ssize_t, and signed long), which is only 32
+# bit wide on 32 bit platforms, to allow indexing datasets with >=2**31 points per axis
+# on disk, as long as you don't load them in memory all at once.
+#
+# Note that hsize_t is unsigned, which can lead to integer underflows.
+#
+# The definition of hsize_t has changed over time in the libhdf5 headers, as noted
+# below.
+# C99 dictates that long long is *at least* 64 bit wide, while uint64_t is *exactly* 64
+# bit wide. So the two definitions are de facto equivalent.
+# The below definition is inconsequential, as it is overridden by whatever version of
+# hdf5.h is installed thanks to the 'cdef extern from "hdf5.h"' block.
+#
+# h5py gets this wildly wrong, as it defines hsize_t as long long, which is signed, and
+# confusingly also defines hssize_t as signed long long - which is an alias for
+# long long:
+# https://github.com/h5py/h5py/blob/eaa9d93cc7620f3e7d8cff6f3a739b7d9834aef7/h5py/api_types_hdf5.pxd#L21-L22
+
+cdef extern from "hdf5.h":
+    # ctypedef unsigned long long hsize_t  # hdf5 <=1.12
+    ctypedef uint64_t hsize_t  # hdf5 >=1.14
+
+cpdef hsize_t stop2count(hsize_t start, hsize_t stop, hsize_t step) noexcept nogil
+cpdef hsize_t count2stop(hsize_t start, hsize_t count, hsize_t step) noexcept nogil
+cpdef Py_ssize_t ceil_a_over_b(Py_ssize_t a, Py_ssize_t b) noexcept nogil
+cpdef Py_ssize_t smallest_step_after(
+    Py_ssize_t x, Py_ssize_t a, Py_ssize_t m
+) noexcept nogil

versioned_hdf5/cytools.py

@@ -0,0 +1,73 @@
+from __future__ import annotations
+
+import cython
+from cython import Py_ssize_t
+
+if cython.compiled:  # pragma: nocover
+    # This is repeated here from the header in order to silence mypy, which now treats
+    # hsize_t as Any, without breaking cython.
+    from cython.cimports.versioned_hdf5.cytools import hsize_t
+
+# numpy equivalent dtype for hsize_t
+from numpy import uint64 as np_hsize_t  # noqa: F401
+
+
+@cython.ccall
+@cython.nogil
+@cython.exceptval(check=False)
+def stop2count(start: hsize_t, stop: hsize_t, step: hsize_t) -> hsize_t:
+    """Given a start:stop:step slice or range, return the number of elements yielded.
+
+    This is functionally identical to::
+
+        len(range(start, stop, step))
+
+    Doesn't assume that stop >= start. Assumes that step >= 1.
+    """
+    # Note that hsize_t is unsigned so stop - start could underflow.
+    if stop <= start:
+        return 0
+    return (stop - start - 1) // step + 1
+
+
+@cython.ccall
+@cython.nogil
+@cython.exceptval(check=False)
+def count2stop(start: hsize_t, count: hsize_t, step: hsize_t) -> hsize_t:
+    """Inverse of stop2count.
+
+    When count == 0 or when step>1, multiple stops can yield the same count.
+    This function returns the smallest stop >= start.
+    """
+    if count == 0:
+        return start
+    return start + (count - 1) * step + 1
+
+
+@cython.ccall
+@cython.nogil
+@cython.exceptval(check=False)
+def ceil_a_over_b(a: Py_ssize_t, b: Py_ssize_t) -> Py_ssize_t:
+    """Returns ceil(a/b). Assumes a >= 0 and b > 0.
+
+    Note
+    ----
+    This module is compiled with the cython.cdivision flag. This causes behaviour to
+    change if a and b have opposite signs and you try debugging the module in pure
+    python, without compiling it. This function blindly assumes that a and b are always
+    the same sign.
+    """
+    return a // b + (a % b > 0)
+
+
+@cython.ccall
+@cython.nogil
+@cython.exceptval(check=False)
+def smallest_step_after(x: Py_ssize_t, a: Py_ssize_t, m: Py_ssize_t) -> Py_ssize_t:
+    """Find the smallest integer y >= x where y = a + k*m for whole k's
+    Assumes 0 <= a <= x and m >= 1.
+
+    a                  x    y
+    | <-- m --> | <-- m --> |
+    """
+    return a + ceil_a_over_b(x - a, m) * m

versioned_hdf5/cytools.pyx

@@ -0,0 +1 @@
+cytools.py

versioned_hdf5/meson.build

@@ -33,7 +33,6 @@ py.extension_module(
     subdir: 'versioned_hdf5',
     dependencies: compiled_deps,
     include_directories: [npy_include_path],
-    cython_args: ['--cplus'],
     override_options : ['cython_language=cpp'],
 )
 
@@ -61,11 +60,21 @@ cython_args = [
 # meson build && pushd build && meson compile && popd && firefox $(find build/ -name "*.html")
 
 py.extension_module(
-    'subchunk_map',
+    'cytools',
     # FIXME This is a symlink. Can't find a way to convince Meson to compile
     # with Cython pure-python .py files
     # (read: https://cython.readthedocs.io/en/latest/src/tutorial/pure.html)
+    'cytools.pyx',
+    install: true,
+    subdir: 'versioned_hdf5',
+    dependencies: compiled_deps,
+    cython_args: cython_args,
+)
+
+py.extension_module(
+    'subchunk_map',
     'subchunk_map.pyx',
     install: true,
     subdir: 'versioned_hdf5',
+    dependencies: compiled_deps,
 )

versioned_hdf5/slicetools.pyx

@@ -1,3 +1,4 @@
+# distutils: language=c++
 import sys
 from functools import lru_cache
 
@@ -18,6 +19,9 @@ from libc.stddef cimport ptrdiff_t, size_t
 from libc.stdio cimport FILE, fclose
 from libcpp.vector cimport vector
 
+from versioned_hdf5.cytools import np_hsize_t
+from versioned_hdf5.cytools cimport count2stop, hsize_t, stop2count
+
 
 cdef FILE* fmemopen(void* buf, size_t size, const char* mode):
     raise NotImplementedError("fmemopen is not available on Windows")
@@ -32,10 +36,6 @@ cdef extern from "hdf5.h":
     ctypedef int hbool_t
     ctypedef int herr_t
     ctypedef int htri_t
-    ctypedef long long hsize_t  # as per C99, uint64 or wider
-    ctypedef signed long long hssize_t  # as per C99, int64 or wider
-    ctypedef signed long long haddr_t
-    ctypedef long int off_t
 
     cdef herr_t H5Dread(
         hid_t dset_id,
@@ -105,11 +105,6 @@ cdef extern from "hdf5.h":
 
 np.import_array()
 
-# Numpy equivalents of Cython types
-np_hsize_t = np.ulonglong
-np_hssize_t = np.longlong
-np_haddr_t = np.longlong
-
 NP_GE_200 = np.lib.NumpyVersion(np.__version__) >= "2.0.0"
 
 
@@ -495,33 +490,6 @@ cdef np.ndarray _preproc_many_slices_idx(obj: ArrayLike, hsize_t ndim, bint fast
         return arr
 
 
-@cython.cdivision(True)
-cpdef hsize_t stop2count(hsize_t start, hsize_t stop, hsize_t step) noexcept nogil:
-    """Given a start:stop:step slice or range, return the number of elements yielded.
-
-    This is functionally identical to::
-
-        len(range(start, stop, step))
-
-    Doesn't assume that stop >= start. Assumes that step >= 1.
-    """
-    # Note that hsize_t is unsigned so stop - start could underflow.
-    if stop <= start:
-        return 0
-    return (stop - start - 1) // step + 1
-
-
-cpdef hsize_t count2stop(hsize_t start, hsize_t count, hsize_t step) noexcept nogil:
-    """Inverse of stop2count.
-
-    When count == 0 or when step>1, multiple stops can yield the same count.
-    This function returns the smallest stop >= start.
-    """
-    if count == 0:
-        return start
-    return start + (count - 1) * step + 1
-
-
 @cython.boundscheck(False)
 @cython.wraparound(False)
 @cython.infer_types(True)

versioned_hdf5/subchunk_map.py

@@ -1,7 +1,8 @@
 from __future__ import annotations
 
 import itertools
-from typing import TYPE_CHECKING, Any, Iterable, Iterator
+from collections.abc import Iterable, Iterator
+from typing import TYPE_CHECKING, Any
 
 import cython
 import numpy as np
@@ -17,6 +18,8 @@
 )
 from numpy.typing import NDArray
 
+from .cytools import ceil_a_over_b, smallest_step_after
+
 if TYPE_CHECKING:
     # TODO import from typing and remove quotes (requires Python 3.10)
     # TODO use type <name> = ... (requires Python 3.12)
@@ -26,23 +29,11 @@
         "slice | NDArray[np.intp] | NDArray[np.bool] | tuple[()] | Py_ssize_t"
     )
 
-
-@cython.cfunc
-def _ceiling(a: Py_ssize_t, b: Py_ssize_t) -> Py_ssize_t:
-    """Returns ceil(a/b)"""
-    return -(-a // b)
-
-
-@cython.cfunc
-def _smallest(x: Py_ssize_t, a: Py_ssize_t, m: Py_ssize_t) -> Py_ssize_t:
-    """Find the smallest integer y >= x where y = a + k*m for whole k's
-    Assumes 0 <= a <= x and m >= 1.
-
-    a                  x    y
-    | <-- m --> | <-- m --> |
-    """
-    n: Py_ssize_t = _ceiling(x - a, m)
-    return a + n * m
+if cython.compiled:  # pragma: nocover
+    from cython.cimports.versioned_hdf5.cytools import (  # type: ignore
+        ceil_a_over_b,
+        smallest_step_after,
+    )
 
 
 @cython.cfunc
@@ -69,12 +60,12 @@ def _subindex_chunk_slice(
     """
     start: Py_ssize_t = max(c_start, s_start)
     # Get the smallest lcm multiple of common that is >= start
-    start = _smallest(start, s_start % s_step, s_step)
+    start = smallest_step_after(start, s_start % s_step, s_step)
     # Finally, we need to shift start so that it is relative to index
     start = (start - s_start) // s_step
 
     stop: Py_ssize_t = min(c_stop, s_stop)
-    stop = _ceiling(stop - s_start, s_step) if stop > s_start else 0
+    stop = ceil_a_over_b(stop - s_start, s_step) if stop > s_start else 0
 
     return slice(int(start), int(stop), 1)
 
@@ -96,7 +87,7 @@ def _subindex_slice_chunk(
     """
     start: Py_ssize_t = max(s_start, c_start)
     # Get the smallest step multiple of common that is >= start
-    start = _smallest(start, s_start % s_step, s_step)
+    start = smallest_step_after(start, s_start % s_step, s_step)
     # Finally, we need to shift start so that it is relative to index
     start -= c_start
 

versioned_hdf5/tests/test_cytools.py

@@ -0,0 +1,43 @@
+import math
+
+from hypothesis import given
+from hypothesis import strategies as st
+
+from ..cytools import ceil_a_over_b, count2stop, smallest_step_after, stop2count
+
+
+def free_slices_st(size: int):
+    """Hypothesis draw of a slice object to slice an array of up to size elements"""
+    start_st = st.integers(0, size)
+    stop_st = st.integers(0, size)
+    # only non-negative steps are allowed
+    step_st = st.integers(1, size)
+    return st.builds(slice, start_st, stop_st, step_st)
+
+
+@given(free_slices_st(5))
+def test_stop2count_count2stop(s):
+    count = stop2count(s.start, s.stop, s.step)
+    assert count == len(range(s.start, s.stop, s.step))
+
+    stop = count2stop(s.start, count, s.step)
+    # When count == 0 or when step>1, multiple stops yield the same count,
+    # so stop won't necessarily be equal to s.stop
+    assert count == len(range(s.start, stop, s.step))
+
+
+@given(st.integers(0, 10), st.integers(1, 10))
+def test_ceil_a_over_b(a, b):
+    expect = math.ceil(a / b)
+    actual = ceil_a_over_b(a, b)
+    assert actual == expect
+
+
+@given(st.data(), st.integers(0, 10), st.integers(1, 4))
+def test_smallest_step_after(data, x, m):
+    a = data.draw(st.integers(0, x))
+    expect = a
+    while expect < x:
+        expect += m
+    actual = smallest_step_after(x, a, m)
+    assert actual == expect

versioned_hdf5/tests/test_slicetools.py

@@ -8,7 +8,8 @@
 from hypothesis import strategies as st
 from numpy.testing import assert_equal
 
-from ..slicetools import count2stop, read_many_slices, spaceid_to_slice, stop2count
+from ..cytools import count2stop
+from ..slicetools import read_many_slices, spaceid_to_slice
 
 max_examples = 10_000
 
@@ -46,26 +47,6 @@ def test_spaceid_to_slice(h5file):
                     assert_equal(a[s.raw], a[sel], f"{(start, count, stride, block)}")
 
 
-def free_slices_st(size: int) -> Any:
-    """Hypothesis draw of a slice object to slice an array of up to size elements"""
-    start_st = st.integers(0, size)
-    stop_st = st.integers(0, size)
-    # only non-negative steps are allowed
-    step_st = st.integers(1, size)
-    return st.builds(slice, start_st, stop_st, step_st)
-
-
-@given(free_slices_st(5))
-def test_stop2count_count2stop(s):
-    count = stop2count(s.start, s.stop, s.step)
-    assert count == len(range(s.start, s.stop, s.step))
-
-    stop = count2stop(s.start, count, s.step)
-    # When count == 0 or when step>1, multiple stops yield the same count,
-    # so stop won't necessarily be equal to s.stop
-    assert count == len(range(s.start, stop, s.step))
-
-
 @st.composite
 def bound_slices_st(draw, arr_size: int, view_size: int) -> slice:
     """Hypothesis draw of a slice object to slice an array of <arr_size> points