Repository navigation
Description
Overview
External contributor @HRSAstro (pyuvimage, Discussion #13, comment item 1) asks for an optional oversample=q on sparse_terms_from_chunks. It would accumulate two fine grids in the same streaming pass: the precision operator at sub-pixel lags, K(Δ) = Σ w cos(2π u·Δ), and the noise-weighted dirty image at sub-pixel positions, D(x) = Re Σ c exp(2πi u·x). Those two grids are enough to form the sparse-likelihood terms of uv-analytic components (points, small Gaussians), so pyuvimage could drop its own pass. This issue covers item B1 of the 2026-10-07 technical review: accumulate and persist the grids. The column-terms helper (B2) is out of scope.
Dependency satisfied: Depends-on: sparse_noise_map_pooling_option is merged as PyAutoArray#619 (merge 2c5cb69). The fine grids inherit its equal-sigma check and pool_noise_map=True pooling.
Coordination note (human-authorised, 2026-10-07): worktree_check_conflict reports claims from community-pages (README only, PyAutoArray#618 / PyAutoGalaxy#649) and scribbler-wave2-radial-panels-regrid (PyAutoGalaxy#641 already MERGED, awaiting close-out). This task proceeds as coordination on disjoint files, the same arrangement sparse-noise-map-pooling used. It does not touch README.md or any radial-panel/regrid plotting code.
Plan
- Lift the type-1 NUFFT core of the existing precision-operator builder into a reusable helper. The existing operator must stay bit-for-bit unchanged, pinned by the existing parity tests.
- Add an optional even
oversample=q(NUFFT path only) to the streamed sparse-terms builder. With it set, the builder accumulates a fine precision-operator grid (full native shape plus a 25 % lag pad, q times finer) and a fine dirty-image grid (twice the field, q times finer). - Accumulate both grids in place across chunks, so memory stays at the two grids plus the NUFFT work grid. Log a memory estimate when
oversampleis set. - Record the grids and
oversample/oversample_padonSparseTerms. Adding two terms refuses mismatched oversampling, or exactly one side carrying fine grids. - Pass the option through
Interferometer.from_streamandapply_sparse_operator_from_chunks. - PyAutoGalaxy: persist the fine grids in the dataset FITS behind an opt-in writer keyword (off by default, about 840 MB at 400 px, q=8). Append the new scalars, and load older files as
None. - Tests: fine K and D agree with today's operator and dirty image at native lags/centres; chunked equals one-shot; the phase centre shifts D but not K;
__add__refusals; odd q and DFT raise; FITS round trip including legacy files.
Tier: judge — merge mode: human /prm
Detailed implementation plan
Affected Repositories
- PyAutoArray (primary)
- PyAutoGalaxy (FITS persistence; depends on the PyAutoArray PR)
Branch Survey
| Repository | Current Branch | Dirty? |
|---|---|---|
| ./array/PyAutoArray | main | clean |
| ./galaxy/PyAutoGalaxy | main | clean |
Suggested branch: feature/sparse-terms-oversampled-fine-grids
Implementation Steps
autoarray/inversion/inversion/interferometer/inversion_interferometer_util.py:- New
_type1_real_grid_from(uv_wavelengths, values, n_modes, pixel_scale_radians, *, eps, chunk_size, centred, zero_nyquist, out=None). It returnsRe nufft2d1(-x, +y, values, n_modes, eps, +1), withx = 2π u Δ,y = 2π v Δ, chunked over visibilities, ifftshifted to wraparound order unlesscentred, Nyquist zeroed ifzero_nyquist. Withoutit adds in place (quadrant-wise for the shift), with no per-chunk full-grid allocation.nufft_precision_operator_via_nufft_fromcalls it withcentred=False, zero_nyquist=True: same operations, same order, bit-for-bit. sparse_terms_from_chunks(..., oversample: Optional[int] = None, oversample_pad: float = 0.25). q must be a positive even integer and the transformer aTransformerNUFFT(exc.InversionException/ValueErrorotherwise, checked before any work).- Fine K:
(2(Ny+p)q, 2(Nx+p)q),(Ny, Nx)= full native shape,p = ceil(oversample_pad * max(Ny, Nx)), valuesw = 1/σ_r², spacingΔ/q, wraparound order like W~,zero_nyquist=False. - Fine D:
(2Ny q, 2Nx q)centred (origin at[Ny q, Nx q], native orientation), from the phase-shiftedc = d_r/σ_r² + i d_i/σ_i². Built as the core onconj(c), which flips both axes to the native orientation. - Both preallocated once and
+=per chunk. OneSparseTermsbuilt at the end viadataclasses.replace. One memory-estimate log line.
- Fine K:
SparseTerms: appendprecision_operator_fine,dirty_image_fine,oversample,oversample_pad(Optional,None).__add__addsoversample,oversample_padto the provenance loop (numeric compare), refuses when exactly one side carries fine grids, checks fine shapes, and sums the fine grids.
- New
autoarray/dataset/interferometer/dataset.py: passoversample/oversample_padthroughInterferometer.from_streamandapply_sparse_operator_from_chunks.- PyAutoGalaxy
autogalaxy/interferometer/model/analysis.py:SPARSE_TERMS_SCALARS_ORDERgetsoversample,oversample_padappended (append-only, NaN = unrecorded), plus headerOVERSAMP. HDUsPRECISION_OPERATOR_FINE/DIRTY_IMAGE_FINEare written only when present andinclude_fine_grids=True(default False). Loaderautogalaxy/aggregator/interferometer/interferometer.pyreads them via_has_hdu; pre-change files giveNone. - Docs: the memory/cost note (400 px, q=8: K 512 MB + D 328 MB plus ~4 GB nufftax work grid; use large chunks, ≳1e6 visibilities, with oversample).
Tests (B1 list)
(i) fine K at lags (iq, jq) == nufft_precision_operator[i, j] within the extent (mixed rtol/atol); (ii) fine D at native pixel centres == dirty_image_native; (iii) chunked == one-shot for both grids; (iv) phase_centre shifts D, not K; (v) __add__ refuses mismatched q and fine/no-fine; odd q / DFT raise; (vi) PyAutoGalaxy FITS round trip, including a legacy file loading as None. W~ bit-for-bit pinned by the existing NUFFT-vs-brute parity tests.
Key Files
autoarray/inversion/inversion/interferometer/inversion_interferometer_util.pyautoarray/dataset/interferometer/dataset.pytest_autoarray/inversion/inversion/interferometer/(sparse terms tests)autogalaxy/interferometer/model/analysis.pyautogalaxy/aggregator/interferometer/interferometer.pytest_autogalaxy/interferometer/model/test_analysis_interferometer.py
Deferred follow-up (not in scope): B2 point_column_terms_from(...) helper (FFT Gaussian smoothing + quintic spline).
Original Prompt
Click to expand starting prompt
SparseTerms oversample=q: accumulate fine precision-operator and dirty-image grids for analytic uv-plane components (streaming P6)
Type: feature
Target: PyAutoArray
Repos:
- PyAutoArray
- PyAutoGalaxy
Themes: - interferometer
- sparse-operator
- community
Autonomy: supervised
Priority: normal
Status: draft
Filed: 2026-10-07
Difficulty: medium
Consequence: judge
Depends-on: sparse_noise_map_pooling_option (merge first; same repo claim)
Witness: withsparse_terms_from_chunks(..., oversample=8)(NUFFT path), (i)precision_operator_finesampled at lags (i q, j q) equalsnufft_precision_operator[i, j]within the masked extent (mixed rtol/atol periiu:529-531), and (ii)dirty_image_finesampled at native pixel centres equalsdirty_image_native; todayoversampledoes not exist.
Review-minutes: 20
Unattended: ready
Source: GitHub Discussion https://github.com/orgs/PyAutoLabs/discussions/13 (external contributor @HRSAstro, pyuvimage), comment Streaming visibilities for memory efficiency .github#13 (comment) item 1; technical review 2026-10-07 (Item B1, accept-with-change).
Request (verbatim)
1. Analytic components on the streamed path
sparse_profile_terms_fromcovers light profiles evaluated on the image grid, but a component that is analytic in the uv-plane (a point source, or a small Gaussian for a marginally resolved source) needs its terms at sub-pixel positions. All three terms of such a column are values of two image-plane functions that can be accumulated in the same streaming pass:
- K(Δ) = Σ w cos(2π u·Δ), i.e. the precision operator at an arbitrary lag
- D(x) = Re Σ (d_re/σ_re² + i d_im/σ_im²) exp(2πi u·x), i.e. the dirty image at an arbitrary position
For a unit point at p, the cross-term with the mesh is Mᵀ k_p with k_p[i] = K(x_i − p), the point-point term is K(p − q) and the data term is D(p). [...]
In pyuvimage we accumulate both on grids 8× finer than the image, with one extra type-1 NUFFT per batch for each, and read them with a quintic spline. [...]
- The precision-operator grid needs lags beyond one field width (we pad by a quarter field), and the dirty-image grid needs twice the field. Otherwise Gaussian-widened columns near the edge pick up the wrap-around.
[...] An optionaloversample=qinsparse_terms_from_chunksthat accumulates these two fine grids would be enough for us to drop our own pass. A helper that returns the column terms for given positions and widths would make it general. Our implementation is inpyuvimage/streamed_points.pyif it's a useful reference.
What
iiu = autoarray/inversion/inversion/interferometer/inversion_interferometer_util.py (origin/main lines). nufft_precision_operator_via_nufft_from (iiu:454-624) already builds W~ as one type-1 NUFFT onto a (2Ny, 2Nx) grid over the masked extent (iiu:2220), spacing hard-wired to the image pixel (iiu:580), Nyquist zeroed (iiu:621-622). A fine grid needs a small refactor, not a flag. The maths in the request is right; K uses w only, so it inherits the equal-sigma assumption (hence the dependency).
Scope: B1 only. Deferred follow-up (not in scope): B2 point_column_terms_from(terms, positions, widths, mapping_matrix=None, *, spline_order=5) helper (FFT-space Gaussian smoothing, quintic spline; not JAX-traceable) — file separately if we want point components in PyAutoLens's own sparse fits.
Plan
- Lift the W~ core loop into
_type1_real_grid_from(uv, values, n_modes, pixel_scale_radians, eps, chunk_size, centred, zero_nyquist); W~ calls it unchanged (existing tests pin it). - API:
sparse_terms_from_chunks(..., oversample: Optional[int] = None, oversample_pad: float = 0.25), pass-through onfrom_stream/apply_sparse_operator_from_chunks. Require even q (native centres land on fine grid points) and the NUFFT transformer, else raise. - Grids:
precision_operator_fineshape(2(Nk+pad)q)^2, pad = ceil(0.25*max(Nk)), Nk = full native shape (points may sit anywhere in the mask), wraparound order like W~,zero_nyquist=False(spline reads across lag N).dirty_image_fineshape(2Ny q, 2Nx q)centred (origin at [Ny q, Nx q]) from the phase-shiftedc = d_re/s_re^2 + i d_im/s_im^2exactly asiiu:2233-2236. SparseTerms(iiu:1880, frozen dataclass): appendprecision_operator_fine,dirty_image_fine,oversample,oversample_pad(all Optional, None).__add__(iiu:1951): addoversample/oversample_padto the provenance loop (iiu:1977-1984) and refuse when exactly one side carries fine grids (pyuvimage silently drops them; refusing is our pattern).- Accumulation: preallocate the two fine arrays and
+=in place across chunks; build oneSparseTermsat the end (do NOT+a fresh SparseTerms per chunk,iiu:2272, for 100s of MB arrays). - Memory/cost (document; consider a memory-estimate log line): 400x400, q=8, 25 % pad → K fine 8000^2 f64 = 512 MB, D fine 6400^2 = 328 MB, plus nufftax 2x-upsampled work grid 16000^2 c128 = 4.1 GB (K) / 2.6 GB (D): peak ~5-6 GB; q=4 ~1.5 GB; contributor's ~112 px at q=8: 40 MB + 26 MB, <0.5 GB work. Two extra type-1 NUFFTs per chunk with FFT cost independent of chunk length (seconds per chunk at 400/q8): document "use large chunks (>= ~1e6 vis) with oversample". JAX: each (K, n_modes, eps) compiles once; ragged last chunk adds signatures — pad to a fixed bucket with zero weights only if measured to matter.
- PyAutoGalaxy persistence (
autogalaxy/interferometer/model/analysis.py:40-66, writer ~:69, loaderautogalaxy/aggregator/interferometer/interferometer.py:62-111): HDUsPRECISION_OPERATOR_FINE,DIRTY_IMAGE_FINEonly when present; appendoversample,oversample_padtoSPARSE_TERMS_SCALARS_ORDER(append-only, NaN = unrecorded) and headerOVERSAMP; loader uses_has_hdu. Decision: fine grids are NOT written intodataset.fitsby default (~840 MB per search output at 400/q8) — config flag to opt in; consideraa.SparseTerms.output_to_fits/from_fitsas the cleaner home.
Tests: (i)/(ii) the witness; (iii) chunked == one-shot for both grids; (iv) phase_centre shifts D, not K; (v) __add__ refuses mismatched q and fine/no-fine; (vi) PyAutoGalaxy FITS round trip incl. an old file without the new HDUs/scalars loading as None; odd q / DFT transformer raise.
Risks: memory and per-chunk FFT cost at large grids; extra JAX compiles per chunk shape; FITS bloat if the default flips.
Activity
Jammy2211 commentedon Oct 7, 2026
Library PRs created
Dependency:
sparse_noise_map_pooling_optionis MERGED as PyAutoArray#619 (2c5cb69). The fine grids inherit its equal-sigma check andpool_noise_mappooling.Coordination (human-authorised): proceeds alongside
community-pages(README only) andscribbler-wave2-radial-panels-regrid(PyAutoGalaxy#641 merged, awaiting close-out) on disjoint files. No README or radial-panel/regrid code was touched.Validation: the witness ran red first. Full
test_autoarray2014 passed / 4 xfailed; fulltest_autogalaxy1357 passed.W~is pinned bit for bit across the refactor. Smokeautolens_workspace/.../datacube/modeling_array_free.pypasses. Heart YELLOW (manifest drift; release validation stale), acknowledged.Measured (CPU, 100² native, 2×1e5-vis chunks, warm): None 2.1 s / 2016 MB peak RSS; q=4 5.2 s / 2159 MB; q=8 5.0 s / 2533 MB.
API: additive only (
oversample,oversample_pad, four optionalSparseTermsfields, opt-ininclude_fine_gridson the Galaxy writer). No workspace migration is needed. A future workspace demo could showoversamplefor point components, alongside the deferred B2 helper.Tier: judge, so merge goes through a human
/prm.Jammy2211 commentedon Oct 7, 2026
Shipped
Both PRs merged (Array first, as planned):
What landed:
sparse_terms_from_chunks(..., oversample=q, oversample_pad=0.25)(evenq, NUFFT transformer) accumulates two fine grids in the same streaming pass —precision_operator_fine(K at arbitrary lags, padded by a quarter field) anddirty_image_fine(D at arbitrary positions, twice the field) — carried onSparseTerms, with__add__refusing mismatchedoversampleor fine/no-fine operands. PyAutoGalaxy persists the fine grids as optional HDUs (opt-in; not written intodataset.fitsby default) and older files load with them asNone.Measured cost (100×100 native, 2 chunks × 1e5 visibilities):
Release: merged, not yet released — pending the next PyAutoArray / PyAutoGalaxy release.
Deferred: the B2 helper (
point_column_terms_from— column terms for given positions/widths via FFT-space Gaussian smoothing + quintic spline) is not in this change; it will be filed separately if wanted.Requested in Discussion https://github.com/orgs/PyAutoLabs/discussions/13 (comment PyAutoLabs/.github#13 (comment), item 1) — thanks @HRSAstro.