From c7848921ffa4a140a470a66781f16166f9ecb6c3 Mon Sep 17 00:00:00 2001 From: Fabian Isensee Date: Tue, 3 Nov 2026 16:20:25 +0100 Subject: [PATCH] Keep the float dtype in mean(), as NumPy does mean() divides the sum by the count, and that division can promote a float32 (or complex64) result to float64: - With NumPy 1.x, a float32 scalar divided by a Python int becomes float64 (value-based casting), so every full a.mean() of float32 returned float64; so did axis reductions over more than 65535 elements (an array over such an int). - With where=, the count is an int64 sum, so a.mean(where=...) of float32 returned float64 with any NumPy. complex64 with where= came back complex128 as well, because the sum of its where() already was. NumPy returns the input's float dtype (or a floating dtype=). The quotient is now cast to that dtype; integer and bool inputs keep their float64 mean. std() and var() take their mean from here too. Test: test_mean_dtype_matches_numpy (float32/float64/complex64/int32/ bool, axis, keepdims, where) has 16 of 40 cases failing without this on NumPy 1.26, and 8 on NumPy 2.5. Co-Authored-By: Claude Opus 5.5 --- src/blosc2/lazyexpr.py | 6 ++++++ tests/ndarray/test_reductions.py | 20 ++++++++++++++++++++ 2 files changed, 26 insertions(+) diff --git a/src/blosc2/lazyexpr.py b/src/blosc2/lazyexpr.py index 732c1141e..a709b111b 100644 --- a/src/blosc2/lazyexpr.py +++ b/src/blosc2/lazyexpr.py @@ -4133,6 +4133,12 @@ def mean( if np.isscalar(num_elements) and num_elements == 0: raise ValueError("mean of an empty array is not defined") out = total_sum / num_elements + # NumPy's dtype for a floating `dtype` or input: NumPy 1.x makes a float32 scalar over a + # Python int float64 (and an array over an int above 65535), and the int64 counts of + # `where` make it float64 with any NumPy + mean_dtype = np.dtype(self.dtype if dtype is None else dtype) + if mean_dtype.kind in "fc" and isinstance(out, (np.ndarray, np.generic)): + out = out.astype(mean_dtype, copy=False) out2 = kwargs.pop("out", None) if out2 is not None: out2[:] = out diff --git a/tests/ndarray/test_reductions.py b/tests/ndarray/test_reductions.py index 89982a550..65a9eeada 100644 --- a/tests/ndarray/test_reductions.py +++ b/tests/ndarray/test_reductions.py @@ -977,3 +977,23 @@ def test_reductions_where_ndarray_and_lazyexpr(): np.testing.assert_allclose(a.std(axis=0, where=mask), np.std(na, axis=0, where=nmask)) np.testing.assert_allclose(a.min(axis=0, where=mask), np.min(na, axis=0, where=nmask, initial=np.inf)) np.testing.assert_allclose(a.max(axis=0, where=mask), np.max(na, axis=0, where=nmask, initial=-np.inf)) + + +@pytest.mark.parametrize("dtype", [np.float32, np.float64, np.complex64, np.int32, np.bool_]) +@pytest.mark.parametrize("axis", [None, 0]) +@pytest.mark.parametrize("keepdims", [False, True]) +@pytest.mark.parametrize("use_where", [False, True]) +def test_mean_dtype_matches_numpy(dtype, axis, keepdims, use_where): + # A float32 mean stays float32: NumPy 1.x makes a float32 scalar over a Python int float64 + # (and an array over a count above 65535), and the int64 counts of where= do so in any version + rng = np.random.default_rng(0) + na = (10 * rng.normal(size=(70_000, 3))).astype(dtype) + nmask = rng.random(na.shape) < 0.7 + a, mask = blosc2.asarray(na), blosc2.asarray(nmask) + kwargs = {"axis": axis, "keepdims": keepdims, **({"where": nmask} if use_where else {})} + res = a.mean(**{**kwargs, "where": mask} if use_where else kwargs) + assert np.asarray(res).dtype == np.mean(na, **kwargs).dtype + # Values against a float64 reference: these means are close to 0 + expected = np.mean(na.astype(np.result_type(na.dtype, np.float64)), **kwargs) + assert np.shape(res) == np.shape(expected) + np.testing.assert_allclose(res, expected, rtol=1e-5, atol=1e-5)