Array-API fallback for Combiner.sigma_clipping when the namespace is not numpy

CHANGES.rst

@@ -9,6 +9,13 @@ New Features
 - Allow ``ImageFileCollection.ccds`` to override the collection's FITS
   extension per call with ``ccd_kwargs["hdu"]`` while preserving ``ext=`` as
   a header filter. [#960]
+- ``Combiner.sigma_clipping`` now clips data in a non-NumPy array namespace
+  with an implementation written in terms of the array API standard that
+  reproduces ``astropy.stats.sigma_clip``'s result up to floating-point
+  rounding of the reductions: a value lying exactly on a bound can be
+  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]
 
 Other Changes and Additions
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^
@@ -134,6 +141,11 @@ Bug Fixes
   ``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]
+- Passing ``axis``, ``copy`` or ``maxiters`` to ``Combiner.sigma_clipping``
+  no longer raises ``TypeError``. [#1001]
+- Correct the ``Combiner.sigma_clipping`` docstring, which said the
+  default ``func`` was ``'median'``; the runtime default has always been
+  ``'mean'``. [#1001]
 
 2.5.1 (2025-07-05)
 ------------------

ccdproc/_nanfuncs.py

@@ -20,10 +20,11 @@
 """
 
 import operator
+from functools import partial
 
 import array_api_compat
 
-__all__ = ["median", "nanmean", "nanmedian", "nanstd", "nansum"]
+__all__ = ["median", "nanmad", "nanmean", "nanmedian", "nanstd", "nansum"]
 
 
 def _promote_to_real(x, xp, device):
@@ -411,3 +412,38 @@ def median(x, /, *, axis=0, xp=None):
     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))
+
+
+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, 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.
+    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)``.
+        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`.
+
+    Returns
+    -------
+    array
+        ``median(|x - median(x)|)`` along ``axis``, with that axis removed.
+        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)
+    if median is None:
+        median = partial(nanmedian, xp=xp)
+    center = xp.expand_dims(median(x, axis=axis), axis=axis)
+    return median(xp.abs(x - center), axis=axis)