From fec5ed2106b8669c1afb9c8b15621a1fdfddb1d3 Mon Sep 17 00:00:00 2001
From: Matt Craig <mattwcraig@gmail.com>
Date: Wed, 26 Aug 2026 20:30:43 -0500
Subject: [PATCH 01/11] Add array-API fallback for sigma_func

sigma_func wraps astropy.stats.median_absolute_deviation, which is
numpy-only, so on every other array namespace the data was converted to
numpy (and on array-api-strict the call failed outright). This is one of
the two remaining astropy.stats call sites behind #929.

Keep the numpy path exactly as it was -- same astropy call, same masked
CCDData handling, no new work before the is_numpy_namespace check -- and
route every other namespace 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 from _nanfuncs
(nanmedian for ignore_nan=True, median otherwise). The fallback promotes
integer and boolean input to the namespace's default real floating
dtype, flattens for axis=None (background_deviation_box), reduces over
tuples of axes by permuting them last and merging them, validates axes,
and excludes the masked pixels of a CCDData. sigma_func stays the same
function object because median_combine tests uncertainty_func identity.

The new tests exercise the fallback directly on every backend (the
strict job uploads no coverage) against astropy over the axis and
ignore_nan grid, including int/bool/float32 input and all-NaN slices,
and check the public entry point against astropy.stats.mad_std, the
namespace/device of the result, the fallback branch on numpy via a
patched is_numpy_namespace, and the CCDData mask. In test_combiner.py,
test_combiner_with_scaling now builds its reference stack with xp.stack
instead of xp.asarray on a tuple of arrays, which array-api-strict
rejects; the test previously failed earlier, in sigma_func, so this
never surfaced.

Verified: numpy 616 passed; array-api-strict 5 failed (the three
Combiner.sigma_clipping tests, #936, #983 -- down from 13), 0 xpassed;
jax 615 passed; dask 609 passed.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01XGArZMnLPC1G7YP62vN18V
---
 ccdproc/core.py                | 147 ++++++++++++++++++++++++++++--
 ccdproc/tests/test_ccdproc.py  | 161 +++++++++++++++++++++++++++++++++
 ccdproc/tests/test_combiner.py |   2 +-
 3 files changed, 300 insertions(+), 10 deletions(-)

diff --git a/ccdproc/core.py b/ccdproc/core.py
index 9b70e154..951bb2c3 100644
--- a/ccdproc/core.py
+++ b/ccdproc/core.py
@@ -5,6 +5,7 @@
 import logging
 import math
 import numbers
+import operator
 import warnings
 
 import array_api_compat
@@ -30,6 +31,7 @@
     _wrap_ccddata_for_array_api,
 )
 from ._nanfuncs import median as _nanfuncs_median
+from ._nanfuncs import nanmedian as _nanfuncs_nanmedian
 from .log_meta import log_to_metadata
 from .utils.slices import slice_from_string
 
@@ -286,6 +288,103 @@ def _median_fallback(array, axis, xp=None):
         return _nanfuncs_median(array, axis=axis, xp=xp)
 
 
+def _mad_fallback(data, axis, ignore_nan, xp=None, mask=None):
+    """
+    Median absolute deviation written purely in terms of the array API.
+
+    This is the non-numpy branch of `sigma_func`;
+    `astropy.stats.median_absolute_deviation` is numpy-only, so any other
+    namespace computes the deviation here, on the device of ``data``.
+
+    Parameters
+    ----------
+    data : array
+        Array whose median absolute deviation is to be calculated. Integer
+        and boolean input is promoted to the namespace's default real
+        floating dtype.
+
+    axis : int, tuple of int or None
+        Axis or axes along which the deviation is computed. ``None``
+        flattens ``data`` first; a tuple reduces over all the listed axes.
+        Negative values count from the last axis.
+
+    ignore_nan : bool
+        If `True`, NaNs are ignored; otherwise a NaN anywhere in a reduced
+        slice makes that slice's result NaN, as `numpy.median` does.
+
+    xp : array namespace, optional
+        Namespace to use. If not provided, it is determined from ``data``.
+
+    mask : array or None, optional
+        Boolean mask with the shape of ``data``. Masked pixels are excluded
+        from the statistics (which forces ``ignore_nan`` on), as astropy
+        does for a masked `~astropy.nddata.CCDData` when ``ignore_nan`` is
+        set.
+
+    Returns
+    -------
+    mad : array
+        Median absolute deviation of ``data`` along ``axis``, with the
+        reduced axes removed (0-d when ``axis`` is `None`), in the namespace
+        and on the device of ``data``.
+
+    Raises
+    ------
+    ValueError
+        If ``axis`` is out of bounds or lists an axis twice.
+
+    Notes
+    -----
+    Both medians are the sort-based fallbacks from `ccdproc._nanfuncs`, so
+    the cost is O(n log n) along the reduced axes rather than the O(n) of
+    `numpy.median`. Slices that are entirely NaN yield NaN without the
+    ``RuntimeWarning`` numpy emits.
+    """
+    xp = xp or array_api_compat.array_namespace(data)
+    device = array_api_compat.device(data)
+
+    # The _nanfuncs medians promote internally, but ``data - center`` below
+    # needs the promoted input too: subtracting a float from an integer array
+    # is not allowed by the standard and raises on array-api-strict.
+    if not xp.isdtype(data.dtype, "real floating"):
+        info = xp.__array_namespace_info__()
+        data = xp.astype(data, info.default_dtypes(device=device)["real floating"])
+
+    if mask is not None:
+        nan = xp.asarray(xp.nan, dtype=data.dtype, device=device)
+        mask = xp.astype(xp.asarray(mask, device=device), xp.bool)
+        data = xp.where(mask, nan, data)
+        ignore_nan = True
+
+    if axis is None:
+        data = xp.reshape(data, (-1,))
+        axis = 0
+    elif isinstance(axis, tuple):
+        ndim = data.ndim
+        axes = []
+        for ax in axis:
+            ax = operator.index(ax)
+            if not -ndim <= ax < ndim:
+                raise ValueError(
+                    f"axis {ax} is out of bounds for array of dimension {ndim}"
+                )
+            axes.append(ax % ndim)
+        if len(set(axes)) != len(axes):
+            raise ValueError(f"duplicate value in 'axis': {axis}")
+        # Move the reduced axes to the end and merge them into one, so that
+        # the single-axis medians below reduce over all of them at once.
+        kept = [ax for ax in range(ndim) if ax not in axes]
+        data = xp.permute_dims(data, tuple(kept + axes))
+        data = xp.reshape(
+            data, tuple(data.shape[ax] for ax in range(len(kept))) + (-1,)
+        )
+        axis = -1
+
+    med = _nanfuncs_nanmedian if ignore_nan else _nanfuncs_median
+    center = med(data, axis=axis, xp=xp)
+    return med(xp.abs(data - xp.expand_dims(center, axis=axis)), axis=axis, xp=xp)
+
+
 @log_to_metadata
 def ccd_process(
     ccd,
@@ -1385,7 +1484,7 @@ def sigma_func(arr, axis=None, ignore_nan=False):
 
     Parameters
     ----------
-    arr : `~astropy.nddata.CCDData` or `numpy.ndarray`
+    arr : `~astropy.nddata.CCDData` or array
         Array whose deviation is to be calculated.
 
     axis : int, tuple of ints or None, optional
@@ -1395,20 +1494,50 @@ def sigma_func(arr, axis=None, ignore_nan=False):
         it counts from the last to the first axis.
         Default is ``None``.
 
+    ignore_nan : bool, optional
+        If `True`, NaNs are ignored when computing the medians; otherwise a
+        NaN in a reduced slice makes that slice's result NaN.
+        Default is ``False``.
+
     Returns
     -------
-    uncertainty : float
-        uncertainty of array estimated from median absolute deviation.
+    uncertainty : array
+        Uncertainty of the array estimated from the median absolute
+        deviation, in the array namespace and on the device of the input.
+        It is 0-d when ``axis`` is ``None``.
+
+    Notes
+    -----
+    For numpy input the deviation is computed by
+    `astropy.stats.median_absolute_deviation`, exactly as before. For every
+    other array namespace it is computed by a fallback written purely in
+    terms of the array API standard, on the device of the input. That
+    fallback uses sort-based medians, so it costs O(n log n) along the
+    reduced axes rather than O(n); it promotes integer and boolean input to
+    the namespace's default real floating dtype, and it yields NaN silently
+    for slices that are entirely NaN, where numpy would warn.
+
+    A masked `~astropy.nddata.CCDData` is handed to astropy as is on numpy.
+    The fallback excludes the masked pixels from the statistics (which
+    implies ``ignore_nan``), as astropy does for a single integer ``axis``
+    with ``ignore_nan=True``.
     """
     if isinstance(arr, CCDData):
-        xp = array_api_compat.array_namespace(arr.data)
+        data = arr.data
+        mask = arr.mask
     else:
-        xp = array_api_compat.array_namespace(arr)
+        data = arr
+        mask = None
+    xp = array_api_compat.array_namespace(data)
+
+    if array_api_compat.is_numpy_namespace(xp):
+        # Pass ``arr`` rather than ``data``: astropy honours a CCDData mask.
+        return xp.asarray(
+            stats.median_absolute_deviation(arr, axis=axis, ignore_nan=ignore_nan)
+            * 1.482602218505602
+        )
 
-    return xp.asarray(
-        stats.median_absolute_deviation(arr, axis=axis, ignore_nan=ignore_nan)
-        * 1.482602218505602
-    )
+    return _mad_fallback(data, axis, ignore_nan, xp=xp, mask=mask) * 1.482602218505602
 
 
 def setbox(x, y, mbox, xmax, ymax):
diff --git a/ccdproc/tests/test_ccdproc.py b/ccdproc/tests/test_ccdproc.py
index 8f51a54e..2844ade0 100644
--- a/ccdproc/tests/test_ccdproc.py
+++ b/ccdproc/tests/test_ccdproc.py
@@ -16,17 +16,23 @@
     StdDevUncertainty,
     VarianceUncertainty,
 )
+from astropy.stats import mad_std, median_absolute_deviation
 from astropy.units.quantity import Quantity
 from astropy.utils.exceptions import AstropyUserWarning
 from astropy.wcs import WCS
 from numpy import array as np_array
+from numpy import asarray as np_asarray
+from numpy import float32 as np_float32
+from numpy import int64 as np_int64
 from numpy import mgrid as np_mgrid
+from numpy import nan as np_nan
 from numpy import random as np_random
 
 from ccdproc.conftest import testing_array_device as xp_device
 from ccdproc.conftest import testing_array_library as xp
 from ccdproc.core import (
     Keyword,
+    _mad_fallback,
     _median_fallback,
     ccd_process,
     cosmicray_lacosmic,
@@ -396,6 +402,161 @@ def __getattr__(self, name):
     assert xp.all(xpx.isclose(result, expected))
 
 
+_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
+_MAD_CLEAN_3D = _MAD_RNG.normal(size=(5, 4, 3))
+
+_MAD_CASES = [
+    # Every axis form sigma_func's callers use, on data with scattered NaNs:
+    # None (background_deviation_box), a single axis (median_combine), a
+    # numpy integer, and tuples in either order with negative entries.
+    *[
+        (_MAD_3D, axis)
+        for axis in [
+            None,
+            0,
+            1,
+            -1,
+            np_int64(1),
+            (0, 1),
+            (0, 2),
+            (2, 0),
+            (-1, 0),
+            (0, 1, 2),
+        ]
+    ],
+    *[(_MAD_RNG.normal(size=(n, 7)), 0) for n in range(1, 7)],  # odd/even lengths
+    (np_array([1.0, 2.0, 3.0, 4.0]), 0),  # 1-D
+    (np_array([1.0, 2.0, 3.0, 4.0]), None),
+    (np_array([[1.0, np_nan], [2.0, np_nan], [3.0, np_nan]]), 0),  # all-NaN column
+    (np_array([np_nan, np_nan, np_nan]), None),  # every value NaN
+    (np_array([[1, 4], [2, 3], [5, 6], [4, 1]]), 0),  # integer input
+    (np_array([[True, False], [False, True], [True, True]]), 0),  # boolean input
+    (np_array([[1.0, 2.0, 3.0], [4.0, 5.0, 7.0]], dtype=np_float32), 1),
+]
+
+
+# The ignore mark lets the astropy reference warn on all-NaN slices where
+# the fallback deliberately does not; test_mad_fallback_all_nan_slice_is_silent
+# owns the fallback's silence.
+@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):
+    # sigma_func only reaches _mad_fallback on non-numpy namespaces, none of
+    # which report coverage, so the fallback is exercised directly here on
+    # every backend.
+    expected_np = np_asarray(
+        median_absolute_deviation(data, axis=axis, ignore_nan=ignore_nan)
+    )
+    expected = xp.asarray(expected_np, device=xp_device)
+
+    result = _mad_fallback(xp.asarray(data, device=xp_device), axis, ignore_nan)
+
+    assert result.shape == expected.shape
+    # Integer and boolean input is promoted to the namespace's default real
+    # dtype, as the _nanfuncs medians do.
+    assert xp.isdtype(result.dtype, "real floating")
+    assert xp.all(xpx.isclose(result, expected, equal_nan=True))
+
+
+def test_mad_fallback_all_nan_slice_is_silent():
+    data = xp.asarray(np_array([[np_nan, np_nan]] * 3), device=xp_device)
+    with warnings.catch_warnings():
+        warnings.simplefilter("error")
+        # The bool() calls force the computation inside the block on lazy
+        # backends, where a warning would otherwise surface at compute time.
+        assert bool(xp.all(xp.isnan(_mad_fallback(data, 0, True))))
+        assert bool(xp.isnan(_mad_fallback(data, None, True)))
+
+
+def test_mad_fallback_rejects_duplicate_axes():
+    data = xp.asarray(_MAD_CLEAN_3D, device=xp_device)
+    with pytest.raises(ValueError, match="duplicate"):
+        _mad_fallback(data, (0, 0), True)
+    # -3 is axis 0 of a 3-D array, so this is a duplicate too.
+    with pytest.raises(ValueError, match="duplicate"):
+        _mad_fallback(data, (0, -3), True)
+    with pytest.raises(ValueError, match="out of bounds"):
+        _mad_fallback(data, (0, 3), True)
+    with pytest.raises(ValueError, match="out of bounds"):
+        _mad_fallback(data, 3, True)
+
+
+@pytest.mark.parametrize(
+    ("data", "axis", "ignore_nan"),
+    [
+        pytest.param(_MAD_CLEAN_3D, None, False, id="all-no_nan"),
+        pytest.param(_MAD_3D, 0, True, id="axis0-ignore_nan"),
+        pytest.param(_MAD_3D, (0, 1), True, id="axes01-ignore_nan"),
+    ],
+)
+def test_sigma_func_matches_mad_std(data, axis, ignore_nan):
+    expected_np = np_asarray(mad_std(data, axis=axis, ignore_nan=ignore_nan))
+    expected = xp.asarray(expected_np, device=xp_device)
+
+    result = sigma_func(
+        xp.asarray(data, device=xp_device), axis=axis, ignore_nan=ignore_nan
+    )
+
+    assert result.shape == expected.shape
+    assert xp.all(xpx.isclose(result, expected, equal_nan=True))
+    if axis is None:
+        # background_deviation_box relies on the 0-d result converting to float.
+        assert result.shape == ()
+        assert float(result) == pytest.approx(float(expected_np))
+
+
+def test_sigma_func_keeps_namespace_and_device():
+    data = xp.asarray(_MAD_3D, device=xp_device)
+    for axis in (None, 0, (0, 2)):
+        result = sigma_func(data, axis=axis, ignore_nan=True)
+        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
+
+
+def test_sigma_func_fallback_branch_on_any_backend(monkeypatch):
+    # On numpy sigma_func hands the data to astropy; every other backend
+    # takes the fallback branch, and none of them reports coverage, so make
+    # numpy take it too and check that the two branches agree.
+    data = xp.asarray(_MAD_3D, device=xp_device)
+    expected = sigma_func(data, axis=0, ignore_nan=True)
+
+    monkeypatch.setattr(array_api_compat, "is_numpy_namespace", lambda _xp: False)
+    result = sigma_func(data, axis=0, ignore_nan=True)
+
+    assert array_api_compat.array_namespace(result) is array_api_compat.array_namespace(
+        data
+    )
+    assert xp.all(xpx.isclose(result, expected, equal_nan=True))
+
+
+def test_sigma_func_ccddata_mask_is_honoured():
+    ccd = ccd_data_func()
+    mask = xp.asarray(RNG(929).random(ccd.shape) > 0.7, device=xp_device)
+    # TODO: Set .mask instead of ._mask when CCDData is array-api compliant
+    ccd._mask = mask
+    nan = xp.asarray(np_nan, dtype=ccd.data.dtype, device=xp_device)
+    nanned = xp.where(mask, nan, ccd.data)
+
+    # A single integer axis with ignore_nan=True is the form median_combine
+    # uses, and the one form for which astropy (the numpy path) also excludes
+    # the masked pixels, so the public function can be checked on every backend.
+    result = sigma_func(ccd, axis=0, ignore_nan=True)
+    expected = sigma_func(nanned, axis=0, ignore_nan=True)
+    assert xp.all(xpx.isclose(result, expected, equal_nan=True))
+
+    # The fallback excludes the masked pixels for every axis and ignore_nan.
+    for axis, ignore_nan in [(None, False), (1, False), ((0, 1), True)]:
+        result = _mad_fallback(ccd.data, axis, ignore_nan, mask=mask)
+        expected = _mad_fallback(nanned, axis, True)
+        assert xp.all(xpx.isclose(result, expected, equal_nan=True))
+
+
 def test_trim_image_fits_section_requires_string():
     ccd_data = ccd_data_func()
     with pytest.raises(TypeError):
diff --git a/ccdproc/tests/test_combiner.py b/ccdproc/tests/test_combiner.py
index 5ff00b6b..002b20ce 100644
--- a/ccdproc/tests/test_combiner.py
+++ b/ccdproc/tests/test_combiner.py
@@ -720,7 +720,7 @@ def test_combiner_with_scaling():
     # Scale by a float
     avg_ccd = combiner.average_combine(scale_to=2.0)
     expected_avg = 2 * xp.mean(
-        xp.asarray((ccd_data.data, ccd_data_lower.data, ccd_data_higher.data))
+        xp.stack((ccd_data.data, ccd_data_lower.data, ccd_data_higher.data))
     )
     assert xp.all(xpx.isclose(xp.mean(avg_ccd.data), expected_avg))
     assert avg_ccd.shape == ccd_data.shape

From dfc5d88432d81ed108128cbe2d270b377d0a5861 Mon Sep 17 00:00:00 2001
From: Matt Craig <mattwcraig@gmail.com>
Date: Wed, 26 Aug 2026 20:30:45 -0500
Subject: [PATCH 02/11] Document the sigma_func fallback and refresh the
 array-API limitations

Add the changelog entry (PR number to be filled in) and bring the
"What limitations should I be aware of?" list in docs/array_api.rst up
to date: it only mentioned the nanmedian fallback, but the combiner also
falls back for nansum/nanmean/nanstd, subtract_overscan for median, and
sigma_func now for the median absolute deviation.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01XGArZMnLPC1G7YP62vN18V
---
 CHANGES.rst        |  5 +++++
 docs/array_api.rst | 30 ++++++++++++++++++++++--------
 2 files changed, 27 insertions(+), 8 deletions(-)

diff --git a/CHANGES.rst b/CHANGES.rst
index 546cb592..30502d3c 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -121,6 +121,11 @@ Bug Fixes
   dtype instead of passing it through unchanged. [#999]
 - Size the overscan model fit in ``subtract_overscan`` with ``shape``
   instead of ``len``, which the array API standard does not provide. [#999]
+- Compute ``sigma_func`` (the ``median_combine`` uncertainty) with a median
+  absolute deviation written purely in terms of the array API when the array
+  namespace is not numpy, instead of converting the data to numpy through
+  ``astropy.stats.median_absolute_deviation``. Numpy input still uses
+  astropy. [#TBD]
 
 2.5.1 (2025-07-05)
 ------------------
diff --git a/docs/array_api.rst b/docs/array_api.rst
index 63e8bdac..1c7d2d62 100644
--- a/docs/array_api.rst
+++ b/docs/array_api.rst
@@ -149,14 +149,28 @@ effect only after they are merged.
 What limitations should I be aware of?
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-+ The ``median`` function is not part of the array API, but most array libraries
-  do provide a NaN-aware ``nanmedian``. When combining images, `ccdproc`_ uses
-  ``nanmedian`` from `bottleneck`_ for `numpy`_ arrays if `bottleneck`_ is
-  installed; otherwise it uses the ``nanmedian`` of the selected array library.
-  If that library has no ``nanmedian`` (for example ``array-api-strict``),
-  `ccdproc`_ falls back to a sort-based implementation written purely in terms
-  of the array API standard. That fallback is correct but slower
-  (O(n log n) along the combination axis rather than O(n)).
++ The NaN-aware reductions ``nansum``, ``nanmean``, ``nanstd`` and
+  ``nanmedian`` are not part of the array API, but most array libraries do
+  provide them. When combining images, `ccdproc`_ uses the versions from
+  `bottleneck`_ for `numpy`_ arrays if `bottleneck`_ is installed; otherwise
+  it uses the ones the selected array library provides. If that library has
+  none (for example ``array-api-strict``), `ccdproc`_ falls back to
+  implementations written purely in terms of the array API standard. They
+  are correct but slower: the sum, mean and standard deviation take a few
+  extra passes over the data, and the median is sort-based (O(n log n)
+  along the combination axis rather than O(n)). The fallbacks promote
+  integer and boolean input to the library's default real floating dtype
+  and yield NaN silently for slices that are entirely NaN.
++ ``median`` is not part of the array API either. ``subtract_overscan``
+  uses the ``median`` of the selected array library when there is one, and
+  otherwise the same sort-based fallback, propagating NaNs as
+  ``numpy.median`` does.
++ ``sigma_func``, the default uncertainty estimate of ``median_combine``,
+  is a median absolute deviation. For `numpy`_ arrays it calls
+  ``astropy.stats.median_absolute_deviation``, which is numpy-only; for
+  every other array library it uses a version written purely in terms of
+  the array API standard, built on the sort-based medians above, so it
+  shares their cost and their silence on all-NaN slices.
 
 Which array library should I use?
 ---------------------------------

From d9001510c0916a57f5601e62f4ac557e412d19a0 Mon Sep 17 00:00:00 2001
From: Matt Craig <mattwcraig@gmail.com>
Date: Wed, 26 Aug 2026 20:30:46 -0500
Subject: [PATCH 03/11] Drop the sigma_func entry from the array-API escape
 baseline

sigma_func no longer converts non-numpy data to numpy, so the escape is
not observed any more. Regenerated with a full-suite dask run
(CCDPROC_ARRAY_LIBRARY=dask CCDPROC_LOG_ARRAY_ESCAPES=1
CCDPROC_WRITE_ESCAPE_BASELINE=1); the only change is this deleted line.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01XGArZMnLPC1G7YP62vN18V
---
 ccdproc/tests/array_escape_baseline.txt | 1 -
 1 file changed, 1 deletion(-)

diff --git a/ccdproc/tests/array_escape_baseline.txt b/ccdproc/tests/array_escape_baseline.txt
index 6b82c205..634a21d9 100644
--- a/ccdproc/tests/array_escape_baseline.txt
+++ b/ccdproc/tests/array_escape_baseline.txt
@@ -23,5 +23,4 @@ core.py      block_reduce                 numpy.asanyarray  BOUNDARY: astropy.nd
 core.py      block_replicate              numpy.asanyarray  BOUNDARY: astropy.nddata.block_replicate is numpy-only
 core.py      cosmicray_median             numpy.asarray     BOUNDARY: scipy.ndimage median_filter/maximum_filter are numpy-only
 core.py      create_deviation             numpy.asarray     BOUNDARY: astropy CCDData.copy()/StdDevUncertainty are numpy-backed
-core.py      sigma_func                   numpy.asanyarray  BOUNDARY: astropy.stats.median_absolute_deviation is numpy-only
 core.py      subtract_overscan            numpy.asanyarray  BOUNDARY: astropy.modeling fit/evaluation (model= path) is numpy-only

From b9f294061522f86b817f9a174eff6e8d00667c63 Mon Sep 17 00:00:00 2001
From: Matt Craig <mattwcraig@gmail.com>
Date: Wed, 26 Aug 2026 20:44:49 -0500
Subject: [PATCH 04/11] Fill in PR number in changelog entry

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01XGArZMnLPC1G7YP62vN18V
---
 CHANGES.rst | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/CHANGES.rst b/CHANGES.rst
index 30502d3c..7d3f0c30 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -125,7 +125,7 @@ Bug Fixes
   absolute deviation written purely in terms of the array API when the array
   namespace is not numpy, instead of converting the data to numpy through
   ``astropy.stats.median_absolute_deviation``. Numpy input still uses
-  astropy. [#TBD]
+  astropy. [#1000]
 
 2.5.1 (2025-07-05)
 ------------------

From dce00b87af2154ba340a600d56919a30c003013f Mon Sep 17 00:00:00 2001
From: Matt Craig <mattwcraig@gmail.com>
Date: Wed, 26 Aug 2026 21:21:22 -0500
Subject: [PATCH 05/11] Honor the CCDData mask in sigma_func on numpy;
 normalize tuple axis for astropy+bottleneck

astropy's median_absolute_deviation only honors the mask of an explicit
numpy.ma.MaskedArray, so hand it one when the CCDData has a mask instead of
relying on numpy.nanmedian noticing the mask of CCDData.__array__'s output
(which only happens on its small-array path and never with bottleneck).

The astropy reference in test_mad_fallback_matches_astropy gets a tuple axis
with negative entries normalized: astropy's bottleneck dispatch transposes
with the tuple as given and raises on a negative entry.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01JnWCg95xE93SbhME52jxGJ
---
 CHANGES.rst                   |  4 ++++
 ccdproc/core.py               | 30 ++++++++++++++----------
 ccdproc/tests/test_ccdproc.py | 44 ++++++++++++++++++++++++-----------
 3 files changed, 53 insertions(+), 25 deletions(-)

diff --git a/CHANGES.rst b/CHANGES.rst
index 7d3f0c30..d690e5b0 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -126,6 +126,10 @@ Bug Fixes
   namespace is not numpy, instead of converting the data to numpy through
   ``astropy.stats.median_absolute_deviation``. Numpy input still uses
   astropy. [#1000]
+- ``sigma_func`` now always excludes the masked pixels of a masked
+  ``CCDData``. Previously the mask was only honoured on numpy for small
+  arrays with a single integer ``axis`` and ``ignore_nan=True``, and never
+  when bottleneck is installed. [#1000]
 
 2.5.1 (2025-07-05)
 ------------------
diff --git a/ccdproc/core.py b/ccdproc/core.py
index 951bb2c3..230f3a51 100644
--- a/ccdproc/core.py
+++ b/ccdproc/core.py
@@ -317,9 +317,8 @@ def _mad_fallback(data, axis, ignore_nan, xp=None, mask=None):
 
     mask : array or None, optional
         Boolean mask with the shape of ``data``. Masked pixels are excluded
-        from the statistics (which forces ``ignore_nan`` on), as astropy
-        does for a masked `~astropy.nddata.CCDData` when ``ignore_nan`` is
-        set.
+        from the statistics (which forces ``ignore_nan`` on), as `sigma_func`
+        arranges on numpy by handing astropy a `numpy.ma.MaskedArray`.
 
     Returns
     -------
@@ -1517,10 +1516,9 @@ def sigma_func(arr, axis=None, ignore_nan=False):
     the namespace's default real floating dtype, and it yields NaN silently
     for slices that are entirely NaN, where numpy would warn.
 
-    A masked `~astropy.nddata.CCDData` is handed to astropy as is on numpy.
-    The fallback excludes the masked pixels from the statistics (which
-    implies ``ignore_nan``), as astropy does for a single integer ``axis``
-    with ``ignore_nan=True``.
+    The masked pixels of a masked `~astropy.nddata.CCDData` are excluded
+    from the statistics on every backend, which implies ``ignore_nan``; a
+    slice that is entirely masked gives NaN.
     """
     if isinstance(arr, CCDData):
         data = arr.data
@@ -1531,11 +1529,19 @@ def sigma_func(arr, axis=None, ignore_nan=False):
     xp = array_api_compat.array_namespace(data)
 
     if array_api_compat.is_numpy_namespace(xp):
-        # Pass ``arr`` rather than ``data``: astropy honours a CCDData mask.
-        return xp.asarray(
-            stats.median_absolute_deviation(arr, axis=axis, ignore_nan=ignore_nan)
-            * 1.482602218505602
-        )
+        if mask is not None:
+            # astropy only honours the mask of a numpy.ma.MaskedArray. Passing
+            # the CCDData itself leaves it to numpy.nanmedian to notice the
+            # mask of the MaskedArray that CCDData.__array__ produces, which
+            # only happens on its small-array path (fewer than 600 elements
+            # along ``axis``) and never once bottleneck is installed.
+            data = np.ma.masked_array(data, mask=mask)
+            ignore_nan = True
+        result = stats.median_absolute_deviation(data, axis=axis, ignore_nan=ignore_nan)
+        if np.ma.isMaskedArray(result):
+            # Entirely masked slices, which the fallback also reports as NaN.
+            result = result.filled(np.nan)
+        return xp.asarray(result * 1.482602218505602)
 
     return _mad_fallback(data, axis, ignore_nan, xp=xp, mask=mask) * 1.482602218505602
 
diff --git a/ccdproc/tests/test_ccdproc.py b/ccdproc/tests/test_ccdproc.py
index 2844ade0..f3511f41 100644
--- a/ccdproc/tests/test_ccdproc.py
+++ b/ccdproc/tests/test_ccdproc.py
@@ -447,8 +447,16 @@ def test_mad_fallback_matches_astropy(data, axis, ignore_nan):
     # sigma_func only reaches _mad_fallback on non-numpy namespaces, none of
     # which report coverage, so the fallback is exercised directly here on
     # every backend.
+    #
+    # The reference gets a tuple axis with its negative entries normalised:
+    # astropy's bottleneck dispatch (astropy.stats.nanfunctions.
+    # _move_tuple_axes_last) transposes with the tuple as given and raises
+    # on a negative entry. The fallback still receives the tuple as written.
+    reference_axis = axis
+    if isinstance(axis, tuple):
+        reference_axis = tuple(ax % data.ndim for ax in axis)
     expected_np = np_asarray(
-        median_absolute_deviation(data, axis=axis, ignore_nan=ignore_nan)
+        median_absolute_deviation(data, axis=reference_axis, ignore_nan=ignore_nan)
     )
     expected = xp.asarray(expected_np, device=xp_device)
 
@@ -535,26 +543,36 @@ def test_sigma_func_fallback_branch_on_any_backend(monkeypatch):
     assert xp.all(xpx.isclose(result, expected, equal_nan=True))
 
 
-def test_sigma_func_ccddata_mask_is_honoured():
+# The ignore mark is for numpy.nanmedian on the entirely masked column of
+# the NaN-filled reference; the masked path itself does not warn.
+@pytest.mark.filterwarnings("ignore:All-NaN slice encountered:RuntimeWarning")
+@pytest.mark.parametrize(
+    ("axis", "ignore_nan"),
+    # axis=0 with ignore_nan=True is the form median_combine uses.
+    [(0, True), (None, False), (1, False), ((0, 1), True)],
+)
+def test_sigma_func_ccddata_mask_is_honoured(axis, ignore_nan):
     ccd = ccd_data_func()
-    mask = xp.asarray(RNG(929).random(ccd.shape) > 0.7, device=xp_device)
+    mask_np = RNG(929).random(ccd.shape) > 0.7
+    mask_np[:, 3] = True  # an entirely masked column gives NaN along axis 0
+    mask = xp.asarray(mask_np, device=xp_device)
     # TODO: Set .mask instead of ._mask when CCDData is array-api compliant
     ccd._mask = mask
     nan = xp.asarray(np_nan, dtype=ccd.data.dtype, device=xp_device)
     nanned = xp.where(mask, nan, ccd.data)
 
-    # A single integer axis with ignore_nan=True is the form median_combine
-    # uses, and the one form for which astropy (the numpy path) also excludes
-    # the masked pixels, so the public function can be checked on every backend.
-    result = sigma_func(ccd, axis=0, ignore_nan=True)
-    expected = sigma_func(nanned, axis=0, ignore_nan=True)
+    # Masked pixels are excluded for every axis and ignore_nan, on numpy
+    # (where astropy is handed a numpy.ma.MaskedArray) as in the fallback,
+    # so the result is that of the NaN-filled data with ignore_nan on.
+    expected = sigma_func(nanned, axis=axis, ignore_nan=True)
+    result = sigma_func(ccd, axis=axis, ignore_nan=ignore_nan)
+    assert result.shape == expected.shape
     assert xp.all(xpx.isclose(result, expected, equal_nan=True))
 
-    # The fallback excludes the masked pixels for every axis and ignore_nan.
-    for axis, ignore_nan in [(None, False), (1, False), ((0, 1), True)]:
-        result = _mad_fallback(ccd.data, axis, ignore_nan, mask=mask)
-        expected = _mad_fallback(nanned, axis, True)
-        assert xp.all(xpx.isclose(result, expected, equal_nan=True))
+    # The public function only reaches the fallback's mask handling on
+    # non-numpy namespaces, so exercise it directly here on every backend.
+    result = _mad_fallback(ccd.data, axis, ignore_nan, mask=mask) * 1.482602218505602
+    assert xp.all(xpx.isclose(result, expected, equal_nan=True))
 
 
 def test_trim_image_fits_section_requires_string():

From 6064b8d645d71dfc3e10142e42f2b4f61b09ff0c Mon Sep 17 00:00:00 2001
From: Matt Craig <mattwcraig@gmail.com>
Date: Sun, 30 Aug 2026 14:59:11 -0500
Subject: [PATCH 06/11] Simplify sigma_func's mask handling and _mad_fallback's
 axis/median logic

Move the CCDData mask to NaN substitution out of _mad_fallback and into
sigma_func, before the numpy/fallback branch split, so both branches
handle a masked array identically: NaN-substituted data with ignore_nan
forced on. The promotion to a real floating dtype happens first, since
xp.where(mask, float_nan, int_data) raises on array-api-strict instead
of promoting silently the way numpy does.

On numpy this replaces the numpy.ma.MaskedArray construction with the
same NaN-filled array astropy already accepts, so numpy now takes its
nanmedian (bottleneck-accelerated, when installed) path instead of
np.ma.median, and an entirely masked slice at axis=None gives NaN
instead of the previous 0.0.

Factor the promote-to-real-dtype step, shared by _nanfuncs._setup,
_mad_fallback and the new mask substitution, into a
ccdproc._nanfuncs._promote_to_real() helper.

In _mad_fallback, replace the hand-rolled tuple-axis validation loop
with numpy's normalize_axis_tuple (numpy 2 only; the next ccdproc
release is expected to require it), which also accepts a list axis; a
bool entry is still rejected explicitly with TypeError, since
normalize_axis_tuple treats it as an int. Also prefer the array
namespace's own median/nanmedian, falling back to the sort-based
ccdproc._nanfuncs implementations only on AttributeError, mirroring
_median_fallback's existing pattern; this drops the O(n log n) sort
cost and the warning-free silence on all-NaN slices for any backend
that provides a native median (dask and jax do; array-api-strict does
not). docs/array_api.rst is updated to match.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01LVxnTLrWKStBdPwcmDNhxA
---
 ccdproc/_nanfuncs.py | 37 +++++++++++++----
 ccdproc/core.py      | 99 ++++++++++++++++++++------------------------
 docs/array_api.rst   |  6 ++-
 3 files changed, 79 insertions(+), 63 deletions(-)

diff --git a/ccdproc/_nanfuncs.py b/ccdproc/_nanfuncs.py
index 75afb378..3f66b2da 100644
--- a/ccdproc/_nanfuncs.py
+++ b/ccdproc/_nanfuncs.py
@@ -26,6 +26,35 @@
 __all__ = ["median", "nanmean", "nanmedian", "nanstd", "nansum"]
 
 
+def _promote_to_real(x, xp, device):
+    """
+    Promote integer and boolean ``x`` to the namespace's default real
+    floating dtype; a real floating ``x``, including float32, passes
+    through unchanged.
+
+    Parameters
+    ----------
+    x : array
+        Input array.
+    xp : array namespace
+        Namespace to use.
+    device : device
+        Device on which to resolve the default real floating dtype.
+
+    Returns
+    -------
+    array
+        ``x``, promoted if necessary.
+    """
+    if xp.isdtype(x.dtype, "real floating"):
+        return x
+    # Promote to the namespace's default real dtype rather than hardcoding
+    # float64: jax without JAX_ENABLE_X64 has no float64 and warns when one
+    # is requested, which pytest's filterwarnings turns into an error.
+    info = xp.__array_namespace_info__()
+    return xp.astype(x, info.default_dtypes(device=device)["real floating"])
+
+
 def _setup(x, axis, xp):
     """
     Validate ``axis``, resolve the namespace and device, promote to float.
@@ -83,13 +112,7 @@ def _setup(x, axis, xp):
     axis = axis % ndim
 
     device = array_api_compat.device(x)
-
-    if not xp.isdtype(x.dtype, "real floating"):
-        # Promote to the namespace's default real dtype rather than hardcoding
-        # float64: jax without JAX_ENABLE_X64 has no float64 and warns when one
-        # is requested, which pytest's filterwarnings turns into an error.
-        info = xp.__array_namespace_info__()
-        x = xp.astype(x, info.default_dtypes(device=device)["real floating"])
+    x = _promote_to_real(x, xp, device)
 
     return x, axis, xp, device
 
diff --git a/ccdproc/core.py b/ccdproc/core.py
index 230f3a51..35f72c21 100644
--- a/ccdproc/core.py
+++ b/ccdproc/core.py
@@ -5,7 +5,6 @@
 import logging
 import math
 import numbers
-import operator
 import warnings
 
 import array_api_compat
@@ -22,6 +21,7 @@
 from astropy.utils import deprecated, deprecated_renamed_argument
 from astropy.wcs.utils import proj_plane_pixel_area
 from numpy import mgrid as np_mgrid
+from numpy.lib.array_utils import normalize_axis_tuple
 from numpy.ma import nomask as np_ma_nomask
 from packaging import version as pkgversion
 from scipy import ndimage
@@ -30,6 +30,7 @@
     _unwrap_ccddata_for_array_api,
     _wrap_ccddata_for_array_api,
 )
+from ._nanfuncs import _promote_to_real
 from ._nanfuncs import median as _nanfuncs_median
 from ._nanfuncs import nanmedian as _nanfuncs_nanmedian
 from .log_meta import log_to_metadata
@@ -288,7 +289,7 @@ def _median_fallback(array, axis, xp=None):
         return _nanfuncs_median(array, axis=axis, xp=xp)
 
 
-def _mad_fallback(data, axis, ignore_nan, xp=None, mask=None):
+def _mad_fallback(data, axis, ignore_nan, xp=None):
     """
     Median absolute deviation written purely in terms of the array API.
 
@@ -303,10 +304,10 @@ def _mad_fallback(data, axis, ignore_nan, xp=None, mask=None):
         and boolean input is promoted to the namespace's default real
         floating dtype.
 
-    axis : int, tuple of int or None
+    axis : int, tuple of int, list of int or None
         Axis or axes along which the deviation is computed. ``None``
-        flattens ``data`` first; a tuple reduces over all the listed axes.
-        Negative values count from the last axis.
+        flattens ``data`` first; a tuple or list reduces over all the
+        listed axes. Negative values count from the last axis.
 
     ignore_nan : bool
         If `True`, NaNs are ignored; otherwise a NaN anywhere in a reduced
@@ -315,11 +316,6 @@ def _mad_fallback(data, axis, ignore_nan, xp=None, mask=None):
     xp : array namespace, optional
         Namespace to use. If not provided, it is determined from ``data``.
 
-    mask : array or None, optional
-        Boolean mask with the shape of ``data``. Masked pixels are excluded
-        from the statistics (which forces ``ignore_nan`` on), as `sigma_func`
-        arranges on numpy by handing astropy a `numpy.ma.MaskedArray`.
-
     Returns
     -------
     mad : array
@@ -330,58 +326,53 @@ def _mad_fallback(data, axis, ignore_nan, xp=None, mask=None):
     Raises
     ------
     ValueError
-        If ``axis`` is out of bounds or lists an axis twice.
+        If a tuple or list ``axis`` is out of bounds for ``data`` or lists
+        an axis more than once (including via a negative alias).
+    TypeError
+        If a tuple or list ``axis`` contains a bool or a non-integer entry.
 
     Notes
     -----
-    Both medians are the sort-based fallbacks from `ccdproc._nanfuncs`, so
-    the cost is O(n log n) along the reduced axes rather than the O(n) of
-    `numpy.median`. Slices that are entirely NaN yield NaN without the
-    ``RuntimeWarning`` numpy emits.
+    The namespace's own ``nanmedian``/``median`` is used when it has one;
+    otherwise the sort-based fallbacks from `ccdproc._nanfuncs` are used
+    instead, which cost O(n log n) along the reduced axes rather than the
+    O(n) of a native median, and yield NaN for an entirely NaN slice
+    without the ``RuntimeWarning`` numpy emits.
     """
     xp = xp or array_api_compat.array_namespace(data)
     device = array_api_compat.device(data)
 
-    # The _nanfuncs medians promote internally, but ``data - center`` below
+    # The medians below promote internally, but ``data - center`` below
     # needs the promoted input too: subtracting a float from an integer array
     # is not allowed by the standard and raises on array-api-strict.
-    if not xp.isdtype(data.dtype, "real floating"):
-        info = xp.__array_namespace_info__()
-        data = xp.astype(data, info.default_dtypes(device=device)["real floating"])
-
-    if mask is not None:
-        nan = xp.asarray(xp.nan, dtype=data.dtype, device=device)
-        mask = xp.astype(xp.asarray(mask, device=device), xp.bool)
-        data = xp.where(mask, nan, data)
-        ignore_nan = True
+    data = _promote_to_real(data, xp, device)
 
     if axis is None:
         data = xp.reshape(data, (-1,))
         axis = 0
-    elif isinstance(axis, tuple):
+    elif isinstance(axis, (tuple, list)):
+        if any(isinstance(ax, bool) for ax in axis):
+            raise TypeError("axis entries must be integers, not bool")
         ndim = data.ndim
-        axes = []
-        for ax in axis:
-            ax = operator.index(ax)
-            if not -ndim <= ax < ndim:
-                raise ValueError(
-                    f"axis {ax} is out of bounds for array of dimension {ndim}"
-                )
-            axes.append(ax % ndim)
-        if len(set(axes)) != len(axes):
-            raise ValueError(f"duplicate value in 'axis': {axis}")
+        axes = normalize_axis_tuple(axis, ndim)
         # Move the reduced axes to the end and merge them into one, so that
         # the single-axis medians below reduce over all of them at once.
         kept = [ax for ax in range(ndim) if ax not in axes]
-        data = xp.permute_dims(data, tuple(kept + axes))
+        data = xp.permute_dims(data, tuple(kept + list(axes)))
         data = xp.reshape(
             data, tuple(data.shape[ax] for ax in range(len(kept))) + (-1,)
         )
         axis = -1
 
-    med = _nanfuncs_nanmedian if ignore_nan else _nanfuncs_median
-    center = med(data, axis=axis, xp=xp)
-    return med(xp.abs(data - xp.expand_dims(center, axis=axis)), axis=axis, xp=xp)
+    def med(d, ax):
+        try:
+            return (xp.nanmedian if ignore_nan else xp.median)(d, axis=ax)
+        except AttributeError:
+            fallback = _nanfuncs_nanmedian if ignore_nan else _nanfuncs_median
+            return fallback(d, axis=ax, xp=xp)
+
+    center = med(data, axis)
+    return med(xp.abs(data - xp.expand_dims(center, axis=axis)), axis)
 
 
 @log_to_metadata
@@ -1510,8 +1501,9 @@ def sigma_func(arr, axis=None, ignore_nan=False):
     For numpy input the deviation is computed by
     `astropy.stats.median_absolute_deviation`, exactly as before. For every
     other array namespace it is computed by a fallback written purely in
-    terms of the array API standard, on the device of the input. That
-    fallback uses sort-based medians, so it costs O(n log n) along the
+    terms of the array API standard, on the device of the input, using the
+    namespace's own ``median``/``nanmedian`` when it provides one and
+    otherwise a sort-based implementation that costs O(n log n) along the
     reduced axes rather than O(n); it promotes integer and boolean input to
     the namespace's default real floating dtype, and it yields NaN silently
     for slices that are entirely NaN, where numpy would warn.
@@ -1528,22 +1520,21 @@ def sigma_func(arr, axis=None, ignore_nan=False):
         mask = None
     xp = array_api_compat.array_namespace(data)
 
+    if mask is not None:
+        device = array_api_compat.device(data)
+        # The promotion has to happen before the ``where`` below: on
+        # array-api-strict, ``xp.where(mask, float_nan, int_data)`` raises
+        # rather than silently promoting the way numpy does.
+        data = _promote_to_real(data, xp, device)
+        nan = xp.asarray(xp.nan, dtype=data.dtype, device=device)
+        data = xp.where(xp.asarray(mask, device=device), nan, data)
+        ignore_nan = True
+
     if array_api_compat.is_numpy_namespace(xp):
-        if mask is not None:
-            # astropy only honours the mask of a numpy.ma.MaskedArray. Passing
-            # the CCDData itself leaves it to numpy.nanmedian to notice the
-            # mask of the MaskedArray that CCDData.__array__ produces, which
-            # only happens on its small-array path (fewer than 600 elements
-            # along ``axis``) and never once bottleneck is installed.
-            data = np.ma.masked_array(data, mask=mask)
-            ignore_nan = True
         result = stats.median_absolute_deviation(data, axis=axis, ignore_nan=ignore_nan)
-        if np.ma.isMaskedArray(result):
-            # Entirely masked slices, which the fallback also reports as NaN.
-            result = result.filled(np.nan)
         return xp.asarray(result * 1.482602218505602)
 
-    return _mad_fallback(data, axis, ignore_nan, xp=xp, mask=mask) * 1.482602218505602
+    return _mad_fallback(data, axis, ignore_nan, xp=xp) * 1.482602218505602
 
 
 def setbox(x, y, mbox, xmax, ymax):
diff --git a/docs/array_api.rst b/docs/array_api.rst
index 1c7d2d62..43160157 100644
--- a/docs/array_api.rst
+++ b/docs/array_api.rst
@@ -169,8 +169,10 @@ What limitations should I be aware of?
   is a median absolute deviation. For `numpy`_ arrays it calls
   ``astropy.stats.median_absolute_deviation``, which is numpy-only; for
   every other array library it uses a version written purely in terms of
-  the array API standard, built on the sort-based medians above, so it
-  shares their cost and their silence on all-NaN slices.
+  the array API standard, preferring the library's own ``median``/
+  ``nanmedian`` and falling back to the sort-based medians above only when
+  the library has neither, so the extra sort cost and the silence on
+  all-NaN slices apply only there.
 
 Which array library should I use?
 ---------------------------------

From 25059129b90e3fcc21bc968c17e77c96fd0a0a55 Mon Sep 17 00:00:00 2001
From: Matt Craig <mattwcraig@gmail.com>
Date: Sun, 30 Aug 2026 14:59:29 -0500
Subject: [PATCH 07/11] Consolidate the _mad_fallback and sigma_func tests

Shrink _MAD_CASES from 23 cases to 6: the odd/even lengths, both
all-NaN cases, the two 1-D cases and the numpy-integer axis case are
already pinned by test_nanfuncs.py::test_matches_numpy, and (0, 1),
(0, 2), (2, 0) and a bool array added nothing over the kept
(-1, 0)-with-negative-entry and int-array cases. test_mad_fallback_
matches_astropy now asserts result.dtype == expected.dtype, which
pins float32 preservation where the previous isdtype check did not.

Reduce the bad-axis coverage to test_mad_fallback_rejects_bad_axis_
tuple, over the tuple forms only; normalize_axis_tuple's own
out-of-bounds handling already exercises the int case through
test_nanfuncs.py::test_bad_axis. The raised message for a duplicate
axis is now numpy's own "repeated axis" rather than the previous
hand-written wording.

Fold test_sigma_func_keeps_namespace_and_device and
test_sigma_func_fallback_branch_on_any_backend into
test_sigma_func_matches_mad_std as a force_fallback parameter, and
drop test_mad_fallback_all_nan_slice_is_silent: once _mad_fallback
prefers a native nanmedian, the silence on an all-NaN slice is a
property of ccdproc._nanfuncs (pinned in test_nanfuncs.py) rather
than of _mad_fallback, and would fail rather than merely duplicate on
a backend with a native, warning nanmedian.

Rename test_sigma_func_ccddata_mask_is_honoured to
..._is_honored, trim its parametrization to the (0, True) and
(None, False) cases that exercise a distinct branch, and add an
all-masked axis=None case that asserts the result is NaN, pinning
the 0.0 bug fixed in the previous commit. Drop the trailing direct
_mad_fallback(..., mask=mask) call, since _mad_fallback no longer
takes a mask argument.

Give each surviving test a one-line leading comment, since the file
uses that style instead of docstrings, and fix "normalised" to
"normalized" in a comment that carries over unchanged.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01LVxnTLrWKStBdPwcmDNhxA
---
 ccdproc/tests/test_ccdproc.py | 185 ++++++++++++++--------------------
 1 file changed, 73 insertions(+), 112 deletions(-)

diff --git a/ccdproc/tests/test_ccdproc.py b/ccdproc/tests/test_ccdproc.py
index f3511f41..22fc64eb 100644
--- a/ccdproc/tests/test_ccdproc.py
+++ b/ccdproc/tests/test_ccdproc.py
@@ -23,9 +23,10 @@
 from numpy import array as np_array
 from numpy import asarray as np_asarray
 from numpy import float32 as np_float32
-from numpy import int64 as np_int64
 from numpy import mgrid as np_mgrid
 from numpy import nan as np_nan
+from numpy import nan_to_num as np_nan_to_num
+from numpy import ones as np_ones
 from numpy import random as np_random
 
 from ccdproc.conftest import testing_array_device as xp_device
@@ -404,54 +405,31 @@ def __getattr__(self, name):
 
 _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
-_MAD_CLEAN_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 = [
-    # Every axis form sigma_func's callers use, on data with scattered NaNs:
-    # None (background_deviation_box), a single axis (median_combine), a
-    # numpy integer, and tuples in either order with negative entries.
-    *[
-        (_MAD_3D, axis)
-        for axis in [
-            None,
-            0,
-            1,
-            -1,
-            np_int64(1),
-            (0, 1),
-            (0, 2),
-            (2, 0),
-            (-1, 0),
-            (0, 1, 2),
-        ]
-    ],
-    *[(_MAD_RNG.normal(size=(n, 7)), 0) for n in range(1, 7)],  # odd/even lengths
-    (np_array([1.0, 2.0, 3.0, 4.0]), 0),  # 1-D
-    (np_array([1.0, 2.0, 3.0, 4.0]), None),
-    (np_array([[1.0, np_nan], [2.0, np_nan], [3.0, np_nan]]), 0),  # all-NaN column
-    (np_array([np_nan, np_nan, np_nan]), None),  # every value NaN
-    (np_array([[1, 4], [2, 3], [5, 6], [4, 1]]), 0),  # integer input
-    (np_array([[True, False], [False, True], [True, True]]), 0),  # boolean input
-    (np_array([[1.0, 2.0, 3.0], [4.0, 5.0, 7.0]], dtype=np_float32), 1),
+    # 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
 ]
 
 
-# The ignore mark lets the astropy reference warn on all-NaN slices where
-# the fallback deliberately does not; test_mad_fallback_all_nan_slice_is_silent
-# owns the fallback's silence.
+# _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). This is a
+# differential check against astropy over the axis forms, not a coverage
+# device.
+#
+# The reference gets a tuple axis with its negative entries normalised:
+# astropy's bottleneck dispatch (astropy.stats.nanfunctions.
+# _move_tuple_axes_last) transposes with the tuple as given and raises
+# on a negative entry. The fallback still receives the tuple as written.
 @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):
-    # sigma_func only reaches _mad_fallback on non-numpy namespaces, none of
-    # which report coverage, so the fallback is exercised directly here on
-    # every backend.
-    #
-    # The reference gets a tuple axis with its negative entries normalised:
-    # astropy's bottleneck dispatch (astropy.stats.nanfunctions.
-    # _move_tuple_axes_last) transposes with the tuple as given and raises
-    # on a negative entry. The fallback still receives the tuple as written.
     reference_axis = axis
     if isinstance(axis, tuple):
         reference_axis = tuple(ax % data.ndim for ax in axis)
@@ -463,116 +441,99 @@ def test_mad_fallback_matches_astropy(data, axis, ignore_nan):
     result = _mad_fallback(xp.asarray(data, device=xp_device), axis, ignore_nan)
 
     assert result.shape == expected.shape
-    # Integer and boolean input is promoted to the namespace's default real
-    # dtype, as the _nanfuncs medians do.
-    assert xp.isdtype(result.dtype, "real floating")
+    assert result.dtype == expected.dtype
     assert xp.all(xpx.isclose(result, expected, equal_nan=True))
 
 
-def test_mad_fallback_all_nan_slice_is_silent():
-    data = xp.asarray(np_array([[np_nan, np_nan]] * 3), device=xp_device)
-    with warnings.catch_warnings():
-        warnings.simplefilter("error")
-        # The bool() calls force the computation inside the block on lazy
-        # backends, where a warning would otherwise surface at compute time.
-        assert bool(xp.all(xp.isnan(_mad_fallback(data, 0, True))))
-        assert bool(xp.isnan(_mad_fallback(data, None, True)))
-
-
-def test_mad_fallback_rejects_duplicate_axes():
-    data = xp.asarray(_MAD_CLEAN_3D, device=xp_device)
-    with pytest.raises(ValueError, match="duplicate"):
-        _mad_fallback(data, (0, 0), True)
-    # -3 is axis 0 of a 3-D array, so this is a duplicate too.
-    with pytest.raises(ValueError, match="duplicate"):
-        _mad_fallback(data, (0, -3), True)
-    with pytest.raises(ValueError, match="out of bounds"):
-        _mad_fallback(data, (0, 3), True)
-    with pytest.raises(ValueError, match="out of bounds"):
-        _mad_fallback(data, 3, True)
+# Duplicate (including a negative alias) and out-of-bounds entries in a
+# tuple axis raise ValueError instead of silently reducing the wrong axes.
+@pytest.mark.parametrize(
+    ("axis", "match"),
+    [
+        pytest.param((0, 0), "repeated axis", id="duplicate"),
+        # -3 is axis 0 of a 3-D array, so this is a duplicate too.
+        pytest.param((0, -3), "repeated axis", id="duplicate-negative-alias"),
+        pytest.param((0, 3), "out of bounds", id="out-of-bounds"),
+    ],
+)
+def test_mad_fallback_rejects_bad_axis_tuple(axis, match):
+    data = xp.asarray(_MAD_3D, device=xp_device)
+    with pytest.raises(ValueError, match=match):
+        _mad_fallback(data, axis, True)
 
 
+# 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).
+@pytest.mark.parametrize("force_fallback", [False, True])
 @pytest.mark.parametrize(
     ("data", "axis", "ignore_nan"),
     [
-        pytest.param(_MAD_CLEAN_3D, None, False, id="all-no_nan"),
+        pytest.param(np_nan_to_num(_MAD_3D), None, False, id="all-no_nan"),
         pytest.param(_MAD_3D, 0, True, id="axis0-ignore_nan"),
         pytest.param(_MAD_3D, (0, 1), True, id="axes01-ignore_nan"),
     ],
 )
-def test_sigma_func_matches_mad_std(data, axis, ignore_nan):
+def test_sigma_func_matches_mad_std(
+    data, axis, ignore_nan, force_fallback, monkeypatch
+):
+    if force_fallback:
+        # On numpy sigma_func hands the data to astropy directly; every
+        # other backend already takes the fallback branch, and none of them
+        # reports coverage, so make numpy take it too.
+        monkeypatch.setattr(array_api_compat, "is_numpy_namespace", lambda _xp: False)
+
     expected_np = np_asarray(mad_std(data, axis=axis, ignore_nan=ignore_nan))
     expected = xp.asarray(expected_np, device=xp_device)
 
-    result = sigma_func(
-        xp.asarray(data, device=xp_device), axis=axis, ignore_nan=ignore_nan
-    )
+    arr = xp.asarray(data, device=xp_device)
+    result = sigma_func(arr, axis=axis, ignore_nan=ignore_nan)
 
     assert result.shape == expected.shape
     assert xp.all(xpx.isclose(result, expected, equal_nan=True))
+    assert array_api_compat.array_namespace(result) is array_api_compat.array_namespace(
+        arr
+    )
+    if xp_device is not None:
+        assert array_api_compat.device(result) == xp_device
     if axis is None:
         # background_deviation_box relies on the 0-d result converting to float.
         assert result.shape == ()
         assert float(result) == pytest.approx(float(expected_np))
 
 
-def test_sigma_func_keeps_namespace_and_device():
-    data = xp.asarray(_MAD_3D, device=xp_device)
-    for axis in (None, 0, (0, 2)):
-        result = sigma_func(data, axis=axis, ignore_nan=True)
-        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
-
-
-def test_sigma_func_fallback_branch_on_any_backend(monkeypatch):
-    # On numpy sigma_func hands the data to astropy; every other backend
-    # takes the fallback branch, and none of them reports coverage, so make
-    # numpy take it too and check that the two branches agree.
-    data = xp.asarray(_MAD_3D, device=xp_device)
-    expected = sigma_func(data, axis=0, ignore_nan=True)
-
-    monkeypatch.setattr(array_api_compat, "is_numpy_namespace", lambda _xp: False)
-    result = sigma_func(data, axis=0, ignore_nan=True)
-
-    assert array_api_compat.array_namespace(result) is array_api_compat.array_namespace(
-        data
-    )
-    assert xp.all(xpx.isclose(result, expected, equal_nan=True))
-
-
-# The ignore mark is for numpy.nanmedian on the entirely masked column of
-# the NaN-filled reference; the masked path itself does not warn.
+# 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.
 @pytest.mark.filterwarnings("ignore:All-NaN slice encountered:RuntimeWarning")
 @pytest.mark.parametrize(
-    ("axis", "ignore_nan"),
-    # axis=0 with ignore_nan=True is the form median_combine uses.
-    [(0, True), (None, False), (1, False), ((0, 1), True)],
+    ("axis", "ignore_nan", "all_masked"),
+    [
+        # axis=0 with ignore_nan=True is the form median_combine uses.
+        pytest.param(0, True, False, id="axis0"),
+        pytest.param(None, False, False, id="axis-none"),
+        pytest.param(None, False, True, id="all-masked-axis-none"),
+    ],
 )
-def test_sigma_func_ccddata_mask_is_honoured(axis, ignore_nan):
+def test_sigma_func_ccddata_mask_is_honored(axis, ignore_nan, all_masked):
     ccd = ccd_data_func()
-    mask_np = RNG(929).random(ccd.shape) > 0.7
-    mask_np[:, 3] = True  # an entirely masked column gives NaN along axis 0
+    if all_masked:
+        mask_np = np_ones(ccd.shape, dtype=bool)
+    else:
+        mask_np = RNG(929).random(ccd.shape) > 0.7
+        mask_np[:, 3] = True  # an entirely masked column gives NaN along axis 0
     mask = xp.asarray(mask_np, device=xp_device)
     # TODO: Set .mask instead of ._mask when CCDData is array-api compliant
     ccd._mask = mask
     nan = xp.asarray(np_nan, dtype=ccd.data.dtype, device=xp_device)
     nanned = xp.where(mask, nan, ccd.data)
 
-    # Masked pixels are excluded for every axis and ignore_nan, on numpy
-    # (where astropy is handed a numpy.ma.MaskedArray) as in the fallback,
-    # so the result is that of the NaN-filled data with ignore_nan on.
     expected = sigma_func(nanned, axis=axis, ignore_nan=True)
     result = sigma_func(ccd, axis=axis, ignore_nan=ignore_nan)
     assert result.shape == expected.shape
     assert xp.all(xpx.isclose(result, expected, equal_nan=True))
-
-    # The public function only reaches the fallback's mask handling on
-    # non-numpy namespaces, so exercise it directly here on every backend.
-    result = _mad_fallback(ccd.data, axis, ignore_nan, mask=mask) * 1.482602218505602
-    assert xp.all(xpx.isclose(result, expected, equal_nan=True))
+    if all_masked:
+        assert bool(xp.isnan(result))
 
 
 def test_trim_image_fits_section_requires_string():

From 5eaa6c5991ab57002563a81d474dd815eb2e6d32 Mon Sep 17 00:00:00 2001
From: Matt Craig <mattwcraig@gmail.com>
Date: Sun, 30 Aug 2026 14:59:35 -0500
Subject: [PATCH 08/11] Correct the sigma_func changelog wording for #1000

The entry about the mask fix said the mask was never honored on numpy
once bottleneck was installed. That was too strong: astropy's
bottleneck dispatch (astropy >= 7) routes only float64 data to
bottleneck, bypassing the mask; other dtypes still go through
numpy.ma even with bottleneck installed. Reword to what is actually
true on every astropy this release supports: not honored for float64
data specifically. Also fix "honoured" to "honored".

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01LVxnTLrWKStBdPwcmDNhxA
---
 CHANGES.rst | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/CHANGES.rst b/CHANGES.rst
index d690e5b0..4510f2a1 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -127,9 +127,9 @@ Bug Fixes
   ``astropy.stats.median_absolute_deviation``. Numpy input still uses
   astropy. [#1000]
 - ``sigma_func`` now always excludes the masked pixels of a masked
-  ``CCDData``. Previously the mask was only honoured on numpy for small
-  arrays with a single integer ``axis`` and ``ignore_nan=True``, and never
-  when bottleneck is installed. [#1000]
+  ``CCDData``. Previously the mask was only honored on numpy for small
+  arrays with a single integer ``axis`` and ``ignore_nan=True``, and, when
+  bottleneck is installed, not for float64 data at all. [#1000]
 
 2.5.1 (2025-07-05)
 ------------------

From dd58f7566e15451a95bbc9bb226fba12198d346b Mon Sep 17 00:00:00 2001
From: Matt Craig <mattwcraig@gmail.com>
Date: Sun, 30 Aug 2026 15:02:15 -0500
Subject: [PATCH 09/11] Update the mask test comment for the unified
 substitution

The endorsed comment wording predates moving the mask->NaN substitution
ahead of the namespace split: numpy no longer goes through numpy.ma, and
the filterwarnings mark is now for numpy.nanmedian warning on the
entirely masked column and on the all-masked input.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01LVxnTLrWKStBdPwcmDNhxA
---
 ccdproc/tests/test_ccdproc.py | 8 +++++---
 1 file changed, 5 insertions(+), 3 deletions(-)

diff --git a/ccdproc/tests/test_ccdproc.py b/ccdproc/tests/test_ccdproc.py
index 22fc64eb..acd14290 100644
--- a/ccdproc/tests/test_ccdproc.py
+++ b/ccdproc/tests/test_ccdproc.py
@@ -502,9 +502,11 @@ def test_sigma_func_matches_mad_std(
         assert float(result) == pytest.approx(float(expected_np))
 
 
-# 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.
+# Masked CCDData pixels are excluded on every backend by the same mask->NaN
+# substitution in sigma_func (which forces ignore_nan on); an all-masked
+# slice gives NaN, not 0.0. The ignore mark: the NaN-filled data reach
+# numpy.nanmedian, which warns on the entirely masked column and on the
+# all-masked input.
 @pytest.mark.filterwarnings("ignore:All-NaN slice encountered:RuntimeWarning")
 @pytest.mark.parametrize(
     ("axis", "ignore_nan", "all_masked"),

From d0b2245d443f2d850d262b42ef123dee76156092 Mon Sep 17 00:00:00 2001
From: Matt Craig <mattwcraig@gmail.com>
Date: Sun, 30 Aug 2026 15:25:51 -0500
Subject: [PATCH 10/11] Require numpy >= 2.0; raise the astropy, reproject and
 astroscrappy floors to match

The oldestdeps CI job (numpy 1.26) fails on the new numpy-2-only
normalize_axis_tuple import; per the review decision on #1000 the fix is
to raise the minimum, not to add a 1.26 shim. Closes #1003.

numpy >= 2.0 forces the rest: astropy 6.0.* caps numpy below 2, so the
oldest astropy becomes 6.1; the reproject 0.9.1 and astroscrappy 1.1.0
wheels are numpy-1 binaries, so the floors move to 0.14 and 1.2, the
oldest releases built against numpy 2. The numpy126 tox factor is gone
(numpy200 is the oldest now, and the bottleneck CI job uses it) and so is
the astroscrappy11 factor, which pinned numpy below 2.

Verified locally in a fresh venv with numpy 2.0.2, astropy 6.1.7,
reproject 0.14.0 and astroscrappy 1.2.0: test_ccdproc.py and
test_nanfuncs.py pass (219 tests).

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01LVxnTLrWKStBdPwcmDNhxA
---
 .github/workflows/ci_tests.yml |  2 +-
 CHANGES.rst                    |  4 ++++
 pyproject.toml                 |  8 ++++----
 tox.ini                        | 15 +++++----------
 4 files changed, 14 insertions(+), 15 deletions(-)

diff --git a/.github/workflows/ci_tests.yml b/.github/workflows/ci_tests.yml
index 2fbbd2ce..57f4346f 100644
--- a/.github/workflows/ci_tests.yml
+++ b/.github/workflows/ci_tests.yml
@@ -56,7 +56,7 @@ jobs:
           - name: 'ubuntu-py312-bottleneck'
             os: ubuntu-latest
             python: '3.12'
-            tox_env: 'py312-test-alldeps-numpy126-bottleneck'
+            tox_env: 'py312-test-alldeps-numpy200-bottleneck'
 
           - name: 'macos-py312-dask'
             os: macos-latest
diff --git a/CHANGES.rst b/CHANGES.rst
index 4510f2a1..205f78d2 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -22,6 +22,10 @@ Other Changes and Additions
   without having to name the interpreter. [#986]
 - ``combine`` no longer accepts an array as its ``array_package`` argument;
   pass the array namespace or module instead, as for ``Combiner``. [#997]
+- The minimum supported numpy is now 2.0 (``normalize_axis_tuple`` is
+  imported from its numpy 2 location with no 1.26 fallback), which raises
+  the minimum astropy to 6.1, reproject to 0.14 and astroscrappy to 1.2,
+  the oldest releases that work with numpy 2. [#1000]
 
 Bug Fixes
 ^^^^^^^^^
diff --git a/pyproject.toml b/pyproject.toml
index 7f5a5ba4..189d7eb2 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -17,10 +17,10 @@ authors = [
 dependencies = [
     "array_api_compat>=1.12.0",
     "array_api_extra>=0.7.0",
-    "astropy>=6.0.1",
-    "astroscrappy>=1.1.0",
-    "numpy>=1.26",
-    "reproject>=0.9.1",
+    "astropy>=6.1",
+    "astroscrappy>=1.2",
+    "numpy>=2.0",
+    "reproject>=0.14",
     "scipy",
 ]
 
diff --git a/tox.ini b/tox.ini
index 7b66af19..9fb76a3c 100644
--- a/tox.ini
+++ b/tox.ini
@@ -52,7 +52,6 @@ description =
     devdeps: with the latest developer version of key dependencies
     oldestdeps: with the oldest supported version of key dependencies
     cov: and test coverage
-    numpy126: with numpy 1.26.*
     numpy200: with numpy 2.0.*
     numpy210: with numpy 2.1.*
     bottleneck: with bottleneck
@@ -66,13 +65,9 @@ description =
 deps =
     cov: coverage
 
-    numpy126: numpy==1.26.*  # currently oldest support numpy version
-    numpy200: numpy==2.0.*
+    numpy200: numpy==2.0.*  # currently oldest supported numpy version
     numpy210: numpy==2.1.*
 
-    astroscrappy11: astroscrappy==1.1.*
-    astroscrappy11: numpy<2.0
-
     bottleneck: bottleneck>=1.3.2
 
     devdeps: astropy>=0.0.dev0
@@ -80,9 +75,9 @@ deps =
 
     # Remember to transfer any changes here to setup.cfg also. Only listing
     # packages which are constrained in the setup.cfg
-    oldestdeps: numpy==1.26.*
-    oldestdeps: astropy==6.0.*
-    oldestdeps: reproject==0.9.1
+    oldestdeps: numpy==2.0.*
+    oldestdeps: astropy==6.1.*
+    oldestdeps: reproject==0.14.0
 
     dask: dask
     # The strict env intentionally installs only the test extras plus
@@ -97,7 +92,7 @@ commands =
     cov: pytest --pyargs ccdproc {toxinidir}/docs --cov ccdproc --cov-config={toxinidir}/pyproject.toml {posargs}
     cov: coverage xml -o {toxinidir}/coverage.xml
     # install astroscrappy after numpy
-    oldestdeps: python -m pip install astroscrappy==1.1.0
+    oldestdeps: python -m pip install astroscrappy==1.2.0
     # Do not care about warnings on the oldest builds
     oldestdeps: pytest --pyargs ccdproc {toxinidir}/docs -W ignore {posargs}
 

From 3487e78f2afe4f370cdc88f4372742425537d560 Mon Sep 17 00:00:00 2001
From: Matt Craig <mattwcraig@gmail.com>
Date: Sun, 30 Aug 2026 15:33:27 -0500
Subject: [PATCH 11/11] Cover _mad_fallback's bool-axis rejection and
 sort-median branch

codecov flagged the four lines the simplification left uncovered: the
TypeError for a bool entry in a tuple axis, and the except branch of
med(), which only runs naturally on a namespace with no native
nanmedian/median (array-api-strict, which uploads no coverage). Hide the
native functions behind a delegating proxy namespace, as
test_median_fallback_without_native_median already does.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01LVxnTLrWKStBdPwcmDNhxA
---
 ccdproc/tests/test_ccdproc.py | 42 ++++++++++++++++++++++++++++-------
 1 file changed, 34 insertions(+), 8 deletions(-)

diff --git a/ccdproc/tests/test_ccdproc.py b/ccdproc/tests/test_ccdproc.py
index acd14290..27ef1ab5 100644
--- a/ccdproc/tests/test_ccdproc.py
+++ b/ccdproc/tests/test_ccdproc.py
@@ -445,23 +445,49 @@ def test_mad_fallback_matches_astropy(data, axis, ignore_nan):
     assert xp.all(xpx.isclose(result, expected, equal_nan=True))
 
 
-# Duplicate (including a negative alias) and out-of-bounds entries in a
-# tuple axis raise ValueError instead of silently reducing the wrong axes.
+# Duplicate (including a negative alias), out-of-bounds and bool entries in
+# a tuple axis raise instead of silently reducing the wrong axes.
 @pytest.mark.parametrize(
-    ("axis", "match"),
+    ("axis", "error", "match"),
     [
-        pytest.param((0, 0), "repeated axis", id="duplicate"),
+        pytest.param((0, 0), ValueError, "repeated axis", id="duplicate"),
         # -3 is axis 0 of a 3-D array, so this is a duplicate too.
-        pytest.param((0, -3), "repeated axis", id="duplicate-negative-alias"),
-        pytest.param((0, 3), "out of bounds", id="out-of-bounds"),
+        pytest.param(
+            (0, -3), ValueError, "repeated axis", id="duplicate-negative-alias"
+        ),
+        pytest.param((0, 3), ValueError, "out of bounds", id="out-of-bounds"),
+        # normalize_axis_tuple would silently treat True as 1.
+        pytest.param((0, True), TypeError, "not bool", id="bool-entry"),
     ],
 )
-def test_mad_fallback_rejects_bad_axis_tuple(axis, match):
+def test_mad_fallback_rejects_bad_axis_tuple(axis, error, match):
     data = xp.asarray(_MAD_3D, device=xp_device)
-    with pytest.raises(ValueError, match=match):
+    with pytest.raises(error, match=match):
         _mad_fallback(data, axis, True)
 
 
+# The except branch of _mad_fallback's med() only runs naturally on
+# namespaces without a native nanmedian/median (array-api-strict is the only
+# such backend in CI, and it does not report coverage), so hide the native
+# functions behind a proxy namespace that otherwise delegates to the backend
+# under test, as test_median_fallback_without_native_median does.
+@pytest.mark.parametrize("ignore_nan", [False, True])
+def test_mad_fallback_without_native_medians(ignore_nan):
+    class _NamespaceWithoutMedians(types.ModuleType):
+        def __getattr__(self, name):
+            if name in ("median", "nanmedian"):
+                raise AttributeError(name)
+            return getattr(xp, name)
+
+    proxy = _NamespaceWithoutMedians("xp_without_medians")
+    data = xp.asarray(_MAD_3D, device=xp_device)
+
+    result = _mad_fallback(data, 0, ignore_nan, xp=proxy)
+
+    expected = _mad_fallback(data, 0, ignore_nan)
+    assert xp.all(xpx.isclose(result, expected, equal_nan=True))
+
+
 # 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).