diff --git a/CHANGES.rst b/CHANGES.rst index 205f78d2..7317bd3b 100644 --- a/CHANGES.rst +++ b/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) ------------------ diff --git a/ccdproc/_nanfuncs.py b/ccdproc/_nanfuncs.py index 3f66b2da..e1f7a201 100644 --- a/ccdproc/_nanfuncs.py +++ b/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) diff --git a/ccdproc/combiner.py b/ccdproc/combiner.py index 3422c3ca..6086b4b5 100644 --- a/ccdproc/combiner.py +++ b/ccdproc/combiner.py @@ -2,6 +2,9 @@ """This module implements the combiner class.""" +import itertools +import operator +import warnings from copy import deepcopy from functools import partial @@ -19,7 +22,7 @@ from astropy.stats import sigma_clip from astropy.utils import deprecated_renamed_argument -from ._nanfuncs import nanmean, nanmedian, nanstd, nansum +from ._nanfuncs import _setup, nanmad, nanmean, nanmedian, nanstd, nansum from .core import _namespace_dtype, _native_numpy, _to_numpy, sigma_func __all__ = ["Combiner", "combine"] @@ -85,6 +88,232 @@ def _default_std(xp=None): return partial(nanstd, xp=xp) +def _nanmadstd(x, /, *, axis=0, xp): + """ + NaN-aware median absolute deviation along ``axis``, scaled to a + standard deviation. + + This is the ``'mad_std'`` deviation option of `Combiner.sigma_clipping` + when the array namespace is not numpy. It uses the same tiered median + as `Combiner.median_combine` (`bottleneck` for numpy, the namespace's + ``nanmedian`` if it has one, otherwise the fallback in + `ccdproc._nanfuncs`), so NaNs are ignored on every backend. + + Parameters + ---------- + x : array + Input array, in the namespace ``xp``. + axis : int, optional + Axis along which to compute the deviation. Default is ``0``. + xp : array namespace + Namespace of ``x``. + + Returns + ------- + array + ``1.4826 * median(|x - median(x)|)`` along ``axis``, with that axis + removed. + """ + return 1.482602218505602 * nanmad(x, axis=axis, xp=xp, median=_default_median(xp)) + + +# The string options of Combiner.sigma_clipping, mapped to the functions +# that produce the tiered (bottleneck -> native -> fallback) implementation +# for a namespace. The names match astropy.stats.sigma_clip's cenfunc and +# stdfunc options so that the two code paths accept the same strings. +_SIGMA_CLIP_CENFUNCS = {"median": _default_median, "mean": _default_average} +_SIGMA_CLIP_STDFUNCS = { + "std": _default_std, + "mad_std": lambda xp: partial(_nanmadstd, xp=xp), +} + + +def _resolve(func, options, kind, xp): + """ + Turn a ``cenfunc``/``stdfunc`` argument into a callable for ``xp``. + + Parameters + ---------- + func : str or callable + A key of ``options`` or a callable ``f(data, axis=axis)``. + options : dict + Map from option name to a function of ``xp`` returning the callable. + kind : str + Name of the argument, used in error messages. + xp : array namespace + Namespace the callable will operate in. + + Returns + ------- + callable + The resolved function. + + Raises + ------ + ValueError + If ``func`` is a string that is not one of ``options`` (this also + catches a non-string, non-callable ``func``, a corner astropy + itself does not guard against either). + """ + if callable(func): + return func + try: + return options[func](xp) + except KeyError: + raise ValueError( + f"{kind} must be one of {sorted(options)} or a callable, got {func!r}" + ) from None + + +def _sigma_clip_mask( + data, + *, + sigma_lower=3, + sigma_upper=3, + axis=0, + maxiters=1, + cenfunc="median", + stdfunc="std", + xp=None, +): + """ + Iterative sigma clipping written purely in terms of the array API. + + Reproduces the mask of ``astropy.stats.sigma_clip``'s compiled path + (the one its string ``cenfunc`` and ``stdfunc`` options take) for any + array-API namespace, without converting the data to numpy, up to + floating-point rounding of the reductions: a value lying exactly on a + bound can be classified differently from astropy. + `Combiner.sigma_clipping` uses this when the namespace is not numpy; + numpy data go to astropy directly. + + Parameters + ---------- + data : array + Data to clip. + sigma_lower, sigma_upper : float or None, optional + 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``. + 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'``. + stdfunc : {'std', 'mad_std'} or callable, optional + Statistic for the deviation, same requirements as ``cenfunc``. + Default is ``'std'``. + xp : array namespace, optional + Namespace of ``data``; resolved from ``data`` when ``None``. + + Returns + ------- + array of bool + Mask of the clipped values, ``True`` where a value is rejected, in + the namespace and on the device of ``data``. + + Raises + ------ + NotImplementedError + If ``axis`` is not a single integer. + ValueError + If ``axis`` is out of bounds, ``maxiters`` is not positive, or a + string ``cenfunc``/``stdfunc`` is not one of the options. + + Notes + ----- + The returned mask is the bounds of the *last* iteration applied to the + original data, ``~isfinite(data) | (data < lower) | (data > upper)``, + rather than the union of the values rejected in each iteration: an + earlier iteration's rejection can be undone if the bounds widen. + Non-finite values are always masked. A slice whose values are all + rejected before the last iteration gets NaN bounds and so keeps only + its non-finite entries masked. This is what astropy's compiled path + does on every astropy version. Astropy's python loop, taken when + ``cenfunc`` or ``stdfunc`` is a callable, did the same with + ``copy=False`` up to astropy 8.0 but from 8.1 masks the union of every + iteration's rejections, so for a callable the numpy path of + `Combiner.sigma_clipping` can differ from this one in those corners. + + The ``'median'``, ``'mean'``, ``'std'`` and ``'mad_std'`` options use + the same NaN-aware reductions as the ``Combiner`` combination methods. + Warnings that those reductions may raise on all-NaN slices are + suppressed here, as astropy suppresses them. Non-floating input, and + floating input narrower than the namespace default, is promoted to the + namespace's default real floating dtype before computing, but only + when both ``cenfunc`` and ``stdfunc`` are strings: that is when + astropy's compiled path (float64 internally) is the reproduction + target. A callable ``cenfunc`` or ``stdfunc`` instead matches + astropy's python loop, which keeps the data's own dtype, so no + promotion happens then either. + + An integer ``maxiters`` runs exactly that many iterations, unlike + astropy, which returns as soon as an iteration rejects nothing, and + never synchronises with the host, so on a lazy backend such as dask the + whole clip stays one graph. ``maxiters=None`` has to test after each + 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) + + if isinstance(cenfunc, str) and isinstance(stdfunc, str): + # B1: widen (never narrow) to the namespace's default real floating + # dtype. _setup above already promotes non-floating input to that + # dtype; a floating input keeps its own dtype there, but astropy's + # *compiled* path -- taken only when both cenfunc and stdfunc are + # strings -- always computes in float64, so float32 input must be + # widened too, or its mask can differ from astropy's near a bound. + # A callable cenfunc or stdfunc instead takes astropy's python + # loop, which applies the callable to the data's own dtype with no + # promotion (verified: cenfunc=np.nanmean/stdfunc=np.nanstd on + # float32 input returns float32 bounds), so this branch leaves a + # callable's input dtype alone to match. xp.result_type picks the + # wider of the two, so this is a no-op once data is already at or + # above the default width (jax without JAX_ENABLE_X64 has no + # float64, so its own default is the ceiling). + info = xp.__array_namespace_info__() + default = info.default_dtypes(device=device)["real floating"] + data = xp.astype(data, xp.result_type(data.dtype, default)) + + if maxiters is not None: + maxiters = operator.index(maxiters) + if maxiters < 1: + raise ValueError("maxiters must be None or a positive integer.") + + # astropy: ``sigma_lower or sigma`` with ``sigma=3``, so None and 0 both + # mean 3. Python floats, because a strict namespace rejects numpy scalars + # in arithmetic with its arrays. + sigma_lower = float(sigma_lower or 3) + sigma_upper = float(sigma_upper or 3) + + center_func = _resolve(cenfunc, _SIGMA_CLIP_CENFUNCS, "cenfunc", xp) + std_func = _resolve(stdfunc, _SIGMA_CLIP_STDFUNCS, "stdfunc", xp) + + nan = xp.asarray(xp.nan, dtype=data.dtype, device=device) + invalid = ~xp.isfinite(data) + filtered = xp.where(invalid, nan, data) + + for _ in range(maxiters) if maxiters else itertools.count(): + with warnings.catch_warnings(): + # All-NaN slices make numpy's nan-functions warn; astropy + # silences the same warnings in its _compute_bounds. + warnings.simplefilter("ignore", RuntimeWarning) + center = xp.expand_dims(center_func(filtered, axis=axis), axis=axis) + deviation = xp.expand_dims(std_func(filtered, axis=axis), axis=axis) + lower = center - deviation * sigma_lower + upper = center + deviation * sigma_upper + rejected = (filtered < lower) | (filtered > upper) + if maxiters is None and not bool(xp.any(rejected)): + break + filtered = xp.where(rejected, nan, filtered) + + return invalid | (data < lower) | (data > upper) + + class Combiner: """ A class for combining CCDData objects. @@ -442,13 +671,14 @@ def sigma_clipping( ---------- low_thresh : positive float or None, optional Threshold for rejecting pixels that deviate below the baseline - value. If negative value, then will be convert to a positive - value. If None, no rejection will be done based on low_thresh. + value. If None, it is treated as 3, as + :func:`~astropy.stats.sigma_clip` does. Default is 3. high_thresh : positive float or None, optional Threshold for rejecting pixels that deviate above the baseline - value. If None, no rejection will be done based on high_thresh. + value. If None, it is treated as 3, as + :func:`~astropy.stats.sigma_clip` does. Default is 3. func : {'median', 'mean'} or callable, optional @@ -457,7 +687,7 @@ def sigma_clipping( function/object and the ``axis`` keyword is used, then it must be able to ignore NaNs (e.g., `numpy.nanmean`) and it must have an ``axis`` keyword to return an array with axis dimension(s) - removed. The default is ``'median'``. + removed. The default is ``'mean'``. dev_func : {'std', 'mad_std'} or callable, optional The statistic or callable function/object used to compute the @@ -468,28 +698,94 @@ def sigma_clipping( removed. The default is ``'std'``. kwd - Any remaining keyword arguments are passed to astropy's - :func:`~astropy.stats.sigma_clip` function. + ``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 + 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 + astropy's :func:`~astropy.stats.sigma_clip` when the data are + numpy arrays; for any other array namespace ``copy`` is + ignored and those other keywords raise `TypeError`. + + Notes + ----- + When the data are numpy arrays the clipping is done by + :func:`~astropy.stats.sigma_clip`. For any other array namespace it + is done by an implementation written in terms of the array API + standard that reproduces the result of astropy's compiled path + (the one the string options take) up to floating-point rounding of + the reductions: a value lying exactly on a bound can be classified + differently from astropy. The mask is the bounds of the last + iteration applied to the data, and non-finite values are always + masked. The string options for ``func`` and ``dev_func`` use the + same NaN-aware reductions as the combination methods on every + backend. + + Pixels that are already masked are not excluded from the statistics + used for the clipping, but they stay masked. """ # Remove in 3.0 _ = kwd.pop("use_astropy", True) - self._data_arr_mask = ( - self._data_arr_mask - | sigma_clip( + xp = self._xp + # Pop rather than get: these were also forwarded through **kwd, + # which made passing any of them a "multiple values" TypeError. + axis = kwd.pop("axis", 0) + copy = kwd.pop("copy", False) + maxiters = kwd.pop("maxiters", 1) + + # This method always requests the mask itself (``masked=True`` + # below, or the mask that _sigma_clip_mask always returns), so + # ``masked`` collides with that and ``return_bounds`` would return a + # tuple this method cannot consume; reject both on every namespace + # rather than let them reach astropy as an accidental "multiple + # values" TypeError or an AttributeError on the tuple. + disallowed = sorted({"masked", "return_bounds"} & kwd.keys()) + if disallowed: + raise TypeError( + f"sigma_clipping does not accept {disallowed}: it always " + "requests the mask itself from astropy.stats.sigma_clip, " + "on every array namespace." + ) + + if array_api_compat.is_numpy_namespace(xp): + clipped = sigma_clip( self._data_arr, sigma_lower=low_thresh, sigma_upper=high_thresh, - axis=kwd.get("axis", 0), - copy=kwd.get("copy", False), - maxiters=kwd.get("maxiters", 1), + axis=axis, + copy=copy, + maxiters=maxiters, cenfunc=func, stdfunc=dev_func, masked=True, **kwd, ).mask - ) + else: + if kwd: + raise TypeError( + "sigma_clipping got unexpected keyword argument(s) " + f"{sorted(kwd)}. These are options of " + "astropy.stats.sigma_clip (such as grow) and are only " + "available when the array namespace is numpy." + ) + clipped = _sigma_clip_mask( + self._data_arr, + sigma_lower=low_thresh, + sigma_upper=high_thresh, + axis=axis, + maxiters=maxiters, + cenfunc=func, + stdfunc=dev_func, + xp=xp, + ) + + self._data_arr_mask = self._data_arr_mask | clipped def _get_scaled_data(self, scale_arg): if scale_arg is not None: diff --git a/ccdproc/core.py b/ccdproc/core.py index 35f72c21..846d86ad 100644 --- a/ccdproc/core.py +++ b/ccdproc/core.py @@ -32,6 +32,7 @@ ) from ._nanfuncs import _promote_to_real from ._nanfuncs import median as _nanfuncs_median +from ._nanfuncs import nanmad as _nanfuncs_nanmad from ._nanfuncs import nanmedian as _nanfuncs_nanmedian from .log_meta import log_to_metadata from .utils.slices import slice_from_string @@ -364,15 +365,16 @@ def _mad_fallback(data, axis, ignore_nan, xp=None): ) axis = -1 - def med(d, ax): + def med(d, axis): try: - return (xp.nanmedian if ignore_nan else xp.median)(d, axis=ax) + return (xp.nanmedian if ignore_nan else xp.median)(d, axis=axis) except AttributeError: fallback = _nanfuncs_nanmedian if ignore_nan else _nanfuncs_median - return fallback(d, axis=ax, xp=xp) + return fallback(d, axis=axis, xp=xp) - center = med(data, axis) - return med(xp.abs(data - xp.expand_dims(center, axis=axis)), axis) + # The deviation itself is _nanfuncs.nanmad, the same computation + # Combiner._nanmadstd uses; only the median tier passed in differs. + return _nanfuncs_nanmad(data, axis=axis, xp=xp, median=med) @log_to_metadata diff --git a/ccdproc/tests/array_escape_baseline.txt b/ccdproc/tests/array_escape_baseline.txt index 634a21d9..ad10f648 100644 --- a/ccdproc/tests/array_escape_baseline.txt +++ b/ccdproc/tests/array_escape_baseline.txt @@ -14,7 +14,6 @@ # that will never leave. Verify/adjust the seeded tags by hand. # combiner.py combine numpy.asarray BOUNDARY: astropy CCDData mask/uncertainty attributes are numpy-backed -combiner.py sigma_clipping numpy.asanyarray BOUNDARY: astropy.stats.sigma_clip is numpy-only core.py _cosmicray_median_array numpy.asarray BOUNDARY: scipy.ndimage median_filter/maximum_filter are numpy-only core.py _to_numpy numpy.asarray BOUNDARY: deliberate host copy for numpy-only consumers (combine() writes output_file through astropy.io.fits) core.py background_deviation_filter numpy.asarray BOUNDARY: scipy.ndimage.generic_filter is numpy-only diff --git a/ccdproc/tests/test_combiner.py b/ccdproc/tests/test_combiner.py index 002b20ce..bc832f87 100644 --- a/ccdproc/tests/test_combiner.py +++ b/ccdproc/tests/test_combiner.py @@ -1,6 +1,7 @@ # Licensed under a 3-clause BSD style license - see LICENSE.rst import math import types +import warnings from functools import partial import array_api_compat @@ -10,9 +11,12 @@ import pytest from astropy.nddata import CCDData from astropy.stats import median_absolute_deviation as mad +from astropy.stats import sigma_clip from astropy.utils.data import get_pkg_data_filename +from astropy.utils.exceptions import AstropyDeprecationWarning from numpy.testing import assert_allclose +import ccdproc.combiner as combiner_module from ccdproc import create_deviation from ccdproc._nanfuncs import nanmean, nanmedian, nanstd, nansum from ccdproc.combiner import ( @@ -23,6 +27,7 @@ _default_median, _default_std, _default_sum, + _sigma_clip_mask, combine, sigma_func, ) @@ -538,6 +543,10 @@ def test_combiner_minmax_min(): def test_combiner_sigmaclip_high(): + # Five frames at 0 / -10 / +10 and a sixth at +1000. With a median + # center and mad_std deviation the sixth frame is above the 3-sigma + # upper bound at every pixel (median 5, bound ~49), so it is masked + # everywhere and the other five nowhere; low_thresh=None means 3 too. ccd_list = [ CCDData(xp.zeros((10, 10)), unit=u.adu), CCDData(xp.zeros((10, 10)) - 10, unit=u.adu), @@ -548,12 +557,17 @@ def test_combiner_sigmaclip_high(): ] c = Combiner(ccd_list) - # using mad for more robust statistics vs. std - c.sigma_clipping(high_thresh=3, low_thresh=None, func="median", dev_func=mad) - assert c._data_arr_mask[5].all() + c.sigma_clipping(high_thresh=3, low_thresh=None, func="median", dev_func="mad_std") + assert xp.all(c._data_arr_mask[5, ...]) + assert not xp.any(c._data_arr_mask[:5, ...]) def test_combiner_sigmaclip_single_pix(): + # Six frames at 0 / -10 / +10, then pixel (5, 5) of the fifth frame is + # set to 25 while the other frames stay within +-10 there. With a median + # center and mad_std deviation the bound at that pixel is ~19.7, so only + # that one value is rejected: not the rest of its frame and not the + # other frames' values at (5, 5). ccd_list = [ CCDData(xp.zeros((10, 10)), unit=u.adu), CCDData(xp.zeros((10, 10)) - 10, unit=u.adu), @@ -570,11 +584,19 @@ def test_combiner_sigmaclip_single_pix(): combo._data_arr = xpx.at(combo._data_arr)[2, 5, 5].set(5) combo._data_arr = xpx.at(combo._data_arr)[3, 5, 5].set(-5) combo._data_arr = xpx.at(combo._data_arr)[4, 5, 5].set(25) - combo.sigma_clipping(high_thresh=3, low_thresh=None, func="median", dev_func=mad) - assert combo._data_arr_mask[4, 5, 5] + combo.sigma_clipping( + high_thresh=3, low_thresh=None, func="median", dev_func="mad_std" + ) + assert bool(combo._data_arr_mask[4, 5, 5]) + assert not xp.any(combo._data_arr_mask[:4, ...]) + assert not xp.any(combo._data_arr_mask[5, ...]) def test_combiner_sigmaclip_low(): + # Five frames at 0 / -10 / +10 and a sixth at -1000. With a median + # center and mad_std deviation the sixth frame is below the 3-sigma + # lower bound at every pixel (median -5, bound ~-49), so it is masked + # everywhere and the other five nowhere; high_thresh=None means 3 too. ccd_list = [ CCDData(xp.zeros((10, 10)), unit=u.adu), CCDData(xp.zeros((10, 10)) - 10, unit=u.adu), @@ -585,9 +607,9 @@ def test_combiner_sigmaclip_low(): ] c = Combiner(ccd_list) - # using mad for more robust statistics vs. std - c.sigma_clipping(high_thresh=None, low_thresh=3, func="median", dev_func=mad) - assert c._data_arr_mask[5].all() + c.sigma_clipping(high_thresh=None, low_thresh=3, func="median", dev_func="mad_std") + assert xp.all(c._data_arr_mask[5, ...]) + assert not xp.any(c._data_arr_mask[:5, ...]) # test that the median combination works and returns a ccddata object @@ -1628,3 +1650,361 @@ def test_combine_array_package_dask_module(tmp_path): result = combine(files, array_package=dask, unit="adu") assert array_api_compat.is_dask_array(result.data) + + +# Sigma clipping off numpy: Combiner.sigma_clipping hands numpy data to +# astropy.stats.sigma_clip and everything else to the array-API +# implementation _sigma_clip_mask. The tests below run the array-API +# implementation on every backend, numpy included, so that it is covered by +# the numpy coverage job, and check it against astropy's mask. + + +def _sigma_clip_datasets(): + """ + Numpy data sets that exercise the corners of astropy's mask; see the + comment on each array below for the corner it targets. + """ + rng = np.random.default_rng(929) + # The plain case: planted outliers among otherwise normal values. + normal = rng.normal(size=(8, 5, 4)) + normal[0, 0, 0] = 10.0 + normal[3, 1, 2] = -25.0 + normal[7, 4, 3] = 6.0 + normal[2, 2, 2] = -7.0 + + # NaNs and infs that must always be masked, a NaN next to finite values + # that must not drag the whole column with it, and an all-NaN column. + nan_inf = normal.copy() + nan_inf[1, 0, 1] = np.nan + nan_inf[5, 3, 3] = np.nan + nan_inf[2, 4, 0] = np.inf + nan_inf[6, 1, 1] = -np.inf + nan_inf[:, 1, 2] = np.nan # an all-NaN column + + # One value barely above a constant slice: with std it is the only + # thing inside the bounds; with mad_std the deviation is 0 and the + # bounds collapse onto the center, so whether the 5.0s are rejected + # comes down to strict "<"/">" comparisons, as astropy uses. + zero_std = np.full((6, 3, 3), 5.0) + zero_std[2, 1, 1] = 5.0 + 1e-12 + + # One column that mean/mad_std clips entirely in the first iteration + # (its MAD is tiny while its mean is far off), and one whose second + # iteration has zero spread, next to two ordinary columns. Non-round + # values, so that no value sits exactly on a bound (see the rounding + # caveat in the _sigma_clip_mask Notes for why that would matter). + collapse = rng.normal(size=(6, 2, 2)) + collapse[:, 0, 0] = [50.0, 0.1, -0.2, 0.05, 1.3, 2.7] + collapse[:, 1, 1] = [0.0, 0.0, 0.0, 0.0, 0.0, 7.0] + + # Integer input; the helper must promote it to the namespace default + # float before subtracting a float center (a strict namespace rejects + # int - float), and the mask must still match astropy's. + ints = np.array([[1, 50], [2, 51], [3, 52], [100, 53]]) + + # Same corners as normal, in float32: pins the promotion in + # _sigma_clip_mask that keeps a float32 clip matching astropy's, which + # always computes in float64. + float32 = normal.astype(np.float32) + + return { + "normal": normal, + "nan_inf": nan_inf, + "zero_std": zero_std, + "collapse": collapse, + "int": ints, + "float32": float32, + } + + +def _sigma_clip_reference(np_data, **kwargs): + """ + Mask of ``astropy.stats.sigma_clip``'s final bounds applied to ``np_data``. + + Parameters + ---------- + np_data : numpy.ndarray + Data to clip, in numpy (the values the backend under test sees). + **kwargs + Passed to `astropy.stats.sigma_clip`: ``sigma_lower``, + ``sigma_upper``, ``maxiters``, ``cenfunc``, ``stdfunc``, ``axis``. + + Returns + ------- + numpy.ndarray of bool + True where ``np_data`` is non-finite or outside the bounds of the + last iteration. For string ``cenfunc``/``stdfunc`` this is also + asserted to equal astropy's own mask, which is what numpy data get + from ``Combiner.sigma_clipping``. + + Notes + ----- + The bounds are what astropy's two code paths agree on. Its compiled + path (string ``cenfunc`` and ``stdfunc``) masks the data outside the + bounds of the last iteration. Its python loop (any callable) did the + same with ``copy=False``, which is what ``Combiner.sigma_clipping`` + passes, up to astropy 8.0; from 8.1 (astropy#19858) it masks the union + of every iteration's rejections whatever ``copy`` is. The two differ + only when the bounds widen between iterations or a slice is clipped + entirely, which the mad_std cases of the data sets provoke. + ``_sigma_clip_mask`` follows the compiled path, so the reference is + built from the bounds. For the compiled path astropy's own mask is + checked against it too: that is what numpy data get from + ``Combiner.sigma_clipping``, and the other backends must agree with it. + """ + _, 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)) + lower = np.reshape(lower, shape) + upper = np.reshape(upper, shape) + with np.errstate(invalid="ignore"): + expected = ~np.isfinite(np_data) | (np_data < lower) | (np_data > upper) + + if isinstance(kwargs.get("cenfunc", "median"), str) and isinstance( + kwargs.get("stdfunc", "std"), str + ): + from_astropy = np.ma.getmaskarray( + sigma_clip(np_data.copy(), masked=True, copy=False, **kwargs) + ) + assert np.array_equal(from_astropy, expected) + return expected + + +# This is a differential test: the only specification of _sigma_clip_mask is +# "astropy.stats.sigma_clip's bounds applied to the data" (see +# _sigma_clip_reference for why the bounds rather than astropy's mask), and +# astropy's behaviour has corners that no hand-written expectation would pin +# down (see the Notes of _sigma_clip_mask: final-bounds-not-union, None/0 +# thresholds meaning 3, fully clipped slices, always-masked non-finite +# values, and the floating-point rounding caveat). The (cenfunc, stdfunc) +# pairs cover the string dispatch table plus a callable pair; (None, 0) +# exercises the None/0 -> 3 threshold rule. Wrapping the _sigma_clip_mask +# call in simplefilter("error") folds in silence on the all-NaN column of +# the nan_inf/float32 datasets (skipped on dask: a lazy backend only runs +# numpy's nan-functions on its chunks when the graph is computed, outside +# this test's control). No value in the datasets sits exactly on a bound +# (see the rounding caveat above) except where a corner intentionally +# provokes it. 6 datasets x 4 (cenfunc, stdfunc) pairs x 4 threshold pairs +# x 3 maxiters = 288 cases per backend. +@pytest.mark.filterwarnings("ignore::astropy.utils.exceptions.AstropyUserWarning") +@pytest.mark.filterwarnings("ignore:invalid value encountered:RuntimeWarning") +@pytest.mark.filterwarnings("ignore:Mean of empty slice:RuntimeWarning") +@pytest.mark.filterwarnings("ignore:Degrees of freedom <= 0:RuntimeWarning") +@pytest.mark.parametrize("maxiters", [1, 3, None]) +@pytest.mark.parametrize( + ("sigma_lower", "sigma_upper"), [(None, 0), (2, 2.5), (1.5, 1.5), (3, 1)] +) +@pytest.mark.parametrize( + ("cenfunc", "stdfunc"), + [ + ("median", "std"), + ("mean", "std"), + ("mean", "mad_std"), + ("callable", "callable"), + ], +) +@pytest.mark.parametrize("dataset", sorted(_sigma_clip_datasets())) +def test_sigma_clip_mask_matches_astropy( + dataset, cenfunc, stdfunc, sigma_lower, sigma_upper, maxiters +): + data = xp.asarray(_sigma_clip_datasets()[dataset], device=xp_device) + # Build the reference from the backend's own values so that a backend + # with a different default float width still sees identical data. + np_data = _to_numpy(data) + + if cenfunc == "callable": + xp_cenfunc, ref_cenfunc = partial(nanmean, xp=xp), np.nanmean + else: + xp_cenfunc = ref_cenfunc = cenfunc + if stdfunc == "callable": + xp_stdfunc, ref_stdfunc = partial(nanstd, xp=xp), np.nanstd + else: + xp_stdfunc = ref_stdfunc = stdfunc + + with warnings.catch_warnings(): + if not array_api_compat.is_dask_namespace(xp): + warnings.simplefilter("error") + result = _sigma_clip_mask( + data, + sigma_lower=sigma_lower, + sigma_upper=sigma_upper, + axis=0, + maxiters=maxiters, + cenfunc=xp_cenfunc, + stdfunc=xp_stdfunc, + xp=xp, + ) + expected = _sigma_clip_reference( + np_data, + sigma_lower=sigma_lower, + sigma_upper=sigma_upper, + axis=0, + maxiters=maxiters, + cenfunc=ref_cenfunc, + stdfunc=ref_stdfunc, + ) + + assert array_api_compat.array_namespace(result) is array_api_compat.array_namespace( + data + ) + assert array_api_compat.device(result) == array_api_compat.device(data) + assert result.dtype == xp.bool + assert result.shape == data.shape + assert bool(xp.all(result == xp.asarray(expected, device=xp_device))) + + +def test_sigma_clip_mask_argument_handling(): + """ + Argument handling of _sigma_clip_mask that the differential grid does + not reach: numpy-scalar thresholds are converted before they meet the + data (a strict namespace rejects them in arithmetic), ``xp=None`` + resolves the namespace from the data, a negative axis counts from the + end, and ``maxiters=0`` and unknown ``cenfunc``/``stdfunc`` strings + raise. Axis validation itself is `_nanfuncs._setup`'s and is covered + by ``test_nanfuncs.py::test_bad_axis``. + """ + data = xp.asarray(_sigma_clip_datasets()["normal"], device=xp_device) + + expected = _sigma_clip_mask(data, sigma_lower=2, sigma_upper=2.5, xp=xp) + result = _sigma_clip_mask( + data, sigma_lower=np.int64(2), sigma_upper=np.float64(2.5), xp=xp + ) + assert bool(xp.all(result == expected)) + + # No xp= here: the helper resolves the namespace from the data. + assert bool(xp.all(_sigma_clip_mask(data) == _sigma_clip_mask(data, xp=xp))) + + # A negative axis is normalised like numpy's. + expected = _sigma_clip_mask(data, axis=2, xp=xp) + assert bool(xp.all(_sigma_clip_mask(data, axis=-1, xp=xp) == expected)) + + with pytest.raises(ValueError, match="maxiters"): + _sigma_clip_mask(data, maxiters=0, xp=xp) + + with pytest.raises(ValueError, match="cenfunc must be one of"): + _sigma_clip_mask(data, cenfunc="mode", xp=xp) + with pytest.raises(ValueError, match="stdfunc must be one of"): + _sigma_clip_mask(data, stdfunc="var", xp=xp) + + +def _sigma_clip_ccd_list(): + """ + The ``normal`` set of `_sigma_clip_datasets` as a list of `CCDData`, + one per slice along axis 0, as arrays of the backend under test. + """ + data = _sigma_clip_datasets()["normal"] + return [CCDData(xp.asarray(image, device=xp_device), unit=u.adu) for image in data] + + +@pytest.mark.parametrize("force_fallback", [False, True]) +def test_sigma_clipping_dispatch(monkeypatch, force_fallback): + """ + ``Combiner.sigma_clipping`` sends numpy data to astropy and everything + else to ``_sigma_clip_mask``, with ``axis``/``maxiters``/``copy`` + reaching astropy correctly, an astropy-only kwarg such as ``grow`` + raising off the numpy path, the ``use_astropy`` deprecation warning + raised on both paths, and a pre-existing mask preserved either way. + ``force_fallback`` forces the array-API branch on every backend, numpy + included, so the numpy coverage job sees it too. + """ + if force_fallback: + # This also flips _default_median/_default_std (combiner.py:31, + # :76), which read the same predicate, so on numpy this clips with + # array_api_compat.numpy.nanmedian/nanstd instead of bottleneck; + # only the dispatch is under test here, not those reductions. + monkeypatch.setattr(array_api_compat, "is_numpy_namespace", lambda _xp: False) + + calls, real = [], combiner_module.sigma_clip + monkeypatch.setattr( + combiner_module, + "sigma_clip", + lambda *a, **k: calls.append(k) or real(*a, **k), + ) + + c = Combiner(_sigma_clip_ccd_list()) + # A pixel that sigma clipping would not reject on its own. + c._data_arr_mask = xpx.at(c._data_arr_mask)[4, 2, 3].set(True) + with pytest.warns(AstropyDeprecationWarning, match="use_astropy"): + c.sigma_clipping( + low_thresh=2, + high_thresh=2.5, + func="median", + dev_func="mad_std", + axis=0, + maxiters=2, + copy=False, + use_astropy=True, + ) + + # Reads the (possibly patched) predicate sigma_clipping itself just + # used, so this reflects which path the call above actually took. + uses_astropy = array_api_compat.is_numpy_namespace(xp) + assert (len(calls) == 1) is uses_astropy + if uses_astropy: + # numpy data must keep going to astropy, with the same arguments + # as before, so that its compiled fast path is used. + assert calls[0]["masked"] is True + assert calls[0]["maxiters"] == 2 + assert calls[0]["copy"] is False + assert calls[0]["axis"] == 0 + + expected = _sigma_clip_reference( + _to_numpy(c._data_arr), + sigma_lower=2, + sigma_upper=2.5, + axis=0, + maxiters=2, + cenfunc="median", + stdfunc="mad_std", + ) + expected[4, 2, 3] = True # the pre-existing mask stays set + assert array_api_compat.array_namespace( + c._data_arr_mask + ) is array_api_compat.array_namespace(c._data_arr) + assert array_api_compat.device(c._data_arr_mask) == array_api_compat.device( + c._data_arr + ) + assert bool(xp.all(c._data_arr_mask == xp.asarray(expected, device=xp_device))) + + # masked and return_bounds are rejected on every namespace: the method + # always requests the mask itself from astropy. + with pytest.raises(TypeError, match="masked"): + c.sigma_clipping(masked=True, return_bounds=True) + + if uses_astropy: + # Forwarded to astropy as before. + c.sigma_clipping(grow=1) + else: + with pytest.raises(TypeError, match="grow"): + c.sigma_clipping(grow=1) + + +def test_combine_sigma_clip_on_any_backend(): + # combine() forwards the namespace's mean and std as callables to + # sigma_clipping; check that the clipping happens on every backend. + ccd_list = _sigma_clip_ccd_list() + result = combine( + ccd_list, + method="average", + sigma_clip=True, + sigma_clip_low_thresh=2, + sigma_clip_high_thresh=2, + ) + np_data = np.stack([_to_numpy(ccd.data) for ccd in ccd_list]) + mask = _sigma_clip_reference( + np_data, + sigma_lower=2, + sigma_upper=2, + axis=0, + maxiters=1, + cenfunc=np.mean, + stdfunc=np.std, + ) + assert mask.any() + 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)) diff --git a/ccdproc/tests/test_nanfuncs.py b/ccdproc/tests/test_nanfuncs.py index cd5f58c3..5536c7b7 100644 --- a/ccdproc/tests/test_nanfuncs.py +++ b/ccdproc/tests/test_nanfuncs.py @@ -4,8 +4,9 @@ import array_api_extra as xpx import numpy as np import pytest +from astropy.stats import median_absolute_deviation -from ccdproc._nanfuncs import median, nanmean, nanmedian, nanstd, nansum +from ccdproc._nanfuncs import median, nanmad, nanmean, nanmedian, nanstd, nansum from ccdproc.conftest import testing_array_device as xp_device from ccdproc.conftest import testing_array_library as xp @@ -105,6 +106,18 @@ 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]) +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) + expected = xp.asarray( + median_absolute_deviation(_some_nan, axis=axis, ignore_nan=True), + device=xp_device, + ) + assert result.shape == expected.shape + assert bool(xp.all(xpx.isclose(result, expected, equal_nan=True))) + + @pytest.mark.parametrize("func", [nansum, nanmean, nanstd, nanmedian, median]) @pytest.mark.parametrize( ("axis", "error"), diff --git a/docs/array_api.rst b/docs/array_api.rst index 43160157..ba435c41 100644 --- a/docs/array_api.rst +++ b/docs/array_api.rst @@ -173,6 +173,18 @@ What limitations should I be aware of? ``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. ++ ``Combiner.sigma_clipping`` uses ``astropy.stats.sigma_clip`` only for + `numpy`_ arrays. For any other array library it uses an implementation + written in terms of the array API standard that reproduces astropy's + result up to floating-point rounding of the reductions: a value lying + exactly on a bound can be classified differently from astropy. + Astropy-only options such as ``grow`` are not available there and raise + ``TypeError``; ``masked`` and ``return_bounds`` are not accepted by + ``sigma_clipping`` at all, on any array library, because the wrapper + always asks astropy for the mask itself. On a lazy library such as + `dask`_, prefer an integer ``maxiters``: ``maxiters=None`` has to + compute the data after every iteration to find out whether anything + else was rejected. Which array library should I use? --------------------------------- diff --git a/docs/image_combination.rst b/docs/image_combination.rst index 67e79ebb..434ad4ea 100644 --- a/docs/image_combination.rst +++ b/docs/image_combination.rst @@ -68,20 +68,25 @@ the list of images. The `~ccdproc.combiner.Combiner.sigma_clipping` method is very flexible: you can specify both the function for calculating the central value and the function -for calculating the deviation. The default is to use the mean (ignoring any -masked pixels) for the central value and the standard deviation (again -ignoring any masked values) for the deviation. +for calculating the deviation. The default is to use the mean for the +central value and the standard deviation for the deviation; pixels that +are already masked are not excluded from these statistics but stay masked. You can mask pixels more than 5 standard deviations above or 2 standard deviations below the median with - >>> combiner.sigma_clipping(low_thresh=2, high_thresh=5, func=np.ma.median) + >>> combiner.sigma_clipping(low_thresh=2, high_thresh=5, func="median") .. note:: - Numpy masked median can be very slow in exactly the situation typically - encountered in reducing ccd data: a cube of data in which one dimension - (in the case the number of frames in the combiner) is much smaller than - the number of pixels. + Prefer the string options ``func="median"`` and ``dev_func="mad_std"`` + to passing a function: they work with every array library, and for + NumPy data they use the compiled fast path of + :func:`~astropy.stats.sigma_clip`. A NaN-aware function such as + ``np.nanmedian`` also works for NumPy data; ``np.ma.median`` can be + very slow in exactly the situation typically encountered in reducing + ccd data: a cube of data in which one dimension (in this case the + number of frames in the combiner) is much smaller than the number of + pixels. Extrema clipping @@ -106,7 +111,7 @@ rejected, loop in the code calling the clipping method: >>> old_n_masked = 0 # dummy value to make loop execute at least once >>> new_n_masked = combiner.mask.sum() >>> while (new_n_masked > old_n_masked): - ... combiner.sigma_clipping(func=np.ma.median) + ... combiner.sigma_clipping(func="median") ... old_n_masked = new_n_masked ... new_n_masked = combiner.mask.sum()