From 850a4d6afd4ad8f1b3460112ef112feb555a34f8 Mon Sep 17 00:00:00 2001
From: Matt Craig <mattwcraig@gmail.com>
Date: Sun, 30 Aug 2026 19:13:03 -0500
Subject: [PATCH 1/3] Hoist None/tuple axis handling into _nanfuncs._setup

Move the flatten-for-None and permute-and-merge steps from
core._mad_fallback into the shared _nanfuncs._setup, as agreed in the
review of #1000. _setup now accepts a single integer, None, or a
tuple/list of axes, and returns a restore callable that maps a
full-shape array back to the caller's layout (the identity for a
single integer axis).

This gives every reduction fallback in _nanfuncs (nansum, nanmean,
nanstd, nanmedian, median, nanmad) tuple-axis and axis=None support
for free, shrinks _mad_fallback to the median tier plus a delegation,
and lifts the single-integer-axis restriction of
Combiner.sigma_clipping's array-API path, which now accepts the same
axis forms astropy.stats.sigma_clip does (the restore callable hands
the clip mask back in the data's own shape).

Bare-bool and non-integer axes now raise TypeError instead of
NotImplementedError, matching the TypeError a bool entry in a tuple
already raised.

Closes #1004

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01LVxnTLrWKStBdPwcmDNhxA
---
 CHANGES.rst                    |   5 +
 ccdproc/_nanfuncs.py           | 185 ++++++++++++++++++++++-----------
 ccdproc/combiner.py            |  39 ++++---
 ccdproc/core.py                |  40 ++-----
 ccdproc/tests/test_combiner.py |  77 +++++++++++++-
 ccdproc/tests/test_nanfuncs.py |  23 +++-
 6 files changed, 250 insertions(+), 119 deletions(-)

diff --git a/CHANGES.rst b/CHANGES.rst
index 2790d045..61144249 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -16,6 +16,11 @@ New Features
   classified differently from astropy (NumPy data still use astropy);
   ``'median'``/``'mean'``/``'std'``/``'mad_std'`` use the namespace's
   NaN-aware reductions or ccdproc's fallbacks. [#1001]
+- ``Combiner.sigma_clipping`` outside NumPy now accepts ``axis=None`` and a
+  tuple of axes, as ``astropy.stats.sigma_clip`` does on the NumPy path, and
+  the reduction fallbacks in ``ccdproc._nanfuncs`` gained the same axis
+  forms; the None/tuple axis handling formerly in ``_mad_fallback`` moved
+  into the shared ``_nanfuncs._setup``. [#1006]
 
 Other Changes and Additions
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^
diff --git a/ccdproc/_nanfuncs.py b/ccdproc/_nanfuncs.py
index e1f7a201..d258f7eb 100644
--- a/ccdproc/_nanfuncs.py
+++ b/ccdproc/_nanfuncs.py
@@ -24,6 +24,10 @@
 
 import array_api_compat
 
+# Host-side axis normalisation for tuple axes: operates on python ints
+# only, never on array data, so it does not tie the fallbacks to numpy.
+from numpy.lib.array_utils import normalize_axis_tuple
+
 __all__ = ["median", "nanmad", "nanmean", "nanmedian", "nanstd", "nansum"]
 
 
@@ -58,16 +62,24 @@ def _promote_to_real(x, xp, device):
 
 def _setup(x, axis, xp):
     """
-    Validate ``axis``, resolve the namespace and device, promote to float.
+    Normalise ``axis``, resolve the namespace and device, promote to float.
+
+    ``axis`` may be a single integer, ``None`` or a tuple/list of integers.
+    ``None`` flattens ``x`` so the caller reduces over everything; a tuple
+    or list moves the listed axes to the end and merges them into one, so
+    the caller's single-axis reduction reduces over all of them at once.
+    Either way the caller only ever sees a single non-negative integer
+    axis.
 
     Parameters
     ----------
     x : array
         Input array.
-    axis : int
-        Axis along which the caller will reduce. Booleans, ``None`` and
-        tuples of axes are rejected; anything else goes through
-        `operator.index`, so numpy integer scalars are accepted.
+    axis : int, tuple of int, list of int or None
+        Axis or axes along which the caller will reduce. Booleans are
+        rejected -- bool subclasses int, so ``axis=True`` would silently
+        mean axis 1 -- while numpy integer scalars are accepted. Negative
+        values count from the last axis.
     xp : array namespace or None
         Namespace to use. ``None`` resolves it from ``x``.
 
@@ -75,47 +87,82 @@ def _setup(x, axis, xp):
     -------
     x : array
         The input, promoted if necessary to the namespace's default real
-        floating dtype.
+        floating dtype, flattened when ``axis`` is ``None``, and with the
+        listed axes moved to the end and merged into one when ``axis`` is
+        a tuple or list.
     axis : int
-        The axis, normalised to a non-negative integer.
+        The single axis of the returned ``x`` to reduce, normalised to a
+        non-negative integer.
     xp : array namespace
         The resolved namespace.
     device : device
         The device ``x`` lives on.
+    restore : callable
+        Maps an array shaped like the returned ``x`` back to the layout of
+        the input ``x``; the identity for a single integer ``axis``.
+        Reductions remove the reduced axis and never need it;
+        ``ccdproc.combiner._sigma_clip_mask`` keeps the full shape and
+        uses it to hand its mask back in the caller's layout.
 
     Raises
     ------
-    NotImplementedError
-        If ``axis`` is not a single integer.
+    TypeError
+        If ``axis``, or an entry of a tuple/list ``axis``, is a bool or
+        not an integer.
     ValueError
-        If ``axis`` is out of bounds for ``x``.
+        If ``axis``, or an entry of a tuple/list ``axis``, is out of
+        bounds for ``x``, or a tuple/list names an axis more than once
+        (including via a negative alias).
     """
+    if xp is None:
+        xp = array_api_compat.array_namespace(x)
+    device = array_api_compat.device(x)
+    x = _promote_to_real(x, xp, device)
+    ndim = x.ndim
+
+    if axis is None:
+        shape = x.shape
+
+        def restore(a):
+            return xp.reshape(a, shape)
+
+        return xp.reshape(x, (-1,)), 0, xp, device, restore
+
+    if isinstance(axis, tuple | list):
+        # normalize_axis_tuple would silently treat True as 1.
+        if any(isinstance(ax, bool) for ax in axis):
+            raise TypeError("axis entries must be integers, not bool")
+        axes = normalize_axis_tuple(axis, ndim)
+        # Move the reduced axes to the end and merge them into one, so that
+        # a single-axis reduction reduces over all of them at once.
+        kept = tuple(ax for ax in range(ndim) if ax not in axes)
+        order = kept + axes
+        permuted_shape = tuple(x.shape[ax] for ax in order)
+        x = xp.reshape(xp.permute_dims(x, order), permuted_shape[: len(kept)] + (-1,))
+        inverse = tuple(order.index(ax) for ax in range(ndim))
+
+        def restore(a):
+            return xp.permute_dims(xp.reshape(a, permuted_shape), inverse)
+
+        return x, len(kept), xp, device, restore
+
     # bool subclasses int -- axis=True would silently mean axis 1 -- so it is
     # rejected explicitly, while operator.index accepts the numpy integer
     # scalars that isinstance(axis, int) would refuse.
-    if axis is None or isinstance(axis, bool):
-        raise NotImplementedError(
-            "NaN-aware reduction fallbacks support only a single integer axis."
-        )
+    if isinstance(axis, bool):
+        raise TypeError("axis must be an integer, not bool")
     try:
         axis = operator.index(axis)
     except TypeError:
-        raise NotImplementedError(
-            "NaN-aware reduction fallbacks support only a single integer axis."
+        raise TypeError(
+            f"axis must be an integer, a tuple or list of integers, or None, "
+            f"got {axis!r}"
         ) from None
 
-    if xp is None:
-        xp = array_api_compat.array_namespace(x)
-
-    ndim = x.ndim
     if not -ndim <= axis < ndim:
         raise ValueError(f"axis {axis} is out of bounds for array of dimension {ndim}")
-    axis = axis % ndim
-
-    device = array_api_compat.device(x)
-    x = _promote_to_real(x, xp, device)
 
-    return x, axis, xp, device
+    return x, axis % ndim, xp, device, lambda a: a
 
 
 def _sum_and_count(x, axis, xp, device, *, keepdims):
@@ -197,19 +244,20 @@ def nansum(x, /, *, axis=0, xp=None):
     x : array
         Input array. Integer and boolean inputs are promoted to the
         namespace's default real floating dtype.
-    axis : int, optional
-        Axis along which to sum. Default is 0. ``None`` and tuples of axes
-        are not supported.
+    axis : int, tuple of int, list of int or None, optional
+        Axis or axes along which to sum. Default is 0. ``None`` sums over
+        every axis; a tuple or list sums over all the listed axes at once.
     xp : array namespace, optional
         Namespace to use. Defaults to ``array_api_compat.array_namespace(x)``.
 
     Returns
     -------
     array
-        Sum of ``x`` along ``axis``, with that axis removed. Slices that are
-        entirely NaN sum to zero, matching `numpy.nansum`.
+        Sum of ``x`` along ``axis``, with the reduced axes removed (0-d
+        when ``axis`` is ``None``). Slices that are entirely NaN sum to
+        zero, matching `numpy.nansum`.
     """
-    x, axis, xp, device = _setup(x, axis, xp)
+    x, axis, xp, device, _ = _setup(x, axis, xp)
     total, _ = _sum_and_count(x, axis, xp, device, keepdims=False)
     return total
 
@@ -223,21 +271,22 @@ def nanmean(x, /, *, axis=0, xp=None):
     x : array
         Input array. Integer and boolean inputs are promoted to the
         namespace's default real floating dtype.
-    axis : int, optional
-        Axis along which to average. Default is 0. ``None`` and tuples of
-        axes are not supported.
+    axis : int, tuple of int, list of int or None, optional
+        Axis or axes along which to average. Default is 0. ``None``
+        averages over every axis; a tuple or list over all the listed axes.
     xp : array namespace, optional
         Namespace to use. Defaults to ``array_api_compat.array_namespace(x)``.
 
     Returns
     -------
     array
-        Mean of ``x`` along ``axis``, with that axis removed. Slices that
-        are entirely NaN yield NaN silently, matching ``bottleneck.nanmean``
+        Mean of ``x`` along ``axis``, with the reduced axes removed (0-d
+        when ``axis`` is ``None``). Slices that are entirely NaN yield
+        NaN silently, matching ``bottleneck.nanmean``
         (the numpy-backend default); `numpy.nanmean` warns here, but a fully
         masked pixel is a routine input for the combiner, not an anomaly.
     """
-    x, axis, xp, device = _setup(x, axis, xp)
+    x, axis, xp, device, _ = _setup(x, axis, xp)
     total, count = _sum_and_count(x, axis, xp, device, keepdims=False)
     return _safe_divide(total, count, xp, device)
 
@@ -255,16 +304,18 @@ def nanstd(x, /, *, axis=0, xp=None):
     x : array
         Input array. Integer and boolean inputs are promoted to the
         namespace's default real floating dtype.
-    axis : int, optional
-        Axis along which to compute the deviation. Default is 0. ``None``
-        and tuples of axes are not supported.
+    axis : int, tuple of int, list of int or None, optional
+        Axis or axes along which to compute the deviation. Default is 0.
+        ``None`` reduces over every axis; a tuple or list over all the
+        listed axes.
     xp : array namespace, optional
         Namespace to use. Defaults to ``array_api_compat.array_namespace(x)``.
 
     Returns
     -------
     array
-        Standard deviation of ``x`` along ``axis``, with that axis removed.
+        Standard deviation of ``x`` along ``axis``, with the reduced axes
+        removed (0-d when ``axis`` is ``None``).
         Slices that are entirely NaN yield NaN silently, matching
         ``bottleneck.nanstd`` (the numpy-backend default); `numpy.nanstd`
         warns here, but a fully masked pixel is a routine input for the
@@ -279,7 +330,7 @@ def nanstd(x, /, *, axis=0, xp=None):
     single-pass form suffers when the values are large relative to their
     spread, which is not unusual for CCD counts.
     """
-    x, axis, xp, device = _setup(x, axis, xp)
+    x, axis, xp, device, _ = _setup(x, axis, xp)
 
     isnan = xp.isnan(x)
     zero = xp.asarray(0, dtype=x.dtype, device=device)
@@ -314,18 +365,20 @@ def nanmedian(x, /, *, axis=0, xp=None):
     x : array
         Input array. Integer and boolean inputs are promoted to the
         namespace's default real floating dtype.
-    axis : int, optional
-        Axis along which to compute the median. Default is 0. Booleans,
-        ``None`` and tuples of axes are not supported; numpy integer
-        scalars are accepted.
+    axis : int, tuple of int, list of int or None, optional
+        Axis or axes along which to compute the median. Default is 0.
+        ``None`` reduces over every axis and a tuple or list over all the
+        listed axes; booleans are rejected, numpy integer scalars are
+        accepted.
     xp : array namespace, optional
         Namespace to use. Defaults to ``array_api_compat.array_namespace(x)``.
 
     Returns
     -------
     array
-        Median of ``x`` along ``axis``, with that axis removed. Slices that
-        are entirely NaN yield NaN silently, matching
+        Median of ``x`` along ``axis``, with the reduced axes removed (0-d
+        when ``axis`` is ``None``). Slices that are entirely NaN yield NaN
+        silently, matching
         ``bottleneck.nanmedian`` (the numpy-backend default);
         `numpy.nanmedian` warns here, but a fully masked pixel is a routine
         input for the combiner, not an anomaly.
@@ -337,7 +390,7 @@ def nanmedian(x, /, *, axis=0, xp=None):
     or ``bottleneck.nanmedian``. Prefer a native ``nanmedian`` when the
     namespace offers one.
     """
-    x, axis, xp, device = _setup(x, axis, xp)
+    x, axis, xp, device, _ = _setup(x, axis, xp)
     ndim = x.ndim
 
     # Replacing NaNs with +inf keeps them past every real value regardless of
@@ -385,18 +438,20 @@ def median(x, /, *, axis=0, xp=None):
     x : array
         Input array. Integer and boolean inputs are promoted to the
         namespace's default real floating dtype.
-    axis : int, optional
-        Axis along which to compute the median. Default is 0. Booleans,
-        ``None`` and tuples of axes are not supported; numpy integer
-        scalars are accepted.
+    axis : int, tuple of int, list of int or None, optional
+        Axis or axes along which to compute the median. Default is 0.
+        ``None`` reduces over every axis and a tuple or list over all the
+        listed axes; booleans are rejected, numpy integer scalars are
+        accepted.
     xp : array namespace, optional
         Namespace to use. Defaults to ``array_api_compat.array_namespace(x)``.
 
     Returns
     -------
     array
-        Median of ``x`` along ``axis``, with that axis removed. Slices that
-        contain any NaN yield NaN, matching `numpy.median`; this is the
+        Median of ``x`` along ``axis``, with the reduced axes removed (0-d
+        when ``axis`` is ``None``). Slices that contain any NaN yield NaN,
+        matching `numpy.median`; this is the
         difference from `nanmedian`, which ignores NaNs entirely.
 
     Notes
@@ -409,7 +464,7 @@ def median(x, /, *, axis=0, xp=None):
     with a final `where` over whether any NaN is present along ``axis``,
     since `nanmedian` alone would silently drop NaNs instead.
     """
-    x, axis, xp, device = _setup(x, axis, xp)
+    x, axis, xp, device, _ = _setup(x, axis, xp)
     nan = xp.asarray(xp.nan, dtype=x.dtype, device=device)
     return xp.where(xp.any(xp.isnan(x), axis=axis), nan, nanmedian(x, axis=axis, xp=xp))
 
@@ -423,14 +478,17 @@ def nanmad(x, /, *, axis=0, xp=None, median=None):
     x : array
         Input array. Integer and boolean inputs are promoted to the
         namespace's default real floating dtype.
-    axis : int, optional
-        Axis along which to compute the deviation. Default is 0. Booleans,
-        ``None`` and tuples of axes are not supported; numpy integer
-        scalars are accepted.
+    axis : int, tuple of int, list of int or None, optional
+        Axis or axes along which to compute the deviation. Default is 0.
+        ``None`` reduces over every axis and a tuple or list over all the
+        listed axes; booleans are rejected, numpy integer scalars are
+        accepted.
     xp : array namespace, optional
         Namespace to use. Defaults to ``array_api_compat.array_namespace(x)``.
     median : callable, optional
-        Reduction used for both medians, called as ``median(x, axis=axis)``.
+        Reduction used for both medians, called as ``median(x, axis=axis)``,
+        always with a single integer ``axis``: a ``None`` or tuple/list
+        ``axis`` has already been flattened or merged away by `_setup`.
         Default is `nanmedian`. A keyword rather than a module-level tier
         (as `ccdproc.combiner._default_median` provides) so this module has
         no dependency on `ccdproc.combiner`.
@@ -438,11 +496,12 @@ def nanmad(x, /, *, axis=0, xp=None, median=None):
     Returns
     -------
     array
-        ``median(|x - median(x)|)`` along ``axis``, with that axis removed.
+        ``median(|x - median(x)|)`` along ``axis``, with the reduced axes
+        removed (0-d when ``axis`` is ``None``).
         Unscaled: multiply by ``1.482602218505602`` for an estimate of the
         standard deviation, as `astropy.stats.mad_std` does.
     """
-    x, axis, xp, device = _setup(x, axis, xp)
+    x, axis, xp, device, _ = _setup(x, axis, xp)
     if median is None:
         median = partial(nanmedian, xp=xp)
     center = xp.expand_dims(median(x, axis=axis), axis=axis)
diff --git a/ccdproc/combiner.py b/ccdproc/combiner.py
index 6086b4b5..fab3cf34 100644
--- a/ccdproc/combiner.py
+++ b/ccdproc/combiner.py
@@ -195,15 +195,17 @@ def _sigma_clip_mask(
         Number of deviations below/above the center beyond which a value is
         clipped. As in astropy, ``None`` and ``0`` mean ``3``. Default is
         ``3``.
-    axis : int, optional
-        Axis along which the center and deviation are computed. Only a single
-        integer axis is supported. Default is ``0``.
+    axis : int, tuple of int, list of int or None, optional
+        Axis or axes along which the center and deviation are computed.
+        ``None`` computes them over the whole array. Default is ``0``.
     maxiters : int or None, optional
         Number of clipping iterations. ``None`` iterates until no more values
         are rejected. Default is ``1``.
     cenfunc : {'median', 'mean'} or callable, optional
         Statistic for the center. A callable must accept ``(data, axis=axis)``
-        and ignore NaNs. Default is ``'median'``.
+        and ignore NaNs; it is always called with a single integer axis (a
+        ``None`` or tuple/list ``axis`` is flattened or merged away first).
+        Default is ``'median'``.
     stdfunc : {'std', 'mad_std'} or callable, optional
         Statistic for the deviation, same requirements as ``cenfunc``.
         Default is ``'std'``.
@@ -213,16 +215,17 @@ def _sigma_clip_mask(
     Returns
     -------
     array of bool
-        Mask of the clipped values, ``True`` where a value is rejected, in
-        the namespace and on the device of ``data``.
+        Mask of the clipped values, ``True`` where a value is rejected,
+        shaped like ``data``, in its namespace and on its device.
 
     Raises
     ------
-    NotImplementedError
-        If ``axis`` is not a single integer.
+    TypeError
+        If ``axis`` is a bool or is not made of integers.
     ValueError
-        If ``axis`` is out of bounds, ``maxiters`` is not positive, or a
-        string ``cenfunc``/``stdfunc`` is not one of the options.
+        If ``axis`` is out of bounds or repeats an axis, ``maxiters`` is
+        not positive, or a string ``cenfunc``/``stdfunc`` is not one of
+        the options.
 
     Notes
     -----
@@ -258,7 +261,7 @@ def _sigma_clip_mask(
     iteration whether anything was rejected, which forces a compute per
     iteration on such backends; prefer an integer there.
     """
-    data, axis, xp, device = _setup(data, axis, xp)
+    data, axis, xp, device, restore = _setup(data, axis, xp)
 
     if isinstance(cenfunc, str) and isinstance(stdfunc, str):
         # B1: widen (never narrow) to the namespace's default real floating
@@ -311,7 +314,9 @@ def _sigma_clip_mask(
             break
         filtered = xp.where(rejected, nan, filtered)
 
-    return invalid | (data < lower) | (data > upper)
+    # ``restore`` maps the mask back to the caller's layout when _setup
+    # flattened (``axis=None``) or merged (tuple ``axis``) the data.
+    return restore(invalid | (data < lower) | (data > upper))
 
 
 class Combiner:
@@ -699,11 +704,11 @@ def sigma_clipping(
 
         kwd
             ``axis`` (default ``0``) and ``maxiters`` (default ``1``) are
-            honoured for every array namespace. Outside numpy ``axis`` must
-            be a single integer: ``None`` and a tuple of axes, both of
-            which :func:`~astropy.stats.sigma_clip` accepts, raise
-            `NotImplementedError` there. ``masked`` and ``return_bounds``
-            are never accepted, on any array namespace -- this method
+            honoured for every array namespace; like
+            :func:`~astropy.stats.sigma_clip`, ``axis`` may be a single
+            integer, ``None`` or a tuple of axes on every one. ``masked``
+            and ``return_bounds`` are never accepted, on any array
+            namespace -- this method
             always asks astropy for the mask itself -- and raise
             `TypeError`. ``copy`` (default `False`) and any other
             astropy-only keyword argument, such as ``grow``, are passed to
diff --git a/ccdproc/core.py b/ccdproc/core.py
index ad5e7e85..13cd601e 100644
--- a/ccdproc/core.py
+++ b/ccdproc/core.py
@@ -21,7 +21,6 @@
 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
@@ -341,29 +340,6 @@ def _mad_fallback(data, axis, ignore_nan, xp=None):
     without the ``RuntimeWarning`` numpy emits.
     """
     xp = xp or array_api_compat.array_namespace(data)
-    device = array_api_compat.device(data)
-
-    # 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.
-    data = _promote_to_real(data, xp, device)
-
-    if axis is None:
-        data = xp.reshape(data, (-1,))
-        axis = 0
-    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 = 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 + list(axes)))
-        data = xp.reshape(
-            data, tuple(data.shape[ax] for ax in range(len(kept))) + (-1,)
-        )
-        axis = -1
 
     def med(d, axis):
         try:
@@ -374,6 +350,10 @@ def med(d, axis):
 
     # The deviation itself is _nanfuncs.nanmad, the same computation
     # Combiner._nanmadstd uses; only the median tier passed in differs.
+    # nanmad's _setup owns the axis handling (``None`` flattens, a tuple
+    # or list is merged into a single trailing axis) and the promotion to
+    # the namespace's default real floating dtype, so ``med`` only ever
+    # sees a single integer axis.
     return _nanfuncs_nanmad(data, axis=axis, xp=xp, median=med)
 
 
@@ -1034,8 +1014,7 @@ def subtract_dark(
     """
     if ccd.shape != master.shape:
         err_str = (
-            f"operands could not be subtracted with "
-            f"shapes {ccd.shape} {master.shape}"
+            f"operands could not be subtracted with shapes {ccd.shape} {master.shape}"
         )
         raise ValueError(err_str)
 
@@ -1749,7 +1728,7 @@ def rebin(ccd, newshape):
     else:
         # check to see that the two arrays are going to be the same length
         if len(ccd.shape) != len(newshape):
-            raise ValueError("newshape does not have the same dimensions as " "ccd.")
+            raise ValueError("newshape does not have the same dimensions as ccd.")
 
         slices = [
             slice(0, old, old / new)
@@ -2083,8 +2062,7 @@ def cosmicray_lacosmic(
             s = "s" if len(bad_args) > 1 else ""
             bads = ", ".join(bad_args)
             raise TypeError(
-                f"The argument{s} {bads} only valid for astroscrappy "
-                "1.1.0 or higher."
+                f"The argument{s} {bads} only valid for astroscrappy 1.1.0 or higher."
             )
 
         if pssl != 0:
@@ -2772,13 +2750,13 @@ def value(self, value):
             self._value = value
         elif isinstance(value, str):
             if self.unit is not None:
-                raise ValueError("keyword with a unit cannot have a " "string value.")
+                raise ValueError("keyword with a unit cannot have a string value.")
             else:
                 self._value = value
         else:
             if self.unit is None:
                 raise ValueError(
-                    "no unit provided. Set value with " "an astropy.units.Quantity."
+                    "no unit provided. Set value with an astropy.units.Quantity."
                 )
             self._value = value * self.unit
 
diff --git a/ccdproc/tests/test_combiner.py b/ccdproc/tests/test_combiner.py
index bc832f87..c5b82af1 100644
--- a/ccdproc/tests/test_combiner.py
+++ b/ccdproc/tests/test_combiner.py
@@ -1755,10 +1755,17 @@ def _sigma_clip_reference(np_data, **kwargs):
     _, lower, upper = sigma_clip(
         np_data.copy(), masked=False, return_bounds=True, **kwargs
     )
-    # The compiled path drops the clipped axis from the bounds while the
-    # python loop keeps it with length one; either way, make them broadcast.
-    axis = kwargs.get("axis", 0) % np_data.ndim
-    shape = tuple(1 if dim == axis else n for dim, n in enumerate(np_data.shape))
+    # The compiled path drops the clipped axes from the bounds while the
+    # python loop keeps them with length one; either way, make them
+    # broadcast. ``axis`` may be an int, a tuple of ints or None here.
+    axis = kwargs.get("axis", 0)
+    if axis is None:
+        axes = tuple(range(np_data.ndim))
+    elif isinstance(axis, tuple):
+        axes = tuple(ax % np_data.ndim for ax in axis)
+    else:
+        axes = (axis % np_data.ndim,)
+    shape = tuple(1 if dim in axes else n for dim, n in enumerate(np_data.shape))
     lower = np.reshape(lower, shape)
     upper = np.reshape(upper, shape)
     with np.errstate(invalid="ignore"):
@@ -1891,6 +1898,38 @@ def test_sigma_clip_mask_argument_handling():
         _sigma_clip_mask(data, stdfunc="var", xp=xp)
 
 
+# Axis forms beyond a single integer: _nanfuncs._setup flattens the data
+# for ``axis=None`` and merges a tuple of axes into one, and the mask must
+# come back in the shape and layout of the input on every backend.
+@pytest.mark.parametrize("axis", [None, (0, 1), (1, 2), (0, -1)], ids=str)
+def test_sigma_clip_mask_axis_forms(axis):
+    data = xp.asarray(_sigma_clip_datasets()["normal"], device=xp_device)
+    np_data = _to_numpy(data)
+
+    result = _sigma_clip_mask(
+        data, sigma_lower=2, sigma_upper=2, axis=axis, maxiters=2, xp=xp
+    )
+
+    # The reference gets a tuple axis with its negative entries normalised:
+    # astropy's bottleneck dispatch transposes with the tuple as given and
+    # raises on a negative entry. The helper receives the tuple as written.
+    ref_axis = axis
+    if isinstance(axis, tuple):
+        ref_axis = tuple(ax % np_data.ndim for ax in axis)
+    expected = _sigma_clip_reference(
+        np_data,
+        sigma_lower=2,
+        sigma_upper=2,
+        axis=ref_axis,
+        maxiters=2,
+        cenfunc="median",
+        stdfunc="std",
+    )
+
+    assert result.shape == data.shape
+    assert bool(xp.all(result == xp.asarray(expected, device=xp_device)))
+
+
 def _sigma_clip_ccd_list():
     """
     The ``normal`` set of `_sigma_clip_datasets` as a list of `CCDData`,
@@ -2008,3 +2047,33 @@ def test_combine_sigma_clip_on_any_backend():
     expected = np.ma.average(np.ma.array(np_data, mask=mask), axis=0)
     assert_allclose(_to_numpy(result.data), expected)
     assert not np.allclose(_to_numpy(result.data), np_data.mean(axis=0))
+
+
+@pytest.mark.filterwarnings("ignore::astropy.utils.exceptions.AstropyUserWarning")
+@pytest.mark.parametrize("axis", [None, (1, 2)], ids=str)
+def test_sigma_clipping_axis_forms_any_backend(axis):
+    # astropy's sigma_clip accepts axis=None and a tuple of axes;
+    # Combiner.sigma_clipping must honour them off the numpy path too,
+    # with the mask coming back in the data's own shape.
+    c = Combiner(_sigma_clip_ccd_list())
+    c.sigma_clipping(
+        low_thresh=2,
+        high_thresh=2,
+        func="median",
+        dev_func="std",
+        axis=axis,
+        maxiters=2,
+    )
+
+    expected = _sigma_clip_reference(
+        _to_numpy(c._data_arr),
+        sigma_lower=2,
+        sigma_upper=2,
+        axis=axis,
+        maxiters=2,
+        cenfunc="median",
+        stdfunc="std",
+    )
+
+    assert c._data_arr_mask.shape == c._data_arr.shape
+    assert bool(xp.all(c._data_arr_mask == xp.asarray(expected, device=xp_device)))
diff --git a/ccdproc/tests/test_nanfuncs.py b/ccdproc/tests/test_nanfuncs.py
index 5536c7b7..7d2136aa 100644
--- a/ccdproc/tests/test_nanfuncs.py
+++ b/ccdproc/tests/test_nanfuncs.py
@@ -38,6 +38,11 @@
     (_some_nan, 1),  # a non-zero axis
     (_some_nan, -1),  # a negative axis
     (_some_nan, np.int64(1)),  # a numpy integer axis
+    (_some_nan, None),  # reduce over everything
+    (_some_nan, (0, 1)),  # a tuple covering every axis
+    (_rng.normal(size=(5, 4, 3)), (0, 2)),  # a tuple of axes
+    (_rng.normal(size=(5, 4, 3)), (-1, 0)),  # a negative entry in a tuple
+    (_rng.normal(size=(5, 4, 3)), (1,)),  # a single-entry tuple
     (np.array([[1.0, np.nan], [2.0, np.nan], [3.0, np.nan]]), 0),  # all-NaN column
     (np.array([[1.0, np.nan], [np.nan, np.nan]]), 0),  # single non-NaN in a slice
     (np.array([np.nan, np.nan, np.nan]), 0),  # every value NaN
@@ -106,7 +111,7 @@ def test_nansum_all_nan_slice_is_zero():
     assert xp.all(xpx.isclose(result, xp.asarray([3.0, 0.0], device=xp_device)))
 
 
-@pytest.mark.parametrize("axis", [0, 1, -1])
+@pytest.mark.parametrize("axis", [0, 1, -1, (0, 1), None], ids=str)
 def test_nanmad_matches_astropy(axis):
     """``nanmad`` reproduces ``median_absolute_deviation(ignore_nan=True)``."""
     result = nanmad(xp.asarray(_some_nan, device=xp_device), axis=axis)
@@ -118,15 +123,25 @@ def test_nanmad_matches_astropy(axis):
     assert bool(xp.all(xpx.isclose(result, expected, equal_nan=True)))
 
 
+@pytest.mark.parametrize("func", [nansum, nanmean, nanstd, nanmedian, median])
+def test_list_axis_matches_tuple(func):
+    """A list of axes means the same as a tuple (numpy itself rejects it)."""
+    data = xp.asarray(_rng.normal(size=(4, 3, 2)), device=xp_device)
+    assert bool(xp.all(func(data, axis=[0, 2]) == func(data, axis=(0, 2))))
+
+
 @pytest.mark.parametrize("func", [nansum, nanmean, nanstd, nanmedian, median])
 @pytest.mark.parametrize(
     ("axis", "error"),
     [
-        (None, NotImplementedError),
-        (True, NotImplementedError),  # bool subclasses int; reject it anyway
-        ((0, 1), NotImplementedError),
+        (True, TypeError),  # bool subclasses int; reject it anyway
+        (1.5, TypeError),
         (2, ValueError),
         (-3, ValueError),
+        ((0, 0), ValueError),  # repeated axis
+        ((0, -2), ValueError),  # repeated via a negative alias
+        ((0, 2), ValueError),  # out-of-bounds entry
+        ((0, True), TypeError),  # bool entry in a tuple
     ],
 )
 def test_bad_axis(func, axis, error):

From c801431d10031280443b86397a326b80a376e7cc Mon Sep 17 00:00:00 2001
From: Matt Craig <mattwcraig@gmail.com>
Date: Mon, 31 Aug 2026 08:50:15 -0500
Subject: [PATCH 2/3] Address review on #1006: fix two _setup edge cases,
 simplify, document

Bug fixes, both found in review:

- Reject np.bool_ axes explicitly, as scalars and as tuple/list entries:
  on the oldest supported numpy (2.0) operator.index still converts
  np.True_ to 1 with only a DeprecationWarning, so axis=(0, np.True_)
  silently reduced axes (0, 1) there. Newer numpy raises TypeError
  itself; the guard makes the behaviour and message uniform.
- Spell out the merged-axis length in _setup's tuple merge instead of
  reshaping with -1, which is ambiguous (and raises) when a kept axis
  has size 0; numpy's reductions return an empty result there. Found
  independently by two reviewers.

Simplifications from review, all behaviour-preserving:

- Fold the axis=None branch into the tuple branch (None is
  tuple(range(ndim))).
- Replace the hand-rolled scalar bounds check with
  normalize_axis_tuple, whose AxisError is a ValueError with an
  identical message; the tuple branch already surfaced AxisError.
- Normalize tuple axes once inside test_combiner's
  _sigma_clip_reference instead of at its call site.
- Shrink test_ccdproc's bad-axis grid to a single delegation check;
  the messages are now pinned in test_nanfuncs::test_bad_axis.

Docs and tests: move _setup's axis-handling discussion into a Notes
section, add brief comments to the permute-and-merge block, explain the
list-equals-tuple test, document that sigma_clipping takes a list of
axes too, note the NotImplementedError->TypeError change in CHANGES,
switch the diff's British spellings to American, and drop the unrelated
string-concatenation cleanups from core.py. New differential rows cover
axis=(), numpy ints in a tuple, and the two empty-array shapes (where
numpy.median itself trips over the same -1 merge; the test substitutes
np.nanmedian as the reference there).

Full runs of the three affected files: numpy 697, dask 701, jax full
suite 968, array-api-strict 689 + 12 pre-existing xfails; ruff and
black clean.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01LVxnTLrWKStBdPwcmDNhxA
---
 CHANGES.rst                    | 10 +++--
 ccdproc/_nanfuncs.py           | 81 +++++++++++++++++++++-------------
 ccdproc/combiner.py            |  4 +-
 ccdproc/core.py                | 12 ++---
 ccdproc/tests/test_ccdproc.py  | 27 ++++--------
 ccdproc/tests/test_combiner.py | 19 ++++----
 ccdproc/tests/test_nanfuncs.py | 62 ++++++++++++++++++++------
 7 files changed, 132 insertions(+), 83 deletions(-)

diff --git a/CHANGES.rst b/CHANGES.rst
index 61144249..26397d55 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -17,10 +17,12 @@ New Features
   ``'median'``/``'mean'``/``'std'``/``'mad_std'`` use the namespace's
   NaN-aware reductions or ccdproc's fallbacks. [#1001]
 - ``Combiner.sigma_clipping`` outside NumPy now accepts ``axis=None`` and a
-  tuple of axes, as ``astropy.stats.sigma_clip`` does on the NumPy path, and
-  the reduction fallbacks in ``ccdproc._nanfuncs`` gained the same axis
-  forms; the None/tuple axis handling formerly in ``_mad_fallback`` moved
-  into the shared ``_nanfuncs._setup``. [#1006]
+  tuple or list of axes, as ``astropy.stats.sigma_clip`` does on the NumPy
+  path, and the reduction fallbacks in ``ccdproc._nanfuncs`` gained the same
+  axis forms; the None/tuple axis handling formerly in ``_mad_fallback``
+  moved into the shared ``_nanfuncs._setup``, and a bool or otherwise
+  non-integer ``axis`` now raises ``TypeError`` rather than
+  ``NotImplementedError``. [#1006]
 
 Other Changes and Additions
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^
diff --git a/ccdproc/_nanfuncs.py b/ccdproc/_nanfuncs.py
index d258f7eb..7dd4bbb4 100644
--- a/ccdproc/_nanfuncs.py
+++ b/ccdproc/_nanfuncs.py
@@ -19,13 +19,16 @@
 keeps the five functions consistent with each other.
 """
 
+import math
 import operator
 from functools import partial
 
 import array_api_compat
 
-# Host-side axis normalisation for tuple axes: operates on python ints
-# only, never on array data, so it does not tie the fallbacks to numpy.
+# Host-side axis handling: normalize_axis_tuple operates on python ints
+# only, never on array data, and np.bool_ appears only in the guards that
+# reject boolean axes, so neither ties the fallbacks to numpy.
+import numpy as np
 from numpy.lib.array_utils import normalize_axis_tuple
 
 __all__ = ["median", "nanmad", "nanmean", "nanmedian", "nanstd", "nansum"]
@@ -62,14 +65,7 @@ def _promote_to_real(x, xp, device):
 
 def _setup(x, axis, xp):
     """
-    Normalise ``axis``, resolve the namespace and device, promote to float.
-
-    ``axis`` may be a single integer, ``None`` or a tuple/list of integers.
-    ``None`` flattens ``x`` so the caller reduces over everything; a tuple
-    or list moves the listed axes to the end and merges them into one, so
-    the caller's single-axis reduction reduces over all of them at once.
-    Either way the caller only ever sees a single non-negative integer
-    axis.
+    Normalize ``axis``, resolve the namespace and device, promote to float.
 
     Parameters
     ----------
@@ -91,7 +87,7 @@ def _setup(x, axis, xp):
         listed axes moved to the end and merged into one when ``axis`` is
         a tuple or list.
     axis : int
-        The single axis of the returned ``x`` to reduce, normalised to a
+        The single axis of the returned ``x`` to reduce, normalized to a
         non-negative integer.
     xp : array namespace
         The resolved namespace.
@@ -113,6 +109,15 @@ def _setup(x, axis, xp):
         If ``axis``, or an entry of a tuple/list ``axis``, is out of
         bounds for ``x``, or a tuple/list names an axis more than once
         (including via a negative alias).
+
+    Notes
+    -----
+    ``axis`` may be a single integer, ``None`` or a tuple/list of integers.
+    ``None`` flattens ``x`` so the caller reduces over everything; a tuple
+    or list moves the listed axes to the end and merges them into one, so
+    the caller's single-axis reduction reduces over all of them at once.
+    Either way the caller only ever sees a single non-negative integer
+    axis.
     """
     if xp is None:
         xp = array_api_compat.array_namespace(x)
@@ -121,24 +126,37 @@ def _setup(x, axis, xp):
     ndim = x.ndim
 
     if axis is None:
-        shape = x.shape
-
-        def restore(a):
-            return xp.reshape(a, shape)
-
-        return xp.reshape(x, (-1,)), 0, xp, device, restore
+        # Reducing over everything is the same as naming every axis.
+        axis = tuple(range(ndim))
 
     if isinstance(axis, tuple | list):
-        # normalize_axis_tuple would silently treat True as 1.
-        if any(isinstance(ax, bool) for ax in axis):
+        # normalize_axis_tuple would treat a bool as an axis: operator.index
+        # turns True into 1, and on the oldest supported numpy (2.0) it
+        # still accepts np.bool_ too, with only a DeprecationWarning.
+        if any(isinstance(ax, bool | np.bool_) for ax in axis):
             raise TypeError("axis entries must be integers, not bool")
+        # Host-side validation and normalization in one call: entries go
+        # through operator.index, negatives are wrapped mod ndim,
+        # out-of-bounds raises AxisError, and a duplicate (even via a
+        # negative alias) raises ValueError.
         axes = normalize_axis_tuple(axis, ndim)
-        # Move the reduced axes to the end and merge them into one, so that
-        # a single-axis reduction reduces over all of them at once.
+        # Move the reduced axes to the end and merge them into one trailing
+        # axis, so that a single-axis reduction reduces over all of them at
+        # once. How the merge interleaves elements is irrelevant: every
+        # reduction here is order-insensitive within the reduced set.
         kept = tuple(ax for ax in range(ndim) if ax not in axes)
         order = kept + axes
         permuted_shape = tuple(x.shape[ax] for ax in order)
-        x = xp.reshape(xp.permute_dims(x, order), permuted_shape[: len(kept)] + (-1,))
+        # The merged length is spelled out because reshape cannot infer it
+        # from -1 when a kept axis has size 0 (total size 0 is ambiguous);
+        # numpy returns an empty result there, and so does this.
+        merged = math.prod(permuted_shape[len(kept) :])
+        x = xp.reshape(
+            xp.permute_dims(x, order), permuted_shape[: len(kept)] + (merged,)
+        )
+        # ``inverse`` undoes ``order``; ``restore`` maps a full-shape array
+        # in the permuted-merged layout back to the caller's layout by
+        # un-merging (reshape) and un-permuting.
         inverse = tuple(order.index(ax) for ax in range(ndim))
 
         def restore(a):
@@ -146,10 +164,11 @@ def restore(a):
 
         return x, len(kept), xp, device, restore
 
-    # bool subclasses int -- axis=True would silently mean axis 1 -- so it is
-    # rejected explicitly, while operator.index accepts the numpy integer
-    # scalars that isinstance(axis, int) would refuse.
-    if isinstance(axis, bool):
+    # bool subclasses int -- axis=True would silently mean axis 1 -- and on
+    # numpy 2.0 operator.index still accepts np.bool_ as well, so both are
+    # rejected explicitly, while numpy integer scalars (which
+    # isinstance(axis, int) would refuse) are accepted.
+    if isinstance(axis, bool | np.bool_):
         raise TypeError("axis must be an integer, not bool")
     try:
         axis = operator.index(axis)
@@ -159,10 +178,10 @@ def restore(a):
             f"got {axis!r}"
         ) from None
 
-    if not -ndim <= axis < ndim:
-        raise ValueError(f"axis {axis} is out of bounds for array of dimension {ndim}")
-
-    return x, axis % ndim, xp, device, lambda a: a
+    # normalize_axis_tuple wraps a negative axis and raises AxisError -- a
+    # ValueError subclass with numpy's own message -- when it is out of
+    # bounds, exactly as the tuple branch above does for entries.
+    return x, normalize_axis_tuple(axis, ndim)[0], xp, device, lambda a: a
 
 
 def _sum_and_count(x, axis, xp, device, *, keepdims):
@@ -174,7 +193,7 @@ def _sum_and_count(x, axis, xp, device, *, keepdims):
     x : array
         Input array, already promoted to a real floating dtype.
     axis : int
-        Axis to reduce, already normalised to a non-negative integer.
+        Axis to reduce, already normalized to a non-negative integer.
     xp : array namespace
         Namespace to use.
     device : device
diff --git a/ccdproc/combiner.py b/ccdproc/combiner.py
index fab3cf34..e0c64824 100644
--- a/ccdproc/combiner.py
+++ b/ccdproc/combiner.py
@@ -704,9 +704,9 @@ def sigma_clipping(
 
         kwd
             ``axis`` (default ``0``) and ``maxiters`` (default ``1``) are
-            honoured for every array namespace; like
+            honored for every array namespace; like
             :func:`~astropy.stats.sigma_clip`, ``axis`` may be a single
-            integer, ``None`` or a tuple of axes on every one. ``masked``
+            integer, ``None`` or a tuple or list of axes on every one. ``masked``
             and ``return_bounds`` are never accepted, on any array
             namespace -- this method
             always asks astropy for the mask itself -- and raise
diff --git a/ccdproc/core.py b/ccdproc/core.py
index 13cd601e..33b08b3c 100644
--- a/ccdproc/core.py
+++ b/ccdproc/core.py
@@ -1014,7 +1014,8 @@ def subtract_dark(
     """
     if ccd.shape != master.shape:
         err_str = (
-            f"operands could not be subtracted with shapes {ccd.shape} {master.shape}"
+            f"operands could not be subtracted with "
+            f"shapes {ccd.shape} {master.shape}"
         )
         raise ValueError(err_str)
 
@@ -1728,7 +1729,7 @@ def rebin(ccd, newshape):
     else:
         # check to see that the two arrays are going to be the same length
         if len(ccd.shape) != len(newshape):
-            raise ValueError("newshape does not have the same dimensions as ccd.")
+            raise ValueError("newshape does not have the same dimensions as " "ccd.")
 
         slices = [
             slice(0, old, old / new)
@@ -2062,7 +2063,8 @@ def cosmicray_lacosmic(
             s = "s" if len(bad_args) > 1 else ""
             bads = ", ".join(bad_args)
             raise TypeError(
-                f"The argument{s} {bads} only valid for astroscrappy 1.1.0 or higher."
+                f"The argument{s} {bads} only valid for astroscrappy "
+                "1.1.0 or higher."
             )
 
         if pssl != 0:
@@ -2750,13 +2752,13 @@ def value(self, value):
             self._value = value
         elif isinstance(value, str):
             if self.unit is not None:
-                raise ValueError("keyword with a unit cannot have a string value.")
+                raise ValueError("keyword with a unit cannot have a " "string value.")
             else:
                 self._value = value
         else:
             if self.unit is None:
                 raise ValueError(
-                    "no unit provided. Set value with an astropy.units.Quantity."
+                    "no unit provided. Set value with " "an astropy.units.Quantity."
                 )
             self._value = value * self.unit
 
diff --git a/ccdproc/tests/test_ccdproc.py b/ccdproc/tests/test_ccdproc.py
index 37f20f49..4af18e37 100644
--- a/ccdproc/tests/test_ccdproc.py
+++ b/ccdproc/tests/test_ccdproc.py
@@ -445,25 +445,16 @@ 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), out-of-bounds and bool entries in
-# a tuple axis raise instead of silently reducing the wrong axes.
-@pytest.mark.parametrize(
-    ("axis", "error", "match"),
-    [
-        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), 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, error, match):
+def test_mad_fallback_bad_axis_propagates():
+    """
+    ``_mad_fallback`` delegates axis validation to ``_nanfuncs._setup``
+    (via ``nanmad``) and lets its error out unchanged; the full grid of
+    bad axis forms, messages included, lives in
+    ``test_nanfuncs.py::test_bad_axis``.
+    """
     data = xp.asarray(_MAD_3D, device=xp_device)
-    with pytest.raises(error, match=match):
-        _mad_fallback(data, axis, True)
+    with pytest.raises(ValueError, match="repeated axis"):
+        _mad_fallback(data, (0, 0), True)
 
 
 # The except branch of _mad_fallback's med() only runs naturally on
diff --git a/ccdproc/tests/test_combiner.py b/ccdproc/tests/test_combiner.py
index c5b82af1..bc0b3403 100644
--- a/ccdproc/tests/test_combiner.py
+++ b/ccdproc/tests/test_combiner.py
@@ -1752,17 +1752,22 @@ def _sigma_clip_reference(np_data, **kwargs):
     checked against it too: that is what numpy data get from
     ``Combiner.sigma_clipping``, and the other backends must agree with it.
     """
+    # Normalize a tuple axis once, up front: astropy's bottleneck dispatch
+    # cannot take negative tuple entries, and the bounds shape below needs
+    # the non-negative values too.
+    axis = kwargs.get("axis", 0)
+    if isinstance(axis, tuple):
+        kwargs["axis"] = axis = tuple(ax % np_data.ndim for ax in axis)
     _, lower, upper = sigma_clip(
         np_data.copy(), masked=False, return_bounds=True, **kwargs
     )
     # The compiled path drops the clipped axes from the bounds while the
     # python loop keeps them with length one; either way, make them
     # broadcast. ``axis`` may be an int, a tuple of ints or None here.
-    axis = kwargs.get("axis", 0)
     if axis is None:
         axes = tuple(range(np_data.ndim))
     elif isinstance(axis, tuple):
-        axes = tuple(ax % np_data.ndim for ax in axis)
+        axes = axis
     else:
         axes = (axis % np_data.ndim,)
     shape = tuple(1 if dim in axes else n for dim, n in enumerate(np_data.shape))
@@ -1910,17 +1915,11 @@ def test_sigma_clip_mask_axis_forms(axis):
         data, sigma_lower=2, sigma_upper=2, axis=axis, maxiters=2, xp=xp
     )
 
-    # The reference gets a tuple axis with its negative entries normalised:
-    # astropy's bottleneck dispatch transposes with the tuple as given and
-    # raises on a negative entry. The helper receives the tuple as written.
-    ref_axis = axis
-    if isinstance(axis, tuple):
-        ref_axis = tuple(ax % np_data.ndim for ax in axis)
     expected = _sigma_clip_reference(
         np_data,
         sigma_lower=2,
         sigma_upper=2,
-        axis=ref_axis,
+        axis=axis,
         maxiters=2,
         cenfunc="median",
         stdfunc="std",
@@ -2053,7 +2052,7 @@ def test_combine_sigma_clip_on_any_backend():
 @pytest.mark.parametrize("axis", [None, (1, 2)], ids=str)
 def test_sigma_clipping_axis_forms_any_backend(axis):
     # astropy's sigma_clip accepts axis=None and a tuple of axes;
-    # Combiner.sigma_clipping must honour them off the numpy path too,
+    # Combiner.sigma_clipping must honor them off the numpy path too,
     # with the mask coming back in the data's own shape.
     c = Combiner(_sigma_clip_ccd_list())
     c.sigma_clipping(
diff --git a/ccdproc/tests/test_nanfuncs.py b/ccdproc/tests/test_nanfuncs.py
index 7d2136aa..fc9c03ff 100644
--- a/ccdproc/tests/test_nanfuncs.py
+++ b/ccdproc/tests/test_nanfuncs.py
@@ -43,6 +43,10 @@
     (_rng.normal(size=(5, 4, 3)), (0, 2)),  # a tuple of axes
     (_rng.normal(size=(5, 4, 3)), (-1, 0)),  # a negative entry in a tuple
     (_rng.normal(size=(5, 4, 3)), (1,)),  # a single-entry tuple
+    (_some_nan, ()),  # an empty tuple reduces over no axes (elementwise)
+    (_some_nan, (np.int64(0), np.int64(1))),  # numpy integers in a tuple
+    (np.zeros((0, 3, 4)), (1, 2)),  # size-0 kept axis: empty result, no error
+    (np.ones((2, 0)), (1,)),  # reducing an axis of size 0
     (np.array([[1.0, np.nan], [2.0, np.nan], [3.0, np.nan]]), 0),  # all-NaN column
     (np.array([[1.0, np.nan], [np.nan, np.nan]]), 0),  # single non-NaN in a slice
     (np.array([np.nan, np.nan, np.nan]), 0),  # every value NaN
@@ -58,10 +62,17 @@
 @pytest.mark.filterwarnings("ignore:Mean of empty slice:RuntimeWarning")
 @pytest.mark.filterwarnings("ignore:Degrees of freedom <= 0:RuntimeWarning")
 @pytest.mark.filterwarnings("ignore:All-NaN slice encountered:RuntimeWarning")
+@pytest.mark.filterwarnings("ignore:invalid value encountered:RuntimeWarning")
 @pytest.mark.parametrize(("func", "reference"), _FUNCS)
 @pytest.mark.parametrize(("data", "axis"), _DATA)
 def test_matches_numpy(func, reference, data, axis):
     """The fallback reproduces its numpy counterpart, in shape and value."""
+    if reference is np.median and isinstance(axis, tuple) and data.size == 0:
+        # numpy.median's own tuple-axis merge reshapes with -1 and trips
+        # over the size-0 case (numpy/numpy _function_base_impl merge) --
+        # the very failure _setup avoids by spelling the merged length
+        # out. On NaN-free input nanmedian is an exact stand-in.
+        reference = np.nanmedian
     converted = xp.asarray(data, device=xp_device)
     if data is _ill_conditioned and bool(xp.all(converted == converted[0])):
         # A float32-default backend (jax without JAX_ENABLE_X64) collapses
@@ -96,10 +107,14 @@ def test_no_warning_on_all_nan_slice(func):
     data = xp.asarray(
         np.array([[1.0, np.nan], [2.0, np.nan], [3.0, np.nan]]), device=xp_device
     )
+    # A slice of size zero is just as routine to stay silent on: several
+    # numpy counterparts warn there too (np.median even divides 0 by 0).
+    empty = xp.asarray(np.ones((2, 0)), device=xp_device)
 
     with warnings.catch_warnings():
         warnings.simplefilter("error")
         func(data, axis=0)
+        func(empty, axis=1)
 
 
 def test_nansum_all_nan_slice_is_zero():
@@ -125,26 +140,47 @@ def test_nanmad_matches_astropy(axis):
 
 @pytest.mark.parametrize("func", [nansum, nanmean, nanstd, nanmedian, median])
 def test_list_axis_matches_tuple(func):
-    """A list of axes means the same as a tuple (numpy itself rejects it)."""
+    """
+    A list of axes means exactly what the equivalent tuple means.
+
+    Worth pinning because the two dispatch paths disagree upstream:
+    numpy's reductions reject a list axis outright while
+    ``astropy.stats.sigma_clip`` accepts one, and
+    ``Combiner.sigma_clipping`` forwards ``axis`` verbatim to whichever
+    path the namespace selects. If the fallbacks copied numpy's rejection,
+    a list would work for numpy data and raise for every other backend, so
+    ``_setup`` treats it like a tuple.
+    """
     data = xp.asarray(_rng.normal(size=(4, 3, 2)), device=xp_device)
     assert bool(xp.all(func(data, axis=[0, 2]) == func(data, axis=(0, 2))))
 
 
 @pytest.mark.parametrize("func", [nansum, nanmean, nanstd, nanmedian, median])
 @pytest.mark.parametrize(
-    ("axis", "error"),
+    ("axis", "error", "match"),
     [
-        (True, TypeError),  # bool subclasses int; reject it anyway
-        (1.5, TypeError),
-        (2, ValueError),
-        (-3, ValueError),
-        ((0, 0), ValueError),  # repeated axis
-        ((0, -2), ValueError),  # repeated via a negative alias
-        ((0, 2), ValueError),  # out-of-bounds entry
-        ((0, True), TypeError),  # bool entry in a tuple
+        (True, TypeError, "not bool"),  # bool subclasses int; reject it anyway
+        # np.bool_ is not a python bool: numpy 2.0 converts it to 1 with
+        # only a DeprecationWarning, so without an explicit guard these
+        # would silently reduce the wrong axes there.
+        (np.True_, TypeError, "not bool"),
+        ((0, np.True_), TypeError, "not bool"),
+        (1.5, TypeError, "must be an integer"),
+        (2, ValueError, "out of bounds"),
+        (-3, ValueError, "out of bounds"),
+        ((0, 0), ValueError, "repeated axis"),  # repeated axis
+        ((0, -2), ValueError, "repeated axis"),  # repeated via a negative alias
+        ((0, 2), ValueError, "out of bounds"),  # out-of-bounds entry
+        ((0, True), TypeError, "not bool"),  # bool entry in a tuple
     ],
 )
-def test_bad_axis(func, axis, error):
-    """Axes the fallbacks cannot handle are rejected instead of silently wrong."""
-    with pytest.raises(error):
+def test_bad_axis(func, axis, error, match):
+    """
+    Axes the fallbacks cannot handle are rejected instead of silently
+    wrong, with the message pinned: ``_mad_fallback`` and the combiner
+    lean on ``_setup`` for axis validation, so
+    ``test_ccdproc.py::test_mad_fallback_bad_axis_propagates`` checks only
+    the delegation and relies on this grid for the message contract.
+    """
+    with pytest.raises(error, match=match):
         func(xp.asarray(np.ones((2, 2)), device=xp_device), axis=axis)

From 9d18599fd2af69115d0a0ac0e373db0606fa6e92 Mon Sep 17 00:00:00 2001
From: Matt Craig <mattwcraig@gmail.com>
Date: Mon, 31 Aug 2026 08:53:11 -0500
Subject: [PATCH 3/3] Template the repeated x/axis/xp docstring entries in
 _nanfuncs

Requested in review on #1006: the six public reductions repeated the
same three Parameters entries with incidental wording drift. They are
now written once in _COMMON_PARAMS and filled into each docstring's
{params} placeholder by the _fill_doc decorator (the scipy doccer
pattern), with only the axis action phrase substituted per function.
Function-specific behaviour -- what an all-NaN slice yields, NaN
propagation -- stays inline in each Returns section, and nanmad keeps
its extra ``median`` entry after the template.

The filled docstrings render identically to the hand-written ones
(help(), Sphinx, and numpydoc all read __doc__, which is complete
after import); the trade-off is that raw source shows the {params}
placeholder. python -OO strips docstrings, so the decorator fills only
when __doc__ exists. The unified axis entry also carries the
booleans-rejected/numpy-scalars-accepted note everywhere, which
previously appeared in only three of the six.

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

diff --git a/ccdproc/_nanfuncs.py b/ccdproc/_nanfuncs.py
index 7dd4bbb4..2c1e45d1 100644
--- a/ccdproc/_nanfuncs.py
+++ b/ccdproc/_nanfuncs.py
@@ -21,6 +21,7 @@
 
 import math
 import operator
+import textwrap
 from functools import partial
 
 import array_api_compat
@@ -33,6 +34,41 @@
 
 __all__ = ["median", "nanmad", "nanmean", "nanmedian", "nanstd", "nansum"]
 
+# The ``x``/``axis``/``xp`` parameters mean the same thing for every public
+# function here, so their docstring entries are written once and filled into
+# each docstring's ``{params}`` placeholder by ``_fill_doc``; only the axis
+# action phrase differs. Function-specific behaviour (what an all-NaN slice
+# yields, NaN propagation, ...) stays inline in each Returns section.
+_COMMON_PARAMS = """\
+x : array
+    Input array. Integer and boolean inputs are promoted to the
+    namespace's default real floating dtype.
+axis : int, tuple of int, list of int or None, optional
+    Axis or axes along which {action}. Default is 0. ``None`` reduces
+    over every axis; a tuple or list over all the listed axes at once.
+    Booleans are rejected, numpy integer scalars are accepted.
+xp : array namespace, optional
+    Namespace to use. Defaults to ``array_api_compat.array_namespace(x)``.\
+"""
+
+
+def _fill_doc(**substitutions):
+    """
+    Fill a function docstring's ``{params}`` placeholder with
+    `_COMMON_PARAMS`, applying ``substitutions`` to the template first.
+    """
+
+    def decorator(func):
+        # ``python -OO`` strips docstrings; there is nothing to fill then.
+        if func.__doc__:
+            params = _COMMON_PARAMS.format(**substitutions)
+            func.__doc__ = func.__doc__.format(
+                params=textwrap.indent(params, "    ").lstrip()
+            )
+        return func
+
+    return decorator
+
 
 def _promote_to_real(x, xp, device):
     """
@@ -254,20 +290,14 @@ def _safe_divide(total, count, xp, device):
     return xp.where(count == 0, nan, quotient)
 
 
+@_fill_doc(action="to sum")
 def nansum(x, /, *, axis=0, xp=None):
     """
     Sum along an axis, ignoring NaNs, using only array-API functions.
 
     Parameters
     ----------
-    x : array
-        Input array. Integer and boolean inputs are promoted to the
-        namespace's default real floating dtype.
-    axis : int, tuple of int, list of int or None, optional
-        Axis or axes along which to sum. Default is 0. ``None`` sums over
-        every axis; a tuple or list sums over all the listed axes at once.
-    xp : array namespace, optional
-        Namespace to use. Defaults to ``array_api_compat.array_namespace(x)``.
+    {params}
 
     Returns
     -------
@@ -281,20 +311,14 @@ def nansum(x, /, *, axis=0, xp=None):
     return total
 
 
+@_fill_doc(action="to average")
 def nanmean(x, /, *, axis=0, xp=None):
     """
     Mean along an axis, ignoring NaNs, using only array-API functions.
 
     Parameters
     ----------
-    x : array
-        Input array. Integer and boolean inputs are promoted to the
-        namespace's default real floating dtype.
-    axis : int, tuple of int, list of int or None, optional
-        Axis or axes along which to average. Default is 0. ``None``
-        averages over every axis; a tuple or list over all the listed axes.
-    xp : array namespace, optional
-        Namespace to use. Defaults to ``array_api_compat.array_namespace(x)``.
+    {params}
 
     Returns
     -------
@@ -310,6 +334,7 @@ def nanmean(x, /, *, axis=0, xp=None):
     return _safe_divide(total, count, xp, device)
 
 
+@_fill_doc(action="to compute the deviation")
 def nanstd(x, /, *, axis=0, xp=None):
     """
     Standard deviation along an axis, ignoring NaNs, via array-API functions.
@@ -320,15 +345,7 @@ def nanstd(x, /, *, axis=0, xp=None):
 
     Parameters
     ----------
-    x : array
-        Input array. Integer and boolean inputs are promoted to the
-        namespace's default real floating dtype.
-    axis : int, tuple of int, list of int or None, optional
-        Axis or axes along which to compute the deviation. Default is 0.
-        ``None`` reduces over every axis; a tuple or list over all the
-        listed axes.
-    xp : array namespace, optional
-        Namespace to use. Defaults to ``array_api_compat.array_namespace(x)``.
+    {params}
 
     Returns
     -------
@@ -368,6 +385,7 @@ def nanstd(x, /, *, axis=0, xp=None):
     return xp.squeeze(xp.sqrt(variance), axis=axis)
 
 
+@_fill_doc(action="to compute the median")
 def nanmedian(x, /, *, axis=0, xp=None):
     """
     Median along an axis, ignoring NaNs, using only array-API functions.
@@ -381,16 +399,7 @@ def nanmedian(x, /, *, axis=0, xp=None):
 
     Parameters
     ----------
-    x : array
-        Input array. Integer and boolean inputs are promoted to the
-        namespace's default real floating dtype.
-    axis : int, tuple of int, list of int or None, optional
-        Axis or axes along which to compute the median. Default is 0.
-        ``None`` reduces over every axis and a tuple or list over all the
-        listed axes; booleans are rejected, numpy integer scalars are
-        accepted.
-    xp : array namespace, optional
-        Namespace to use. Defaults to ``array_api_compat.array_namespace(x)``.
+    {params}
 
     Returns
     -------
@@ -448,22 +457,14 @@ def nanmedian(x, /, *, axis=0, xp=None):
     return xp.where(xp.squeeze(n, axis=axis) == 0, nan, result)
 
 
+@_fill_doc(action="to compute the median")
 def median(x, /, *, axis=0, xp=None):
     """
     Median along an axis, using only array-API functions.
 
     Parameters
     ----------
-    x : array
-        Input array. Integer and boolean inputs are promoted to the
-        namespace's default real floating dtype.
-    axis : int, tuple of int, list of int or None, optional
-        Axis or axes along which to compute the median. Default is 0.
-        ``None`` reduces over every axis and a tuple or list over all the
-        listed axes; booleans are rejected, numpy integer scalars are
-        accepted.
-    xp : array namespace, optional
-        Namespace to use. Defaults to ``array_api_compat.array_namespace(x)``.
+    {params}
 
     Returns
     -------
@@ -488,22 +489,14 @@ def median(x, /, *, axis=0, xp=None):
     return xp.where(xp.any(xp.isnan(x), axis=axis), nan, nanmedian(x, axis=axis, xp=xp))
 
 
+@_fill_doc(action="to compute the deviation")
 def nanmad(x, /, *, axis=0, xp=None, median=None):
     """
     Median absolute deviation along ``axis``, ignoring NaNs.
 
     Parameters
     ----------
-    x : array
-        Input array. Integer and boolean inputs are promoted to the
-        namespace's default real floating dtype.
-    axis : int, tuple of int, list of int or None, optional
-        Axis or axes along which to compute the deviation. Default is 0.
-        ``None`` reduces over every axis and a tuple or list over all the
-        listed axes; booleans are rejected, numpy integer scalars are
-        accepted.
-    xp : array namespace, optional
-        Namespace to use. Defaults to ``array_api_compat.array_namespace(x)``.
+    {params}
     median : callable, optional
         Reduction used for both medians, called as ``median(x, axis=axis)``,
         always with a single integer ``axis``: a ``None`` or tuple/list