Array-API-native k×k window filters for cosmicray_median, median_filter, background_deviation_filter and ccdmask

[mwcraig]

Depends on #1009 (rebased on top of it; merge that first).

Part of #971, section 3. The masked-median fix for #984 is a planned follow-up (PR-2) that builds on the NaN handling added here.

What changed and why

scipy.ndimage's window filters are numpy-only, so the four call sites that used them silently copied a jax or dask array to the host, and failed outright for an array on a non-default device. This adds ccdproc/_windowfilters.py, an implementation of the four window reductions ccdproc.core needs, written purely in terms of the array API standard:

helper	replaces
window_rank	ndimage.percentile_filter
window_median	ndimage.median_filter
window_any	ndimage.maximum_filter (on a boolean mask)
window_reduce	ndimage.generic_filter

The approach is deliberately literal: the array is padded once — the standard has no pad, so reflect and nearest are built from flip/concat of edge slices — and each window offset is taken as a shifted slice, with the offsets stacked into a trailing axis, which turns a window reduction into an ordinary reduction over that axis. The rank filters reduce with a new _nanfuncs._nanrank, factored out of nanmedian's sentinel-sort machinery; that is why they reproduce ndimage's rank exactly, including the upper-middle element of an even window.

The stack holds prod(size) copies of the input — 3.8 GiB for an 11x11 window over a 2048x2048 float64 image — so _windowed processes the output in bands of rows sized to a 256 MiB module-constant budget. Bands are cut from the padded array with the window overhang included, so each output pixel sees exactly the window it would have seen unbanded; band_rows stays on the private helpers (a test pins banded == unbanded) and is not exposed on any public function.

Dispatch policy: the numpy path is unchanged

Four wrappers in core.py — _dispatch_median_filter, _dispatch_percentile_filter, _dispatch_maximum_filter, _dispatch_generic_filter — send numpy input to scipy.ndimage exactly as before and everything else to _windowfilters. Concentrating the choice in one place per filter also leaves the seam #984 needs, where numpy input carrying a mask will start taking the native path too.

Call sites converted: _cosmicray_median_array's mbox median, gbox growth and rbox replacement; the public median_filter; background_deviation_filter (which gains an xp argument, matching background_deviation_box); and ccdmask's median and its two percentile filters. ccdmask's byblocks branch is untouched.

median_filter on a non-numpy array now accepts only size and mode. The argument list is bound against ndimage's own signature rather than duplicated, so a repeated or unknown argument still gets ndimage's error message, and footprint/origin/output/cval/axes raise TypeError naming the argument. The numpy path stays a verbatim passthrough.

cosmicray_median's CCDData branch now assigns through nccd._mask, as ccd_process already does: NDDataArray's mask setter runs the value through np.asarray, which was the real source of the "cosmicray_median numpy.asarray" escape the baseline attributed to ndimage.

Documented divergences from ndimage

• Integer input is promoted to the namespace's default real floating dtype; ndimage keeps an integer dtype.
• Only ndimage's reflect and nearest boundary modes exist; anything else raises rather than quietly behaving like reflect.
• A window needing more padding than its axis holds raises, where ndimage re-reflects.
• Cost is O(k2 log k2) per pixel for a k-by-k window, from a sort, against ndimage's O(k**2) selection.
• NaN handling. The rank filters exclude NaNs from a window and take the rank among the values that remain; ndimage sorts NaNs in with the values, above every real number. The one caller this reaches is ccdmask, whose input is a flat ratio that may well contain NaN — it opens by masking the non-finite pixels — so a ratio with NaN in it can give a slightly different mask off numpy. Every other caller filters finite data. The exclusion is deliberate: it is what will let cosmicray_median keep masked pixels out of its median in PR-2. test_windowfilters.py::test_ccdmask_window_filters_exclude_nan_off_numpy pins it in both directions, against an explicit sliding_window_view rank reference off numpy and against ndimage on numpy, asserting first that the two references really disagree.

All of these are in the docs/array_api.rst limitations list.

Verified backends

Dev environment: numpy 2.1.0, scipy 1.14.0, astropy 7.1.0, jax 0.5.0 (JAX_ENABLE_X64=1), dask 2025.7.0, array-api-strict 2.5 (on the non-default device1).

Affected files = test_nanfuncs.py test_windowfilters.py test_blocks.py test_ccdmask.py test_cosmicray.py test_ccdproc.py.

backend	affected files	full pytest ccdproc
numpy	591 passed	1203 passed, 5 skipped
jax	586 passed, 5 skipped	1196 passed, 36 skipped
dask	586 passed, 5 skipped	1197 passed, 16 skipped (run with docs, under the escape ratchet)
array-api-strict	560 passed, 31 xfailed	1166 passed, 10 skipped, 32 xfailed

No failures and no XPASS on any backend. Against the pre-PR branch point, array-api-strict loses 11 xfails, all of them tests that only failed because they went through ndimage.

The escape-baseline ratchet was run the way CI's py312-alldeps-dask-enforce job runs it (CCDPROC_ARRAY_LIBRARY=dask CCDPROC_LOG_ARRAY_ESCAPES=1 CCDPROC_ENFORCE_ESCAPE_BASELINE=1 pytest ccdproc docs). Before the baseline edit it reported exactly three entries as "not hit this run" — _cosmicray_median_array, background_deviation_filter and cosmicray_median; after deleting them it reports OK: no library escapes outside the baseline.

Deviations from the agreed plan

1. Banding is not skipped on dask. The plan said to skip it and let dask's chunking bound memory. That made the dask suite unusable: test_cosmicray_median_rbox took 202 s and test_background_deviation_filter 230 s. The cause is that slicing a chunked band into prod(size) differently-offset windows makes dask realign every one of them — a 21x21 window over a 100x100 image builds a graph of 585,149 tasks, against 912 when the band is first collapsed to a single chunk. So each band is now rechunk(-1)ed (guarded by is_dask_namespace), and banding stays on for dask, since it is the band budget that makes collapsing safe. The dask affected-file run went from 496 s to ~15 s.
2. dask's PerformanceWarning about the chunk-count multiplication is suppressed inside _stack_from_padded, matched by message so dask need not be imported. That multiplication is the algorithm and no caller can act on it, and ccdproc's pytest configuration turns warnings into errors, so leaving it would fail every dask test touching a window filter. Documented in the module docstring.
3. Test-side numpy-isms had to be fixed for the marker removals to mean anything. Dropping the ndimage backend_xfails exposed that add_cosmicrays built its rays by copying the image to numpy — impossible for an array on a non-default device — and that several assertions used .sum(), .mean(), .std() or np.array() on backend arrays. add_cosmicrays now blends a host-side overlay in with a single where; the assertions use xp.* or _to_numpy.
4. One extra core.py fix: nccd.mask | crarr coerced the namespace array to numpy once crarr stopped being a numpy array, so the existing mask is brought into the data's namespace first.
5. One extra marker pruned: test_cosmicray_lacosmic_detects_inconsistent_units began XPASSing on array-api-strict once add_cosmicrays worked (it raises before ever reaching astroscrappy), so its marker is gone, per the "prune on XPASS" note in docs/array_api.rst.
6. The helpers are N-D rather than 2-D-only. The padding was already axis-generic, so generality cost nothing.

Not verified

• The Sphinx docs build was not run locally. sphinx-astropy is not installed in the dev environment. The first CI run caught two nitpicky cross-references to the private _windowfilters module in public docstrings; those are now spelled as literals. docs/array_api.rst was checked with docutils (clean apart from Sphinx-only roles) and the .. _scipy.ndimage: link target the new bullets use was added.
• CuPy. Not installed, so untested, as usual for this repo.
• The XPASS-driven marker removals were confirmed only on macOS, not in CI — docs/array_api.rst asks for CI confirmation before deleting a marker, so that one is worth a second look on the CI logs.

🤖 Generated with Claude Code

https://claude.ai/code/session_01F6f9L1GyMHrKnvLrrgWbWN

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 98.31%. Comparing base (2f866f7) to head (9e0b027).
⚠️ Report is 1 commits behind head on main.

Additional details and impacted files

@@            Coverage Diff             @@
##             main    #1010      +/-   ##
==========================================
+ Coverage   98.15%   98.31%   +0.15%     
==========================================
  Files          10       11       +1     
  Lines        2007     2191     +184     
==========================================
+ Hits         1970     2154     +184     
  Misses         37       37              

Flag	Coverage Δ
dask	97.53% <99.03%> (+0.12%)	⬆️
jax	97.62% <99.03%> (+0.16%)	⬆️
numpy	97.48% <93.71%> (-0.42%)	⬇️

Flags with carried forward coverage won't be shown. Click here to find out more.

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[mwcraig]

Rebased onto main now that #1009 has merged (aa71121 → 9e0b027). The five block-function commits that this branch carried were dropped automatically as already-applied upstream, leaving the nine window-filter commits.

Four files conflicted:

• ccdproc/_nanfuncs.py — Array-API-native block_reduce, block_average and block_replicate off numpy #1009's review dropped the _fill_doc(template, /, ...) generalisation and inlined _blocks' Parameters so that _nanfuncs.py stays unchanged from main. Rather than reintroduce it, _nanfuncs.py is left exactly as main has it (apart from this branch's _nanrank factoring) and _windowfilters._window_doc now fills its own {params} placeholder in six lines of its own. Happy to go back to a shared template argument instead if you'd rather — four functions × a five-parameter block is a better case for it than _blocks was.
• CHANGES.rst — the incoming side still carried the pre-review wording of Array-API-native block_reduce, block_average and block_replicate off numpy #1009's block-functions entry, tagged [#971]. Dropped that as superseded; kept Array-API-native block_reduce, block_average and block_replicate off numpy #1009's two merged entries followed by this branch's window-filter entry, still [#1010].
• docs/array_api.rst — both sides extended the same block-functions bullet. Kept Array-API-native block_reduce, block_average and block_replicate off numpy #1009's full version (the float32/astropy#20360 and jax int32 sentences) and appended this branch's three window-filter bullets after it.
• ccdproc/tests/array_escape_baseline.txt — resolved by content, then regenerated the documented way (CCDPROC_WRITE_ESCAPE_BASELINE=1 with the dask backend over a full-suite run) rather than trusting the hand-merge. Content came back identical, confirming the merge: the three scipy.ndimage entries this branch removes are gone and Array-API-native block_reduce, block_average and block_replicate off numpy #1009's removals stay gone, leaving four. The regeneration also re-narrowed the column padding, which had still been sized for the longest now-deleted name; that is its own commit and changes no entry.

One extra commit, Use the shared helpers from #1009. #1009's review added the array-like guard in _block_dispatch — fall back to numpy when no array library claims the input, but re-raise when it does advertise __array_namespace__. _median_filter_array was written before that landed and resolves its namespace with a bare array_namespace, which regressed ccdproc.median_filter: it documents itself as a passthrough for scipy.ndimage.median_filter and accepted any array-like on every previous release, but a nested list was raising TypeError: list is not a supported array type. It now uses the same guard, so lists and tuples go back to ndimage unchanged. Three tests added. No other refactoring.

Two things I deliberately left alone, both yours to call:

• background_deviation_filter(data, bbox, xp=None) uses a raw caller-supplied xp without passing it through core._namespace_from_module, so a raw module (dask.array rather than the compat namespace) would reach _windowfilters missing xp.bool and the device keyword. It matches its sibling background_deviation_box, which has the same property on main, so fixing one without the other seemed worse than leaving both.
• background_deviation_filter has the same array-like question as median_filter did, but there the pre-existing background_deviation_box behaviour is the stated model, so I did not change it.

tox, all green, on the rebased head:

env	result
codestyle	passed
strict	1240 passed, 39 skipped, 32 xfailed
py312-alldeps-dask-enforce	1269 passed, 42 skipped (baseline ratchet green)
py312-alldeps-jax	1267 passed, 44 skipped
py312-alldeps (numpy, extra)	1269 passed, 42 skipped

— Written by Claude at @mwcraig's direction.

[Copilot]

🟡 Changes recommended

Scalar inputs fail and low-precision rank arithmetic can select incorrect window values.

Get a fresh assessment by requesting another Copilot review.

Pull request overview

Adds Array-API-native window filters while preserving existing SciPy behavior for NumPy inputs.

Changes:

• Implements banded median, percentile, boolean-any, and generic window reductions.
• Dispatches affected core operations by array namespace.
• Expands backend tests, documentation, and escape-baseline coverage.

File summaries

File	Description
ccdproc/_windowfilters.py	Adds native window-filter implementations.
ccdproc/_nanfuncs.py	Adds reusable NaN-aware rank selection.
ccdproc/core.py	Dispatches filters between SciPy and native implementations.
ccdproc/tests/test_windowfilters.py	Tests filtering, padding, banding, and dispatch.
ccdproc/tests/test_nanfuncs.py	Tests NaN-aware rank behavior.
ccdproc/tests/test_cosmicray.py	Enables backend-neutral cosmic-ray tests.
ccdproc/tests/test_ccdproc.py	Removes obsolete backend exclusion.
ccdproc/tests/test_ccdmask.py	Tests dispatcher integration.
ccdproc/tests/array_escape_baseline.txt	Removes resolved NumPy escapes.
docs/array_api.rst	Documents behavior and limitations.
CHANGES.rst	Records the new functionality.

Review details

• Files reviewed: 11/11 changed files
• Comments generated: 3
• Review effort level: Balanced

---

💡 Add a code-review agent skill or configure MCP servers for context-aware, tailored reviews. Learn more in the docs.

[mwcraig]

Half right, and the half that is wrong matters for how it gets fixed.

Wrong: it cannot "stop matching ndimage", because ndimage never takes 2-D float16 at all:

ndimage.median_filter(x16, size=(45,47))  ->  RuntimeError: array type not supported

Right: the arithmetic really is done in s.dtype, so float16 breaks the rank this module documents (int(prod(size) * percentile / 100)). np.float16(2115) is 2116.0, so the median index is 1058 instead of 1057. Against the same call on float64-promoted input:

window_median(x16, (45,47)) vs float64 reference: 2319 / 3600 pixels differ

float32 is safe for any real window (exact integers to 2**24, i.e. k up to 4096), so float16 is the only exposure.

The fix is the same one as the off-by-one at _nanfuncs.py:585: compute the index in xp.float64, not s.dtype, and multiply before dividing, which is what scipy does (rank = int(float(filter_size) * percentile / 100.0)):

index = xp.astype(xp.astype(n, xp.float64) * percentile / 100, n.dtype)

A float16 regression test is worth adding, but as "matches the float64-promoted result", not "matches ndimage".

— Written by Claude at @mwcraig's direction.

[mwcraig]

Valid, and the suggested shape of the fix is right too.

window_median(np.asarray(1.0), 1)            -> IndexError: tuple index out of range
ndimage.median_filter(np.asarray(1.0), size=1) -> array(1.)
window_median(np.arange(5.0), 3)             -> OK (1-D is fine)

The IndexError comes from _default_band_rows's x.shape[0], and _windowed would hit the same thing two lines later. Everything below that does work for 0-D — I checked by hand: _normalize_size(1, 0) is (), _pad_windows is a no-op, _stack_from_padded gives shape (1,) and _nanrank(..., -1) returns 3.0.

So one guard in _windowed is enough:

    if x.ndim == 0 or band_rows >= shape[0]:
        return reduction(_stack_from_padded(_collapse(padded, dask), size, shape, xp))

plus an early return 1 in _default_band_rows for x.ndim == 0. The alternative is to narrow the docstring from "Input array of any rank" to "rank at least 1" — no ccdproc.core caller passes 0-D. Either is fine; the guard is two lines, so I'd do that.

— Written by Claude at @mwcraig's direction.

[mwcraig]

Valid for median_filter; not demonstrated for cosmicray_median.

median_filter is a public entry point with no finiteness check, and the divergence is easy to hit:

x = np.arange(25.0).reshape(5, 5); x[1, 1] = x[2, 2] = np.nan
ccdproc.median_filter(x, 3)                 # ndimage
ccdproc.median_filter(xps.asarray(x), 3)    # window_median
# -> 13 of 25 pixels differ

So the sentence "The one ccdproc function this reaches is ccdmask" is wrong as written and should name median_filter too.

cosmicray_median is a weaker case: it also does no finiteness check, but a NaN in data makes rarr NaN, crarr False there on both backends, and with rbox=0 nothing propagates — a 20x20 frame with one NaN gave identical crmask on numpy and strict. It can still differ through rbox > 0, so "every other caller filters finite data" overstates it; "every other ccdproc caller is normally given finite data" would be accurate.

background_deviation_filter is genuinely unaffected — sigma_func is not an order statistic.

— Written by Claude at @mwcraig's direction.

[mwcraig]

Adversarial review of 9e0b027, focused on science/correctness bugs and on complexity/redundancy/LOC. Two subagent passes, with the top findings run empirically in .tox/strict.

Fix first (both invisible to the current tests):

1. Percentile rank is off by one whenever percentile/100 is not exactly representable (_windowfilters.py:640, _nanfuncs.py:585). Breaks the "reproduces scipy.ndimage exactly" contract for e.g. 29.0 or 57.0 at size 10; the test grid and its _rank_reference share the same regrouping so they cannot see it.
2. cosmicray_median on integer images now mixes int data with a float median off numpy (core.py:2828, 2849): strict raises, jax/dask silently change the output dtype.

Other correctness: _normalize_size rejects numpy ints that ndimage accepts; background_deviation_filter no longer takes plain lists; complex input is silently truncated instead of raising; dask unknown-shape gives an opaque range() error; dask PerformanceWarning is matched by message text; _median_filter_array re-zips args in the wrong positional order. Plus two performance items (window_any builds the full k^2 stack for a running OR; isnan computed twice in _sorted_with_nan_last).

LOC: _windowfilters.py is 755 lines with roughly 158 executable; the _dispatch_* block in core.py is 110 lines holding 8. The inline comments identify about 330 removable lines with no behaviour change, ranked by lines saved over risk: dispatch docstrings (-90), merging the two pad helpers (-50), deleting _window_stack (-43), generalising _fill_doc instead of copying it (-35), _ITEMSIZES (-40, the one item with a tiny observable effect), test dedup (-64), the shared preamble in the three public filters (-18), _collapse (-22), doc repetition (-27).

Details on each line.

[mwcraig]

Science bug: percentile rank off by one. fraction = percentile / 100 loses precision before _nanrank multiplies by n. scipy.ndimage computes n * percentile / 100; this computes n * (percentile / 100). For a 10x10 window (n=100), int(100*29/100.0) is 29 but int(100*(29/100.0)) is 28. Verified in .tox/strict: window_rank(x, 10, 29.0) vs ndimage.percentile_filter(x, 29.0, size=10) on a 24x24 float64 image differs at 432 of 576 pixels (max abs diff 0.16). Same for 57.0 at sizes 10, 20, 100. ccdproc's own 30.9/69.1/50 happen to be safe at realistic sizes, and _PERCENTILES in the tests only covers those, so the grid is blind. Any ccdmask call with non-default ncsig/nlsig gets a backend-dependent result.

Fix: pass percentile through to _nanrank and divide by 100 after multiplying by n.

[mwcraig]

Confirmed, with one correction to the scope.

scipy does the multiply first, in C double (_filters.py, _rank_filter):

rank = int(float(filter_size) * percentile / 100.0)

Reproduced against ndimage.percentile_filter on a 24x24 float64 image:

pct=29.0 size=10  n=100  int(n*p/100)=29  int(n*(p/100))=28  diff 432/576  max 0.16
pct=57.0 size=10  n=100  ->57 vs 56       diff 443/576  max 0.17
pct=57.0 size=20  n=400  ->228 vs 227     diff 360/576  max 0.05
pct=30.9/69.1/50 size=7                    diff 0/576

Correction: "Any ccdmask call with non-default ncsig/nlsig gets a backend-dependent result" is too strong. For the two percentiles ccdmask actually uses, I scanned every (nlsig, ncsig) in 2..63 — 3844 box shapes — and only n = 3000 (50x60 and 60x50) diverges. So ccdmask is close to safe by luck; the bug is in window_rank as a general filter, which is reason enough.

Fix as you describe, and take the dtype part of Copilot's _nanfuncs.py:586 note at the same time — see my reply at _nanfuncs.py:585.

— Written by Claude at @mwcraig's direction.

[mwcraig]

This is where the regrouping bites (see comment at _windowfilters.py:640). Taking percentile here and computing xp.astype(n, s.dtype) * percentile / 100 reproduces ndimage exactly, which is what the docstring promises.

[mwcraig]

Agreed, with one amendment: don't compute it in s.dtype.

Multiplying first fixes the percentile / 100 truncation, but s.dtype is the image dtype, and for float16 the count itself is already wrong before any multiplication — np.float16(2115) is 2116.0 (that is Copilot's point at line 586, which is otherwise mis-stated; ndimage refuses 2-D float16 outright, so there is no parity to lose).

One change covers both:

def _nanrank(x, percentile, axis, xp):
    ...
    # ndimage computes int(float(size) * percentile / 100.0) in C double:
    # multiply before dividing, and never in the image dtype.
    index = xp.astype(xp.astype(n, xp.float64) * percentile / 100, n.dtype)
    index = xp.minimum(index, n - 1)

with window_rank passing percentile straight through instead of fraction = percentile / 100. float64 is mandated by the standard, so every namespace has it; jax without jax_enable_x64 silently lands on float32, which I measured as differing from float64 in 255 of 2,004,002 (n, percentile) pairs — against the current code, which is systematically wrong.

The docstring's fraction parameter becomes percentile in [0, 100], and _rank_reference in the tests has to be regrouped the same way or it will keep hiding this.

— Written by Claude at @mwcraig's direction.

[mwcraig]

Science bug: integer images. ndimage.median_filter preserves int dtype, but window_median promotes to float via _promote_to_real. So off numpy, data - marr mixes int and float: array-api-strict raises TypeError: array_api_strict.int64 and array_api_strict.float64 cannot be type promoted together (confirmed), and jax/dask silently return float ndata where numpy returns int. Raw CCD frames from FITS are commonly int16/uint16, so this is the normal case, not an edge case. The comment above ("make sure that mdata is the same type as data") is now false.

Fix at the right depth: promote data once at the top of _cosmicray_median_array (as _blocks._promote_for_division does), or have _dispatch_median_filter cast the native result back to the input dtype.

[mwcraig]

Confirmed. _cosmicray_median_array(int16_image, ...):

numpy  int16 -> OK, ndata.dtype = int16
strict int16 -> TypeError: array_api_strict.int16 and array_api_strict.float64
                cannot be type promoted together
dask   int16 -> OK, ndata.dtype = float64   (line 2849, rbox > 0)

and the same through the public cosmicray_median. So three backends, three answers: raise, int16, float64.

Of the two fixes you suggest, promoting once at the top of _cosmicray_median_array is the better one:

    data = _blocks._promote_for_division(data, xp)

rarr = (data - marr) / error_image is a true division, so data has to be floating anyway for array-api-strict; casting marr back to int instead would only move the TypeError to that line. It does change the public return dtype for integer input from int16 to float64 on numpy — a behaviour change, but the one that makes every backend agree, and it needs a CHANGES.rst line.

The stale comment at 2848 ("make sure that mdata is the same type as data") should go with it.

— Written by Claude at @mwcraig's direction.

[mwcraig]

Same mixed-kind problem as line 2828: xp.where(crarr, mdata, data) with float mdata and int data.

[mwcraig]

Confirmed, and this is the line where the backends actually produce different values rather than an error. With rbox > 0:

numpy int16 -> ndata.dtype = int16
dask  int16 -> ndata.dtype = float64

Line 2828 raises first on array-api-strict, so on that backend this line is only reachable once 2828 is fixed — which is why both should be fixed together, by promoting data once at the top of _cosmicray_median_array (see my reply at 2828). That makes xp.where(crarr, mdata, data) same-kind on every backend and leaves this line untouched.

— Written by Claude at @mwcraig's direction.

[mwcraig]

isinstance(size, int) rejects numpy integers, which scipy.ndimage accepts. Verified: window_median(x, np.int64(3)) raises TypeError: 'numpy.int64' object is not iterable (falls through to tuple(size) at line 184), and window_median(x, (np.int64(3), np.int64(3))) raises ValueError: size entries must be positive integers, while ndimage.median_filter(x, size=np.int64(3)) succeeds. So ccdmask(ratio, ncmed=np.int64(7)) works on numpy and dies on jax/dask/strict.

_blocks._block_size (same PR stack) already handles this correctly (try: list(block_size) plus float(size).is_integer()). Reuse it or factor one shared helper rather than a second, stricter validator.

[mwcraig]

Confirmed, including the ccdmask consequence:

window_median(x, np.int64(3))         -> TypeError: 'numpy.int64' object is not iterable
window_median(x, (np.int64(3),) * 2)  -> ValueError: size entries must be positive integers
ndimage.median_filter(x, np.int64(3)) -> OK
ccdmask(CCDData(numpy),  ncmed=np.int64(7)) -> OK
ccdmask(CCDData(strict), ncmed=np.int64(7)) -> ValueError: size entries must be positive integers

_blocks._block_size(np.int64(3), 2) returns (3, 3) and _block_size(3.0, 2) returns (3, 3), so that helper is already correct.

One caveat on "reuse it": the two validators are not interchangeable. _block_size broadcasts a length-1 sequence to ndim and accepts integral floats, both to match astropy.nddata; ndimage accepts np.int64(3) but rejects 3.0 (TypeError: 'float' object cannot be interpreted as an integer), and _normalize_size has to keep rejecting 3.0 to stay ndimage-compatible. So share the scalar-vs-sequence detection, not the whole function:

    try:
        sizes = tuple(size)
    except TypeError:
        sizes = (size,) * ndim
    if any(not isinstance(k, numbers.Integral) or k < 1 for k in sizes):
        raise ValueError(...)

numbers.Integral covers int, numpy integers and anything else that registers, and still rejects 3.0.

— Written by Claude at @mwcraig's direction.

[mwcraig]

Test mirrors the implementation. This _rank_reference (and the one in test_nanfuncs.py:197) re-derives _nanrank step for step, including the n * fraction regrouping, which is exactly why the off-by-one at _windowfilters.py:640 slipped through. A genuinely independent reference that borrows ndimage's own padding is 7 lines:

def _rank_reference(data, size, percentile, mode="reflect"):
    def pick(w):
        v = np.sort(w[~np.isnan(w)])
        return np.nan if v.size == 0 else v[min(int(v.size * percentile / 100), v.size - 1)]
    return ndimage.generic_filter(data, pick, size=size, mode=mode)

It also retires _NUMPY_PAD_MODES and the sliding_window_view import. Please also add a percentile like 29.0 or 57.0 at size 10 to _PERCENTILES so the grid covers the regrouping case.

[mwcraig]

Agreed, and the proposed reference checks out. I ran it:

• Against ndimage on finite input over the full _SIZES x _MODES grid at percentiles 0, 29, 30.9, 50, 57, 69.1, 100 — exact everywhere (with maximum_filter as the comparison at 100, as the existing test already does).
• Against the current _rank_reference on the NaN image from test_nan_windows_match_an_explicit_rank_reference and on the inf-carrying flat ratio from _flat_ratio_with_non_finite_pixels (used at line 483) — identical, so it is a drop-in for both call sites. ~np.isnan keeps the infinities, as the current one does.
• With percentile / 100 hoisted out of pick, i.e. the bug at _windowfilters.py:640 transplanted into the reference, it differs at 434/576 pixels for 29.0 at size 10. So this reference would have caught it and the current one structurally cannot.
• Cost is negligible: the nine reference calls in the NaN test take 0.01 s, and the 31x29 / 7x7 ccdmask case 0.02 s.

Agreed on adding 29.0 at size 10 to the grid too.

One correction on test_nanfuncs.py:197. That reference is not at fault in the same way: _nanrank's signature takes fraction, and min(floor(n * fraction), n - 1) is its documented contract, so mirroring it there is mirroring the spec, not the implementation. The regrouping bug lives entirely in window_rank's decision to divide first. If _nanrank grows a percentile argument per your comment at _nanfuncs.py:585, then that reference must change to np.trunc(n * percentile / 100) in that order — and there is no independent oracle available for it (ndimage has no NaN-aware reduction along an axis), so the honest fix there is just to write the arithmetic in ndimage's grouping and say why.

— Written by Claude at @mwcraig's direction.

[mwcraig]

Redundant grid. window_median is a one-line delegate to window_rank(..., 50.0, ...), and test_window_rank_matches_ndimage already parametrises 50.0 over the same _SIZES x _MODES, so these 14 cases are a strict subset of the 70-case rank grid. Keep one median-vs-ndimage case for the public name and drop the rest.

[mwcraig]

Mostly right, but the subset claim needs one fact checked first, and it holds: I compared ndimage.median_filter with ndimage.percentile_filter(..., 50.0) over the whole _SIZES x _MODES grid and they are bit-identical, including the even sizes 4 and (4, 2) where the conventions could plausibly have differed. So with window_median a one-line delegate to window_rank(..., 50.0, ...), the 14 median cases really are covered by the 70-case rank grid.

Two things the median test still does carry, both cheap to keep:

• It exercises the public delegate. If window_median ever passed 0.5 instead of 50.0, the rank grid would not notice.
• It compares against median_filter rather than percentile_filter, so it is the only thing pinning the scipy invariant the delegation rests on.

So keep two cases rather than one — an odd size and an even one, @pytest.mark.parametrize("size", [3, 4], ids=str) at the default mode — and drop the other 12. About -12 cases for no loss of coverage.

— Written by Claude at @mwcraig's direction.

[mwcraig]

Helpers duplicated from test_blocks.py. _as_test_array/_assert_matches have the same bodies as _to_xp/_assert_matches in test_blocks.py:161-178, and test_blocks' _assert_same_namespace_and_device is open-coded here as a whole test (263-279). Move the three into pytest_fixtures.py (or conftest) and import them in both modules. Also test_result_stays_in_the_input_namespace_and_device (254) and test_integer_input_is_promoted_to_a_real_floating_dtype (282) repeat almost the same four-lambda call list; one module-level _CALLS and one merged test asserting namespace, device and dtype in a single pass would do. About -34 lines across these.

[mwcraig]

Three claims here, and they do not all hold up.

_as_test_array vs _to_xp — yes, identical two-line bodies. Sharing saves two lines and costs an import; worth doing only as a side effect of moving the third helper.

_assert_matches — no, the bodies differ, and materially. test_blocks' version asserts result.dtype == expected.dtype and calls xpx.isclose(result, expected); this one casts the reference to float, makes no dtype assertion, and passes equal_nan=True. Both differences are load-bearing here: these filters deliberately promote integer input (so dtype equality with a numpy reference is wrong), and test_window_of_only_nan_is_nan and the whole NaN grid depend on equal_nan. Merging them means a helper with two flags, which is not a saving.

_assert_same_namespace_and_device — agreed, that one is genuinely open-coded at 276-279.

On merging the two parametrised tests: the lists are four lambdas and three, and window_any is excluded from the dtype test deliberately, because it returns bool rather than a real floating dtype. One list with the expected dtype kind attached fixes that and is better than either test alone, since window_any's dtype stops being unpinned:

_CALLS = [
    pytest.param(lambda d: window_median(d, 3), "real floating", id="window_median"),
    pytest.param(lambda d: window_rank(d, 3, 69.1), "real floating", id="window_rank"),
    pytest.param(lambda d: window_any(d, 3), "bool", id="window_any"),
    pytest.param(lambda d: window_reduce(d, 3, xp.std), "real floating", id="window_reduce"),
]

fed integer input, asserting namespace, device and xp.isdtype(result.dtype, kind) in one pass.

Realistic total is about -20, not -34: the shared helper is worth ~5 and the merged test ~15.

— Written by Claude at @mwcraig's direction.

[mwcraig]

Docs repetition. The stack-and-reduce/banding rationale and the dask chunk-count warning are explained here, again at greater length in _windowed's Notes (541-559), pointed to from _collapse's Notes, and restated in the comment at 375-381. Keep the full explanation once in _windowed and cut this to the first paragraph plus a pointer. Similarly the NaN divergence is written out in window_rank Notes, window_median Notes, the core.py banner, median_filter Notes, background_deviation_filter Notes, CHANGES, and a 17-line docs bullet; window_median already says "see there", so its copy can go.

[mwcraig]

Half of this holds; the NaN half does not survive a grep.

The banding/warning half — agreed, with a correction on where it lands. The module docstring's paragraph at 28-34 and the comment at 375-381 are the same seven-line explanation of the dask PerformanceWarning twice over, and the suppression is in _stack_from_padded, not _windowed, so the pointer has to go there. Cut 28-34 to one sentence ("one dask PerformanceWarning is deliberately swallowed in _stack_from_padded; see there") and cut the banding sentence out of 14-24, leaving the stack-and-reduce rationale and the cost, with _windowed keeping the full banding argument. About -10.

The NaN half — the repetition is not where the comment says. I grepped NaN across core.py: the only hits in the two public docstrings are at 472-473, 498, 1629-1630 and 1649-1654, which are nanmedian/sigma_func semantics. Neither median_filter's Notes (2201-2209) nor background_deviation_filter's (2803-2809) mentions the divergence at all — they cover dtype promotion and the mode restriction only. And window_median's Notes is already nothing but a pointer ("Exactly window_rank at percentile=50; see there for how NaNs are treated..."), so there is no copy there to delete.

So the real state is: one full statement in window_rank, a two-line summary in the core.py banner (which is where a reader of the dispatch code needs it), a clause in CHANGES.rst:46, and the docs bullet. That is not over-documented — and per Copilot's comment at docs/array_api.rst:230, median_filter and cosmicray_median reach this too, so those two public Notes arguably need a sentence added rather than removed.

— Written by Claude at @mwcraig's direction.

[mwcraig]

Docs bloat. This 17-line bullet re-argues the NaN divergence from scratch. The user-visible fact fits in two sentences: the median/percentile window filters exclude NaNs from a window and rank among the remaining values, where scipy.ndimage sorts NaNs in with them, so a flat ratio containing NaN can give ccdmask a slightly different mask off numpy; every other caller filters finite data first. Also worth adding here that dask arrays need a fully known shape for the window filters, not just the block functions.

[mwcraig]

Agreed — the bullet runs 214-232, 19 lines, and a reader of a "what differs between array libraries" page needs the fact and the consequence, not the derivation. Draft:

+ The window filters that take an order statistic -- the median and the
  percentile -- **exclude NaNs from a window** and rank among the values
  that remain, where `scipy.ndimage`_ **sorts NaNs in with them**, above
  every real number. Results on data containing NaN therefore differ
  between `numpy`_ and every other array library: ``median_filter`` and
  ``cosmicray_median`` accept such data, and ``ccdmask``'s flat ratio
  routinely contains it, so a ratio with NaN in it can give a slightly
  different mask off `numpy`_. Infinities are ordinary values to both.
  These filters also need a fully known shape, so a `dask`_ array with
  unknown chunk sizes must have ``compute_chunk_sizes()`` called on it
  first.

11 lines, so about -8, and it folds in the dask known-shape note. Two deliberate changes beyond shortening: the "(The exclusion is deliberate: it is what will let cosmicray_median keep masked pixels out of its median.)" parenthesis goes — that is a note to ourselves about #984, not something a user can act on — and "The one ccdproc_ function this reaches is ccdmask" is replaced, because it is not true: Copilot's comment at line 230 is right that the public median_filter takes arbitrary user arrays with no finiteness check, and cosmicray_median does not enforce finite input either.

— Written by Claude at @mwcraig's direction.