Array-API fallback for Combiner.sigma_clipping when the namespace is not numpy - #1001
Conversation
Codecov Report✅ All modified and coverable lines are covered by tests. Additional details and impacted files@@ Coverage Diff @@
## main #1001 +/- ##
==========================================
+ Coverage 97.83% 97.96% +0.12%
==========================================
Files 9 9
Lines 1850 1912 +62
==========================================
+ Hits 1810 1873 +63
+ Misses 40 39 -1
Flags with carried forward coverage won't be shown. Click here to find out more. ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
There was a problem hiding this comment.
Pull request overview
Adds an Array API-native sigma-clipping fallback while preserving Astropy behavior for NumPy.
Changes:
- Implements backend-native clipping and MAD-based deviation.
- Fixes keyword forwarding and adds extensive backend tests.
- Updates documentation, changelog, and escape tracking.
Reviewed changes
Copilot reviewed 6 out of 6 changed files in this pull request and generated 4 comments.
Show a summary per file
| File | Description |
|---|---|
ccdproc/combiner.py |
Implements fallback and dispatch logic. |
ccdproc/tests/test_combiner.py |
Adds cross-backend regression coverage. |
ccdproc/tests/array_escape_baseline.txt |
Removes the resolved NumPy escape. |
docs/array_api.rst |
Documents fallback limitations. |
docs/image_combination.rst |
Updates recommended clipping usage. |
CHANGES.rst |
Records the feature and keyword fix. |
💡 Add a code-review agent skill or configure MCP servers for context-aware, tailored reviews. Learn more in the docs.
| "astropy.stats.sigma_clip (such as grow, masked or " | ||
| "return_bounds) and are only available when the array " | ||
| "namespace is numpy." |
There was a problem hiding this comment.
Checked, and you are right about both: on the numpy path masked=... raises TypeError: sigma_clip() got multiple values for keyword argument 'masked' and return_bounds=True raises AttributeError: 'tuple' object has no attribute 'mask' (both pre-existing; main has the same masked=True, **kwd call). They are fixed by the wrapper, not astropy-only.
Will do: the message and the kwd docstring name grow as the example of astropy-only keywords, and masked/return_bounds are rejected with a clear TypeError on every namespace instead of the accidental errors above.
— Written by Claude at @mwcraig's direction.
| ``axis`` (default ``0``) and ``maxiters`` (default ``1``) are | ||
| honoured for every array namespace. ``copy`` (default `False`) |
There was a problem hiding this comment.
Agreed. Astropy's sigma_clip takes None or a tuple as well (checked through Combiner on numpy: axis=None and axis=(0, 1) both run), and _sigma_clip_mask raises NotImplementedError for those. Will do: this paragraph states that outside numpy axis must be a single integer and that None/tuples raise NotImplementedError.
— Written by Claude at @mwcraig's direction.
| :func:`~astropy.stats.sigma_clip`. A NaN-aware function such as | ||
| ``np.nanmedian`` also works; ``np.ma.median`` can be very slow in |
There was a problem hiding this comment.
Partly. Checked what actually happens when np.nanmedian reaches the fallback: on array-api-strict it fails (np.nanmedian returns an ndarray, which xp.expand_dims then rejects), on jax it works through a host round-trip via __array__, and on dask it dispatches to dask's own nanmedian through __array_function__. So "numpy-only" is not accurate either; what is true is that it is only guaranteed for NumPy data, which is the context of this page (the combiner here holds numpy CCDData).
Will do: "A NaN-aware function such as np.nanmedian also works for NumPy data". The second half of the suggestion, using the namespace's NaN-aware callable on other backends, is exactly what the string options already do (_default_median/_default_std pick bottleneck, the namespace's nanmedian, or ccdproc's fallback), which is why the note says to prefer them; no change there.
— Written by Claude at @mwcraig's direction.
| result, but the astropy-only options (``grow``, ``masked``, | ||
| ``return_bounds``, ...) are not available there and raise ``TypeError``. |
There was a problem hiding this comment.
Agreed on the substance, see the reply on the TypeError message in combiner.py: masked and return_bounds are consumed by the wrapper and passing them on numpy already fails, just less clearly. Will do: this bullet says "the astropy-only options such as grow are not available there and raise TypeError", and masked/return_bounds are documented as not accepted by sigma_clipping at all (the wrapper always asks for the mask).
— Written by Claude at @mwcraig's direction.
mwcraig
left a comment
mwcraig
left a comment
There was a problem hiding this comment.
Verdict: needs small changes. The design and the core algorithm are sound. Two things should change before merge: float32 promotion (B1, inline) and the "reproduces astropy's result" / "merges cleanly" claims (C1 inline, B2 below). The rest is simplification (S1-S4) and test trimming (inline on the test file).
Verified
- numpy path unchanged: the same nine arguments reach
astropy.stats.sigma_clip; thekwd.popchange only affects the case that previously raisedTypeError: got multiple values. - Dispatch (astropy iff numpy), the
TypeErrorfor astropy-only kwargs (raised before any work), the mask OR, and the namespace/device/dtype/shape of the result are correct on numpy, array-api-strict (device1), jax and dask. - dask: an integer
maxitersstays one graph (1 graph execution formaxiters=3vs 6 formaxiters=Noneover five iterations). jax:jax.jitwith staticmaxiters=3traces and runs;maxiters=NoneraisesTracerBoolConversionError, as the docstring says. - Pre-existing-mask handling is unchanged from main and from 2.5.1 (
sigma_clip(self.data_arr.data, ...)never saw the mask either), so the new Note is accurate. - Bound comparisons are strict
</>on both sides, as in astropy's C loop. - The escape-baseline removal is right:
numpy.asanyarrayis now reached only in the numpy branch, and the astropy reference calls in the tests are blamed on test frames.
B2 - PR description: "merge cleanly in either order" - git merge-tree --write-tree sigma-clipping-array-api sigma-func-array-api reports CONFLICT (content) in CHANGES.rst and docs/array_api.rst. Trivial, but whichever lands second needs a rebase; please adjust the description.
C2 - upstream, not this PR: astropy 8.0.1's compiled sigma_clip bus-errors nondeterministically on the (unchanged) numpy path. faulthandler places it in _sigmaclip_fast -> _sigma_clip_fast; it crashed one scan in three on
x = np.array([[0.8, -0.8, 1.1, 0.7, -0.2, 0.5],
[-0.3, 0.2, 1.0, 0.3, 1.2, -0.4],
[-0.1, -0.9, 0.1, 0.4, -2.4, -1.3]], dtype=np.float32)
sigma_clip(x, sigma=2, cenfunc="mean", stdfunc="mad_std", maxiters=3, axis=0, masked=True, copy=False)and replayed 20x in isolation without crashing. Consistent with wirth_median(buffer, count=0) in compute_bounds.c after a slice is fully rejected (there is no count == 0 guard). Pre-existing; the array-API fallback is immune. Worth an upstream report.
Changelog placement: the entry is under Bug Fixes; the fallback is a feature and only the axis/copy/maxiters fix is a bug - consider splitting as #999 did.
docs/image_combination.rst:70-72 (not in the diff): "the mean (ignoring any masked pixels) ... the standard deviation (again ignoring any masked values)" is false on both paths and now contradicts the Note added at combiner.py:723. Worth fixing while this section is being touched.
Test runs (test_combiner.py, this branch)
| backend | passed | failed | skipped |
|---|---|---|---|
| numpy | 661 | 0 | 6 |
| array-api-strict (device1) | 658 | 7 (all median_combine -> sigma_func -> astropy MAD np.asanyarray; #1000's scope) |
2 |
| jax (X64) | 665 | 0 | 2 |
| dask | 665 | 0 | 2 |
test_sigma_clip_mask_matches_astropy: 540 cases per backend - numpy 1.3 s, strict 2.0 s, jax 4.9 s, dask 17.2 s.
Random trials vs astropy.stats.sigma_clip(..., masked=True).mask (1-3-D shapes, every axis, tie-rich integer grids, NaN/inf, all-NaN and constant slices, float64/float32/int32/int64, thresholds incl. None/0, maxiters 1-5/None): numpy 44/5000, strict 29/3000 (device1 subset 11/600), jax 7/600, dask 9/800 mismatches. Every one is either B1 (float32, fixable) or C1 (exact ties, intrinsic) - see the inline comments. Final bounds agree with return_bounds=True in 1982/2000 float64 trials; the 18 exceptions are the C1 ties.
| if not xp.isdtype(data.dtype, "real floating"): | ||
| # The namespace default rather than float64: jax without X64 has no | ||
| # float64 and warns when one is requested. | ||
| info = xp.__array_namespace_info__() | ||
| data = xp.astype(data, info.default_dtypes(device=device)["real floating"]) |
There was a problem hiding this comment.
B1 - float32 is clipped in float32; astropy clips it in float64. This only promotes non-floating dtypes, but astropy's _sigmaclip_fast hands float32 data to a double-precision gufunc, so its bounds are float64. Near a bound the masks differ, on the real off-numpy path:
import numpy as np, array_api_strict as xps
from astropy.stats import sigma_clip
from ccdproc.combiner import _sigma_clip_mask
col = np.array([-0.4, 1.0, 0.4, 1.3, -0.8], dtype=np.float32)[:, None]
sigma_clip(col, sigma=1.5, maxiters=1, cenfunc="median", stdfunc="std", axis=0, masked=True).mask[:, 0]
# [False False False False True]
_sigma_clip_mask(xps.asarray(col), sigma_lower=1.5, sigma_upper=1.5, maxiters=1, xp=xps)[:, 0]
# [False False False False False]
_sigma_clip_mask(xps.asarray(col.astype(np.float64)), sigma_lower=1.5, sigma_upper=1.5, maxiters=1, xp=xps)[:, 0]
# [False False False False True]In 12 000 float32 trials on rounded (tie-rich) data: 8 mismatches, all 8 gone after upcasting; multi-iteration clips amplify it because one flipped value changes the next iteration's statistics. So the same float32 stack gets different masks on numpy and on every other backend.
Fix: promote unconditionally to the namespace default real dtype (drop the if). That is float64 everywhere except jax without X64, where nothing better exists. The differential test never sees this because every dataset in _sigma_clip_datasets is float64 - please add a float32 one.
There was a problem hiding this comment.
Agreed; the strict reproduction matches. Will do: promote unconditionally (right after the _setup call from S1) and add a float32 set to _sigma_clip_datasets.
One refinement to "the namespace default real dtype": a bare astype(default) narrows float64 input on a namespace whose default is float32 (torch; jax without X64 cannot hold float64 anyway). Taking the wider of the two avoids that and still gives float64 everywhere the default is float64:
info = xp.__array_namespace_info__()
default = info.default_dtypes(device=device)["real floating"]
if not xp.isdtype(data.dtype, "real floating"):
data = xp.astype(data, default)
data = xp.astype(data, xp.result_type(data.dtype, default)) # widen, never narrow(array_api_strict.result_type(float32, float64) is float64; the second astype is a no-op when nothing changes.)
— Written by Claude at @mwcraig's direction.
| # Same axis rules as the NaN-aware reduction fallbacks in _nanfuncs. | ||
| if axis is None or isinstance(axis, bool): | ||
| raise NotImplementedError( | ||
| "sigma clipping outside numpy supports only a single integer axis." | ||
| ) | ||
| try: | ||
| axis = operator.index(axis) | ||
| except TypeError: | ||
| raise NotImplementedError( | ||
| "sigma clipping outside numpy supports only a single integer axis." | ||
| ) from None | ||
| ndim = data.ndim | ||
| if not -ndim <= axis < ndim: | ||
| raise ValueError(f"axis {axis} is out of bounds for array of dimension {ndim}") | ||
| axis = axis % ndim |
There was a problem hiding this comment.
S1 - this block is _nanfuncs._setup. _setup(x, axis, xp) (_nanfuncs.py:29-94) applies the same rules - None/bool rejected, operator.index, out-of-bounds ValueError, % ndim, namespace and device resolution, and the same default_dtypes(...)["real floating"] promotion - and its messages contain the same "single integer axis" / "out of bounds" substrings the tests match. data, axis, xp, device = _setup(data, axis, xp) replaces lines 256-294 (~25 lines) and gives B1's fix one home (a keyword on _setup, or one astype after it). test_sigma_clip_mask_bad_axis then shrinks to the negative-axis check; the validation itself is already covered by test_nanfuncs.py::test_bad_axis.
There was a problem hiding this comment.
Agreed, it is _setup with the messages reworded. Will do: data, axis, xp, device = _setup(data, axis, xp) followed by the unconditional promotion from B1, kept in _sigma_clip_mask rather than added as a _setup keyword so that the _nanfuncs reductions keep returning float32 for float32 input. test_sigma_clip_mask_bad_axis shrinks to the negative-axis and maxiters=0 checks (they move into the argument-handling test from the test-file thread).
Side effect worth stating: the NotImplementedError text becomes _setup's "NaN-aware reduction fallbacks support only a single integer axis", which is less pointed for someone calling Combiner.sigma_clipping(axis=None) off numpy. The single-integer restriction gets spelled out in the sigma_clipping docstring (Copilot's comment on the kwd paragraph), so it is discoverable from the method the user actually called.
— Written by Claude at @mwcraig's direction.
| def _resolve_sigma_clip_func(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``. | ||
| TypeError | ||
| If ``func`` is neither a string nor callable. | ||
| """ | ||
| if callable(func): | ||
| return func | ||
| if isinstance(func, str): | ||
| 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 | ||
| raise TypeError(f"{kind} must be a string or a callable, got {type(func).__name__}") | ||
|
|
There was a problem hiding this comment.
S2 - 38 lines that can be ~6. There is no other string->callable dispatch in ccdproc for this to share, so the general helper buys nothing:
def _resolve(func, options, kind, xp):
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_default_mad_std can go too ("mad_std": lambda xp: partial(_nanmadstd, xp=xp) in the table). The non-str/non-callable TypeError branch is a corner astropy itself does not guard; test_sigma_clip_mask_bad_string loses one assertion.
There was a problem hiding this comment.
Will do, as written: _resolve replaces _resolve_sigma_clip_func, _default_mad_std becomes the lambda xp: partial(_nanmadstd, xp=xp) table entry, and the TypeError assertion in test_sigma_clip_mask_bad_string goes. A non-string non-callable then gets "stdfunc must be one of ['mad_std', 'std'] or a callable, got 3", which says everything the TypeError did.
— Written by Claude at @mwcraig's direction.
| ``1.4826 * median(|x - median(x)|)`` along ``axis``, with that axis | ||
| removed. | ||
| """ | ||
| # TODO: use _nanfuncs.<name> once PR A (sigma_func fallback) merges |
There was a problem hiding this comment.
S3 - duplicate of #1000's core._mad_fallback. The companion PR adds _mad_fallback(data, axis, ignore_nan, xp=, mask=) in core.py with the same promotion and the same 1.482602218505602; after both merge one of the two should go. Trade-off to settle there: _mad_fallback always uses the sort-based _nanfuncs medians, _nanmadstd uses the tiered _default_median (native nanmedian on jax/dask). The better direction is probably for _mad_fallback to adopt _default_median and for this to become 1.4826... * _mad_fallback(x, axis, ignore_nan=True, xp=xp) (combiner already imports from core). Either way, name #1000 here rather than "PR A".
There was a problem hiding this comment.
Will do: the TODO names #1000.
On the consolidation, one constraint your sketch skips: _default_median lives in combiner.py, and combiner imports from core, so _mad_fallback cannot import it without a cycle. Either the _default_* tier moves out of combiner.py or _mad_fallback takes the median as a parameter, e.g. _mad_fallback(x, axis, ignore_nan=True, xp=xp, median=_default_median(xp)), with _nanmadstd reduced to 1.482602218505602 * _mad_fallback(...). I prefer the parameter (smaller diff, _mad_fallback keeps its sort-based default for sigma_func). Since it changes #1000's function, it belongs in whichever of the two lands second, in the rebase B2 forces anyway.
— Written by Claude at @mwcraig's direction.
| iteration = 0 | ||
| while maxiters is None or iteration < maxiters: | ||
| iteration += 1 | ||
| 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) |
There was a problem hiding this comment.
C4 - no early exit for an integer maxiters (by design, so dask stays one graph): maxiters=10 does ten rounds of reductions on eager backends even when the clip converged after two, whereas astropy's C loop returns as soon as nothing changes. Results are identical; the docstring already says "runs exactly that many iterations" - consider adding "unlike astropy, which stops early once nothing changes" so nobody expects the early return.
S4 (optional): for _ in (range(maxiters) if maxiters else itertools.count()): drops the manual iteration counter.
There was a problem hiding this comment.
Will do both: the Notes get "unlike astropy, which returns as soon as an iteration rejects nothing", and the loop becomes for _ in (range(maxiters) if maxiters else itertools.count()): (maxiters is already validated to be None or >= 1 just above, so the truthiness test is safe).
— Written by Claude at @mwcraig's direction.
| (O(n log n) along the combination axis rather than O(n)). | ||
| + ``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 |
There was a problem hiding this comment.
"reproduces astropy's result" - see the C1 comment on _sigma_clip_mask: qualify with the rounding caveat (exact ties; float32 until B1 is fixed).
There was a problem hiding this comment.
Will do: "reproduces astropy's result" becomes "reproduces astropy's result up to floating-point rounding of the reductions (a value lying exactly on a bound can be classified differently)". The float32 half of the caveat goes away with the B1 fix, so it will not be mentioned here.
— Written by Claude at @mwcraig's direction.
| instead of ``len``, which the array API standard does not provide. [#999] | ||
| - ``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`` (NumPy data still use astropy); |
There was a problem hiding this comment.
Same caveat as docs/array_api.rst (C1 comment on _sigma_clip_mask): "reproduces" is true up to floating-point rounding of the reductions.
There was a problem hiding this comment.
Will do: same qualification as in docs/array_api.rst, in the feature entry (see the changelog-placement reply on the PR for the split).
— Written by Claude at @mwcraig's direction.
| def _mad_dev_func(x, axis=0): | ||
| """``median_absolute_deviation`` for whichever backend is under test.""" | ||
| if array_api_compat.is_numpy_namespace(xp): | ||
| return mad(x, axis=axis) | ||
| center = xp.expand_dims(nanmedian(x, axis=axis, xp=xp), axis=axis) | ||
| return nanmedian(xp.abs(x - center), axis=axis, xp=xp) | ||
|
|
||
|
|
||
| @pytest.mark.parametrize( | ||
| "dev_func", [_mad_dev_func, "mad_std"], ids=["callable", "mad_std"] | ||
| ) | ||
| def test_combiner_sigmaclip_high(dev_func): |
There was a problem hiding this comment.
Test LOC: 501 added lines can be ~270 with the same set of distinct branches covered. Starting here: _mad_dev_func and the [callable, "mad_std"] parametrize of the three old tests add ~30 lines to cover a callable dev_func through Combiner, which combine(sigma_clip=True) (callables xp.mean/xp.std) and the grid already cover - switch the three tests to dev_func="mad_std" and drop the helper.
| test | ~LOC | unique path | proposal |
|---|---|---|---|
_mad_dev_func + parametrize (:546-620) |
+30 | callable dev_func via Combiner |
dev_func="mad_std" only (-25) |
_sigma_clip_datasets (:1669) |
58 | data | keep; trim the docstring that repeats the inline comments; add a float32 set (-10) |
_sigma_clip_reference (:1729) |
38 | astropy 8.1 union-vs-bounds | keep |
grid comment + test_sigma_clip_mask_matches_astropy (:1769) |
77 | differential | keep; 5-line comment; zip (cenfunc, stdfunc) into 3 pairs; (None, 0) instead of (3, 3); simplefilter("error") around the _sigma_clip_mask call (skip on dask). 540 -> 180 cases, dask 17 s -> ~6 s (-20) |
test_sigma_clip_mask_is_silent_on_nan |
21 | warning silence | folded into the grid (-21) |
test_sigma_clip_mask_none_threshold_means_three |
20 | None/0 -> 3; xp=None |
(None, 0) in the grid; one xp=None call in the argument-handling test (-18) |
test_sigma_clip_mask_final_bounds_not_union |
33 | none (see separate comment) | drop (-33) |
accepts_numpy_scalar_thresholds / bad_axis / bad_string |
36 | strict scalar conversion; axis normalisation; bad strings | one test_sigma_clip_mask_argument_handling (~18 lines); axis validation lives in test_nanfuncs::test_bad_axis once S1 is done (-18) |
test_nanmadstd_matches_astropy |
12 | axis != 0 for the helper | keep |
uses_astropy_only_for_numpy / accepts_axis_and_maxiters / rejects_astropy_only_kwargs / use_astropy_still_deprecated / keeps_existing_mask / fallback_branch_on_any_backend |
114 | spy iff numpy; axis/maxiters regression; TypeError off numpy; deprecation; mask OR; forced fallback | one parametrized test, ~45 lines (-69) |
test_combine_sigma_clip_on_any_backend |
24 | combine(sigma_clip=True) |
keep |
Sketch of the consolidated dispatch test (every branch retained):
@pytest.mark.parametrize("force_fallback", [False, True])
def test_sigma_clipping_dispatch(monkeypatch, force_fallback):
if force_fallback:
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())
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)
uses_astropy = array_api_compat.is_numpy_namespace(xp)
assert (len(calls) == 1) is uses_astropy
if uses_astropy:
assert calls[0]["masked"] is True and calls[0]["maxiters"] == 2 and calls[0]["copy"] is False
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
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)))
if uses_astropy:
c.sigma_clipping(grow=1)
else:
with pytest.raises(TypeError, match="grow"):
c.sigma_clipping(grow=1)Grid parametrization:
@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", "mad_std"), ("callable", "callable")])
@pytest.mark.parametrize("dataset", sorted(_sigma_clip_datasets())) # incl. a float32 set
...
with warnings.catch_warnings():
if not array_api_compat.is_dask_namespace(xp):
warnings.simplefilter("error")
result = _sigma_clip_mask(...)Each table entry (median, mean, std, mad_std) and the callable path are still each exercised on every dataset x threshold x maxiters; (None, 0) exercises None/0 -> 3 on both sides of the reference.
There was a problem hiding this comment.
Mostly agreed; row by row, with one pushback on the grid pairs.
_mad_dev_func+ parametrize: drop.combine(sigma_clip=True)forwardsxp.mean/xp.stdas callables (combiner.py:1376-1379) and the grid runs_sigma_clip_maskwith callables, so a callabledev_functhroughCombineris covered twice already. The three tests usedev_func="mad_std"._sigma_clip_datasets: keep, trim the docstring, add the float32 set (B1).- Grid:
(None, 0)in place of(3, 3), yes.simplefilter("error")around the_sigma_clip_maskcall with the dask skip, yes: checked that 0 of 300 cases raise on numpy, jax and array-api-strict with the pairs below, sotest_sigma_clip_mask_is_silent_on_nanfolds in. Pairs: pushback on three.("mean", "std")is the defaultfunc/dev_funcpair ofCombiner.sigma_clipping, and with your three pairs it would appear in no differential test at all (the consolidated dispatch test usesmedian/mad_std,combineuses callables). Four pairs,("median", "std"),("mean", "std"),("mean", "mad_std"),("callable", "callable"): 6 datasets x 4 x 4 x 3 = 288 cases (from 540; your 180 is 216 once the float32 set is in). The union-vs-final corner survives the trimming: it shows up in everymad_stdpairing (84 cases today, see the other thread). none_threshold_means_three: drop;(None, 0)covers the semantics against astropy and thexp=Nonecall moves to the argument-handling test.final_bounds_not_union: drop (other thread).accepts_numpy_scalar_thresholds/bad_axis/bad_string-> onetest_sigma_clip_mask_argument_handling(numpy scalars,xp=None, negative axis,maxiters=0, bad strings); axis validation istest_nanfuncs::test_bad_axis's job once S1 lands.- Dispatch tests -> your
test_sigma_clipping_dispatchsketch, with two additions: keepcalls[0]["axis"] == 0in the numpy branch (the regression this PR fixes is precisely those three keywords reaching astropy once), and the C3 comment in theforce_fallbackbranch. Theuse_astropydeprecation warning is raised by the decorator, sopytest.warnsin the sketch is right on both branches. test_nanmadstd_matches_astropy,_sigma_clip_reference,test_combine_sigma_clip_on_any_backend: keep.
I have not timed the dask run after the cut; 540 -> 288 cases should land it near half of 17 s rather than 6 s.
— Written by Claude at @mwcraig's direction.
| assert bool(xp.all(result == reference)) | ||
|
|
||
|
|
||
| def test_sigma_clip_mask_final_bounds_not_union(): |
There was a problem hiding this comment.
Drop this one (-33 lines): the grid with maxiters=3 on collapse/normal is compared against a bounds-derived reference, which is exactly the final-bounds-not-union property, and the worked example already lives in the _sigma_clip_mask Notes.
There was a problem hiding this comment.
Dropping it, after checking the premise rather than assuming it: over the 360 grid cases with maxiters 3 or None, the union-of-iterations mask differs from the final-bounds mask in 84 of them (every mad_std pairing, on all five datasets, 1-14 elements each), and each of those is asserted against the bounds-derived reference, so the property is exercised. With std the bounds never widen on these datasets, which is all the hand-built column adds, and that worked example is already in the _sigma_clip_mask Notes.
— Written by Claude at @mwcraig's direction.
| assert bool(c._data_arr_mask[3, 1, 2]) | ||
|
|
||
|
|
||
| def test_sigma_clipping_fallback_branch_on_any_backend(monkeypatch): |
There was a problem hiding this comment.
C3: monkeypatching array_api_compat.is_numpy_namespace globally also flips _default_median/_default_std (combiner.py:31, 76), so on numpy this runs array_api_compat.numpy.nanmedian rather than bottleneck. Harmless for what is asserted, but it is not the reductions numpy users get - worth a one-line comment, or patch only the dispatch predicate if the consolidated test above keeps this branch.
There was a problem hiding this comment.
Correct: _default_median/_default_std (combiner.py:31, :76) read the same array_api_compat.is_numpy_namespace, so on numpy the forced fallback clips with array_api_compat.numpy.nanmedian/nanstd instead of bottleneck. Patching only the dispatch would need the dispatch to go through a named predicate (both sites look the function up on the same module object at call time), which is more machinery than the test warrants. Will do: a one-line comment in the force_fallback branch of the consolidated test saying the reductions differ from what numpy users get and that only the dispatch is under test.
— Written by Claude at @mwcraig's direction.
|
Replies to the points in the body of review 5061416614 that have no inline thread:
— Written by Claude at @mwcraig's direction. |
| return partial(nanstd, xp=xp) | ||
|
|
||
|
|
||
| def _nanmadstd(x, /, *, axis=0, xp): |
There was a problem hiding this comment.
Shouldn't this be in nanfuncs?
There was a problem hiding this comment.
The computation, yes; the tiering, no. The MAD itself is the kind of pure array-API reduction _nanfuncs holds, and _nanfuncs is the one module both combiner and core can import (it imports only array_api_compat; combiner imports core, so #1000's core._mad_fallback can never reach combiner._default_median, the cycle from the S3 thread). What does not fit there is the median choice: _nanmadstd picks it with _default_median(xp) (bottleneck / native nanmedian / fallback), while every _nanfuncs function is the fallback tier itself, single axis, no dispatch.
So, in this PR: the pure part moves to _nanfuncs with the median as a parameter, and only the one-line tiered wrapper stays here.
# _nanfuncs.py
def nanmad(x, /, *, axis=0, xp=None, median=None):
"""
Median absolute deviation along ``axis``, ignoring NaNs.
``median`` is the reduction used for both medians, ``median(x, axis=)``;
the default is `nanmedian`.
"""
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)
# combiner.py
def _nanmadstd(x, /, *, axis=0, xp):
return 1.482602218505602 * nanmad(x, axis=axis, xp=xp, median=_default_median(xp))test_nanmadstd_matches_astropy becomes a test_nanfuncs.py check of nanmad against median_absolute_deviation(ignore_nan=True); the _default_median wiring is already exercised by every mad_std case of the grid.
_mad_fallback then keeps its None/tuple axis and mask handling and calls nanmad(..., median=partial(_nanfuncs.median, xp=xp)) when ignore_nan is false, dropping its own promotion and both median calls. That edits #1000's function, so it stays where the S3 thread put it: whichever PR lands second, now one step shorter.
— Written by Claude at @mwcraig's direction.
| 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'``. |
There was a problem hiding this comment.
Why is the default changed here?
There was a problem hiding this comment.
It isn't; only the docstring changes. The signature is byte-identical on main (ccdproc/combiner.py:430-432) and on this branch (:661-663):
def sigma_clipping(
self, low_thresh=3, high_thresh=3, func="mean", dev_func="std", **kwd
):The runtime default has been the mean since the first commit of the class (0835aee, 2014: func=ma.mean). 5d2abc0 (2022-08-18, "Always use astropy's sigma clipping and default to string-based functions") turned that into func='mean' and, in the same commit, added the sentence "The default is 'median'" two lines below it; git log -S"func='median'" and -S'func="median"' on the file find no commit, ever. combine() resolves its own defaults to xp.mean/xp.std (:1098-1101 on main) as well. So this line makes the docstring say what the code has always done; the dev_func sentence ("The default is 'std'") was already right. I can add a changelog line if you want the correction called out.
— Written by Claude at @mwcraig's direction.
| @pytest.mark.parametrize( | ||
| "dev_func", [_mad_dev_func, "mad_std"], ids=["callable", "mad_std"] | ||
| ) | ||
| def test_combiner_sigmaclip_high(dev_func): |
There was a problem hiding this comment.
This needs a comment describing what it does
There was a problem hiding this comment.
Will do. Per the thread above, the [callable, "mad_std"] parametrize and _mad_dev_func go and the test takes dev_func="mad_std", so the comment goes on the plain test, replacing "using mad for more robust statistics vs. std":
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.(test_combiner_sigmaclip_low gets the mirror-image comment.)
— Written by Claude at @mwcraig's direction.
| @pytest.mark.parametrize( | ||
| "dev_func", [_mad_dev_func, "mad_std"], ids=["callable", "mad_std"] | ||
| ) | ||
| def test_combiner_sigmaclip_single_pix(dev_func): |
There was a problem hiding this comment.
Will do; same shape as :557 once the parametrize is gone. The existing "add a single pixel in another array..." comment is replaced by:
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).— Written by Claude at @mwcraig's direction.
|
|
||
|
|
||
| def _sigma_clip_reference(np_data, **kwargs): | ||
| """ |
There was a problem hiding this comment.
Make this a real docstring
There was a problem hiding this comment.
Taking this as "numpydoc, like the rest": summary, Parameters, Returns, with the why-bounds explanation demoted to Notes.
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 current paragraph: compiled path vs python loop, copy=False,
astropy 8.1 / astropy#19858, why the bounds rather than the mask>
"""If you meant the # This is a differential test block above test_sigma_clip_mask_matches_astropy instead: that becomes the test's docstring, at the 5-line length agreed in the :557 thread.
— Written by Claude at @mwcraig's direction.
| assert bool(xp.all(result == expected)) | ||
|
|
||
|
|
||
| def test_sigma_clip_mask_bad_axis(): |
There was a problem hiding this comment.
Will do. Per the :557 thread this test merges with accepts_numpy_scalar_thresholds and bad_string into one test_sigma_clip_mask_argument_handling, so the docstring goes on that:
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``.
"""— Written by Claude at @mwcraig's direction.
| assert bool(xp.all(xpx.isclose(result, expected, equal_nan=True))) | ||
|
|
||
|
|
||
| def _sigma_clip_ccd_list(): |
There was a problem hiding this comment.
Will do:
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.
"""— Written by Claude at @mwcraig's direction.
Fixes a real bug the review found: _sigma_clip_mask only promoted non-floating input, so float32 data was clipped in float32 while astropy's compiled path always computes in float64, producing a different mask near a bound. It now widens (never narrows) to the namespace default whenever both cenfunc and stdfunc are strings, matching astropy's own promotion (a callable takes astropy's python loop instead, which keeps the data's own dtype, so this leaves a callable's input alone too). masked and return_bounds now raise a clear TypeError on every namespace, including numpy, instead of the accidental "multiple values"/AttributeError they triggered before. Refactors, as sketched in review: - _sigma_clip_mask's axis/dtype/device setup is now _nanfuncs._setup, same as every other NaN-aware fallback. - _resolve_sigma_clip_func/_default_mad_std collapse into a single _resolve() helper. - The pure MAD computation moves to _nanfuncs.nanmad(median=...); _nanmadstd becomes a one-line tiered wrapper around it. The TODO to dedup with astropy#1000's _mad_fallback stays, updated for nanmad. - The manual iteration counter becomes a for loop over range(maxiters) or itertools.count() when maxiters is None. Docs/docstrings: add the floating-point rounding caveat, spell out that axis outside numpy must be a single integer, document that masked/return_bounds are never accepted, and drop the stale "negative value" sentence (that's astropy#1002's question). Tests: drop the _mad_dev_func parametrize (redundant with the combine() callable coverage and the grid); add a float32 dataset; collapse the (cenfunc, stdfunc) x sigma-pair grid to four pairs with (None, 0) replacing a redundant (3, 3), wrapped in simplefilter("error") (skipped on dask) so it also covers silence on NaN; drop three tests whose properties the grid now covers; merge the argument-handling and dispatch tests each into one parametrized test; retarget the _nanmadstd differential test at _nanfuncs.nanmad. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01LVxnTLrWKStBdPwcmDNhxA
Split the changelog entry as astropy#999 did: the array-API sigma-clipping fallback stays under New Features (with the rounding caveat), the axis/copy/maxiters TypeError fix moves under Bug Fixes, and a new line records the sigma_clipping docstring's default-value correction (it wrongly said 'median'; the runtime default has always been 'mean'). docs/array_api.rst: state the rounding caveat, note that grow is an example of the astropy-only options that raise TypeError off numpy, and that masked/return_bounds are never accepted by sigma_clipping at all, on any array library. docs/image_combination.rst: the sigma_clipping paragraph no longer implies masked pixels are excluded from the statistics (they are not, on either code path); the np.nanmedian callable example is qualified as only guaranteed for NumPy data. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01LVxnTLrWKStBdPwcmDNhxA
Combiner.sigma_clipping wrapped astropy.stats.sigma_clip unconditionally, which converts the data to numpy and fails outright on array-api-strict (ccdproc#929). NumPy users must not be slowed down, so numpy data still go to astropy exactly as before (masked=True, all extra keyword arguments forwarded, compiled fast path intact). Every other namespace now uses _sigma_clip_mask, an iterative clip written purely in terms of the array API standard that reproduces astropy's mask: the last iteration's bounds applied to the original data, non-finite values always masked, None/0 thresholds meaning 3 as in astropy. The 'median'/'mean'/'std'/'mad_std' string options map to the same tiered NaN-aware reductions the combination methods use (new private _nanmadstd for 'mad_std'); astropy-only keyword arguments raise TypeError off numpy. An integer maxiters runs without any host synchronisation, so the clip stays a single graph on dask. Also fixes a latent bug: axis, copy and maxiters were read with kwd.get and then forwarded again through **kwd, so passing any of them raised "TypeError: got multiple values for keyword argument". They are popped now. Docstring fixes: func defaults to 'mean' (not 'median'); None thresholds are treated as 3 rather than "no rejection"; kwd semantics per namespace. Docs: image_combination.rst recommends func="median" over np.ma.median; array_api.rst lists the limitations off numpy. The sigma_clipping escape baseline entry is deleted (confirmed by a full-suite dask regeneration). Verified: the design was checked against astropy.stats.sigma_clip(..., masked=True, copy=False).mask over 864 structured cases (numpy and array-api-strict on device1), 2400 random numpy trials, and 243 cases each on jax and dask, with 0 mismatches. The new differential test (5 data sets x 3 cenfunc x 3 stdfunc x 4 thresholds x 3 maxiters) runs on every backend. Full suite: numpy 1123 passed; array-api-strict 10 failed (down from 13; all remaining are ccdproc#929 sigma_func, ccdproc#936 and ccdproc#983), 0 xpassed; jax 1122 passed; dask 1116 passed. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01XGArZMnLPC1G7YP62vN18V
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01XGArZMnLPC1G7YP62vN18V
Record in the test file why _sigma_clip_mask is tested by comparison against astropy over a full cross product rather than a curated list (the corners of astropy's behaviour only appear for particular combinations of center, deviation, thresholds and maxiters), what each data set is for, and what the grid costs. Comments only; no behaviour change. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01XGArZMnLPC1G7YP62vN18V
astropy 8.1 (astropy#19858) changed the python loop that sigma_clip takes for a callable cenfunc/stdfunc: with masked=True it now masks the union of every iteration's rejections regardless of copy, where it used to apply the final bounds to the data for copy=False. Its compiled path (string options) still applies the final bounds, which is what _sigma_clip_mask reproduces. Derive the reference mask from return_bounds=True, which both astropy paths and versions agree on, and keep checking astropy's own mask against it on the compiled path. Docstrings note the corner where astropy >= 8.1's numpy path can differ for a callable. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01JnWCg95xE93SbhME52jxGJ
Fixes a real bug the review found: _sigma_clip_mask only promoted non-floating input, so float32 data was clipped in float32 while astropy's compiled path always computes in float64, producing a different mask near a bound. It now widens (never narrows) to the namespace default whenever both cenfunc and stdfunc are strings, matching astropy's own promotion (a callable takes astropy's python loop instead, which keeps the data's own dtype, so this leaves a callable's input alone too). masked and return_bounds now raise a clear TypeError on every namespace, including numpy, instead of the accidental "multiple values"/AttributeError they triggered before. Refactors, as sketched in review: - _sigma_clip_mask's axis/dtype/device setup is now _nanfuncs._setup, same as every other NaN-aware fallback. - _resolve_sigma_clip_func/_default_mad_std collapse into a single _resolve() helper. - The pure MAD computation moves to _nanfuncs.nanmad(median=...); _nanmadstd becomes a one-line tiered wrapper around it. The TODO to dedup with astropy#1000's _mad_fallback stays, updated for nanmad. - The manual iteration counter becomes a for loop over range(maxiters) or itertools.count() when maxiters is None. Docs/docstrings: add the floating-point rounding caveat, spell out that axis outside numpy must be a single integer, document that masked/return_bounds are never accepted, and drop the stale "negative value" sentence (that's astropy#1002's question). Tests: drop the _mad_dev_func parametrize (redundant with the combine() callable coverage and the grid); add a float32 dataset; collapse the (cenfunc, stdfunc) x sigma-pair grid to four pairs with (None, 0) replacing a redundant (3, 3), wrapped in simplefilter("error") (skipped on dask) so it also covers silence on NaN; drop three tests whose properties the grid now covers; merge the argument-handling and dispatch tests each into one parametrized test; retarget the _nanmadstd differential test at _nanfuncs.nanmad. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01LVxnTLrWKStBdPwcmDNhxA
Split the changelog entry as astropy#999 did: the array-API sigma-clipping fallback stays under New Features (with the rounding caveat), the axis/copy/maxiters TypeError fix moves under Bug Fixes, and a new line records the sigma_clipping docstring's default-value correction (it wrongly said 'median'; the runtime default has always been 'mean'). docs/array_api.rst: state the rounding caveat, note that grow is an example of the astropy-only options that raise TypeError off numpy, and that masked/return_bounds are never accepted by sigma_clipping at all, on any array library. docs/image_combination.rst: the sigma_clipping paragraph no longer implies masked pixels are excluded from the statistics (they are not, on either code path); the np.nanmedian callable example is qualified as only guaranteed for NumPy data. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01LVxnTLrWKStBdPwcmDNhxA
The consolidation agreed in the S3 review thread, landing in the PR that merged second: _mad_fallback delegates its center/abs/median tail to _nanfuncs.nanmad, the same computation Combiner's 'mad_std' option uses, passing its native-first median tier as the median parameter. The MAD arithmetic now lives in exactly one place. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01LVxnTLrWKStBdPwcmDNhxA
396fd55 to
c1700b9
Compare
codecov flagged the new TypeError as the one uncovered patch line: no test called sigma_clipping with masked or return_bounds. Assert the rejection in the dispatch test, on both its branches, since it applies on every namespace by design. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01LVxnTLrWKStBdPwcmDNhxA
Part of #929 (the
Combiner.sigma_clippinghalf). Companion PR: #1000 (thesigma_funchalf), now merged; this branch is rebased ontomainwith it (the predicted trivial conflicts were the adjacent bullets inCHANGES.rst/docs/array_api.rst).As agreed in the S3 review thread, the
_mad_fallback/_nanmadstddedup lands here, in the PR that merged second: the MAD computation lives once, in_nanfuncs.nanmad, and bothcore._mad_fallbackandCombiner's'mad_std'option call it with their own median tier. The rebase also picks up #1000's dependency floors (numpy ≥ 2.0, astropy ≥ 6.1), which fixes the oldest-deps job's one red cell here: astropy 6.0's python sigma-clip loop promoted float32 to float64 unconditionally (sigma_clipping.py:483), a corner 6.1 fixed, so the float32×callable grid case now agrees.Combiner.sigma_clippingwrappedastropy.stats.sigma_clipunconditionally, which converts the data to numpy and fails outright onarray-api-strict(3 of the 13 remaining strict failures).Policy (numpy is untouched): when the namespace is numpy, the data still go to
astropy.stats.sigma_clipexactly as before —masked=True, all extra keyword arguments forwarded, compiled fast path intact. Every other namespace uses a new private_sigma_clip_mask, an iterative clip written purely in terms of the array API that reproduces astropy's mask: the last iteration's bounds applied to the original data, non-finite values always masked,None/0thresholds meaning 3 as in astropy. The'median'/'mean'/'std'/'mad_std'string options map to the same tiered (bottleneck → native → fallback) NaN-aware reductions the combination methods use (new private_nanmadstdfor'mad_std'); astropy-only keyword arguments (grow,masked,return_bounds, …) raiseTypeErroroff numpy. An integermaxitersruns with no host synchronisation, so on dask the clip stays one graph.Latent bug fixed:
axis,copyandmaxiterswere read withkwd.getand then forwarded again through**kwd, so passing any of them raisedTypeError: got multiple values for keyword argument. They are popped now.Docstring fixes:
funcdefaults to'mean', not'median';Nonethresholds are treated as 3 (as astropy does) rather than "no rejection" — the docstring still described the pre-2.4 behaviour. TheNone/0/negative threshold semantics are a behaviour question, not part of this PR — see #1002.Tests: the three
test_combiner_sigmaclip_*tests are parametrized over a backend-generic MAD callable and"mad_std"and no longer use numpy-only.all(). A differential test runs_sigma_clip_maskon every backend againstastropy.stats.sigma_clip(..., masked=True, copy=False).maskover 5 data sets (outliers, NaN/inf incl. an all-NaN column, zero spread, a column that collapses mid-loop, ints) × 3 cenfuncs × 3 stdfuncs × 4 threshold pairs ×maxiters ∈ {1, 3, None}; plus tests for silence on all-NaN slices,None/0→ 3, final-bounds-not-union semantics, numpy scalar thresholds, bad axes/strings,_nanmadstdvsmad_std, a spy proving astropy is called iff numpy, theaxis/maxitersregression, astropy-only kwargs, theuse_astropydeprecation, existing-mask preservation, andcombine(sigma_clip=True)(previously not exercised at all). Design was prototype-verified against astropy over 864 structured cases (numpy + strictdevice1), 2400 random numpy trials and 243 cases each on jax and dask, 0 mismatches.Docs: changelog entry; one bullet in
docs/array_api.rst;docs/image_combination.rstnow recommendsfunc="median"overnp.ma.median. Thecombiner.py sigma_clippingline leaves the escape baseline (confirmed by a full-suite dask regeneration; the diff is that one line).Verified locally
sigma_functests for the companion PR, Array API: units/Quantity handling with non-numpy arrays #936, Consider marray as a uniform masked-array representation across array backends #983), 0 xpassed🤖 Generated with Claude Code
https://claude.ai/code/session_01XGArZMnLPC1G7YP62vN18V