🍾 Array-API fallback for sigma_func (MAD) when the namespace is not numpy

[mwcraig]

Part of #929 (the sigma_func half). Companion PR: #1001 (the Combiner.sigma_clipping half), built off the same main.

sigma_func — the default uncertainty_func of median_combine, also used by background_deviation_box/background_deviation_filter — wrapped astropy.stats.median_absolute_deviation, which is numpy-only: on any other array namespace the data was converted to numpy, and on array-api-strict the call failed outright (8 of the 13 remaining strict failures).

Policy: when the namespace is numpy, sigma_func still calls astropy for the deviation, and the numpy path is unchanged for plain arrays and unmasked CCDData. A masked CCDData now has its mask honored on every backend — masked pixels are excluded from the statistics via a mask→NaN substitution made before the namespace split, and an all-masked slice gives NaN — as the second CHANGES bullet says. Every other namespace goes to a new private core._mad_fallback that computes the median absolute deviation purely in terms of the array API, on the input's device, using the sort-based medians in _nanfuncs (nanmedian for ignore_nan=True, median otherwise). The fallback promotes int/bool input to the namespace default float, flattens for axis=None, reduces over tuples of axes (permute last + merge), and validates axes. sigma_func stays the same function object because median_combine tests uncertainty_func is sigma_func.

Also: test_combiner_with_scaling built its reference with xp.asarray((a, b, c)) on a tuple of arrays, which array-api-strict rejects; it now uses xp.stack. The test never got that far before because it failed first in sigma_func.

Tests: the fallback is exercised directly on every backend (the strict job uploads no coverage) against astropy over axis ∈ {None, int, numpy int, tuples incl. negative} × ignore_nan, with int/bool/float32 input and all-NaN slices; the public entry point is checked against astropy.stats.mad_std, for namespace/device of the result, for the fallback branch on numpy (patched is_numpy_namespace), and for the CCDData mask. Design was prototype-verified with zero mismatches vs astropy on numpy, strict (device1), jax and dask.

Docs: changelog entry; the "What limitations should I be aware of?" list in docs/array_api.rst is refreshed (it only mentioned the nanmedian fallback, but nansum/nanmean/nanstd, median and now the MAD also fall back). The core.py sigma_func line leaves the escape baseline (confirmed by a full-suite dask regeneration; the diff is that one line).

Dependencies: the minimum numpy is now 2.0 (normalize_axis_tuple comes from its numpy 2 location, numpy.lib.array_utils, with no 1.26 shim — the oldestdeps job caught it), which forces astropy ≥ 6.1 (6.0 caps numpy below 2) and reproject ≥ 0.14 / astroscrappy ≥ 1.2 (the oldest numpy-2 binaries). The oldestdeps env and the bottleneck CI job (numpy126 → numpy200) move with it. Closes #1003.

Verified locally

• numpy: 616 passed
• array-api-strict: 5 failed (down from 13: the three sigma_clipping tests for the companion PR, Array API: units/Quantity handling with non-numpy arrays #936, Consider marray as a uniform masked-array representation across array backends #983), 0 xpassed
• jax (X64): 615 passed
• dask: 609 passed; escape-baseline enforce run clean

🤖 Generated with Claude Code

https://claude.ai/code/session_01XGArZMnLPC1G7YP62vN18V

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 97.83%. Comparing base (57fd985) to head (3487e78).
⚠️ Report is 1 commits behind head on main.

Additional details and impacted files

@@            Coverage Diff             @@
##             main    #1000      +/-   ##
==========================================
+ Coverage   97.78%   97.83%   +0.05%     
==========================================
  Files           9        9              
  Lines        1808     1850      +42     
==========================================
+ Hits         1768     1810      +42     
  Misses         40       40              

Flag	Coverage Δ
dask	97.01% <100.00%> (+0.06%)	⬆️
jax	97.18% <100.00%> (+0.06%)	⬆️
numpy	97.72% <100.00%> (+0.05%)	⬆️

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.

[Copilot]

Pull request overview

Adds an Array-API-native MAD fallback for sigma_func, avoiding NumPy conversion on non-NumPy backends.

Changes:

• Implements _mad_fallback with axis, mask, NaN, dtype, namespace, and device handling.
• Expands MAD tests across backends and fixes strict-array scaling setup.
• Updates Array API documentation, changelog, and escape baseline.

Reviewed changes

Copilot reviewed 6 out of 6 changed files in this pull request and generated 2 comments.

Show a summary per file

File	Description
ccdproc/core.py	Adds and integrates the MAD fallback.
ccdproc/tests/test_ccdproc.py	Adds fallback and public API tests.
ccdproc/tests/test_combiner.py	Uses xp.stack for backend compatibility.
ccdproc/tests/array_escape_baseline.txt	Removes the resolved NumPy escape.
docs/array_api.rst	Documents reduction fallbacks and limitations.
CHANGES.rst	Records the new behavior.

Suppressed comments (1)

ccdproc/core.py:1539

• Do not override the caller's ignore_nan choice merely because the CCDData has a mask. Astropy's explicit MaskedArray path already excludes masked pixels and only additionally masks ordinary NaNs when ignore_nan=True; forcing this flag makes unmasked NaNs get ignored even when the caller requested propagation.

            ignore_nan = True

---

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

[mwcraig]

Verdict: needs small changes — one low-severity bug, otherwise merge-ready in substance.

What I verified: the numpy path is unchanged for plain arrays and unmasked CCDData (checked against the old np.asarray(median_absolute_deviation(...) * 1.4826...) formula over the axis × ignore_nan grid). _mad_fallback matches astropy.stats.median_absolute_deviation on 300 seeded random cases per backend (ndim 1–4, NaN fractions 0/0.1/0.5, every axis form) on numpy, array-api-strict on device1, jax (X64) and dask, plus the mask variants (all-True, all-False, fully masked column, int mask, numpy mask with foreign data). On dask the result is a lazy dask.array (no eager compute); on jax it traces under jax.jit; float32 is preserved and the result stays on the input device. The removal of the core.py sigma_func numpy.asanyarray baseline line is justified: an enforce run on dask over all sigma_func callers logs no sigma_func site.

Two notes that do not anchor to a diff line:

• The PR description says the numpy path keeps the "same masked-CCDData handling". It does not — the numpy masked path changed deliberately (np.ma.masked_array + forced ignore_nan), and that is a fix (old behaviour for a masked CCDData with a NaN at axis=0, ignore_nan=False was [0., nan]; new is [1.48, 2.97]). The second #1000 CHANGES bullet describes it correctly; the PR text should not claim the numpy path is untouched.
• Possible follow-up rather than for this PR: the None/tuple axis support (permute the reduced axes last + merge) belongs in _nanfuncs._setup, where it would apply to all five fallbacks at once. _mad_fallback would shrink to ~6 lines and _median_fallback/the combiner fallbacks would gain tuple axes for free.

Merge with #1001: only textual conflicts, in CHANGES.rst and docs/array_api.rst (adjacent bullet insertions); no code overlap (#1001 does not touch core.py, this PR does not touch combiner.py; the baseline file auto-merges).

Test runs on this branch (test_ccdproc.py + test_combiner.py):

• numpy: 231 passed, 6 skipped
• array-api-strict (device1): 219 passed, 4 failed (test_combiner_sigmaclip_{high,low,single_pix} → #1001; test_unit_mismatch_behaves_as_expected → #936), 2 skipped, 12 xfailed, 0 xpassed
• jax (X64): 235 passed, 2 skipped
• dask: 235 passed, 2 skipped

[mwcraig]

Bug (low severity): an entirely masked CCDData with axis=None returns 0.0 here, not NaN. astropy's median_absolute_deviation does result.item() for masked input when axis is None, and np.ma.masked.item() is 0.0, so by the time np.ma.isMaskedArray(result) runs there is nothing left to .filled(np.nan). This contradicts the docstring above ("a slice that is entirely masked gives NaN"), the same call at axis=0 (NaN) and the fallback on every other backend (NaN).

import numpy as np; from astropy.nddata import CCDData
from ccdproc.core import sigma_func, _mad_fallback
d = np.arange(6.).reshape(3, 2); m = np.ones((3, 2), bool)
sigma_func(CCDData(d, unit="adu", mask=m))            # array(0.)
sigma_func(CCDData(d, unit="adu", mask=m), axis=0)    # array([nan, nan])
_mad_fallback(d, None, False, mask=m)                 # array(nan)

The simplification suggested on the np.ma.masked_array line fixes this as a side effect.

[mwcraig]

Reproduced, including the axis=0 vs axis=None asymmetry. The .filled(np.nan) guard never sees the axis=None case because astropy has already turned np.ma.masked into 0.0 with .item().

Will do: fix it via the mask→NaN substitution in sigma_func proposed on the np.ma.masked_array line, so numpy goes through astropy's nanmedian path and an entirely masked input gives NaN for every axis, and drop the np.ma.isMaskedArray special case. test_sigma_func_ccddata_mask_is_honoured gets an all-masked axis=None case so this cannot come back.

— Written by Claude at @mwcraig's direction.

[mwcraig]

Simplification (≈ −15 lines, also fixes the axis=None bug): do the mask→NaN substitution once in sigma_func, before the numpy/fallback branch, and drop the mask handling from both branches:

if mask is not None:
    nan = xp.asarray(xp.nan, dtype=data.dtype, device=device)
    data = xp.where(xp.asarray(mask, device=device), nan, data)
    ignore_nan = True

Then _mad_fallback loses its mask= parameter (and the xp.where block plus its docstring), and this branch loses the np.ma.masked_array construction, the np.ma.isMaskedArray/.filled(np.nan) special case and the comment explaining astropy's mask quirks. On numpy astropy then takes its nanmedian (bottleneck when installed) path instead of np.ma.median, which is markedly slower, and an all-NaN slice at axis=None gives NaN. The two branches become structurally identical, and the direct _mad_fallback(..., mask=mask) call in test_sigma_func_ccddata_mask_is_honoured goes away. (Integer CCDData data promotes to float in the where, which astropy would have done anyway.)

[mwcraig]

Will do, with two adjustments found while checking it:

1. The where needs the promotion first. On array-api-strict, xp.where(mask, float64_nan, int64_data) raises TypeError: array_api_strict.float64 and array_api_strict.int64 cannot be type promoted together (numpy promotes silently, so an integer CCDData would pass on numpy and fail elsewhere). So the promotion helper from the other thread runs in sigma_func when mask is not None, before the where; unmasked numpy input still reaches astropy untouched.
2. Warnings: the numpy masked path moves from np.ma.median (silent on an entirely masked slice) to np.nanmedian, which emits RuntimeWarning: All-NaN slice encountered, once per such slice. That is numpy's normal behaviour for sigma_func(arr, ignore_nan=True) on NaN-filled data, so it is acceptable, but the comment above the filterwarnings mark on test_sigma_func_ccddata_mask_is_honoured ("the masked path itself does not warn") becomes false and will be rewritten.

Along with that: _mad_fallback loses mask=, the numpy branch loses the np.ma construction and the .filled(np.nan) special case, the direct _mad_fallback(..., mask=mask) call in the test goes, and the test gains an all-masked axis=None case for the 0.0 bug.

— Written by Claude at @mwcraig's direction.

[mwcraig]

Simplification (≈ −10 lines): this is numpy's normalize_axis_tuple, which is pure-Python integer logic (touches no arrays) and numpy is a hard dependency:

axes = list(normalize_axis_tuple(axis, ndim))

It raises AxisError (a ValueError subclass) for out-of-bounds and ValueError("repeated axis") for duplicates, accepts numpy ints, and accepts lists — which closes a small inconsistency: sigma_func(a, axis=[0, 1]) works on numpy (astropy accepts it) but raises NotImplementedError from _nanfuncs._setup on every other backend because only isinstance(axis, tuple) is tested here. With it, an int axis can be treated as a 1-tuple and the elif split disappears (permute/reshape of a single axis is a no-op view), and the operator import goes. Caveats: with the numpy>=1.26 pin the import location differs (numpy.lib.array_utils on 2.x, numpy.core.numeric on 1.26), so a 3-line try/except import; and normalize_axis_tuple(True, n) silently gives (1,), so keep the bool rejection if that matters.

[mwcraig]

Will do. Checked the details:

• Import: numpy.lib.array_utils on 2.x; numpy.core.numeric still resolves on numpy 2.5 but with a DeprecationWarning (which pytest's filterwarnings = error would turn into a failure), so it is try: from numpy.lib.array_utils import normalize_axis_tuple / except ImportError: from numpy.core.numeric import ... for 1.26.
• AxisError is a ValueError (and IndexError) subclass; (0, -3) gives ValueError: repeated axis. test_mad_fallback_rejects_duplicate_axes keeps its match= strings adjusted.
• The inconsistency is real: sigma_func(a, axis=[0, 1]) works on numpy and raises NotImplementedError on strict, jax and dask; axis=True is silently axis 1 on numpy (astropy accepts it) and raises on the others. The bool rejection stays on the fallback, as a TypeError.

One quibble with the comment: with a single int axis the permute_dims is a transpose, not a no-op, unless that axis is already last. It is a view on numpy and lazy on dask, and the sort has to move data along that axis regardless, so it is not worth a special case; the elif split and the operator import go.

— Written by Claude at @mwcraig's direction.

[mwcraig]

Let's assume that the next ccdproc release (or at least the one with the array api stuff) will require numpy 2, so we can drop the 1.26 stuff. Open an issue to remind us to bump the requirement

[mwcraig]

These three lines duplicate _nanfuncs._setup (_nanfuncs.py:87-92) verbatim. Worth a tiny _promote_to_real(x, xp, device) helper in _nanfuncs used by both (or xpx.default_dtype, if the pinned array_api_extra>=0.7.0 has it — the installed 0.11.1 does).

[mwcraig]

Will do: a _promote_to_real(x, xp, device) in _nanfuncs, used by _setup, _mad_fallback, and the mask substitution that moves into sigma_func (which needs the promotion before its where on strict, see that thread). xpx.default_dtype only arrived in array-api-extra 0.8.0 (checked the 0.7.0 and 0.8.0 wheels: absent in 0.7.0), so using it means bumping the >=0.7.0 pin; the helper avoids that. Say if you would rather bump.

— Written by Claude at @mwcraig's direction.

[mwcraig]

Question: _median_fallback a few functions up and combiner._default_median both try the namespace's native median/nanmedian first and only fall back to the sort-based _nanfuncs versions when there is none; this always sorts, even on jax/dask/cupy/torch, which all have nanmedian (O(n), fused under jit on jax). Two full sorts of the stack per call. Is that intentional — to keep all-NaN slices silent, which native nanmedian on numpy-backed dask would not be? If so, a one-line comment here saying so would stop the next reader from "fixing" it; if not, native-first would match the rest of the module and the docs would no longer need to advertise the sort cost for every non-numpy library.

[mwcraig]

Not intentional in the sense of a silence guarantee: the sort path was chosen to have one code path on every backend, and the silence was documented after the fact. Checked what native-first would actually do:

• dask (array_api_compat.dask.array) and jax both have median and nanmedian; array-api-strict has neither, so it would still fall back. dask's need an integer axis (axis=None raises NotImplementedError there), which is exactly what _mad_fallback passes after its None/tuple handling, so the merge-axes step stays.
• The silence buys nothing on dask: Combiner.median_combine on dask already uses native xp.nanmedian through _default_median, which emits All-NaN slice encountered at compute time for the very same slices; jax's nanmedian is silent.

Will do: try xp.nanmedian/xp.median and fall back on AttributeError, mirroring _median_fallback; the docstring and docs/array_api.rst then claim the sort cost only for libraries without a native version; test_mad_fallback_all_nan_slice_is_silent goes, since silence becomes a property of the _nanfuncs medians (pinned in test_nanfuncs.py) rather than of _mad_fallback.

— Written by Claude at @mwcraig's direction.

[mwcraig]

Nit: a list or bool axis raises NotImplementedError (from _nanfuncs._setup), which the Raises section does not mention. Moot if normalize_axis_tuple is adopted for the tuple path.

[mwcraig]

Agree. Will do: with normalize_axis_tuple handling the int/tuple/list forms (see the thread on the axis loop) a list axis works, a bool one is rejected explicitly, and the Raises section will list exactly what is raised (ValueError for out-of-bounds/duplicate, TypeError for a bool or non-integer entry) instead of leaking NotImplementedError from _nanfuncs._setup.

— Written by Claude at @mwcraig's direction.

[mwcraig]

Nit: currently false on numpy for axis=None (see the comment on the .filled(np.nan) line below) — fix the code or the sentence.

[mwcraig]

Will do: the code side. After the mask→NaN substitution moves into sigma_func, the numpy branch computes an entirely masked axis=None input through np.nanmedian and returns NaN, so the sentence becomes true as written.

— Written by Claude at @mwcraig's direction.

[mwcraig]

Nit: "never when bottleneck is installed" is only true for float64 data — astropy's _DtypeDispatch sends only f8 arrays to bottleneck, so a masked float32 CCDData still went through numpy's small-array path. Suggest "not for float64 data when bottleneck is installed". The rest of the sentence (fewer than 600 elements along a single integer axis, ignore_nan=True) checks out against numpy's _nanmedian/_nanmedian_small and CCDData.__array__ returning a MaskedArray.

[mwcraig]

Agree. Checked the installed astropy (8.0.1): stats.nanfunctions._DtypeDispatch.__call__ routes to bottleneck only when dt.kind == "f" and dt.itemsize == 8, and median_absolute_deviation reaches it only on the ignore_nan=True path. One nuance: the f8-only dispatch arrived in astropy 7.0; astropy 6.0.x (still allowed by the >=6.0.1 pin) hands every dtype to bottleneck, so the original sentence was right there and wrong on anything current.

Will do: reword to "…and, when bottleneck is installed, not for float64 data at all", which is correct on every supported astropy.

— Written by Claude at @mwcraig's direction.

[mwcraig]

Assume the minimum Astro is at least 7 for the array api release

[mwcraig]

Test LOC: the new block is 179 lines; about 70 of them (≈ 38%) re-cover paths already pinned elsewhere, mostly in test_nanfuncs.py. What each test uniquely covers:

Test	LOC	Unique path	Redundant with
_MAD_* fixtures + _MAD_CASES (23 cases)	35	data	axis -1, np_int64(1), the 6 odd/even lengths, both all-NaN cases and the two 1-D cases are covered by test_nanfuncs.py::test_matches_numpy (_DATA, lines 36-42); (0,1), (0,2), (2,0) add nothing over (-1,0); bool adds nothing over int; the float32 case only asserts isdtype(..., "real floating"), so it tests nothing float32-specific
test_mad_fallback_matches_astropy (46 params)	32	axis None / int / tuple; promotion; median vs nanmedian	—
test_mad_fallback_all_nan_slice_is_silent	10	nothing new: the silence is _nanfuncs' property, pinned by test_nanfuncs.py::test_no_warning_on_all_nan_slice for both medians; abs(nan - nan) cannot warn	drop
test_mad_fallback_rejects_duplicate_axes	13	tuple duplicate (incl. negative alias), tuple out-of-bounds	int out-of-bounds is _setup's, covered by test_nanfuncs.py::test_bad_axis
test_sigma_func_matches_mad_std (3)	24	public wrapper × constant; 0-d float()	—
test_sigma_func_keeps_namespace_and_device	11	namespace/device of result	same computation as above; two asserts
test_sigma_func_fallback_branch_on_any_backend	16	numpy coverage of the one return _mad_fallback(...) * c line	same test with a monkeypatch
test_sigma_func_ccddata_mask_is_honoured (4)	32	numpy masked path (axis None / axis given); fallback mask path	(1, False) and ((0,1), True) hit no branch not hit by (0, True) / (None, False)

Consolidated sketch (~110 lines including imports, ≈ −70; −75 if the mask substitution moves into sigma_func):

_MAD_3D = ...  # as now
_MAD_CASES = [
    *[(_MAD_3D, ax) for ax in (None, 0, (-1, 0), (0, 1, 2))],  # flatten, int, neg+unsorted tuple, all axes
    (np_array([[1, 4], [2, 3], [5, 6], [4, 1]]), 0),            # promotion
    (<float32 2x3 array>, 1),                                   # dtype preserved
]

@pytest.mark.filterwarnings("ignore:All-NaN slice encountered:RuntimeWarning")
@pytest.mark.parametrize("ignore_nan", [False, True])
@pytest.mark.parametrize(("data", "axis"), _MAD_CASES)
def test_mad_fallback_matches_astropy(data, axis, ignore_nan):
    expected = xp.asarray(np_asarray(median_absolute_deviation(data, axis=axis, ignore_nan=ignore_nan)), device=xp_device)
    result = _mad_fallback(xp.asarray(data, device=xp_device), axis, ignore_nan)
    assert result.shape == expected.shape
    assert result.dtype == expected.dtype   # promotion to float64 AND float32 preservation in one line
    assert xp.all(xpx.isclose(result, expected, equal_nan=True))

@pytest.mark.parametrize(("axis", "match"), [((0, 0), "duplicate"), ((0, -3), "duplicate"), ((0, 3), "out of bounds")])
def test_mad_fallback_rejects_bad_axis_tuple(axis, match): ...

@pytest.mark.parametrize("force_fallback", [False, True])
@pytest.mark.parametrize(("data", "axis", "ignore_nan"), [...the same 3 params...])
def test_sigma_func_matches_mad_std(data, axis, ignore_nan, force_fallback, monkeypatch):
    if force_fallback:
        monkeypatch.setattr(array_api_compat, "is_numpy_namespace", lambda _xp: False)
    ...existing body...
    assert array_api_compat.array_namespace(result) is array_api_compat.array_namespace(data)
    if xp_device is not None:
        assert array_api_compat.device(result) == xp_device

@pytest.mark.filterwarnings("ignore:All-NaN slice encountered:RuntimeWarning")
@pytest.mark.parametrize(("axis", "ignore_nan"), [(0, True), (None, False)])
def test_sigma_func_ccddata_mask_is_honoured(axis, ignore_nan): ...  # as now, two params fewer

The negative-tuple normalisation of the astropy reference only matters when bottleneck is installed, which it is not in any tox env; keeping the two lines is fine but the comment could say so.

[mwcraig]

Mostly agree; will do the consolidation along the lines sketched: one parametrised test_mad_fallback_matches_astropy over the six cases with result.dtype == expected.dtype (which does pin float32 preservation, unlike the current isdtype check), the bad-axis test reduced to the tuple cases (the int case becomes numpy's own AxisError once normalize_axis_tuple is in), the namespace/device asserts and the force_fallback monkeypatch folded into test_sigma_func_matches_mad_std, and the mask test down to (0, True) and (None, False) plus an all-masked axis=None case for the 0.0 bug.

Two corrections:

• bottleneck is in CI: ci_tests.yml runs py312-test-alldeps-numpy126-bottleneck and py312-alldeps-dask-bottleneck (the tox bottleneck factor). _MAD_3D is float64, so on those jobs astropy's _DtypeDispatch sends the reference through _apply_bottleneck → _move_tuple_axes_last, which builds other_axes with i not in axis and then calls np.transpose(array, (1, 2, -1, 0)) for (-1, 0), a four-entry order for a 3-D array (ValueError: axes don't match array). The normalisation and its comment stay as they are; the comment already says why.
• test_mad_fallback_all_nan_slice_is_silent goes, but for the reason in the native-first thread: once _mad_fallback prefers a native nanmedian, silence on dask is no longer its property at all, so the test would start failing there rather than merely being redundant.

— Written by Claude at @mwcraig's direction.

[mwcraig]

On the two notes in the review body of review 5061397724 that have no diff line:

• PR description. Agreed, "same masked-CCDData handling" is wrong; the numpy masked path changed on purpose (and changes again with the mask→NaN substitution being adopted from the inline threads). The description will be edited to say that the numpy path is unchanged for plain arrays and unmasked CCDData, and that a masked CCDData now has its mask honoured on every backend, as the second #1000 CHANGES bullet says.
• None/tuple axis support in _nanfuncs._setup. Agreed as a follow-up, not for this PR: with normalize_axis_tuple adopted here first, moving the permute-last-and-merge step into _setup is a small change that gives _median_fallback and the combiner fallbacks tuple axes for free and shrinks _mad_fallback to a few lines. To be filed as an issue once this PR is in.
• Merge order. Noted: Array-API fallback for Combiner.sigma_clipping when the namespace is not numpy #1001 rebases after this one; the CHANGES.rst/docs/array_api.rst conflicts are adjacent bullets only.

The Copilot review's suppressed comment on line 1539 (do not force ignore_nan for a masked CCDData) is the same point as its inline comment on line 356 and is answered there.

— Written by Claude at @mwcraig's direction.

[mwcraig]

American spelling throughout please

[mwcraig]

Will do. Grepping the added lines of the whole diff (356 lines, five files) for the usual -our/-ise/-yse/-re forms finds four, all fixed in the next push:

• ccdproc/core.py:1533 — honours → honors (this comment)
• ccdproc/tests/test_ccdproc.py:451 — normalised → normalized
• ccdproc/tests/test_ccdproc.py:554 — test_sigma_func_ccddata_mask_is_honoured → ..._honored
• CHANGES.rst:130 — honoured → honored

Nothing in docs/array_api.rst or test_combiner.py. The subject of b8b0e9f ("Honour ... normalise ...") has the same two words and gets reworded when the branch is rewritten for the test consolidation.

— Written by Claude at @mwcraig's direction.

[mwcraig]

American spelling please

[mwcraig]

Renamed to test_sigma_func_ccddata_mask_is_honored; the other three in the diff (honours at core.py:1533, normalised at test_ccdproc.py:451, honoured at CHANGES.rst:130) are listed in the core.py thread.

— Written by Claude at @mwcraig's direction.

[mwcraig]

Couldn't we add one of the backends to coverage instead of this?

[mwcraig]

Two of them already are, and this test's premise is wrong. Since c2ce2e2 (2026-08-23, before this PR was opened) ci_tests.yml runs py313-jax-cov and py312-alldeps-dask-enforce-cov and uploads both with flags: jax / flags: dask alongside the numpy upload; there is no codecov.yml, so codecov merges the three uploads per commit (that merged report is what codecov/patch is green on here). Only the strict job is uncovered: it sits in the continue-on-error job without a -cov factor, and adding one would not help while it fails by design, because the upload step runs only after a successful tox.

On jax and dask the real branch, return _mad_fallback(...) * 1.482602218505602, is executed by test_sigma_func_matches_mad_std, by test_sigma_func_ccddata_mask_is_honoured (including the mask= path) and by every median_combine test in test_combiner.py, so the monkeypatch buys nothing. Changes:

• this test goes, and with it the force_fallback parameter promised in the thread at line 405 (superseded);
• the two "none of which report coverage" comments (test_mad_fallback_matches_astropy and the tail of the mask test) go, and so does the trailing direct _mad_fallback(ccd.data, ..., mask=mask) call in the mask test;
• the direct-call test_mad_fallback_matches_astropy stays, because it is the differential check against astropy over the axis forms, not a coverage device; its comment will say that;
• the PR description gets the same correction (it says the strict job uploads no coverage, which is true but was the wrong reason).

No CI change needed.

— Written by Claude at @mwcraig's direction.

[mwcraig]

Every test needs a comment explaining what it does...

[mwcraig]

Will do. The file uses a leading # comment rather than a docstring (0 docstrings across its 58 tests), so each test in the block gets a one-line comment saying what it checks and, where it is not obvious, why. For the four tests that survive the consolidation:

• test_mad_fallback_matches_astropy — _mad_fallback agrees with astropy.stats.median_absolute_deviation on every backend for each axis form sigma_func's callers use, and handles dtype like the _nanfuncs medians (int promoted, float32 kept).
• test_mad_fallback_rejects_bad_axis_tuple — duplicate (including a negative alias) and out-of-bounds entries in a tuple axis raise ValueError instead of silently reducing the wrong axes.
• test_sigma_func_matches_mad_std — the public entry point matches astropy.stats.mad_std, stays in the input's namespace and device, and for axis=None gives a 0-d result that converts to float (what background_deviation_box relies on).
• test_sigma_func_ccddata_mask_is_honored — masked CCDData pixels are excluded on every backend (numpy via numpy.ma, the rest via the fallback), equivalent to NaN-filling with ignore_nan on; an all-masked slice gives NaN, not 0.0.

The existing multi-line comments that explain a non-obvious step (the bottleneck axis normalization, the bool() forcing on lazy backends) stay below those.

— Written by Claude at @mwcraig's direction.

[mwcraig]

Can this be simplified?

[mwcraig]

Yes; this is the block the consolidation agreed in the thread at line 405 shrinks. After it:

_MAD_RNG = np_random.default_rng(929)
_MAD_3D = _MAD_RNG.normal(size=(5, 4, 3))
_MAD_3D[[0, 1, 2, 4], [1, 2, 0, 3], [0, 2, 1, 1]] = np_nan  # scattered NaNs

_MAD_CASES = [
    # axis=None (background_deviation_box), an int (median_combine), a tuple
    # with a negative, unsorted entry, and every axis.
    *[(_MAD_3D, axis) for axis in (None, 0, (-1, 0), (0, 1, 2))],
    (np_array([[1, 4], [2, 3], [5, 6], [4, 1]]), 0),  # int -> default float
    (np_array([[1.0, 2.0, 3.0], [4.0, 5.0, 7.0]], dtype=np_float32), 1),  # float32 kept
]

35 lines → 11, 23 cases → 6, with _MAD_CLEAN_3D gone: the bad-axis test raises before it computes anything so it can use _MAD_3D, and the one ignore_nan=False case of test_sigma_func_matches_mad_std uses np_nan_to_num(_MAD_3D). The dropped cases are covered in test_nanfuncs.py (-1, np_int64, the six odd/even lengths, both all-NaN inputs, the 1-D inputs) or add nothing over a kept case ((0, 1)/(0, 2)/(2, 0) vs (-1, 0), bool vs int), and result.dtype == expected.dtype in the test replaces the isdtype check so the float32 case pins something.

— Written by Claude at @mwcraig's direction.