Include a primary beam correction for interferometer inversions

[HRSAstro]

Primary beam in interferometer inversions

As far as I can tell, autoarray has no primary-beam handling at present: the interferometer transform maps the image straight to visibilities, so an inversion reconstructs the apparent sky (primary beam × true sky). In pyuvimage we currently correct for the beam afterwards, and we have added an option to include it in the model. However, this correction is even more important for gravitational lens modelling, where not including the primary beam response would affect conservation of surface brightness.

What it changes. With P the primary beam on the image grid, the response becomes A' = T P M, so

• F' = Mᵀ P W̃ P M
• D' = Mᵀ P D̃ (D̃ the weighted dirty image)

The w-tilde kernel, the dirty image and dᵀN⁻¹d are unchanged, because P does not enter the dirty beam. On the sparse path this is two element-wise multiplications by P around the existing convolution, so it costs nothing extra. If P is treated as constant across a mesh pixel (P_m), this reduces to F' = P_m F P_m and D' = P_m D. That is equivalent to a regularisation matrix P_m⁻¹ H P_m⁻¹ on the apparent sky, and that is how we implemented it: a wrapper around any regularization (PrimaryBeamPrior). The evidence then compares the two treatments directly.

What we see. On a mock with emission out to P ≈ 0.3, the prior on the true sky gave a lower overall error and a total flux close to the truth, and the evidence preferred it, though faint emission near the beam edge came out biased low. So it seems a useful option rather than a new default.

Mosaics. The stronger case for doing this upstream is probably mosaics, where each pointing j has its own beam and kernel:

• F = Σ_j Mᵀ P_j W̃_j P_j M
• D = Σ_j Mᵀ P_j D̃_j

This would need one w-tilde kernel per pointing, which sparse_terms_from_chunks could accumulate per pointing in the same pass.

Questions

1. Would you want this in autoarray, and if so, where should P live: on the dataset, the mapper, or the sparse operator?
2. A Gaussian beam (FWHM ≈ 1.13 λ/D) covers ALMA single pointings (D = 12m for ALMA). For heterogeneous arrays (all antennas the same size), this is easy, but it's complicated for heterogeneous arrays. Should it also accept a user-supplied beam image?
3. Is mosaic support of interest? Happy to share the pyuvimage implementation and mock if helpful.

[Jammy2211]

@HRSAstro thank you, the maths checks out for a single pointing and we would like this in autoarray. Our preferred shape: the primary beam lives on the Interferometer dataset as an optional image on the real-space mask, and is applied inside both the transformer (dense/NUFFT path) and the sparse operator, so that light profiles, linear light profiles and the mapper all see the same P.

Two things we found while checking. The regularisation-wrapper form is only exact when P is constant across a mesh pixel, and for a lensed source the mesh lives in the source plane, so evaluating P on the mesh grid is the wrong plane; applying P on the image grid avoids that. And the proposal covers the mapper only, but profile images and model visibilities go through the transformer and need P too, or a fit mixes true-sky and apparent-sky components.

We will take a user-supplied beam image as the core API with a Gaussian helper (FWHM ≈ 1.13 λ/D) on top, and defer mosaics until there is a concrete target.

To guide the implementation, could you share a working Python snippet or small example from your build (ideally the P ≈ 0.3 mock and the script that produced the evidence comparison), and say whether both runs there used the same regularisation coefficient?

[HRSAstro]

upstream_snippets.py

I attached a Python file that includes the primary beam methods. Note that this code also contains snippets of the matrix-free solution implementation (#39).

[HRSAstro]

I don't have a direct comparison of the regularisation with/without PB. But some old tests on mock data showed the regularisation coefficient converged to different values with the PB in the solve, but it depended on the type of regularisation - for Matern solve the regularisation went down, but for adaptive it went slightly higher. The evidence favoured PB for Matern but was about the same for adaptive. The mocks were not gravitational lenses, though.

[Jammy2211]

@HRSAstro thanks for the file and the regularisation notes. We have opened PyAutoLabs/PyAutoArray#640 and the implementation is now up as PyAutoLabs/PyAutoArray#644: the beam lives on the Interferometer dataset (user-supplied image on the real-space mask, plus a Gaussian helper with FWHM 1.13 λ/D), is multiplied in on the image grid inside both transformers and the sparse w-tilde operator so profiles, linear profiles and mappers all see the same P, and the dirty image and kernel stay unchanged as in your derivation. Your P⁻¹HP⁻¹ wrapper remains expressible as a custom regularisation. Your observation that the preferred coefficient shifts with the beam for Matérn but not adaptive priors is useful; we will check that on the lensed mock.