astropy / astropy/ccdproc

combine(): uncertainty and mask always allocated at full output size regardless of mem_limit, and each tile re-reads whole input files

Open
#1,012 0 comments 0 reactions 0 assignees View on GitHub
Dominant language
Python
Stars
93
Forks
92
Avg merge
14h 44m
Merged PRs (30d)
30

Description

This issue was written and opened by Claude Code at the direction of Matt Craig (@mwcraig).

Environment: ccdproc 2.4.3, astropy 7.1.0, numpy 2.3.5.

`combine()` (`ccdproc/combiner.py`) already tiles the pixel array by rows/columns to respect
`mem_limit`, and already honors an explicit `dtype=` for the combined **data** accumulator
(line 897: `ccd.data = ccd.data.astype(dtype)`; defaults to float64 at line 891 when
`dtype=None`). Two things keep the peak from ever coming down as far as `mem_limit` implies:

1. Regardless of `dtype`/`mem_limit`, if the sample image lacks uncertainty/mask, `combine()`
unconditionally attaches full-size arrays to the output *before* tiling starts:
`ccd.uncertainty = StdDevUncertainty(np.zeros_like(ccd.data))` (line 902) and
`ccd.mask = np.zeros_like(ccd.data, dtype=bool)` (line 907). There is no kwarg to skip
either. At 4096x4096, `dtype=float64` costs 17 bytes/pixel (285 MB) for these three arrays
alone; `dtype=float32` costs 9 bytes/pixel (151 MB) — still paid no matter how small
`mem_limit` is.
2. Even though chunk count is derived from `mem_limit` (line 922), each tile still does
`imgccd = CCDData.read(image, **ccdkwargs)` for *every* input file (lines 952, 987) —
`no_chunks * n_images` full-file reads — instead of reading only the band via `hdu.section`.
Chunking shrinks the retained accumulator, not the transient per-file read.

Confirmed on the versions above:
```python
>>> out = ccdproc.combine(imgs, method='average', mem_limit=1e4, dtype=np.float32)
>>> out.data.dtype, out.mask.dtype, out.uncertainty.array.dtype
(dtype('float32'), dtype('bool'), dtype('float32'))
```
(mask/uncertainty attached even though none of the 5 inputs had either, `mem_limit` far below
one image's size.)

Reproducer sketch:
```python
import numpy as np, tracemalloc
from astropy.nddata import CCDData
import astropy.units as u
import ccdproc

imgs = [CCDData(np.random.default_rng(i).random((4096, 4096)).astype(np.float32), unit=u.adu)
for i in range(20)]
tracemalloc.start()
combined = ccdproc.combine(imgs, method='average', mem_limit=5e7, dtype=np.float32,
sigma_clip=True, sigma_clip_func='median',
sigma_clip_dev_func='mad_std')
print(tracemalloc.get_traced_memory()[1] / 1e6, "MB peak")
```

Measurements (4096x4096, 20 float32 flats, tracemalloc peak, excluding ~235 MB process
baseline), comparing `combine()` against a ~40-line hand-written band loop over
`ccdproc.Combiner` that reads bands with `hdu.section`:

| Variant | Peak | Time |
|---------------------------------------------------------------------------------|--------:|------:|
| `combine(mem_limit=1e8)`, float64 accumulator, callable clip funcs | 403 MB | 38 s |
| `combine(mem_limit=1e8)`, `dtype=float32`, string clip funcs, mask/uncertainty dropped right after return | 230 MB | 21 s |
| Hand-written band loop, float32, `mem_limit=1e8`, average | 212 MB | 13 s |
| Hand-written band loop, median | 258 MB | 38 s |
| Hand-written band loop, average, `mem_limit=5e7` | 140 MB | 13 s |

The band loop's output matched `combine()` exactly (max abs diff 0.0); its peak scales down
with `mem_limit` because it never attaches mask/uncertainty and reads only the needed band per
file.

Suggested improvements: make uncertainty/mask opt-in on the output (only allocate when
`combine_uncertainty_function` is given or an input already has one; make mask opt-in
separately); read each tile via `hdu.section` instead of `CCDData.read()`-ing the whole file
per chunk.

Contributor guide

Open the contributing guide

Research direction

Read ccdproc/combiner.py around lines 891-987, then run the provided combine() reproducer with tracemalloc to establish current allocations and reads. Trace how mem_limit, uncertainty, mask, and each tile's CCDData.read() interact. Done means preserving the reported output while avoiding unconditional full-size arrays and reducing per-tile input reads to the requested bands.

Written by the indexing model from the issue text.

Assessment

Tech stack
numpy, python
Domain
data, performance
Issue type
Bug
Difficulty
4/5
Estimated time
3-5 days
Activity status
Active
Clarity
Mostly clear
Newbie friendliness
45/100

Get new issues in your inbox

A short digest of beginner-friendly GitHub issues.