Repository navigation
Faster NDArray reads, writes, fills, serialization and expressions; fixes for crashes and wrong results - #744
Closed
FabianIsensee wants to merge 19 commits into
Closed
FabianIsensee wants to merge 19 commits into
FabianIsensee wants to merge 19 commits into
Conversation
NDArray.__setitem__ calls mark_indexes_stale() after every write. Loading the store resolves the urlpath and stats every cached persistent path, and on NFS each stat() of a file with buffered writes first pushes them to the server. Return early when no cached store has indexes and the array's vlmeta has no index store. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
NDArray.get_slice_numpy() decompressed with the GIL held, so threads reading separate arrays, or separate handles of one file, ran one at a time; on NFS this included the wait for the server. Release the GIL around b2nd_get_slice_cbuffer(). The GIL was also what kept two threads from using one super-chunk at the same time. Replace the per-NDArray read_lock with one reentrant lock per super-chunk, shared by the SChunk wrapper, the NDArray, its views and vlmeta, and hold it in every entry point into C. Keep the GIL when a Python postfilter, codec or filter could run during the read, and on wasm32. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01SwMZj7Wi7Sc8HtQpELyNXB
Hold the operands' handle locks while the miniexpr prefilter of a lazy expression reads their super-chunks from C. Count the reads that run without the GIL, so that exit and the registration of Python codecs and filters wait for them, and free handle locks that other threads held at a fork in the child. Lock the metalayer counters too. Skip the new tests on wasm32, shorten them, retry the GIL check on busy machines, and add tests for lazy expressions and fork. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01SwMZj7Wi7Sc8HtQpELyNXB
Move _may_have_indexes next to the store helpers, iterate over a snapshot of the persistent store cache, and use plain helpers in the two tests. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
The miniexpr prefilter read the chunks of every in-memory operand from sc->data, and the fast matmul prefilter did so for every operand. Frames have no chunk table there: in-memory arrays created with contiguous=True, arrays from ndarray_from_cframe() (copy=False, the default) and, for matmul, any on-disk array. Any lazy expression or matmul with such an operand dereferenced a NULL table and crashed the interpreter. Read frame-backed in-memory operands through blosc2_schunk_get_chunk(), as on-disk ones already are, and keep frame-backed operands off the fast matmul path: in its Python checks, and in _set_pref_matmul() itself, which now raises ValueError for them before changing any state. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01SwMZj6jNaqsesMKnBde8Lr
The general indexing path (ndindex) raised, or returned a wrong shape, for keys that NumPy indexes fine: 554 of 18000 random keys in a comparison with NumPy. Fix the classes found: - slice bounds beyond the index range, and bool bounds or steps, next to ..., None or NumPy integers raised IndexError or TypeError; normalize them as NumPy reads them (clamped, bools as ints) before ndindex sees them; - 0-d integer arrays index as integers (AttributeError before); - scalar bools are 0-d boolean indices, as in NumPy: their dimension goes among adjacent advanced indices (integers included), else in front, and next to array indices they broadcast (True: no effect, False: an empty result or IndexError). The result shape comes from NumPy on a stride-0 stand-in for the array, so it costs no memory; - None next to (adjacent) array indices is supported; - 0-d arrays take NumPy scalar bools (AttributeError before) and raise IndexError for integer lists (ValueError before). Advanced indices separated by slices, ... or None still raise NotImplementedError (NumPy moves their dimensions to the front there). After this, none of the 18000 keys gives another result than NumPy; the remaining differences are that limitation and which error a key that is invalid for several reasons reports first. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01SwMZj6jNaqsesMKnBde8Lr
Keys made of integers (Python or NumPy, not bools), unit-step slices, one Ellipsis and, for reads, None went through the general path: 8-30 us of Python per call (ndindex as soon as a NumPy integer, ... or None appears), more than decompressing a small selection. NDArray.get_basic_numpy() classifies such keys in Cython, resolves them against the shape and reads the result under one acquisition of the handle lock, through _get_slice_cbuffer(), the GIL-releasing read that get_slice_numpy() now shares with it. get_basic_start_stop() resolves the same keys (without None) for __setitem__. Any other key, or one that does not fit the shape, gives None and takes the unchanged general path, which raises the same exceptions as before. NumPy integers without __index__ (timedelta64) are left to it too, where they raise IndexError as in NumPy. SWMR: a read resolves the key against the cached shape and relies on the refresh inside the C read; if b2nd_array_t.last_tick shows that this refresh changed the shape (or the key did not fit, or the read failed), it refreshes, resolves and reads again, so a key is never resolved against a stale shape, at one file check per on-disk read instead of two. A write refreshes under the handle lock before resolving. Keys with None now return arrays that own their data (they were views of a temporary). Measured (Ryzen 9 7950X): a[3, 5, 7] of a 300x512x512 float32 array 34 -> 25 us in memory and 48 -> 35 us on disk, a[500] of a 1000x1000 array 21.5 -> 13.9 us, a[np.int64(7), np.int32(9)] 37.5 -> 7.3 us, iterating over 2000 rows 64 -> 44 ms; a[np.int64(i), j] = v 22 us less. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01SwMZj6jNaqsesMKnBde8Lr
Creating a small array cost ~30 us, mostly Python bookkeeping: - compute_chunks_blocks() and _check_ndarray_kwargs() converted CParams and Storage with dataclasses.asdict(), which deep-copies; a shallow dict (vars() for Storage, whose fields are plain values) is enough; - every SChunk wrapper, so every NDArray, built its cparams and dparams eagerly under the handle lock; they are now built on first use; - create_cparams_from_kwargs() looked up Filter.BYTEDELTA.value (a slow enum descriptor) twice per filter slot; now once per call; - is_inside_new_expr() and is_inside_ne_evaluate(), which arange, linspace, eye, reshape, fromiter and lazy functions such as clip run on every call, used inspect.stack(), which reads the source of every frame; walking the frames gives the same answer in under 1 us. __setitem__ also wrapped plain NumPy values in a SimpleProxy, whose constructor runs compute_chunks_blocks(); other array types keep the wrapping, NumPy arrays give the same outcome without it (checked for contiguous, strided, Fortran-order, 0-d, broadcast and converted values). Measured (Ryzen 9 7950X): zeros((100, 100)) 31 -> 20 us, eye() 311 -> 179 us, reshape() 437 -> 276 us, clip(a, 0, 1) 84 -> 12 us; a 3x3 write into a 100x100 array ~12 us less (asdict() ~5 us, SimpleProxy ~7 us). Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01SwMZj6jNaqsesMKnBde8Lr
- update_special() built its special chunk in a temporary SChunk with default cparams, so on an NDArray's super-chunk the chunk header had another block size than the array (e.g. 20000 instead of 800 bytes) and every later read of the array failed: b2nd reads blocks selectively and checks their number per chunk. The new SChunk._special_chunk() makes zero, NaN, uninitialized and repeated-value chunks with the super-chunk's own cparams. - update_data() and insert_data() with copy=False did not free their compressed chunk when the C call failed (nor on a compression error), so it leaked wherever c-blosc2 did not free it itself. They now free it whenever the call fails. This needs the next c-blosc2 release (or the companion c-blosc2 PR) for its caller-keeps-the-chunk-on-failure rule; with an older c-blosc2, a failed frame update can free the chunk twice. - decompress_chunk() without dst fetched the whole compressed chunk only to read its size from the header, then read it again to decompress it; it now gets the lazy chunk (header and block offsets) for the size. astype(float64) of a 300x512x512 float32 array on disk: 385 -> 366 ms. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01SwMZj6jNaqsesMKnBde8Lr
a[region] = scalar built np.full() of the whole region, and a[region] = ndarray decompressed all of the source first: an uncompressed copy of the region in memory (a fill of a large on-disk array could exceed RAM), then one set_slice. Scalars: every chunk the region fully covers gets one shared special (repeated-value) chunk from SChunk._special_chunk(), and the rest is written in boxes of bounded size (32 MB, or one chunk) from one buffer sized to the largest box. Edge chunks get a special chunk only for zero, as writing them leaves zeros in their padding (which a resize brings to light). Chunks go last to first: a frame compacts its payload behind an updated chunk, and that tail then holds the tiny special chunks already written. Lossy codecs and filters, and dictionaries, keep compressing the values. NDArray values: read box by box along the destination's chunk boundaries into reused buffers and written with one set_slice per box; when both arrays encode chunks the same way (dtype, chunks, blocks, lossless cparams) and the region starts at a chunk boundary, covered chunks without padding are copied compressed. A value that may share storage with the destination (same super-chunk or same file) is still read in full first. Both hold the handle locks (HandleLocks with the source) for the whole write, so that other threads see it as atomic as one set_slice. Regions smaller than a chunk, and arrays whose frame takes the cross-process lock (locking=True or BLOSC_LOCKING, destination or source), keep the single set_slice, which holds that lock from the shape refresh to the last chunk. A shared flat chunk index helper now also serves detect_aligned_chunks(), and __setitem__ marks indexes stale through one helper. Measured (300x512x512 float32, default chunks, Ryzen 9 7950X): a[:] = 0.0 41 -> 0.2 ms with no extra memory (was 303 MB); the first a[:] = 1 into an existing on-disk array 588 -> 29 ms; a[:] = b with the same encoding 142 -> 18 ms and 306 -> 28 MB peak. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01SwMZj6jNaqsesMKnBde8Lr
asarray() of an NDArray, and so astype(), sliced every chunk with __getitem__ and cast it into fresh arrays: chunk-sized allocations per chunk, which glibc trims and page-faults back in each time. When source and destination have the same shape, chunks and blocks (the default for astype), chunks without padding are now decompressed into one reused buffer, cast item by item into another and compressed from there; edge chunks keep the per-chunk path (their padding may hold anything, and casting it could warn). Chunks encoded the same way in both arrays are copied compressed (e.g. asarray(a, copy=True, urlpath=...)). Measured (300x512x512 float32, Ryzen 9 7950X, this change toggled): astype(np.uint8) 110 -> 56 ms, with nnU-Net-like chunks 286 -> 267 ms; astype(np.float64) unchanged (compression-bound). Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01SwMZj6jNaqsesMKnBde8Lr
NDArray.to_cframe() and SChunk.to_cframe() got the frame from b2nd_to_cframe() / blosc2_schunk_to_buffer() in a C buffer and copied it into a new bytes object. They now call blosc2_schunk_to_buffer_alloc() with a callback that creates the bytes object of the final size, and the frame is written into it directly: no second copy of the whole frame. SChunk.to_cframe() still releases the GIL only under set_releasegil(True); the callback takes it back to create the bytes. Needs blosc2_schunk_to_buffer_alloc() from the next c-blosc2 release (or the companion c-blosc2 PR). Measured with it (300x512x512 float32 in memory, default chunks): to_cframe() 54 -> 15 ms and peak memory 456 -> 230 MB. Adds round-trip tests of to_cframe() / from_cframe() for every storage kind and for special, updated and resized chunks. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01SwMZj6jNaqsesMKnBde8Lr
- expr[:] and expr[()] with miniexpr ran the prefilter into a temporary uncompressed NDArray and read that back into NumPy. A clevel=0 chunk is a header followed by the raw prefilter output, so when the chunks tile the flat array (blocks C-contiguous in their chunk, chunks in the array, only the last one partial) each chunk is now produced in its final place in the NumPy result, last to first, its header landing in bytes not computed yet. - Full std()/var() of float32/float64 (no axis, no where, same dtype) are two single-pass reductions: the mean, then the mean of (x - m) * (x - m) with m inlined as a scalar, instead of a broadcast (1, ..., 1) operand that kept the second pass off the fast paths. Other cases, such as dtype=np.float64 for float32 data, keep the general path. std() and var() share one implementation. - Chunked reductions (axis=..., argmin, cumulative_sum, ...) and chunked evaluation decompress the chunks of operands whose blocks are C-contiguous runs of their chunk straight into reused chunk-sized buffers, compute flat chunk indices without NumPy and skip ndindex for whole arrays. The unused `aligned` plumbing is gone. - where(cond, x, y) with a scalar x or y is fused before the operands are specialized, so the scalar is inlined and miniexpr is used (it always fell back to numexpr). - No dtype inference for the bare-array LazyExpr that the NDArray reductions build, and one regex for all constructor checks. Measured (300x512x512 float32 in memory, 28 threads, Ryzen 9 7950X): (a + b)[:] 120 -> 62 ms, a.std() 68 -> 21 ms, (a.mean(), a.std()) 78 -> 31 ms, a.sum(axis=0) 28 -> 22 ms, where(a > 100, a, 0).compute() 184 -> 153 ms; 100x100: a.std() 0.61 -> 0.38 ms. std/var may differ from before in the last ulps (same accuracy as NumPy). Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01SwMZj6jNaqsesMKnBde8Lr
The bundled miniexpr (3f4db93) returns wrong results, or crashes, for several classes of expressions, all reproducible on main: reductions with a constant subexpression ((o * -2).sum() is 0.0 for an array of ones); two-argument float functions and true division of integer or bool operands computed in integer arithmetic (bool / bool died of SIGFPE); constants typed by the result dtype ((a / 3) > 2 on integers, int32 < 2.5); % and remainder() with C sign rules; mixed-dtype typing (garbage values, and int16 ((i + i) * abs(i))[:] corrupted the heap); bool + bool added as numbers; grammar differences (-x ** 2, ** associativity, chained comparisons). The new miniexpr_guard module reproduces miniexpr's typing of an expression and flags these cases, replacing the check for integer arctan2 only; fast_eval() and reduce_slices() then evaluate them with numexpr, which matches NumPy (1.5-4x slower for those expressions, e.g. (a * -2).sum() 10 ms but wrong -> 39 ms; others unchanged, the check is a cached lookup). It agreed with NumPy on 10.5k random expressions, but is not a proof. check_dtype() also gave integer or bool result dtypes to hypot, copysign and nextafter of integer subexpressions, to true division involving bools ((a < b) / b, bool / bool) and to the bool //, % and ** that NumPy computes in int8, so their results were truncated on every evaluation path. It now asks NumPy (the operator or ufunc on empty operands, with a scalar next to an array kept as it is, so that it promotes as a Python scalar); int32 / int32 keeps its float32. Array operands' dtypes are read directly, which pays for that probe. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01SwMZj6jNaqsesMKnBde8Lr
LazyExpr._collect_flat_indices_from_bool_ndarray(), used for 1-D masks, took every chunk as chunk_len items, so the padding of the last chunk could yield positions past the end of the array: blosc2.full(1000, True, chunks=(300,)) gave 1200 positions. Clip each chunk to the items inside the array. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01SwMZj6jNaqsesMKnBde8Lr
C-Blosc2 uses BLOSC_NTHREADS, when set to a number it accepts, for every compression and decompression context, whatever nthreads the parameters ask for. blosc2.nthreads, and so the CParams/DParams defaults and numexpr's thread count, ignored it: with BLOSC_NTHREADS=4, Python reported 28 threads and gave numexpr 28 while C used 4. Start blosc2.nthreads from what C-Blosc2 gives a compression context (the new blosc2_ext.context_nthreads(), via blosc2_ctx_get_cparams()), and say in set_nthreads()'s docs that the variable wins at the C level. Importing blosc2 also overrode NUMEXPR_NUM_THREADS, which set_nthreads() recommends setting; it now leaves numexpr alone then, and otherwise gives it at most numexpr.MAX_THREADS (BLOSC_NTHREADS can go beyond that, which numexpr refuses with an error). Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01SwMZj6jNaqsesMKnBde8Lr
List the faster indexing, writes, casts, serialization and expressions, the bug fixes, and the compatibility notes (the next C-Blosc2 release is required; None keys return arrays owning their data; std/var rounding) in the 4.14.2 section. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01SwMZj6jNaqsesMKnBde8Lr
Drop this commit before merging: once the c-blosc2 changes are merged and released, point BLOSC2_BUNDLED_VERSION back to Blosc/c-blosc2 at a release tag that contains them and bump BLOSC2_MIN_VERSION, since this branch uses blosc2_schunk_to_buffer_alloc() and the documented ownership of chunks after failed writes. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01SwMZj6jNaqsesMKnBde8Lr
This was referenced Oct 30, 2026
Member
|
Thanks for this contribution, appreciated. As this is really big, I think we really should split this. Our proposal (well, my agent and mine) would be to close this and open 6 new PRs: Sub-PR 1: Bug fixes & ASan/UBSan/LSan cleanups (Merge first)
Sub-PR 2: Skip byte shuffle for 1-byte types (Fast merge)
Sub-PR 3: Frame header-only reads, NFS handle adoption & slice locking
Sub-PR 4: One-pass frame copy &
|
Author
|
First of all thank you for your openness to PRs like this! Things move fast nowadays and it really is cool to see blosc2 moving fast as well :-) |
Member
|
Thanks for your understanding. Yeah, using Claude for sending the new PRs is fine. |
This was referenced Nov 3, 2026
Merged
Closed
This was referenced Nov 3, 2026
Open
Author
|
Closing in favour of the split, as proposed in #744 (comment):
#742 and #743 stay open as they are. The c-blosc2 side replaces Blosc/c-blosc2#833 with Blosc/c-blosc2#835 to Blosc/c-blosc2#840. As you suggested, two parts are not resubmitted until #741 has landed:
The miniexpr bug reports linked in the description above still apply. (The split and this comment were made with Claude Code on my behalf.) |
This was referenced Nov 3, 2026
FabianIsensee
added a commit
to FabianIsensee/python-blosc2
that referenced
this pull request
Nov 3, 2026
…ctions check_dtype() gave integer or bool result dtypes to hypot, copysign and nextafter of integer subexpressions, to true division involving bools ((a < b) / b, bool / bool) and to the bool //, % and ** that NumPy computes in int8, so their results were truncated on every evaluation path: hypot(a + 1, b) of int8 arrays gave [2 3 4 5] for NumPy's [2.83 3.61 4.47 5.39], (a < b) / b gave [0 0 0 0] for [0.5 0 0 0]. It now asks NumPy for these operators and functions: the operator or ufunc on empty operands, with a scalar next to an array kept as it is, so that it promotes as a Python scalar; int32 / int32 keeps its float32. Array operands' dtypes are read directly instead of through blosc2.result_type(), which more than pays for the probe: building (a + b) * 2 - f takes 62 instead of 70 us. The new test checks the dtype on every path and the values where the typing follows NumPy (compute() and partial slices). expr[:] and reductions of integer true division and two-argument float functions of two integer arrays still compute in integer arithmetic inside the bundled miniexpr, which is for miniexpr to fix. Taken from Blosc#744 (55660b4), without its miniexpr guard. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
FabianIsensee
added a commit
to FabianIsensee/python-blosc2
that referenced
this pull request
Nov 3, 2026
…rs in chunked paths - fast_eval()'s NumPy-output path wrote the miniexpr results of each chunk straight into the result only for 1-D arrays. It now does so for any array whose chunks are consecutive runs of the C-order result: blocks C-contiguous in their chunk, chunks C-contiguous in the array, only the last one partial (_chunks_tile_flat()), as with the default partitions. An N-D `out` must be C-contiguous for that; a 1-D one may still be strided. Other arrays keep the uncompressed temporary NDArray. - Chunked evaluation and reductions (axis=..., argmin, cumulative_sum, ...) decompress the chunks of operands whose blocks are C-contiguous runs of their chunk straight into reused chunk-sized buffers (_chunk_buffers()), instead of allocating per chunk; the unused `aligned` plumbing is gone. blocks_contiguous() moves to module level in ndarray.py for that. - One compiled regex for all constructor checks. Measured on a Ryzen 9 7950X (300x512x512 float32 in memory, 28 threads), on top of the std/var and where() change: (a + b)[:] 116 -> 71 ms, a.sum(axis=0) 29.3 -> 26.3 ms. Taken from the lazy-expression commit of Blosc#744 (b28fb0e), which was deferred until this branch lands; its own NumPy-output mechanism is replaced by extending the one of this branch. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
FabianIsensee
added a commit
to FabianIsensee/python-blosc2
that referenced
this pull request
Nov 3, 2026
…rs in chunked paths - fast_eval()'s NumPy-output path wrote the miniexpr results of each chunk straight into the result only for 1-D arrays. It now does so for any array whose chunks are consecutive runs of the C-order result: blocks C-contiguous in their chunk, chunks C-contiguous in the array, only the last one partial (_chunks_tile_flat()), as with the default partitions. An N-D `out` must be C-contiguous for that; a 1-D one may still be strided. Other arrays keep the uncompressed temporary NDArray. - Chunked evaluation and reductions (axis=..., argmin, cumulative_sum, ...) decompress the chunks of operands whose blocks are C-contiguous runs of their chunk straight into reused chunk-sized buffers (_chunk_buffers()), instead of allocating per chunk; the unused `aligned` plumbing is gone. blocks_contiguous() moves to module level in ndarray.py for that. - One compiled regex for all constructor checks. Measured on a Ryzen 9 7950X (300x512x512 float32 in memory, 28 threads), on top of the std/var and where() change: (a + b)[:] 116 -> 71 ms, a.sum(axis=0) 29.3 -> 26.3 ms. Taken from the lazy-expression commit of Blosc#744 (b28fb0e), which was deferred until this branch lands; its own NumPy-output mechanism is replaced by extending the one of this branch. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
FrancescAlted
pushed a commit
that referenced
this pull request
Nov 3, 2026
Two common cases fell off miniexpr, and without numexpr they now run chunk by chunk through evaluate_chunk(): - Full std()/var() of float32/float64 (no axis, no where, same dtype) reduced (x - mean) ** 2 with the mean as a broadcast (1, ..., 1) operand. They are now two single-pass reductions: the mean, then the mean of (x - m) * (x - m) with m inlined as a scalar. The sign of m is folded into the operator, as miniexpr reductions mis-evaluate a negative literal on the right of a binary operator. Other cases keep the general path; std() and var() share one implementation. - where(cond, x, y) with a scalar x or y: specialize_miniexpr_inputs() inlined the scalar and dropped its operand before the where() was fused in, so the fused expression referenced a missing operand, miniexpr failed to compile it and the evaluation fell back. The where() is now fused before specializing. - The bare-array LazyExpr that the NDArray reductions build takes the array's dtype instead of inferring it. Measured on a Ryzen 9 7950X (300x512x512 float32 in memory, 28 threads), against this branch at 7db7b09: a.std() 470 -> 53 ms, (a.mean(), a.std()) 498 -> 78 ms, where(a > 100, a, 0).compute() 690 -> 469 ms; 100x100 a.std() 0.73 -> 0.52 ms. std/var may differ from before in the last ulps. Taken from the lazy-expression commit of #744 (b28fb0e), which was deferred until this branch lands. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Hi Francesc,
this is the python-blosc2 side of Blosc/c-blosc2#833 (it needs that one:
blosc2_schunk_to_buffer_alloc()and the ownership rule for failed writes, so this PR temporarily builds against that branch; the last commit switches back to a release tag once it is merged).It contains #742 and #743: their original commits are included unchanged (merged at the bottom of this branch), because several changes here build on #743's per-handle locks. If you merge those two first, this diff shrinks accordingly; otherwise this PR supersedes them and you can close them when merging this one. Same approach: performance improvements with identical results, plus bugs we found on the way, each with a test that fails without the fix.
Performance (Ryzen 9 7950X, 28 threads; medians, interleaved against main + #742/#743)
a[3, 5, 7](memory / disk)blosc2.zeros((100, 100))a[:], 315 MB float32 (memory / disk)a[:] = 0.0a[:] = b(NDArray)to_cframe()(a + b)[:]a.std()/a.mean(), a.std()astype(np.uint8)On the Xeon node the bulk gains are larger (e.g.
a[:]1.8x,(a + b)[:]3x,std5.3x, writing a whole on-disk array 2x on NFS).What changed:
...,None) in__getitem__/__setitem__; less work when creating small arrays (lazy cparams/dparams, nodataclasses.asdict(), noinspect.stack()).__setitem__with a scalar fills fully covered chunks with special chunks instead of buildingnp.full(region); with an NDArray value it goes chunk by chunk and copies compressed chunks when both arrays encode them the same way. Arrays opened withlocking=Truekeep the single atomicset_slicepath.asarray/astypeof NDArrays cast whole chunks in reused buffers.expr[:]evaluates straight into the NumPy result (no temporary NDArray); fullstd/varare two single-pass reductions; chunked reductions read chunks into reused buffers.to_cframe()writes the frame straight into the returnedbytes.decompress_chunk()reads only the chunk header for the size.Bug fixes
contiguous=True, orndarray_from_cframe(cf)with the defaultcopy=False) and the fastmatmulpath with on-disk or frame-backed operands dereferenced a NULL chunk table.hypot/copysign/nextafterof integer subexpressions and true division of bools were inferred as integers and truncated on every path.update_specialon an NDArray's super-chunk wrote chunks with the wrong block size (every later read failed);insert_data(copy=False)leaked on failure; bool-mask positions ran past the array length.__getitem__'s general path raised or returned wrong shapes where NumPy returns a result (e.g. huge slice bounds, bool steps next to.../None/NumPy ints,(None, array), 0-d integer arrays, scalar bools next to arrays).BLOSC_NTHREADSdecides the real thread count in C butblosc2.nthreadsignored it;NUMEXPR_NUM_THREADSwas overridden at import.a[region] = scalarbuilt an uncompressednp.fullof the whole region (a fill of a large on-disk array could exceed RAM), anda[region] = other_ndarraydecompressed the whole source first; both now work chunk by chunk.update_data/insert_data(copy=False)calls freed their chunk again after c-blosc2 had (double free) or leaked it, depending on where the C call failed; fixed together with c-blosc2's ownership rule (needs Faster slice reads/writes, serialization and NFS opens; fixes for memory-safety bugs c-blosc2#833).miniexpr: bugs we could only route around
These are in the bundled miniexpr (pinned 3f4db93), reproduced on main;
src/blosc2/miniexpr_guard.pydetects affected expressions and sends them to numexpr (correct, 1.5-4x slower for those expressions only). Reports with repros, source locations and suggested fixes are in miniexpr_bug_reports.md:(o * -2).sum()-> 0.0 on a ones array/on integer/bool operands computed in integer arithmetic ((a / b)[:]truncates,bool / boolSIGFPE,hypot,logaddexp, ...); constants typed by the result dtype ((a / 3) > 2wrong for ints,int32 < 2.5compares with 2)%/remainderfollow C sign rules (also for floats)((i + i) * abs(i))[:]aborts in glibc)-x ** 2,**associativity and chained comparisons parse differentlyThe guard reproduces miniexpr's typing in Python; it agrees with NumPy on 10.5k random expressions but is not a proof: two shapes from that test still crash (sums of
where(...) * (b - 0)with uint8/int8,where(c > 3, a, b) - (b + a)with uint8/bool/float32), and narrow-integer reductions disagree with NumPy's wrap-around (uint8(a ** 2).sum()). A simpler, safe alternative would be using miniexpr only for same-dtype float operands; happy to switch if you prefer that until miniexpr is fixed.Found, not fixed (questions for you)
update_data& co. allocate a fresh chunk-sized output per chunk, and main's read path happened to warm glibc's heap for that size, which the faster read path no longer does. Planned follow-up: reuse the output buffer (python side, or with anupdate_buffercall in c-blosc2).NDArray.schunk.lockingreports False even when the array was opened withlocking=True.locking=True, a read can still fail withnchunk exceeds the number of chunkswhen another process shrinks the array: the slice is computed from a shape read under an earlier lock (c-blosc2 only exports the exclusiveblosc2_schunk_lock(); a shared variant would let us hold one lock across refresh, slice computation and read).pickle.dumps/copy.copy/copy.deepcopyraise TypeError for NDArray, SChunk and LazyExpr (nothing inside python-blosc2 needs them; user code with process pools has to pass paths orto_cframe()bytes).expr[...],expr[:, :]and partial slices are evaluated by numexpr,expr[:]andcompute()by miniexpr, so they differ in the last bits (transcendental functions, float32 arithmetic with float literals, int32 + float32 computed in float32 by numexpr). Routing them all through the direct miniexpr path would make(a + b)[...]2x faster; we left it because it changes those spellings' results.Behaviour changes
Nonereturn arrays that own their data (was a view of a temporary).std/varresults may differ from before in the last ulps (same accuracy as NumPy).arctan2cases with mixed int/float operands moved from numexpr to miniexpr (correct values, ±1 ulp;.sum()now float64, matchingexpr.dtype).blosc2.nthreadsreflectsBLOSC_NTHREADS.Tests
11071 passed (2 failures: the
b2viewhttp tests that need aiohttp, also on main); new tests for every fix; randomized differential tests against main (+#742/#743): getitem keys (~150k vs NumPy), setitem regions, fills, expressions, reductions, cframes (identical apart fromstdrounding of constant arrays).AI disclosure: I used Claude Code (Anthropic) for the investigation, the code, the tests and a first draft of this description. I reviewed and tested everything.
Best,
Fabian