[FabianIsensee]

Hi Francesc,

this is a PR into #741, as you suggested in #741 (comment). It fixes the two slowdowns I reported in #741 (comment). It is one commit on top of improved-miniexpr at 7db7b09, taken from the deferred lazy-expression commit of #744.

Why: two common cases fell off miniexpr, and since numexpr is gone they now run chunk by chunk through evaluate_chunk().

What:

• 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, because miniexpr reductions mis-evaluate a negative literal on the right of a binary operator ((o0 - -0.5)).
  • Other cases (axis, where, other dtypes) keep the general path. std() and var() now 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. The fused expression then referenced a missing operand, miniexpr failed to compile it, and the evaluation fell back. The where() is now fused before specializing.
• Small: 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, default chunks, 28 threads. Medians of two interleaved rounds; all builds against c-blosc2 main (e5086bf).

	main	#741 (7db7b09)	#741 + this PR
a.std()	75 ms	470 ms	53 ms
a.mean(), a.std()	88 ms	498 ms	78 ms
blosc2.where(a > 100, a, 0).compute()	485 ms	690 ms	469 ms
a.std(), 100x100	0.62 ms	0.73 ms	0.52 ms
(a + b)[:]	118 ms	112 ms	111 ms
a.sum(axis=0)	29.5 ms	29.4 ms	29.7 ms

std/var may differ from before in the last ulps.

Tests:

• New: test_where_scalar_uses_miniexpr checks that miniexpr compiles the fused expression; its 8 scalar cases fail without this change.
• New: test_std_var_full_reduction (float32/float64, ddof, keepdims, means of 0 and +-1e4, with and without miniexpr), ..._wider_dtype and ..._nonfinite.
• pytest tests: 10707 passed. The 2 failures are the b2view http tests that need aiohttp, which also fail without this change and on main.
• ruff passes.

Not in this PR: the rest of that commit (expr[:] straight into the NumPy result: (a + b)[:] about 1.6x; chunked reductions into reused buffers) and the result-dtype fix of check_dtype() (hypot, copysign, nextafter of integers; true division of bools). I'll send those as separate PRs into #741 if you like.

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

Best,
Fabian

🤖 Generated with Claude Code

[FrancescAlted]

Thanks for the optimization work. I reviewed this PR and reproduced two correctness regressions in the newly enabled scalar where() miniexpr path, compared with the PR base (7db7b09). Could you take a look?

1. Float32 comparisons can select the wrong branch

import blosc2
import numpy as np

a = blosc2.asarray(
    np.full(12, 0.1, dtype=np.float32),
    chunks=(6,), blocks=(3,),
)
print(blosc2.where(a > 0.1, 1, 0)[:])

The base and NumPy return all zeros; this PR returns all ones. blosc2.where(a > 0.1, a, 0)[:] is affected too: the base returns zeros, while the PR returns 0.1.

This appears to expose a difference in scalar-promotion semantics: comparing the stored float32 value against a higher-precision 0.1 literal changes the condition. The newly fused/specialized expression in fast_eval() now reaches that native path instead of falling back.

2. Large uint64 branch scalars silently lose precision

a = blosc2.asarray(
    np.arange(12, dtype=np.uint64),
    chunks=(6,), blocks=(3,),
)
expr = blosc2.where(a > 5, np.uint64(2**63 + 1), a)
print(expr[:])
print(expr.compute()[:])

The selected elements should equal 9223372036854775809, as they do on the base. This PR returns 9223372036854775808 through both evaluation methods. The scalar appears to lose integer precision along the native specialization/compilation path.

Both issues concern the changed dispatch around src/blosc2/lazyexpr.py:1893–1910. Ideally the native path should preserve the comparison semantics and exact integer constants; otherwise affected cases should retain fallback. Regression tests for these two examples would help.

Validation: using the blosc2 conda environment and isolated copies of the PR/base lazyexpr.py with the same installed native extension, the two changed test files passed under the default selection (1,417 passed, 7,072 deselected), including all 110 newly added cases. The examples above fail on PR head d153173 and pass on its base. I did not rebuild the native dependencies for this review.

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

[FabianIsensee]

Oh my bad, I will take a look immediately

[FabianIsensee]

Thanks for catching these, Francesc. Both come from miniexpr typing the whole fused where(cond, x, y) by its output dtype:

1. With an int64 output (where(a > 0.1, 1, 0)), the literal 0.1 turned into 0. With a float32 output, miniexpr compared the float32 array against the double 0.1. NumPy rounds 0.1 to float32 first. Arithmetic in the condition is also computed in the output dtype, so where(a + b == c, 1, 0) with float32 operands was wrong as well.
2. miniexpr parses literals into C doubles, so 2**63 + 1 became 2**63. Inside where(), 64-bit integer operands and branch values also go through doubles.

Before this PR, array branches were already fused, and they had the same problems. For example, where(a > 0.1, a, b) returns 0.1 on 7db7b09. The scalar case only seemed right because it fell back.

Fix (new commit on top, python-only):

• where() is fused only when computing the condition in the output dtype gives NumPy's result. That means one of:

  • every operand of the condition has the output dtype;
  • the condition only compares operands that cast safely to the output dtype, against literals that their dtypes hold exactly.

  It is never fused with 64-bit integers. Anything else falls back to the chunked evaluation. Common patterns stay on miniexpr: where(a > 100, a, 0), where(mask, a, 0), where(labels == 3, a, 0), where(a > 100, 1.0, 0.0), where(a + b > 100, a, 0).

• Expressions with integer literals beyond 2**53, written or inlined, are not sent to miniexpr. This also fixes a + np.uint64(2**63 + 1), which was wrong before this PR too.

Tests:

• Your two examples.
• A NumPy comparison grid: 7 condition dtypes × 7 literals × >/== × 5 output dtypes × scalar/array/mixed branches, through [:] and compute().
• Arithmetic conditions, large integer constants, and checks that the patterns above still compile with miniexpr.
• Against the PR head without the fix, 38 of the 70 new and touched tests fail. A random differential run (7,900 where() expressions, mixed dtypes) has 0 mismatches with the fix, compared with 34 per 1,300 before.
• pytest tests: 10767 passed. The 2 failures are the b2view http tests that need aiohttp, as on main. With NumPy 1.26 the new tests pass, and the NEP 50 grids are skipped like the existing ones.

Not changed here, but found on the way (all also on 7db7b09): reductions have the same literal problem ((a > 0.1).sum() with float32 a counts 12 instead of 0; (a == 2**53).sum() on int64). Plain comparisons against a NumPy float64 scalar use float32 (a > np.float64(0.1)). These need fixing in miniexpr itself.

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

[FrancescAlted]

LGTM. Thanks for noting the new issues; we will address them separately.

[FrancescAlted]

Update: the latest commits implementing plans/migrate-typed-expression-evaluator.md address the correctness issues reported above, including those listed as “Not changed here.” Ordinary expressions now use typed evaluation by default, preserving intermediate operation dtypes independently of output/reduction dtypes, strong NumPy scalar types, and exact integers.

Reported issue	Current behavior
where() output dtype changes condition comparisons/arithmetic	Intermediate operations retain their own dtypes
Large integer literals and 64-bit branches lose precision through doubles	Typed execution preserves exact integers
(a > 0.1).sum() and large-integer comparison reductions	Comparison typing is preserved independently of the accumulator
a > np.float64(0.1) incorrectly uses float32 semantics	Strong NumPy scalar types survive expression fusion

Validation against the unset/default evaluator with NumPy 2.5.4: 75 relevant regression tests passed, covering the where() comparison grid, arithmetic conditions, large integers, scalar promotion, and reductions. The migration plan also records successful cross-platform, NumPy 1.26, actual Pyodide, and installed-wheel CI qualification.

This fixes ordinary-expression behavior by replacing the uniform-dtype execution route, not by repairing the old uniform evaluator. The temporary legacy escape hatch and explicit compatibility DSL kernels remain separate contracts.