Repository navigation
Introduce Menudet 1.0 draft for portable array and table computation - #741
FrancescAlted wants to merge 266 commits into
Conversation
|
Hi Francesc, as you suggested, I left the miniexpr guard and the lazy-expression speedups of #744 out of the split until this PR is in. To see what will still be needed after it, I built its current head (c29b074) and 1. Most of the miniexpr wrong results are still there. The guard's 418 test cases compare with NumPy: 270 fail on Partial slices such as import numpy as np, blosc2
o = blosc2.asarray(np.ones((37, 40), dtype=np.int64))
blosc2.lazyexpr("o * -2", {"o": o}).sum() # 0 NumPy: -2960
f = blosc2.asarray(np.ones((37, 40)))
blosc2.lazyexpr("f - -0.5", {"f": f}).max() # 1.0 NumPy: 1.5
a = blosc2.asarray(np.arange(1, 13, dtype=np.int32))
b = blosc2.asarray(np.array([5, 3, 7, 2] * 3, dtype=np.int32))
blosc2.lazyexpr("(a + b) / b", {"a": a, "b": b})[:][:4] # [1. 1. 1. 3.] NumPy: [1.2 1.667 1.429 3.]
blosc2.lazyexpr("hypot(a, b)", {"a": a, "b": b})[:][:4] # [5. 3. 7. 4.] NumPy: [5.099 3.606 7.616 4.472]
blosc2.lazyexpr("a < 2.5", {"a": a})[:][:3] # [T F F] NumPy: [T T F]
blosc2.lazyexpr("(a / 3) > 2", {"a": a})[:][5:9] # [F F F T] NumPy: [F T T T]
x = blosc2.asarray(np.array([-7.0, -5.0, 4.0, 5.0]))
blosc2.lazyexpr("x % 3", {"x": x})[:] # [-1. -2. 1. 2.] NumPy: [2. 1. 1. 2.] (also remainder() of ints)
y = blosc2.asarray(np.array([1.0, 2.0, 3.0]))
blosc2.lazyexpr("-y ** 2", {"y": y})[:] # [1. 4. 9.] NumPy: [-1. -4. -9.]
blosc2.lazyexpr("2 ** 3 ** 2 * y", {"y": y})[:] # [64. 128. 192.] NumPy: [512. 1024. 1536.]
i = blosc2.asarray((np.arange(1000, dtype=np.int16) % 300) - 150)
blosc2.lazyexpr("(i + i) * abs(i)", {"i": i})[:][:3] # different garbage on every run (main: "double free or corruption")
c = blosc2.asarray(np.array([True, False, True, False]))
blosc2.lazyexpr("c / c", {"c": c})[:] # SIGFPE, the interpreter dies (main too)All of these give the same results on Since this PR drops numexpr, the guard as written (send such expressions to numexpr) has nowhere to send them any more. I think these are better fixed in miniexpr than worked around in python-blosc2, so I am not resubmitting the guard. If you would rather have a Python-side guard on top of this PR (e.g. declining these signatures in the native admission), let me know. 2.
In 3. The lazy-expression speedups. On top of the split PRs (#746-#749), the deferred commit (b28fb0e) speeds up these cases:
None of this is in this PR. Its single-pass AI disclosure: the measurements and this comment were prepared with Claude Code; I reviewed them. Best, |
|
I could also create a PR in miniexpr to solve the bugs we found, once this is merged |
|
@FabianIsensee Ok, with the recent change, this should be very close to be ready. It would be nice if you can test this and suggest improvements; PRs to this PR are welcome too. |
|
On it. I presume you mean the speed-ups? What do we do with the miniexpr bugs? Shall I work on that as well (or rather let claude loose on it)? |
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>
|
Both miniexpr bugs and speed-ups are welcome in this PR. The new miniexpr in this PR is meant to mimic NumPy as close as possible, so any divergence from the NumPy artihmetic can be discussed (proposed PRs welcome) here too. |
Fusing where(cond, x, y) with scalar x/y (previous commit) exposed miniexpr to cases it gets wrong, as Francesc found in review: - where(a > 0.1, 1, 0) with a float32 `a` of 0.1 selected 1. miniexpr types the literals of the fused expression by its output dtype and computes in it: with an int64 output 0.1 became 0, with a float32 output it compared against the double 0.1, where NumPy rounds it to float32 first. Arithmetic in the condition also ran in the output dtype (where(a + b == c, 1, 0) with float32 operands). - where(a > 5, np.uint64(2**63 + 1), a) gave 2**63: miniexpr parses literals into C doubles. Inside where() it also passes 64-bit integer operands and branch values through doubles. where() is now fused only when computing the condition in the output dtype changes nothing: all its operands have the output dtype, or it only compares operands that cast safely to it against literals their dtypes hold exactly; never with 64-bit integers. Otherwise it falls back to the chunked evaluation. This also fixes array branches, which were fused before and had the same problems (where(a > 0.1, a, b)). Masks, label maps, thresholds and same-dtype arithmetic stay on miniexpr; where(a > 100, a, 0).compute() keeps its speed. Expressions with integer literals beyond 2**53, written or inlined scalars, are no longer sent to miniexpr (also a + np.uint64(2**63 + 1)). Tests: the review examples, a NumPy differential grid over condition dtypes, literals, comparisons, output dtypes and scalar/array branches, arithmetic conditions, large integer constants, and checks that the common mixed-dtype where() patterns still compile with miniexpr. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Two CI failures of the previous commits: - NumPy 1.26: a.std(ddof=1) and a.var(ddof=1) of float32 returned float64. The full std/var returns a float32 scalar, and under NumPy 1.x's value-based casting a float32 scalar times the Python int count of the n / (n - ddof) rescaling becomes float64 (with keepdims, a float32 array times a count above 65535 does too). The rescaled variance is now cast back to its own dtype; with NumPy 2 this changes nothing. This fixes the 24 test_std_var_full_reduction[...-1-...- float32] failures. - WebAssembly: test_large_integer_constants_exact and test_where_literal_condition_matches_numpy added 2**53 to a plain np.arange(), which is int32 on wasm32 and raised OverflowError. They build int64 inputs explicitly now. Not changed here: on NumPy 1.26, mean() of float32 already returns float64 for full reductions and for axis reductions over more than 65535 elements, and mean(where=) does on NumPy 2 too (also on 7db7b09). Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
…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 #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>
Keep std/var and where() with scalars on the miniexpr path
Write expr[:] into the NumPy result for N-D arrays; reuse chunk buffers in chunked paths
Supersedes #735.
Menudet is a little language for portable computation on arrays and tables, with Python-like syntax and native execution through miniexpr. Validated kernels can be saved and loaded without reconstructing Python functions, while preserving defined numeric semantics and logical computation groups. Menudet 1.0 remains a draft specification.
Main changes
@blosc2.menudetauthoring with annotated or inferred dtypes, and native execution without Python reconstruction or execution fallback.deserialize="full"for legacy Python-specific recipes and reject unsupported saves before destination writes.Performance and validation
For
x + (y - 1)over 100K float64 NumPy operands on macOS arm64, with 12 configured threads and fresh NumPy outputs, preparation-cold medians from 151 samples were:Preparation/compiler caches were cleared before each cold sample; inputs were prebuilt and CPU/OS caches were not flushed. Measurements used a rebuilt native extension containing the pinned uniform-cast optimization. Compressed-output LazyExpr is excluded from this comparison.