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

[mwcraig]

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:

>>> 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:

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.

[mwcraig]

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

More data on point 2 (each tile calls CCDData.read on every input file). The cost is in time as well as memory, and it grows as the square of the number of images: the number of chunks is proportional to n_images (no_chunks = memory_factor * size_of_an_img * n_images / mem_limit), and every chunk opens all n_images files.

Environment: ccdproc 2.4.3, astropy 7.1.0, numpy 2.3.5; the loop is unchanged on main (0c21068). Inputs: 4096x4096 float32 single-HDU FITS files given as paths, mem_limit=5e7, dtype='float32', no clipping, macOS, warm page cache.

	10 images	40 images
chunks (average / median)	79 / 118	315 / 472
CCDData.read calls (average / median)	790 / 1,180	12,600 / 18,880
combine, average	16 s	223 s
combine, median (trivial combine_uncertainty_function)	23 s	363 s
files opened once, bands read with hdu.section, Combiner per band: average	1.0 s	3.9 s
same, median	1.7 s	4.8 s

Four times the images takes 14 to 16 times as long. The read itself is lazy (the data comes back as a view of a memmap, and only the tile is copied), so nearly all of that time is the fixed cost of CCDData.read: about 16 ms per call here for parsing the header and building the WCS, paid no_chunks * n_images times. The output of the open-once loop is bit-identical to combine's in all four cases, and its peak RSS is no higher.

A smaller change than reading with hdu.section also removes the N² term: read each file once before the tile loop and reuse the CCDData objects. That took 19 s (average) and 29 s (median) for 40 images. But the data of each is then a live memmap, and the pages touched stay mapped for the whole combine: peak RSS was 3.1 GB for 40 images against about 0.5 GB now, which defeats mem_limit for anyone whose limit is on RSS (a cgroup, say). Keeping the HDUs open and reading each band with .section (or memmap=False) does not have that problem.