combine(): uncertainty and mask always allocated at full output size regardless of mem_limit, and each tile re-reads whole input files
- 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
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