Tracking: remaining work for the array-API migration

This is the canonical checklist for finishing the array-API (array-library-agnostic) migration. It supersedes #909 as the place to track status; #909 and its children (#910–#913) remain as the CuPy-specific findings.

**State as of 2026-08-25 (post-#989/#990):** the August batch (#972–#980, plus #985 CI coverage), #987 (NaN-aware sum/mean/std/median fallbacks, closing #986), #988 (combiner internals: `xp.any`/`xp.count_nonzero`, namespace bool mask dtypes, `nbytes`-free memory sizing), #989 (array-API `median` fallback for `subtract_overscan`, wrapping `nanmedian` with numpy's NaN-propagating semantics) and #990 (test-only hygiene batch addressing #969/#970 and the test-body numpy-isms) are all merged. #992 (combine-method `mask=` bucket, 63 → 54) and #993 (`gain_correct`/`flat_correct` `device=`, 54 → 41) are merged too. #995 (namespace dtypes: `ccd_process` bool, `transform_image` mask, `combine()`/`ImageFileCollection` FITS byte order) and #994 (`clip_extrema` rank-comparison rewrite) are merged too (2026-08-25, `6724c8e`): strict is now **31** on `main`, and every section-1 bucket that was on the strict failure list is closed. What is left on strict is 18 test-body fixes (section 2) and 13 upstream-blocked tests (section 3). CI exercises numpy, jax, dask (including dask+bottleneck and a dask escape-baseline gate); `array-api-strict` runs as a soft (`continue-on-error`) job. CuPy and torch are not in CI. User docs live in `docs/array_api.rst`. 36 `backend_xfail` markers remain in the test suite. **Later on 2026-08-25, three PRs are open on top of `6724c8e`:** #996 (narrows the two stale `test_cosmicray.py` markers, 5 XPASS → 0), #997 (#982 and the correlated add/subtract leak in the wrapper) and #998 (the section-2 test bodies, strict 31 → 15). Merged together locally they give strict **15 failed, 496 passed, 36 skipped, 45 xfailed, 0 xpassed**; numpy 557 passed / 35 skipped; jax and dask green apart from xfails. Of the 15, 11 are #929, one each is #936, #983 (`test_rebin_ccddata[True-True]`) and the two new ccdproc items below (`combine(output_file=)` writing a namespace array, `combine(dtype=int)`). **All three merged 2026-08-25 (`main` = `4d7d21b`); strict on `main` verified at exactly those 15.** #999 then took the two ccdproc items plus the `len(oscan)` one-liner; **merged 2026-08-26 (`main` = `be80055`)**. Strict on `main` verified at **13 failed / 503 passed / 36 skipped / 45 xfailed / 0 xpassed**, and every remaining failure is upstream-blocked: 8 × #929 via the default `sigma_func` (`median_absolute_deviation`), 3 × #929 via `Combiner.sigma_clipping` (`sigma_clip`), #936, #983. **Section 1 is complete.** #946 and #943 are closed. **2026-08-26, later:** the #929 half that is fixable on the ccdproc side is now in two PRs off `be80055`: #1000 (`sigma_func` → xp-native MAD when the namespace is not numpy; strict 13 → 5) and #1001 (`Combiner.sigma_clipping` → xp-native clip reproducing `astropy.stats.sigma_clip`; strict 13 → 10). Numpy keeps astropy in both. Merged together locally: strict **2 failed / 1129 passed / 0 xpassed** (#936, #983 only), the two `BOUNDARY` lines leave the escape baseline, and the dask enforce run is clean.

**State as of 2026-08-30 (post-#1005):** #1000, #1001 (with #1003's numpy>=2.0/astropy>=6.1 dependency bump) and #1005 are merged. Strict on `main` (`1bd015c`) is **0 failed / 839 passed / 10 skipped / 46 xfailed / 0 xpassed**, and the `py313-strict` job now runs in the regular CI matrix (the `continue-on-error` carve-out and `strict_status.yml` are removed; a strict failure fails CI). The two astropy-blocked tests carry `backend_xfail` markers citing #936 and #983. The escape baseline is 10 entries, all `BOUNDARY`; regenerating it from a full dask run is byte-identical and the enforce ratchet holds on dask and jax. Full-suite results per backend: numpy 890 passed; dask 879 passed (ratchet enforced); jax 885 passed (ratchet enforced). All 30 remaining `backend_xfail` markers were re-triaged 2026-08-30: every one is a genuine numpy-only boundary (astroscrappy, scipy.ndimage, reproject, astropy modeling/nddata/units).

Current strict result on `main` (1b9b621): **63 failed, 418 passed, 36 skipped, 46 xfailed, 5 xpassed** — down from 85 with #989 + #990. With #993 (`711bb26`) the strict job was 41; #995 and #994 (`6724c8e`) take it to **31 failed, 469 passed, 36 skipped, 43 xfailed, 5 xpassed** with no new failures; the 5 XPASSes are still the stale `test_cosmicray.py` markers. The job is `continue-on-error`, so read the pytest summary line in its log rather than the job conclusion. Every one of the 63 has been traced to a cause (2026-08-25):

| # | cause | where | section |
|---|---|---|---|
| 12 | `average_combine`/`sum_combine`/`median_combine` build the result with `CCDData(..., mask=mask)`; astropy's `NDData.mask` setter forces numpy and cannot convert a `device1` strict array | `combiner.py:726`, `:804`, `:616` | 1 |
| 9 | `gain_correct` calls `xp.asarray(gain_value)` without `device=` ("two different devices") | `core.py:968` | 1 |
| 2 | `flat_correct` calls `xp.asarray(flat_mean)` without `device=` (same class) | `core.py:1057` | 1 |
| 4 | `ccd_process` bad-pixel mask uses builtin `bool` instead of `xp.bool` | `core.py:383` | 1 |
| 4 | `clip_extrema` fancy-index assignment | `combiner.py:406` | 1 |
| 2 | `combine()` passes a numpy dtype (`ccd.data.dtype.type`) to the requested namespace's `asarray` — strict warns, `filterwarnings=error` fails | `combiner.py:1082` | 1 |
| 2 | uncertainty propagation for multiply/divide delegates to astropy's `_propagate_multiply_divide`, whose `np.abs(...)` (`nduncertainty.py:832`) densifies the strict array and the mixed expression then fails | `_ccddata_wrapper_for_array_api.py:305–330` | 1 |
| 1 | `transform_image` applies the user transform to the **bool** mask; strict rejects `10 * bool_array` (numpy silently promotes) | `core.py:1153` | 1 |
| 8 | `astropy.stats` densification: default `sigma_func` (`median_absolute_deviation`, 5) and `sigma_clipping` (`sigma_clip`, 3) | `core.py:1291`, `combiner.py:491` | 3 (#929) |
| 6 | test bodies call astropy `CCDData.multiply(...)` on strict data to build scaled inputs; astropy's `_prepare_then_do_arithmetic` does `np.result_type(ref, operand)` on a strict dtype ("Could not convert Array … to a NumPy dtype") — this is the NDData-arithmetic upstream item already in #940, not a ccdproc code path | `test_combiner.py` (`test_combiner_with_scaling*`, `test_combiner_result_dtype`, `test_combine_overwrite_output`) | 2 |
| 5 | `combine()` given FITS filenames / numpy `CCDData` with no `array_package=` returns numpy data, and the tests then compare it with strict arrays via `xp.all(xpx.isclose(...))`; `array_package` is opt-in by design, so the tests should pass it | `test_combine_average_fitsimages`, `test_combine_numpyndarray`, `test_combine_average_ccddata`, `test_combine_limitedmem_fitsimages`, `test_combine_limitedmem_scale_fitsimages` | 2 |
| 2 | test bodies `ccd.write()` a `device1` strict `CCDData`; `astropy.io.fits` cannot coerce it | `test_combiner_image_file_collection_input`, `test_combine_image_file_collection_input` | 2 |
| 2 | test body `xp.asarray(np_mgrid[...]) / 10.0` — integer strict array divided by a float | `test_ccdproc.py:304` (`test_subtract_overscan_model`) | 2 |
| 1 | test helper `_make_mean_scaler` uses `.mean()` / `np.ma.average` — **this is the "`Combiner.scaling` setter bug" from #990's deferred list; the setter is fine** | `test_combiner.py:52` (`test_3d_combiner_with_scaling`) | 2 |
| 1 | astropy `_arithmetic` decorator leaks `_config_ccd_requires_unit=False` (astropy/astropy#20268) | `test_generator_ccds_without_unit` | 2/3 |
| 1 | astropy `NDData.mask` setter (#983's blocker) | `test_rebin_ccddata[True-True]` | 2/3 |
| 1 | `Quantity` machinery on a `device1` strict array | `test_unit_mismatch_behaves_as_expected` | 3 (#936) |

So: **36 are ccdproc bugs fixable now** (section 1), **17 are test-body fixes** (section 2), and **10 are upstream** (#929 ×8, #936, astropy#20268) — plus the rebin/mask-setter one that is fixable with the `_mask` workaround if we want it. The 5 xpassed are all in `test_cosmicray.py`: the two `test_cosmicray_median_masked_column[masked_array-*]` cases (stale since #979) and the three `gain_apply=False` cases of `test_cosmicray_gain_correct_uncertainty_namespace` (the `gain_apply=True` cases still xfail, so that marker needs narrowing to the `True` parametrizations, not removing).

### 1. Fixable-now bugs in ccdproc (no upstream dependency)
- [x] #965 `Combiner.__init__` builds a nested `xp.asarray([...])`; needs `xp.stack` (44 strict failures — biggest single win)
- [x] #966 `subtract_dark(scale=True)` device-propagation bug
- [x] #967 `rebin()` uses the `.astype()` method instead of `xp.astype`
- [x] #968 `create_deviation` float × bool-mask multiplication (`core.py:459`)
- [x] #963 `background_deviation_box` discards the `xpx.at` result (silent no-op on immutable backends)
- [x] #962 Variance/InverseVariance uncertainty wrappers escape to numpy on jax (only StdDev is array-API-aware)
- [x] #932 dead code in `cosmicray_median` clobbers masked-input handling
- [x] #906 xp-native fallback median for backends without `scipy.ndimage`
- [x] #986 no array-API fallback for NaN-aware mean/sum/std -- fixed in #987, which also merged the #906 median fallback into `ccdproc/_nanfuncs.py`
- [x] `subtract_overscan` calls `xp.median` (`core.py:629`), which is not in the array API standard — fixed in #989 with a `median` fallback in `_nanfuncs.py`
- [x] `Combiner.average_combine`/`sum_combine`/`median_combine` pass `mask=` to the `CCDData(...)` constructor (`combiner.py:726`, `:804`, `:616`), so astropy's numpy-forcing `mask` setter runs on the result. `ccd_process` already sidesteps the same setter by assigning `nccd._mask` (`core.py:381`); do the same here. **12 strict failures — biggest remaining bucket** — **PR #992** (also fixes the same setter sites inside `combine()` and its numpy-only `mask.copy()`; the remaining members of the bucket now stop at #929 (`median_combine`'s default `sigma_func`) or at `test_writeable_after_combine`, which writes a `device1` array to FITS and belongs with the section-2 "write a numpy copy" tests)
- [x] `gain_correct` calls `xp.asarray(gain_value)` without `device=` (`core.py:968`), so the scalar lands on the default device and the multiply fails when the image is elsewhere. Surfaced by #989; **9 strict failures** (`test_gain.py` ×5, `test_gain_correct*` ×3, `test_ccd_process_does_not_change_input`) — **PR #993** (also makes a plain-number gain a `float`, since the standard does not promote `int64` with `float64`)
- [x] `flat_correct` has the same bug for `xp.asarray(flat_mean)` (`core.py:1057`). 2 strict failures (`test_flat_correct`, `test_flat_correct_norm_value`) — **PR #993** (the two tests then needed their flats created on the testing device; fixed in the same PR)
- [x] `ccd_process` builds its bad-pixel mask with `xp.asarray(mask, dtype=bool)` — Python builtin `bool` where `xp.bool` is required (`core.py:383`; same class as #968/#988's fixes). Surfaced by #989; 4 strict failures — **PR #995** (the same four tests then hit `log_meta._replace_array_with_placeholder`, which assumed anything without `__len__` is NDData and raised `AttributeError` on a strict array; fixed in the same PR, plus the test bodies now build a bool mask and compare with `==`)
- [x] `Combiner.clip_extrema` assigns through fancy indexing (`combiner.py:406`, `IndexError: Fancy indexing __setitem__ is not supported`); needs `xpx.at` or a `where`-based rewrite. 4 strict failures — **PR #994** (no scatter at all: `argsort` twice gives each image's per-pixel rank, and `ranks < nlow | ranks >= n - nhigh` is the mask; bit-identical to the old code on 8400 random trials incl. ties/NaNs)
- [x] `combine()` converts file-read data with `xp.asarray(ccd.data, dtype=ccd.data.dtype.type)` (`combiner.py:1082`, and the uncertainty line below it), handing a numpy dtype to the target namespace; strict warns and the suite's `filterwarnings = error` turns it into a failure. Let `xp.asarray` infer the dtype or map it through the namespace. 2 strict failures (`test_combine_ccd_with_uncertainty_and_mask_from_fits`) — **PR #995** (dropping `dtype=` alone is not enough: strict and jax reject big-endian `>f8` outright, so a `_native_numpy` helper converts to native byte order in numpy first; the same three sites in `ImageFileCollection.ccds`/`data` get the same fix, and its mask has to go through `_mask` because the public setter forces numpy). Only `[mean]` clears; `[function]` then stops at `_make_mean_scaler` (section 2)
- [x] `_CCDDataWrapperForArrayAPI` uncertainty: `_propagate_multiply`/`_propagate_divide` call `super()._propagate_multiply_divide`, and astropy's implementation uses `np.abs(...)` (`nduncertainty.py:832`), which densifies the strict operand; the wrapper needs to own that formula the way it owns add/subtract (upstream fix would be `abs()` instead of `np.abs`; worth a note in #940). 2 strict failures (`test_flat_correct_deviation`, `test_flat_correct_data_uncertainty`) — **PR #993** (the mixin now owns `_propagate_multiply_divide`; after the `device=` fix this bucket also sat behind the five `test_gain.py` failures; the two strict xfail markers in `test_ccddata_wrapper_for_array_api.py` covering it are dropped)
- [x] Same class of leak for `add`/`subtract` with `uncertainty_correlation != 0`: the wrapper still delegates to astropy's `_propagate_add_sub`, whose correlation term uses `np.sqrt(this * other)` (`nduncertainty.py:742`), so a correlated add/subtract through the wrapper fails on strict with `TypeError: Expected Array or Python scalar; got numpy.ndarray`. Uncorrelated add/subtract is fine. Not covered by any test and not on the strict failure list; found while covering #993's multiply/divide override. Fix is the same shape as #993 (the mixin owns the formula with `xp.sqrt`) — **PR #997** (`test_wrapped_arithmetic_correlated_uncertainty` now covers `add`/`subtract` too)
- [x] `transform_image` applies `transform_func` to the bool mask directly (`core.py:1153`); numpy promotes `10 * mask` silently, strict raises. Cast the mask to a numeric dtype before transforming (the `mask > 0` afterwards already assumes a numeric result). 1 strict failure (`test_transform_image[True-True]`) — **PR #995** (cast to the data dtype; the test body also discarded the `xpx.at(...).set()` result, same class as #963, and now checks the transformed mask's contents)
- [x] #982 `combine()` does not normalise a raw module passed as `array_package`, so `combine(files, array_package=dask.array)` fails; `Combiner.__init__` was fixed in #976 and the two entry points should agree — **PR #997** (the crash itself went inert with #995, whose `device=`-free conversion a raw `dask.array` accepts, and `combine()` re-derives `xp` from the converted array anyway; the normalisation now matches `Combiner`, an array is still accepted as `array_package`, and regression tests cover raw `numpy`/`dask.array`)
- [x] `combine(dtype=...)`/`Combiner(dtype=...)` hand the user's dtype straight to `xp.astype`/`xp.asarray` (`combiner.py:1104`, `:191`); a builtin Python `int` is a valid NumPy dtype but array-api-strict rejects it (`AttributeError: type object 'int' has no attribute '_np_dtype'`). 1 strict failure (`test_combiner_result_dtype`, whose `dtype=int` case is the point of the test), reachable since #998 cleared the `.multiply()` failure in front of it — **PR #999** (`core._namespace_dtype` resolves the name with `numpy.dtype` and looks it up on the namespace; the namespace's own dtype objects pass through)
- [x] `subtract_overscan` calls `len(oscan)` on the overscan array (`core.py:698`, `xp.arange(len(oscan))`); array-api-strict has no `__len__`. Masked on the strict job because `test_subtract_overscan_model` is xfailed on #933 (with `oscan.shape[0]` in its place the test stops in `astropy.modeling.fitting._convert_input`, so the marker is right), but it is a one-line ccdproc numpy-ism — **PR #999**
- [x] ~~`Combiner.scaling` setter has an array-API bug exposed when test bodies stop using numpy `.mean()`~~ — diagnosed 2026-08-25: not a setter bug. `test_3d_combiner_with_scaling` fails inside the test helper `_make_mean_scaler` (`test_combiner.py:52`, `.mean()`/`np.ma.average`), and `test_combiner_with_scaling` fails earlier on `ccd_data.multiply(3)` in the test body (astropy NDData arithmetic). Both moved to section 2
- [x] `Combiner._get_nan_substituted_data` calls the numpy method `self._data_arr_mask.any()` (`combiner.py:517`) — fixed in #988, along with the further numpy-isms the fix exposed in the same functions (`.sum(axis=0)` mask counts, `len()` on arrays, list passed to `xp.reshape`, int→float promotion before `sqrt`)
- [x] `combine()` memory sizing uses `.nbytes` on data/uncertainty/mask/flags — fixed in #988 via `_array_size_in_bytes` (element count × `finfo`/`iinfo` bit width)
- [x] Builtin `bool` passed where `xp.bool` is required, in the mask paths of `Combiner` and `_CCDDataWrapperForArrayAPI` — fixed in #988 (wrapper mask setters plus the test sites building masks in the strict namespace)
- [x] Combiner: median combine and the uncertainty calculation still drop to plain numpy (see comments around `combiner.py` median/uncertainty paths; #351)

### 2. Test-suite hygiene
- [x] #969 `test_subtract_overscan` calls numpy-only `.copy()` — fixed in #990 (the two `median=True` parametrizations went green with #989)
- [x] #970 raw `np.zeros_like` in `test_ccdproc_logging.py` / `test_rebin.py` — fixed in #990 (`test_rebin_ccddata[True-True]` stays red on astropy's numpy-forcing `CCDData.mask` setter, #983's blocker)
- [x] Numpy `.sum()` / `.mean()` method calls in `test_combiner.py` test bodies — fixed in #990, along with ellipsis indexing, int/float promotion and `xp.asarray(list)`→`xp.stack` cases
- [x] `test__overscan_schange` calls `xp.allclose` — fixed in #990
- [x] Six `test_combiner.py` tests build scaled inputs with astropy `CCDData.multiply(...)` on strict data (`test_combiner_with_scaling`, `test_combiner_with_scaling_uncertainty` ×3, `test_combiner_result_dtype`, `test_combine_overwrite_output`); astropy's `np.result_type` in `_prepare_then_do_arithmetic` cannot take a strict dtype. Build the inputs in the namespace (`CCDData(ccd.data * 3, unit=...)`) instead — this is the NDData-arithmetic limitation already listed in #940 and is what `_CCDDataWrapperForArrayAPI` exists for. 6 strict failures — **PR #998** (the four jax `backend_xfail` markers for the same `.multiply()` calls are dropped, they XPASS; `test_combiner_with_scaling_uncertainty[average|sum]` clear, and the other four stop on the next cause: `test_combine_overwrite_output` on `combine(output_file=)` (section 3), `test_combiner_result_dtype` on `combine(dtype=int)` (section 1), `test_combiner_with_scaling` and `test_combiner_with_scaling_uncertainty[median_combine]` on #929)
- [x] Five `combine()` tests feed FITS filenames or numpy `CCDData` without `array_package=` and then compare the (numpy) result against strict arrays with `xp.all(xpx.isclose(...))` (`test_combine_average_fitsimages`, `test_combine_numpyndarray`, `test_combine_average_ccddata`, `test_combine_limitedmem_fitsimages`, `test_combine_limitedmem_scale_fitsimages`). Pass `array_package=xp` (or compare in numpy). 5 strict failures — **PR #998** (`array_package=xp` plus `Combiner(ccd_list, xp=xp)` for the reference; `test_combine_average_ccddata` converts its `CCDData` list by hand because `array_package` is ignored for `CCDData` input)
- [x] `test_combiner_image_file_collection_input` / `test_combine_image_file_collection_input` write a `device1` strict `CCDData` with `ccd.write()`; `astropy.io.fits` cannot coerce it. Write a numpy copy. 2 strict failures — **PR #998** (`numpy_copy`/`numpy_ccddata` helpers in `pytest_fixtures.py`; also clears `test_writeable_after_combine[average|sum_combine]`, while `[median_combine]` stays on #929)
- [x] `test_subtract_overscan_model` divides an integer strict array by a float (`test_ccdproc.py:304`); cast `np_mgrid` to float first. 2 strict failures — expect these to then land on #933 (`astropy.modeling`) — **PR #998** (they do, via the `len(oscan)` numpy-ism listed in section 1; now `backend_xfail("array-api-strict")` citing #933)
- [x] `_make_mean_scaler` (`test_combiner.py:52`) uses `.mean()` and `np.ma.average` — rewrite in `xp`. 2 strict failures (`test_3d_combiner_with_scaling`, and `test_combine_ccd_with_uncertainty_and_mask_from_fits[function]` once #995 lands) — **PR #998** (both pass)
- [x] #946 device-mismatch noise from inline array creation in strict-job test bodies — after the section-1 `device=` fixes above the remaining device errors are all astropy setters/`Quantity`, not test bodies; re-check once those land and close if nothing is left — re-checked after #998: none of the 15 remaining strict failures is a test-body device error; **#946 closed 2026-08-25**
- [x] `test_generator_ccds_without_unit` fails on the strict job because of an **astropy bug**, not test ordering per se: the `_arithmetic` decorator in `astropy/nddata/ccddata.py` sets the module-global `_config_ccd_requires_unit = False` and restores it with no `try`/`finally`, so *any* exception inside `CCDData.add/subtract/multiply/divide` leaves the unit requirement disabled process-wide. `test_unit_mismatch_behaves_as_expected` raises inside `.subtract()` by design, so the flag leaks on every backend — reproduced on plain numpy with just those two tests. The full numpy run stays green only because intervening tests (e.g. `test_gain`) complete a successful arithmetic pass, which resets the flag; on strict those tests fail too, so nothing repairs it. Filed upstream as astropy/astropy#20268; listed in #940 — **PR #998**: removing the `.multiply()` calls from `test_combiner.py` also removed the incidental successful arithmetic that had been resetting the flag, so the leak surfaced on plain numpy too; the test now resets `astropy.nddata.ccddata._config_ccd_requires_unit = True` itself and passes on every backend regardless of order, and the marker is gone
- [x] #943 jax `ccds()` platform-dependent `ValueError` behaviour (Linux CI only) — same astropy leak as above; the xfail reason at `test_image_collection.py:438` should cite astropy/astropy#20268 rather than "platform-dependent" — **PR #998** drops the marker altogether (see the previous bullet); **#943 closed 2026-08-25**
- [x] Correct the xfail reason added in #958 that cites astroscrappy (astroscrappy is monkeypatched out in that test) — already done on `main` in `d33d7cf`

### 3. Design issues / blocked on upstream (tracked in #940)
- [x] #929 `Combiner.sigma_clipping` and the default `sigma_func` densify via `astropy.stats` — **done: PR #1000 (`sigma_func`) and PR #1001 (`sigma_clipping`), both merged 2026-08-30**, numpy unchanged, every other namespace xp-native; leaves strict at 2 (#936, #983). History: 8 strict failures (5 via `median_absolute_deviation` at `core.py:1291`, 3 via `sigma_clip` at `combiner.py:491`); the largest remaining bucket once section 1 is done. After #998 it is **11**: the three `median_combine` callers `test_combiner_with_scaling`, `test_combiner_with_scaling_uncertainty[median_combine]` and `test_writeable_after_combine[median_combine]` reach the default `sigma_func` once their first cause is cleared. With #999 merged the split is **8 via `sigma_func`** (`core.py:1409`; `test_sigma_func_for_ccddata`, `test_combiner_dtype`, `test_combiner_median`, `test_combiner_with_scaling`, `test_combine_result_uncertainty_and_mask[median_combine-*]`, `test_writeable_after_combine[median_combine]`, `test_combiner_with_scaling_uncertainty[median_combine]`) and **3 via `sigma_clipping`** (`combiner.py:480`; `test_combiner_sigmaclip_{high,low,single_pix}`). The `sigma_func` half is cheap on the ccdproc side (MAD = `median(|x - median(x)|)` on top of the `_nanfuncs.py` median); the `sigma_clipping` half is harder because the method forwards `**kwd` and string/callable `func`/`dev_func` to astropy — options are an xp-native clip or `backend_xfail` markers
- [x] `combine(output_file=...)` writes the namespace result with `ccd.write(output_file, ...)` (`combiner.py:1273`); `astropy.io.fits` needs NumPy, so on strict it fails with `AttributeError: 'Array' object has no attribute 'astype'` (`ccddata.py:399`). 1 strict failure (`test_combine_overwrite_output`, reachable since #998). Same class as the CPU-only policy (#935): the writer needs a NumPy copy — **PR #999** (new `core._to_numpy` host-copy helper, device-generic; the writer gets a NumPy `CCDData` and the namespace result is returned unchanged. This is the first library-side host-copy helper, so #935 can build on it)
- [ ] #936 `Quantity`/units handling with non-numpy arrays — 1 strict failure (`test_unit_mismatch_behaves_as_expected`)
- [ ] #984 `cosmicray_median` does not exclude masked pixels from the median filter itself, so bright masked regions pull the local median (follow-up to #932/#979; runs in numpy, so it also belongs with the CPU-only policy in #935)
- [ ] #983 consider `marray` as a uniform masked-array representation across backends; bounded by astropy's numpy-only `NDData.mask` setter
- [ ] #933 `subtract_overscan` model fitting goes through numpy-only `astropy.modeling`
- [ ] #930 `wcs_project` mixes numpy output from `reproject` with xp operations
- [ ] #940 upstream blockers: `astropy.stats.sigma_clip`, `median_absolute_deviation`, `block_reduce`/`block_replicate`, NDData arithmetic (`np.result_type` in `_prepare_then_do_arithmetic`, `np.abs` in `_propagate_multiply_divide`; the `_CCDDataWrapperForArrayAPI` workaround), the `_config_ccd_requires_unit` leak (astropy/astropy#20268), `astropy.modeling`, `Quantity`, `reproject` (CPU-only), astroscrappy (permanent → document under the CPU-only policy, #935), `__array__(copy=)` warning (#934)

### 4. Tooling / infrastructure
- [ ] #945 escape-baseline keys are fragile across Python versions (PEP 709)
- [ ] #947 refactor `_escape_triage.py`; register pytest hooks explicitly
- [ ] #938 sparse backend as a densification detector (sub-items in that issue)
- [ ] Escape logger only patches `np.asarray`/`np.asanyarray`/`np.ma.asanyarray`; it misses `np.array()` and ufunc-driven coercion, so it is a sampler rather than a complete detector
- [ ] Decide between a parametrized `xp` fixture (SciPy-style, all backends in one run) and the current one-backend-per-`CCDPROC_ARRAY_LIBRARY` approach
- [x] Promote the `array-api-strict` job from soft to required — **done in #1005 (merged 2026-08-30)**: the last ccdproc-owned strict failure (`flat_correct` using the numpy-only `mask.any()` method) is fixed, #936/#983 carry `backend_xfail` markers, `py313-strict` moved into the regular `ci-tests` matrix, and the expected-failures job + `strict_status.yml` are deleted. (Marking the job *required* in branch protection is a separate admin step; `main` currently has no required status checks configured at all.) History: sections 1–2 have landed (2026-08-26); the residue is exactly 13 upstream-blocked tests (11 × #929, #936, #983), of which 11 are fixed by #1000/#1001; the residue after those is #936 and #983, which need `backend_xfail` markers citing the astropy items so the job can be green
- [x] #1004 hoist `None`/tuple axis handling into `_nanfuncs._setup` — **done in #1006 (merged 2026-08-31)**: `_setup` now takes an int, `None` or a tuple/list of axes and returns a `restore` callable mapping full-shape arrays back to the caller's layout, so every `_nanfuncs` reduction gained tuple/`None` axis support for free, `core._mad_fallback` shrank to the median tier plus a delegation to `nanmad`, and `Combiner.sigma_clipping` off numpy accepts the same axis forms as `astropy.stats.sigma_clip` (int, `None`, tuple or list). Review caught two edge cases before merge: `np.bool_` axes were silently treated as axis 1 on numpy 2.0 (now `TypeError`), and the tuple merge reshaped with `-1`, which raises when a kept axis has size 0 (`numpy.median` itself has this bug; the fallbacks now return an empty result like numpy's other reductions). Behaviour note: a bool or otherwise non-integer `axis` now raises `TypeError` rather than `NotImplementedError`
- [ ] Decide whether CuPy and/or torch get a CI job

### 5. Close out
- [x] Re-triage the `backend_xfail` markers after the fixes above and drop the stale ones — **done 2026-08-30 (#1005 dropped the stale `test_flat_correct_masked_flat_with_immutable_array` marker; the 30 remaining markers are all genuine numpy-only boundaries)**. History: the five known-stale ones below are narrowed in **PR #996** (merged; strict XPASS 5 → 0), and #998 drops five jax markers (four `.multiply()` ones in `test_combiner.py`, the `test_generator_ccds_without_unit` one). Known stale before #996: `test_cosmicray.py::test_cosmicray_median_masked_column[masked_array-1.0]` and `[masked_array-None]` XPASS on the strict job, and their reason (a `scipy.ndimage.median_filter` device failure) no longer applies after #979; the `test_cosmicray_gain_correct_uncertainty_namespace` marker (`test_cosmicray.py:154`) XPASSes for the three `gain_apply=False` cases and should be narrowed to the `gain_apply=True` parametrizations
- [ ] Re-run the suite against CuPy and close the parts of #909–#913 that are subsumed by the issues above
- [ ] Confirm `docs/array_api.rst` and `CHANGES.rst` reflect the final set of supported backends

Start by reading the migration checklist in this issue, then inspect the named entry points in ccdproc/core.py, ccdproc/combiner.py, and _ccddata_wrapper_for_array_api.py. Run the py313-strict and backend test jobs to identify an unclaimed remaining item. Done means the selected migration work is addressed without new failures and the strict array-API checks remain clean.

Written by the indexing model from the issue text.