[FabianIsensee]

Hi Francesc,

a third PR into #741, independent of #752 and #753: it is based on improved-miniexpr (7db7b09) and merges cleanly with both, in any order. It is the Python-side part of the miniexpr guard of #744, without the guard.

The bug: check_dtype() gave integer or bool result dtypes to hypot, copysign and nextafter of integer subexpressions, to true division involving bools, and to the bool //, % and ** that NumPy computes in int8. Their results were truncated on every evaluation path, e.g. for int8 arrays:

blosc2.hypot(a + 1, b)[:]   # [2 3 4 5] (int8)      NumPy: [2.83 3.61 4.47 5.39] (float16)
((a < b) / b)[:]            # [0 0 0 0] (int8)      NumPy: [0.5 0. 0. 0.] (float64)

The fix: for these operators and functions, check_dtype() 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 instead of through blosc2.result_type(), which more than pays for the probe: building (a + b) * 2 - f takes 62 instead of 70 µs.

Tests:

• New: test_lazyexpr_result_dtype_matches_numpy (17 expressions x bool/int8/uint16/int32/int64, built from functions and from strings). It checks the dtype on every path, and the values on compute() and on partial slices. 64 of its cases fail without this change.
• pytest tests: 10757 passed. The 2 failures are the b2view http tests that need aiohttp, as on main.
• ruff passes.

Still wrong, but not fixable here: expr[:] and the reductions of integer true division and of two-argument float functions of two integer arrays still compute in integer arithmetic inside the bundled miniexpr: (a < b) / b now has the right dtype, but [:] still gives [0. 0. 0. 0.]. compute() gets them right, because it goes through the native graph. That is for miniexpr to fix; I'd like to send PRs there once this branch is in, as offered above.

AI disclosure: the code, the tests and the measurements were done with Claude Code; I reviewed them.

Best,
Fabian

🤖 Generated with Claude Code

[FrancescAlted]

Thanks for this fix. I reviewed PR head aa351ae against its base (7db7b09) and found one regression in the _operand_dtype() shortcut. Could you take a look?

Non-native-endian operands now fail existing dtype checks

import blosc2
import numpy as np

a = blosc2.asarray(np.arange(1, 7, dtype=np.int32))
b = np.ones(6, dtype=">i4")
print((a & b)[:])

On the base, this returns [1, 0, 1, 0, 1, 0] through both [:] and .compute()[:]. On this PR, constructing a & b raises:

ValueError: Invalid operand type for &: >i4

The previous blosc2.result_type(value) normalized byte order, whereas convert_dtype(value.dtype) preserves it. As a result, the existing comparisons against np.int32/np.int64 in check_dtype() reject non-native-endian integer operands. The other bitwise operators are affected too.

The same distinction also bypasses the intentional int32-division special case: for a = blosc2.asarray(np.arange(1, 7, dtype=">i4")), (a / a).dtype changes from float32 on the base to float64 on this PR.

This concerns src/blosc2/lazyexpr.py:3644–3647. Normalizing byte order for these dtype checks (for example, with dtype.newbyteorder("="), while retaining the actual operand storage unchanged) should preserve the old behavior. A regression test with a native NDArray and a non-native-endian NumPy operand would cover the working example above.

Validation: using the blosc2 conda environment and an isolated PR module with the installed native extension, 1,487 tests passed, 10 skipped, across the lazy-expression, reduction, and NumPy-prefilter test files under the default selection; 900 additional native-endian dtype probes passed. The reproducer above was checked against both PR and base. I did not rebuild native dependencies.

Posted by an AI coding assistant at the maintainer’s request.

[FabianIsensee]

Thanks, Francesc, you are right: _operand_dtype() kept the operand's byte order, and the dtype checks in check_dtype() compare against native dtypes. Fixed in f1fa158: it now returns np.result_type(dtype), which is what blosc2.result_type() computed, without building the array. Your example gives [1 0 1 0 1 0] again, and a / a of a >i4 array is float32 again. test_lazyexpr_non_native_byte_order covers both, with NumPy and NDArray operands.

Looking for more of the same, I found two other changes against the base and removed them:

• arctan2 of complex, datetime or string operands raised NumPy's TypeError when the expression was built; with the probe it was only raised on evaluation. It is raised when built again (test_arctan2_rejects_complex_when_built).
• Structured and void operands go through blosc2.result_type() again, which rejects some of them.

To check for others, I compared the base (7db7b09) and this PR on about 116k expressions: all operators and two-argument functions, and the unary functions, on bool, all int, uint, float and complex dtypes in both byte orders, plus datetime, timedelta, string, void and structured; NDArray, NumPy, LazyExpr and scalar operands; built from functions and from strings; checking the dtype, [:], compute(), a partial slice and sum() against NumPy, with NumPy 2.5 and 1.26.
Every difference left is the intended one, except four places where the corrected dtype runs into a problem that is already there:

• (td / i).sum() of a timedelta64 and an integer array: the base typed td / i as float64 (NumPy counts timedelta64 as an integer type), and its sum() happened to work. Now td / i is timedelta64 as in NumPy, so [:] gives timedeltas, but sum() raises the DTypePromotionError that (td * 2).sum() already raises on main.
• blosc2.pow(bool_array, 2) is int8, like bool_array ** 2 in NumPy, not int64 like np.power(bool_array, 2), because blosc2.pow is **.
• In an expression with a non-bool result, bool + bool is computed as an integer addition, not as NumPy's or: (a + a) + 1 gives 3 where NumPy gives 2, on main as well. (a + a) ** b and hypot(a + a, 300) now have NumPy's (non-bool) dtype, so they go the same way.
• With NumPy 1.26, uint8_array ** 300 is int16, as NumPy 1.26 gives it, and its sum() hits the reductions of overflowing integer powers, which do not wrap like [:] does: on the base, (a + a) ** 9 of a uint8 array sums to 80906496000, where NumPy gives 0.

pytest tests: 10794 passed with NumPy 2.5, 10640 with NumPy 1.26 and no numexpr; the 2 b2view http tests that need aiohttp fail, as on the base. ruff passes.

The comparison also turned up some problems that are on main already, unrelated to this PR. I can open issues for them if you like:

• Expressions of big-endian NDArray operands read the data with the native byte order on the full-array paths: (a + b)[:] with a >i4 NDArray b gives 16777217 for 2 ([:3] is right).
• The bool / bool SIGFPE from my comment on Introduce Menudet 0.1 draft for portable array and table computation #741 also hits integer / bool arrays with a False in them, depending on what ran before; on main, some ** expressions die the same way.
• Comparisons with a negative Python int, e.g. x > -3, often give wrong elements or sums, and the same expression differs from run to run (206 of the 378 such cases on main).
• Reductions of timedelta64 expressions raise DTypePromotionError ((td * 2).sum()).

AI disclosure: the fix, the tests and the comparison were done with Claude Code; I reviewed them.

Best,
Fabian

[FabianIsensee]

One more thing: the WASM job on my fork showed that the new test_lazyexpr_result_dtype_matches_numpy had failed there since the first commit (32 cases). I had not run it under WebAssembly before.

30 of them are miniexpr's: on WebAssembly, compute() and slices also run the bundled miniexpr, which computes integer true division and hypot/nextafter of integer arrays in integer arithmetic (0.0 for NumPy's 0.5); natively that only happens on expr[:]. The other 2 are int32_array / 2: a Python int is int32 on wasm32, so it gets the float32 of int32 / int32 there, while it is float64 on 64-bit platforms. That is the same on the base; tell me if you want the int32 rule to apply only to two arrays, then it would be float64 everywhere.

I did not want to weaken the test, since these should pass once miniexpr is fixed, so b62727a marks exactly these 32 cases as xfail(strict=True) on WebAssembly (_WASM_KNOWN_FAILURES); the test itself is unchanged (014a7b1, which only dropped the value checks on WebAssembly, is undone by it). With strict, a fix in miniexpr shows up as XPASS, and the entries can be removed. The library code is the same as in f1fa158. CI on my fork passes on all platforms (except mac which queues and then gets cancelled, probably due to a lack of runners), WASM included: Tests, Tests (WASM) (32 xfailed).

AI disclosure: this fix was also done with Claude Code; I reviewed it.

[FrancescAlted]

Thanks for the follow-up fixes and WASM investigation. The strict xfails are useful tracking, but I suggest separating the two kinds of WASM failures before merging:

1. The two int32_array / 2 cases: please replace these xfails with ordinary tests asserting the existing platform-dependent dtype policy. On wasm32, the Python integer triggers our intentional int32/int32 → float32 exception, so this is not itself a native numerical bug. The test should still check values. Any change to restrict that exception to array operands and make scalar division consistent across platforms can be considered separately.

2. The other 30 numerical cases: retain narrowly scoped xfail(strict=True) entries, with specific reasons and links to tracked issues. These represent genuinely incorrect numerical results, not an accepted platform policy.

3. Upstream tracking: please open miniexpr issues with minimal reproducers for integer true division and integer-input floating functions, including the WASM interpreter path. Native fixes should cover operation/intermediate typing, not just the final output dtype.

4. Python-Blosc2 follow-up: please track safe NumPy fallback guards until corrected native dependencies are bundled. Guards should cover [:], .compute(), partial slices, and reductions, and must actually reach a correct evaluator. Simply disabling the prefilter may not suffice if chunk evaluation selects another affected native path. Blanket conversion of all inputs to float would also be unsafe for integer overflow, comparisons, and large-integer precision.

5. Please rerun native and WASM CI after adjusting the tests.

Only items 1, 2, and 5 require code/test changes in this PR. The numerical fixes can follow separately, with tracking issues in place before merging; there is no need to expand this PR into a complete native typing fix. Does that split sound reasonable?

I inspected the successful fork WASM CI run and its 32 expected failures; I have not independently run the WASM build locally.

Posted by an AI coding assistant at the maintainer’s request.

[FabianIsensee]

Thanks, Francesc, the split sounds right. Done:

1. The two int32_array / 2 cases are ordinary tests again: where a Python int is int32, as on wasm32, the test expects the float32 of our int32 / int32 rule, elsewhere NumPy's float64, and it checks the values in both cases (bf840e5).
2. The 30 miniexpr cases stay xfail(strict=True), each with its reason and issue: True division of integer operands is computed in integer arithmetic miniexpr#26 (integer true division) and hypot(), atan2(), logaddexp() and nextafter() of integer operands are computed in integer arithmetic miniexpr#27 (hypot(), nextafter() and other float functions of integers).
3. The miniexpr issues have C reproducers. These also fail with SIMD and the JIT disabled, and they include nested expressions such as a / b + 0.25 and hypot(a, b) > 2, where the intermediate nodes need the float type too.
4. The guard is tracked in Keep integer true division, hypot() and nextafter() away from the bundled miniexpr until it is fixed #756, with your requirements.
5. CI on my fork: Tests, Tests (WASM) (30 xfailed).

Fixes for both miniexpr issues are prepared; I can send them as miniexpr PRs whenever you want them.

AI disclosure: the test change and the issues were done with Claude Code; I reviewed them.

[FrancescAlted]

Hi Fabian, thanks for the fixes and the careful regression coverage here. Since this PR was opened, the changes on improved-miniexpr have substantially superseded its main correctness fix:

• fb00ab60 introduced the typed expression evaluator, including NumPy-based operation-result dtype probing in expression_evaluator.binary_result_dtype().
• 71ce1504 made that evaluator the default. check_dtype() now uses that probing for ordinary typed expressions.
• 6dd0dbb9 also preserves zero-dimensional NumPy operands during promotion, which should be retained in any follow-up.

I checked the added tests from this PR against the current branch (fe932dc6), without applying the PR changes: 197 passed, 15 skipped, no failures, using the blosc2 conda environment. This was native validation, not an actual WASM run. The current typed-evaluator tests also cover representative division, hypot, and nextafter cases across full reads, compute(), partial slices, and reductions.

Could you please close this PR as superseded and adapt any still-useful parts into a follow-up against the new typed evaluator, rather than just resolving the merge conflicts? In particular, selected additional regression tests and the operand-dtype performance optimization may still be useful. Please preserve strong/zero-dimensional scalar promotion, avoid duplicating the new dtype-probing machinery, and reassess the old WASM xfails against current behavior instead of carrying them over unchanged. Any remaining legacy-evaluator fixes can be considered separately.

Thanks again for the investigation and tests!

Posted by an AI coding assistant at the maintainer’s request.

[FabianIsensee]

Hey @FrancescAlted,
thanks for the detailed response. I can definitely do that.
Note that I have a bunch of PRs ready for miniexpr. It seems like you are actively working on that at the moment. Do you want them ASAP (it's again speed improvements and bug fixes) or do you prefer having them once your current work is done? They are all ready and only need to be opened. Nicely separated so that's easier to review.
Best,
Fabian

[FrancescAlted]

Hi Fabian. Yeah, miniexpr/blosc2 integration has been a moving target for the past month or so, and I hope to settle down specs very soon. Your work is important to reach stability, so I think that sending feedback now (via PRs or issues) would be a good thing right now. Go ahead and see how it goes.

[FabianIsensee]

A lot of my fixes touched me_compile while you seem to be moving away from that. Is me_compile still relevant or will it be superseded by the Menudet graph engine?

[FrancescAlted]

Yes, you’ve understood the direction correctly. We’re moving ordinary expression evaluation toward the typed Menudet graph engine, rather than the legacy me_compile path.

The goal is portable execution with explicitly defined, NumPy-style numerical semantics: intermediate operations retain their own types instead of inheriting the final output dtype. The same validated graphs can execute natively through interpretation or optional JIT, including standalone and WASM deployments.

me_compile remains relevant for existing paths, but the typed graph engine is the direction for new correctness work. Your reproducers and regression tests remain valuable; ideally, native numerical fixes should target the shared typed semantics and be checked through graph execution.

The best architectural reference is Native numerical graph preparation. For the numerical contract, see NumPy arithmetic semantics and NumPy function semantics.

This is still evolving and not a claim of complete NumPy compatibility (which could be difficcult to achieve in its entirety). Those documents rather describe the supported subset and remaining qualification boundaries.

[FabianIsensee]

Thanks Francesc, that makes the direction clear.

While checking my PRs against the new default, I found that the typed evaluator is much slower than the legacy one for several common expressions (times in ms):

Expression	Legacy evaluator	Legacy evaluator + planned miniexpr PRs	Typed evaluator (new default)
float32 sin(a)[:]	54.4	53.2	44.0
float32 (a > 100).sum()	32.0	24.1	15.6
float32 (a + b)[:]	36.3	36.1	42.3
float64 (a * 2 + 1)[:]	41.5	41.9	56.0
float32 where(a > 100, a, 0)[:]	38.1	26.0	64.0
int32 ((a + b) * abs(a))[:]	42.2	38.9	78.5
bool p, float32 a, b: where(p, a, b)[:]	41.6	41.5	171
float32 where((a > 90) & (a < 110), a, 0)[:]	71.7	43.4	639
float32 a, float64 b: ((a > 100) & (b < 100))[:]	38.1	36.5	476
float32 a, float64 b: ((a > 100) & (b < 100)).sum()	62.3	61.8	481
int32 (a * 3 + 1)[:]	29.6	28.4	431
int32 (a % b)[:]	54.3	53.6	2044
bool (p | q)[:]	17.5	17.0	1851
bool (p & q).sum()	14.3	14.3	1726
int16 a, uint8 b: (a + b)[:]	20.5	19.3	2822
float32 a.std() (300×512×512, default chunks)	56.6	56.1	51.6

Legacy evaluator: yesterday's b49bfc42, where it was the default. Typed evaluator: today's fe932dc6. Both use the bundled miniexpr 687e02b. The middle column is the legacy evaluator with the seven miniexpr PRs listed at the end merged onto 687e02b; they leave the typed evaluator's timings unchanged (within a few percent). Operands: 200×512×512 (52M items) in memory, chunks (10, 512, 512), blocks (1, 128, 512), unless noted. Medians of 3 alternating rounds of 5–7 runs each. AMD Ryzen 9 7950X, Linux, 28 blosc2 threads, NumPy 2.5.4.

What I found so far:

• i16 + u8, integer % and bool |/& run through evaluate_chunk()'s native chunk schedule (BLOSC_ME_JIT_TRACE=1 prints engine=evaluate_chunk), at about 53 ns per element; NumPy needs 0.2 ns for i16 + u8 in memory. Smaller chunks and BLOSC_ME_JIT=0 change nothing, so it is not compile time. compute() is as slow as [:] for these.
• a * 3 + 1 and the masked where() stay on the typed prefilter but are 9–15x slower for [:]; their compute() takes as long as with the legacy evaluator.
• Faster: sin() and reductions of comparisons.
• With BLOSC2_EXPRESSION_EVALUATOR=legacy on fe932dc6, reductions of expressions with comparisons now go to NumPy: (a > 100).sum() takes 96 ms instead of about 30 ms, the mask sum 5.8 s instead of about 60 ms.

Reproducer (5.2M items, a few seconds)

import os, time
import numpy as np
import blosc2

def best(f, n=5):
    f()
    ts = []
    for _ in range(n):
        t = time.perf_counter(); f(); ts.append(time.perf_counter() - t)
    return min(ts) * 1e3

rng = np.random.default_rng(0)
shape = (20, 512, 512)  # 5.2M items
kw = {"chunks": (10, 512, 512), "blocks": (1, 128, 512)}
f32 = blosc2.asarray(rng.normal(100, 30, size=shape).astype(np.float32), **kw)
i32a = blosc2.asarray(rng.integers(-1000, 1000, size=shape, dtype=np.int32), **kw)
i32b = blosc2.asarray(rng.integers(1, 1000, size=shape, dtype=np.int32), **kw)
i16 = blosc2.asarray(rng.integers(-1000, 1000, size=shape, dtype=np.int16), **kw)
u8 = blosc2.asarray(rng.integers(0, 255, size=shape, dtype=np.uint8), **kw)
p = blosc2.asarray(rng.random(shape) > 0.5, **kw)
q = blosc2.asarray(rng.random(shape) > 0.3, **kw)
cases = {
    "i16 + u8": ("a + b", {"a": i16, "b": u8}),
    "i32 a % b": ("a % b", {"a": i32a, "b": i32b}),
    "bool p | q": ("p | q", {"p": p, "q": q}),
    "i32 a * 3 + 1": ("a * 3 + 1", {"a": i32a}),
    "where((a > 90) & (a < 110), a, 0)": ("where((a > 90) & (a < 110), a, 0)", {"a": f32}),
}
print("BLOSC2_EXPRESSION_EVALUATOR =", os.environ.get("BLOSC2_EXPRESSION_EVALUATOR", "(default)"))
for name, (expr, ops) in cases.items():
    getitem = best(lambda: blosc2.lazyexpr(expr, ops)[:])
    compute = best(lambda: blosc2.lazyexpr(expr, ops).compute())
    print(f"{name:36s} [:] {getitem:7.1f} ms   compute() {compute:7.1f} ms")

On fe932dc6:

BLOSC2_EXPRESSION_EVALUATOR = legacy
i16 + u8                             [:]     2.1 ms   compute()     9.6 ms
i32 a % b                            [:]     5.5 ms   compute()    30.3 ms
bool p | q                           [:]     2.0 ms   compute()     9.0 ms
i32 a * 3 + 1                        [:]     2.5 ms   compute()    61.2 ms
where((a > 90) & (a < 110), a, 0)    [:]     6.8 ms   compute()    87.7 ms
BLOSC2_EXPRESSION_EVALUATOR = typed
i16 + u8                             [:]   279.2 ms   compute()   286.4 ms
i32 a % b                            [:]   207.2 ms   compute()   217.2 ms
bool p | q                           [:]   183.9 ms   compute()   190.9 ms
i32 a * 3 + 1                        [:]    40.2 ms   compute()    59.6 ms
where((a > 90) & (a < 110), a, 0)    [:]    59.4 ms   compute()    87.5 ms

On the miniexpr PRs: following your answer, I will open the ones that fix crashes, wrong values or speed in the existing paths without changing their semantics: heap corruption in and/or/not with numeric operands, wrong reduction results (constant subexpressions, comparisons of bools, bool min/max on padded blocks), exact 64-bit integer operations, complex powers (in Menudet kernels too), float32 nextafter(), faster arithmetic loops, and a test-only fix for the WASM stack overflow in test_dsl_artifact on improved-miniexpr. I am holding back the four that would bring NumPy semantics into me_compile (a typing pass for true division and literals, signed/unsigned promotion, floor modulo, Python operator precedence). If you want them for the existing paths, I can open them too. Otherwise I will check their reproducers against graph execution.

AI disclosure: the measurements and this comment were prepared with Claude Code; I reviewed them.

[FrancescAlted]

@FabianIsensee Yes, that kind of feedback and possible enhancements would be great, and it is the kind of help needed at this time. Looking forward to your PRs/issues to miniexpr.

[FabianIsensee]

Closing this as superseded, as you suggested. The useful parts are in #758:

• The NumPy probe this PR added to check_dtype() is what binary_result_dtype() in your typed evaluator does, so there is nothing to port.
• _operand_dtype() is not worth a separate change: blosc2.result_type() builds only a one-element array per operand, so it saved about 1 µs. On fe932dc, blosc2.result_type() of a 52M-element NDArray takes about 1 µs and check_dtype() about 7 µs.
• The tests are in Test NumPy's result dtypes on the typed evaluator (follow-up to #754) #758, adapted to the typed evaluator. They pass there, including the 30 cases that were strict xfails on WebAssembly (also with blosc2.IS_WASM = True set natively), so Test NumPy's result dtypes on the typed evaluator (follow-up to #754) #758 drops those xfails. With BLOSC2_EXPRESSION_EVALUATOR=legacy, the 80 cases the legacy evaluator still gets wrong are strict xfails.
• Test NumPy's result dtypes on the typed evaluator (follow-up to #754) #758 also adds a strict xfail for a bug in the typed evaluator: comparing an int8 array with 300 raises OverflowError, while NumPy 2 compares (all True).

AI disclosure: the checks and this comment were done with Claude Code; I reviewed them.

🤖 Generated with Claude Code