Skip to content

fix: overflow-safe stats without needless copies#4223

Open
ianna wants to merge 5 commits into
mainfrom
ianna/stats_overflow_fix
Open

fix: overflow-safe stats without needless copies#4223
ianna wants to merge 5 commits into
mainfrom
ianna/stats_overflow_fix

Conversation

@ianna

@ianna ianna commented Jul 22, 2026

Copy link
Copy Markdown
Member

Where the overflow actually is

ak.sum's reducer already promotes integer/bool leaves to the platform integer
(int64/uint64) accumulator — see _promote_integer_rank in
src/awkward/_reducers.py and the awkward_reduce_sum<result, input, offsets>
kernel templated separately on accumulator and input dtype. So:

  • ak.sum(int32_array) already accumulates in int64. A plain integer mean
    without weights
    (sum(x) / count) is overflow-safe up to int64 today and
    needs no cast.
  • The real overflow is the element-wise products computed in Python before
    the reduction: x * x, x * weight, x * x * weight in var, std,
    covar, corr, moment. int32 * int32 stays int32 and overflows
    per-element — long before int64 accumulation can help.

So the original #3527 PR's "cast the whole input to float64 up front, always" is
heavier than needed: it copies floating inputs that never overflow, and it casts
in a spot that isn't the only thing that matters.

The two objections on the #3527 PR, addressed

  1. "Binary floats can't represent decimals exactly." True, but that's an
    argument against floating-point in general — it applies equally to x / 2,
    which we already promote. The alternative here is silent integer overflow,
    which is strictly worse than rounding error. NumPy makes exactly this choice
    for mean/var/std.
  2. "The cast triggers device-to-host copies on GPU." A dtype cast via
    ak.values_astype runs on the array's own backend — cupy.ndarray.astype on
    the cuda backend. It stays on-device; there is no host round-trip. The cost on
    GPU is float64 memory (2× an int32 buffer) and lower float64 throughput,
    not a D2H transfer. The solution below only pays that cost for integral
    inputs
    and only once per call.

Solution

Two tiers. Tier 1 is a small, self-contained Python PR. Tier 2 is the fully
copy-free version and needs a reducer tweak -- a separate PR to follow.

Tier 1 — promote only integral leaves, once, on-backend

A shared helper (e.g. in src/awkward/_layout.py):

import awkward as ak
from awkward._nplikes.numpy_like import NumpyMetadata

np = NumpyMetadata.instance()


def _has_integral_leaf(layout) -> bool:
    """True if any NumpyArray leaf is bool/integer (kinds 'b', 'i', 'u')."""
    found = False

    def action(node, **kwargs):
        nonlocal found
        if node.is_numpy and node.dtype.kind in "biu":
            found = True
        return None

    ak._do.recursively_apply(layout, action, return_array=False)
    return found


def promote_integral_to_float64(array):
    """Cast bool/integer data to float64 to avoid overflow in summations,
    matching NumPy's reduction dtype.

    Floating and complex data is returned unchanged (no copy). The cast runs on
    the array's own backend, so on GPU it stays on device — no host transfer.
    """
    if _has_integral_leaf(array.layout):
        return ak.operations.values_astype(array, np.float64, highlevel=True)
    return array

Then, in each statistic that forms products, promote once and reuse the
result across every reduction. For ak.var (src/awkward/operations/ak_var.py),
the change is one line plus reuse:

    x = ctx.wrap(x_layout)
    weight = ctx.wrap(weight_layout, allow_other=True)

    # Overflow guard: integer squares/products overflow before reduction.
    # Promote integral leaves to float64 once (no-op + zero copy for floats).
    x = promote_integral_to_float64(x)
    if weight is not None:
        weight = promote_integral_to_float64(weight)

    ...
    # sumwx, sumwxx (and sumw) now reduce float64 data — one promotion, reused.

The actual scope is smaller than "every stat," because most paths are already
safe:

  • ak_var.py — needs it (sum(x*x)). ak_std.py inherits the fix for
    free: it just calls ak.var and takes sqrt.
  • ak_moment.py — needs it (sum(x**n)).
  • ak_mean.py — needs it only in the weighted branch (sum(x*weight)).
    The unweighted path (sum(x)/count) is already overflow-safe to int64 via
    the reducer's int64 accumulator, and the final division yields float64, so
    it is left copy-free.
  • ak_covar.py / ak_corr.py — need nothing. They center first,
    (x - xmean) * (y - ymean), and xmean is already float64 (it comes from
    ak.mean), so the products are float before any multiply. Their internal
    ak.mean call inherits the weighted-branch fix if integer weights are passed.

So the real change is three files: ak_var.py, ak_moment.py, and the weighted
branch of ak_mean.py — plus the shared helper.

If you additionally want unweighted integer mean to match NumPy's
float64-accumulator behavior (guarding int64 overflow too), use Tier 2 rather
than adding a copy.

Net effect vs. the #3527 PR: floating inputs are untouched (zero copy — the common
case in analysis), integral inputs are promoted exactly once instead of
implicitly re-cast inside each x*x / x*weight, and the cast provably stays on
the GPU backend.

@github-actions github-actions Bot added the type/fix PR title type: fix (set automatically) label Jul 22, 2026
@ikrommyd

Copy link
Copy Markdown
Member

I've been doing exactly the same in #3527 and you were arguing against it since it has to copy the whole array at once. Which is true, when it's an array of integers, a whole array of floats has to be allocated at once of the same size to do so. It's not like numpy's np.sum(..., dtype=np.float64) which does the casting while iterating inside the loop.

@ianna
ianna force-pushed the ianna/stats_overflow_fix branch from 7f3de4a to 3a2ed18 Compare July 22, 2026 18:57
@ianna

ianna commented Jul 22, 2026

Copy link
Copy Markdown
Member Author

I've been doing exactly the same in #3527 and you were arguing against it since it has to copy the whole array at once. Which is true, when it's an array of integers, a whole array of floats has to be allocated at once of the same size to do so. It's not like numpy's np.sum(..., dtype=np.float64) which does the casting while iterating inside the loop.

tier 2 PR to follow

@ianna
ianna requested a review from TaiSakuma July 22, 2026 19:00
@github-actions

Copy link
Copy Markdown

The documentation preview is ready to be viewed at http://preview.awkward-array.org.s3-website.us-east-1.amazonaws.com/PR4223

@codecov

codecov Bot commented Jul 22, 2026

Copy link
Copy Markdown

Codecov Report

✅ All modified and coverable lines are covered by tests.
⚠️ Please upload report for BASE (main@27a23cd). Learn more about missing BASE report.
⚠️ Report is 1 commits behind head on main.
✅ All tests successful. No failed tests found.

Additional details and impacted files
Files with missing lines Coverage Δ
src/awkward/_layout.py 92.34% <100.00%> (ø)
src/awkward/operations/ak_mean.py 94.44% <100.00%> (ø)
src/awkward/operations/ak_moment.py 100.00% <100.00%> (ø)
src/awkward/operations/ak_var.py 89.65% <100.00%> (ø)

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

type/fix PR title type: fix (set automatically)

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants