Skip to content

perf: build nufft_precision_operator_from as a type-1 NUFFT - #541

Merged
Jammy2211 merged 2 commits into
mainfrom
feature/interferometer-preload-nufft-type1
Sep 8, 2026
Merged

Jammy2211 merged 2 commits into
mainfrom
feature/interferometer-preload-nufft-type1

Conversation

@Jammy2211

@Jammy2211 Jammy2211 commented Sep 8, 2026

Copy link
Copy Markdown
Collaborator

Stacked on #540; merge #540 first, GitHub retargets this PR to main.

Summary

nufft_precision_operator_from builds the interferometer sparse-operator preload — the compact (2Ny, 2Nx) array that depends only on the (dy, dx) offset between image pixels. Until now it was built by brute force: an O(N_pix · K) cosine accumulation over every image pixel × visibility pair, which is minutes to hours at real dataset sizes and is why the workspace ships a separate "prepare the preload overnight and cache it to .npy" example.

The same array is the real part of a type-1 (adjoint) NUFFT of the inverse-variance weights w = 1/σ² evaluated on the doubled-extent grid:

P = Re nufft2d1(-x, y, w, (2Nx, 2Ny), eps, +1)   with x = 2π·u·Δ, y = 2π·v·Δ
    → ifftshift → zero the Nyquist row Ny and column Nx → contiguous float64

which costs O(N_vis · nspread² + M log M) for M = 4·Ny·Nx. The sign convention (-x, +y, isign +1) is the one of eight candidate mappings that reproduces the brute force; the derivation and the search are recorded in the function's docstring.

method="nufft" is the new default. The numpy and JAX brute-force builders are untouched and stay as the reference implementations the NUFFT is pinned against.

The array is bit-for-bit the same object (to the pins below), so existing workspace .npy preload caches remain valid — nothing needs regenerating.

Measured (autolens_profiling#229, phase 3)

dataset brute force type-1 NUFFT speed-up
alma 2101 s 7.3 s 289× wall (111× CPU-seconds)
alma_high 22 s

CPU-seconds caveat: the 289× is wall-clock. The brute-force builder is single-threaded, the NUFFT spreader is multi-threaded, so in CPU-seconds the honest figure is 111×; the wall-clock number is what a user experiences on a workstation.

Pins

  • Mixed-tolerance parity against _via_np_from, seeded 16×16 mask at K=300: max deviation 4.4e-14 of peak (rtol=1e-10, atol=1e-10·P[0,0] — mixed because the off-peak entries are many orders of magnitude below the peak, so a pure rtol would pin round-off).
  • Chunked (chunk_size=64, K > 64) == one-shot: 5.3e-14.
  • Wrong-sign control (uv[:, 0] negated) fails the pin at 0.187 of peak — i.e. the pin has real discriminating power.
  • Padding row Ny / column Nx exactly zero, neighbours non-zero; P[i, j] == P[-i, -j].
  • End-to-end: the existing test_interferometer.py sparse-vs-mapping comparisons run through the new default unchanged.

eps and chunking

eps = 1e-12 saturates fp64 at every instrument profiled — the measured max|Δ| is already at the round-off floor, and it is the only value that holds the array pin. chunk_size is not an optimisation, it is a memory ceiling: the spreader's gather buffer is one-shot ~15 GB at K = 5e6, so chunking is mandatory at real scale. When not given it is taken from the transformer's own chunk_size if it is a TransformerNUFFT.

Fallback semantics

Two loud, logged fallbacks to method="numpy" — never silent:

  • PYAUTO_DISABLE_JAX=1 (the test-mode kill switch already honoured in dataset.py), or
  • nufftax is not importable.

Both emit a logger.warning naming the O(N_pix · K) cost being paid. Any other unusable method raises.

API Changes

Additive only — every existing keyword is kept with its existing meaning, so no caller breaks.

  • nufft_precision_operator_from gains method="nufft" (new default), eps=1e-12 and chunk_size=None. use_jax is now only honoured when a brute-force method is selected (method="numpy" + use_jax=True = the JAX brute force); it is ignored under the default method="nufft", which already runs on JAX. method="jax" is the explicit way to ask for the JAX brute force. chunk_k still chunks the brute-force builders.
  • Interferometer.apply_sparse_operator and Interferometer.psf_precision_operator_from gain method="nufft", eps=None and nufft_chunk_size=None (None = take the transformer's).
  • New public builder nufft_precision_operator_via_nufft_from(...).
  • Behaviour change: callers get the NUFFT builder instead of the numpy brute force — same array, ~100–300× faster — whether or not they pass use_jax=True. Existing workspace apply_sparse_operator(use_jax=True) calls therefore land on the fast path with no edit; opting back into a brute force now takes an explicit method="numpy" or method="jax".

See full details below.

Test Plan

Full API Changes (for automation & release notes)

Added

  • autoarray.inversion.inversion.interferometer.inversion_interferometer_util.nufft_precision_operator_via_nufft_from(noise_map_real, uv_wavelengths, shape_masked_pixels_2d, grid_radians_2d, *, eps=1e-12, chunk_size=None) -> np.ndarray — builds the preload as a type-1 NUFFT.

Changed Signature

  • nufft_precision_operator_from(..., *, method="nufft", eps=1e-12, chunk_size=None, chunk_k=2048, show_progress=False, show_memory=False, use_jax=False)method, eps, chunk_size added; all pre-existing keywords unchanged.
  • Interferometer.apply_sparse_operator(nufft_precision_operator=None, batch_size=128, method="nufft", eps=None, nufft_chunk_size=None, chunk_k=2048, show_progress=False, show_memory=False, use_jax=False)
  • Interferometer.psf_precision_operator_from(..., method="nufft", eps=None, nufft_chunk_size=None, ...)

Changed Behaviour

  • The default preload builder is now the type-1 NUFFT rather than the numpy brute force. The array is identical to within 4.4e-14 of peak; cached .npy preloads stay valid.
  • The "N_vis · N_pix ≳ 1e7 crossover" warning in dataset.py was rewritten: the preload is now seconds, so the warning is about the DFT transformer, not the preload.

Migration

  • None required. To opt back into the reference builder: method="numpy" (or method="jax" for the JAX one).
  • Workspace scripts that pass use_jax=True get the NUFFT path unchanged — the kwarg is now a no-op under the default method, so no workspace edit is needed to buy the speed-up.

Closes #539

Generated by the PyAutoLabs agent workflow.

🤖 Generated with Claude Code

https://claude.ai/code/session_018hLF3ZAcz5MmaSJBEcLkvF

The sparse-operator preload is the real part of a type-1 (adjoint) NUFFT of the
inverse-variance weights 1/sigma^2 evaluated on the doubled-extent grid:
`nufft2d1(-x, y, w, (2Nx, 2Ny), eps, +1)`, followed by an `ifftshift` and zeroing
of the Nyquist row `Ny` and column `Nx`. That replaces the brute-force
O(N_pix * K) cosine accumulation with an O(K + N_pix log N_pix) transform.

- New `nufft_precision_operator_via_nufft_from`, and `method="nufft"` (the new
  default) on the `nufft_precision_operator_from` dispatcher, with `eps=1e-12`
  and mandatory visibility chunking taken from the transformer.
- The numpy and JAX brute-force builders are kept untouched as the reference
  implementations (`method="numpy"` / `"jax"`; `use_jax=True` still routes to
  the JAX one).
- Loud, logged fallback to `method="numpy"` when `PYAUTO_DISABLE_JAX=1` is set
  or nufftax is not importable; never silent, anything else raises.
- `dataset.py`: `apply_sparse_operator` and `psf_precision_operator_from` gain
  `method`, `eps` and `nufft_chunk_size`, defaulting to the transformer's own
  `eps` / `chunk_size` when it is a `TransformerNUFFT`.

Pins: mixed-tolerance parity against the numpy builder at 4.4e-14 of peak at
K=300; chunked == one-shot at 5.3e-14; the wrong-sign control fails at 0.187 of
peak. Measured in autolens_profiling#229: alma 2101 s -> 7.3 s wall (289x;
111x in CPU-seconds), alma_high 22 s.

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_018hLF3ZAcz5MmaSJBEcLkvF
@Jammy2211 Jammy2211 added the pending-release PR queued for the next release build label Sep 8, 2026
Base automatically changed from feature/interferometer-apply-operator-rfft2 to main September 8, 2026 17:22
…ce (#539)

Every existing workspace `apply_sparse_operator(use_jax=True)` call was mapped
onto `method="jax"`, silently putting the production path on the brute force;
`use_jax` is now only honoured when a brute-force method is already selected.

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_018hLF3ZAcz5MmaSJBEcLkvF
@Jammy2211
Jammy2211 merged commit 9bd7679 into main Sep 8, 2026
3 checks passed
@Jammy2211
Jammy2211 deleted the feature/interferometer-preload-nufft-type1 branch September 8, 2026 17:28
@Jammy2211

Copy link
Copy Markdown
Collaborator Author

Workspace PR: PyAutoLabs/autogalaxy_workspace#238 — corrects the interferometer preload prose in autogalaxy_workspace ("minutes to hours" -> built as a type-1 NUFFT in seconds). It waits on the next autoarray release per the pending-release chain.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

pending-release PR queued for the next release build

Projects

None yet

Development

Successfully merging this pull request may close these issues.

perf: build nufft_precision_operator_from as a type-1 NUFFT

1 participant