2.1. Fitting with SHAMFI

The simplest method to fit the basis functions in image space is to set up a number of linear equations, where for an image with $s \times t$ pixels, each pixel value P at location x, y can be described by

(7) \begin{equation} \begin{bmatrix} B_{0,0}(x_0,y_0) & \quad \cdots & \quad B_{p_{\text{max}},0}(x_0,y_0) \\ \vdots & \quad \ddots & \quad \vdots \\ B_{0,0}(x_s,y_t) & \quad \dots & \quad B_{p_{\text{max}},0}(x_s,y_t) \end{bmatrix} \begin{bmatrix} C_{0,0} \\ \vdots \\ C_{p_{\text{max}},0} \end{bmatrix} = \begin{bmatrix} P_{0,0} \\ \vdots \\ P_{s,t} \end{bmatrix},\end{equation}

where $C_{p_1,p_2}$ is a scalar coefficient for the basis function $B_{p_1,p_2}$ , for all basis functions where $p_1 + p_2 \leq p_{\text{max}}$ . This reduces the problem to the well-studied linear matrix $\mathbf{A}\mathbf{x} = \mathbf{b}$ . For convenience, we use the python function numpy.linalg.lstsq,Footnote g which solves for $\mathbf{x}$ by minimising the Euclidean two-norm $|| \mathbf{b} - \mathbf{A} \mathbf{x} ||^2$ .

We choose to fit the basis functions in image space, as it is far easier to isolate emission from a single object. Each visibility is the sum of the emission of all sources within the field of view, which is of particular concern for large field-of-view instruments such as the MWA. One can simply crop or edit in image space to remove difficult/unwanted pixels. CLEANed restored images offer an easy data set to fit, as the restoring beam tends to smooth out any sharp model pixelisation effects from the CLEAN fitting procedure. By design, these images are ‘true’ sky emission convolved with a Gaussian kernel known as the ‘restoring beam’, which is an approximation of the synthesised PSF. The synthesised PSF is the image-based manifestation of the incomplete sampling of an interferometer in the Fourier transform space in which it samples data, and is the reason images must be deconvolved. For a more accurate and instrument-independent shapelet fit, one must account for this. Rather than deconvolve the CLEANed image however, we achieve the equivalent by convolving the basis functions themselves by the restoring beam kernel, as outlined in the LOFAR cookbook mentioned in Section 1.

Right-ascension (RA) and declination ( $\delta$ ) are related to x, y via a rotation by the position angle $\phi_{\text{PA}}$ , defined as the increasing angle to East from North. Aside from fitting $C_{p_1,p_2}$ , $\beta_1$ , $\beta_2$ , and $\phi_{\text{PA}}$ are all variables. Minimising $|| \mathbf{b} - \mathbf{A} \mathbf{x} ||^2$ for many values of $\beta_1,\beta_2$ , and $\phi_{\text{PA}}$ is computationally expensive. We instead settle for an initial 2D-Gaussian fit [defined in Equation (10)] to the entire object to determine $\phi_{\text{PA}}$ , and then fit a range of $\beta_1,\beta_2$ in a grid, using the residuals of each fit to determine the optimal parameters. Examples of these steps are given in Sections 4.2 and 6.1.

3.1. Point source generation

WODEN analytically generates a sky model for calibration directly in visibility space via the measurement equation (c.f. Thompson, Moran, & Swenson Reference Thompson, Moran and Swenson2001)

(8) \begin{align} & V(u,v,w) = \\[4pt] & \int \mathcal{B}(l,m) I(l,m) \exp[-2\pi i(ul + vm + w(n-1))] \frac{dldm}{n},\nonumber \end{align}

where V(u, v, w) is the measured visibility at baseline co-ordinates u, v, w, given the sky intensity I(l, m) and instrument beam pattern $\mathcal{B}(l,m)$ , which are functions of the direction cosines l, m, with $n=\sqrt{1-l^2-m^2}$ . Equation (9) is discretised for point sources such that

(9) \begin{align} & V(u_i,v_i,w_i) =\\[5pt] & \sum_j \mathcal{B}(l_j,m_j)I(l_j,m_j) \exp[-2\pi i(u_il_j + v_im_j + w_i(n_j-1))],\nonumber \end{align}

where $u_i,v_i,w_i$ are the visibility co-ordinates of the $i{\text{th}}$ baseline, and $l_j$ , $m_j$ , $n_j$ is the sky position of the $j{\text{th}}$ point source.

This is an embarrassingly parallel problem, well suited to optimisation on GPUs. The basic application of a GPU to a problem such as this is to calculate all iterations over i, j in parallel, and then sum over j after the fact. While this efficiently leverages the parallel computational structure of the GPU, it requires large arrays to store the $i \times j$ outputs, which can quickly use large amounts of GPU memory when j becomes large.

In recent years, the CUDA GPU-language (Nickolls et al. Reference Nickolls, Buck, Garland and Skadron2008) has made progress in making efficient memory-managed read–modify–write operations on a single address, or so-called atomic operations.Footnote h These allow calculations to be made in parallel, and the results written out to a single address with ‘hardware-level’ thread-safe management, with minimal slowdown. We make use of atomicAdd, which adds the output of a thread to a single memory address, to perform the summation over j in Equation (10). This reduces the size of the array needed to store outputs significantly.

We make no attempt to add any instrumental or ionospheric effects, and write WODEN to simply evaluate Equation (10) with $\mathcal{B}(l,m) \equiv 1$ for a given array layout and observation parameters. We use WODEN here simply to test the modelling efficiency and model-generation costs of both the MS CLEAN and shapelet methods.

3.2. Gaussian and shapelet generation

To simulate peeling using both shapelets and MS CLEAN outputs, we add the functionality to simulate Gaussians and shapelets into WODEN. We define a 2D Gaussian in image space as

(10) \begin{equation}G(x,y) = \exp \left( -4\ln(2)\left[\frac{(x - x_0)^2}{\theta_{\text{maj}}^2} + \frac{(y - y_0)^2}{\theta_{\text{min}}^2} \right] \right),\end{equation}

where $x_0,y_0$ is the central point, and $\theta_{\text{maj}},\theta_{\text{min}}$ are the full-width half-maximum (FWHM) major and minor axes. Right-ascension and declination are related to x, y via a rotation by the position angle $\phi_{\text{PA}}$ . We utilise the RTS methodology of inserting a visibility ‘envelope’ $\xi$ into Equation (10) to create a Gaussian or shapelet source:

(11) \begin{multline}V(u_i,v_i,w_i) = \\\sum_j \xi_j(u_i,v_i)I(l_j,m_j) \exp[-2\pi i(u_il_j + v_im_j + w_i(n_j-1))],\end{multline}

where for a Gaussian,

(12) \begin{align}\xi_j = \exp\left( -\frac{\pi^2}{4\ln(2)} \left( k_x^2\theta_{\text{maj}j}^2 + k_y^2\theta_{\text{min}j}^2\right) \right), \end{align}

(13) \begin{align}k_x = \cos(\phi_{PAj})v_i + \sin(\phi_{PAj})u_i, \end{align}

(14) \begin{align}k_y = -\sin(\phi_{PAj})v_i + \cos(\phi_{PAj})u_i, \end{align}

[noting that as visibilities are the Fourier dual of image space, the normalisation factor is inverted from Equation (10)] and for a shapelet,

(15) \begin{align}\xi_j = \sum^{p_k +p_l < p_{\text{max}}}_{k,l} C_{p_k,p_l} \Tilde{B}_{p_k,p_l}(k_x,k_y), \end{align}

(16) \begin{align}k_x = \frac{\pi}{\sqrt{2\ln(2)}} \left[\cos(\phi_{PA})v_{i,j} + \sin(\phi_{PA})u_{i,j} \right], \end{align}

(17) \begin{align}k_y = \frac{\pi}{\sqrt{2\ln(2)}} \left[-\sin(\phi_{PA})v_{i,j} + \cos(\phi_{PA})u_{i,j} \right], \end{align}

where $u_{i,j},v_{i,j}$ are visibility co-ordinates for baseline i, calculated with a phase-centre $RA_j,\delta_j$ , which corresponds to the central position $x_0,y_0$ used to fit the shapelet model in image space. Inserting the visibility envelope $\xi$ in this way causes a convolution in visibility space, which in turn causes a multiplication in image space. Using Equation (10) causes a convolution with a Gaussian kernel that sums to one, meaning the visibility can simply be multiplied by I(l, m) to give the correct flux density. When using Equation (15), one must scale the coefficients $C_{p_k,p_l}$ to ensure the kernel sums to one, to use Equation (11) for both Gaussian and shapelet components. In Equation (15), then we set $C_{p_k,p_l} \equiv C_{p_k,p_l} / I(l_j,m_j)$ , and note that the coefficients output by SHAMFI have this division by integrated flux density applied, hence Equation (15) is implemented in WODEN exactly as it appears.

The shapelet basis function values $\Tilde{B}_{p_k,p_l}(u,v)$ can be calculated by interpolating from one-dimensional (1D) look-up tables of $\Tilde{B}(k_x;\,1)$ [c.f. Equation (1)], and scaling by the appropriate $\beta$ . When using the scalings presented in Equations (13), (14), (16), and (17), the first-order shapelet basis function is equal to a 2D Gaussian when $\beta_1=\theta_{\text{maj}}, \beta_2=\theta_{\text{min}}$ , explicitly $G(x,y,\theta_{\text{maj}},\theta_{\text{min}}) \equiv B_{0,0}(x,y,\theta_{\text{maj}},\theta_{\text{min}})$ . This is useful to test the consistency between simulating Gaussians and shapelets.

For all simulation bench-marking used in this paper, we use a GeForce GTX 1080 Ti NVIDIA GPU, with 12 GB of memory. This allows us to store all time steps in GPU memory, meaning the optimisation over j includes all time steps. The net result of these optimisations is that per frequency channel, we launch a single CUDA kernel instance each for all point sources, Gaussians, and shapelets, all of which are optimised with a consistent strategy. As u, v, w only differ in frequency by a scaling factor, we only calculate them once per time step.

4.1. Very Large Array (VLA) Fornax A Model

The goal of this work is to compare shapelet and MS CLEAN methods to generate models of extended radio galaxies. To fairly test the two methods as implemented in SHAMFI and WSClean, in reality one would start with calibrated visibility data, deconvolve and image this data using WSClean to create MS CLEAN components, and create an image to use when modelling with SHAMFI. This would ensure both methods use the same data set, with the same angular resolution, and both suffer any systematics inherent to the data. As stated in Sections 1 and 3, we choose to test the methods first using simulated data, to reduce the number of systematics. To do this, we need to create simulated ‘observed’ data in visibility space. We choose to simulate observed visibility data using point sources alone, to minimise biasing towards either method. To build a realistic and high-resolution model of a complex and extended radio galaxy, we use a public VLA image of Fornax A at 1.4 GHz, obtained from the NASA/IPAC Extragalactic Database,Footnote i a digitised version of the data presented in Fomalont et al. (Reference Fomalont, Ebneter, van Breugel and Ekers1989). The online FITS file is a cut-out around the galaxy, containing some CLEAN residuals. To create a point source model, we first apply a Gaussian taper around the edge of the image (see Figure 1). Not only does this reduce CLEAN residuals around the image, but also removes the sharp boundary edge, minimising future ringing when imaging from the simulated visibility data. We then move the image from B1950 to J2000 co-ordinates, convert each pixel into a point source [explicitly we convert each pixel into a Dirac delta function, allowing us to run Equation (10) over all pixels], and scale the flux density of all pixels to sum 500 Jy at 180 MHz. We give each pixel a single spectral index of $-0.8$ . As the restoring beam of this CLEAN-ed image is smaller than the resolution of the instruments used in simulation, we make no attempt to remove it from the image.

Figure 1. Upper panel: Fornax A image obtained from NED, which includes CLEAN residuals around the edge of the lobes. Over-plotted in red is a boundary line outside of which a Gaussian taper was applied to create the model used in WODEN. Lower panel: The image after applying the taper and transforming from B1950 to J2000 centred coordinates. Any pixels that were not masked (the grey region) in the lower image were converted into point sources.

We simulate six MWA observations of 30.72 MHz bandwidth and 10 kHz spectral resolution using WODEN. These simulations use the observational parameters of the real data used to image Fornax A in Section 6.1. As the MWA has recently been upgraded with added receiver elements, dubbed Phase II, it can be operated in various configurations, with different angular resolutions and sensitivities to diffuse or compact emission (see Wayth et al. Reference Wayth2018, for details). Three of the observations use the Phase I configuration (maximum baseline $\sim3\,$ km), with the other three the Phase II extended configuration (maximum baseline $\sim$ 5.5 km, hereto referred to as Phase II). The Phase I observations have 2 s time resolution and the Phase II observations 0.5 s resolution to reduce time decorrelation on the longer baselines. We simulate 56 time steps for all observations regardless of the time resolution to reduce computational load. We use WSClean to perform an MS CLEAN on the simulated data, with the results shown in Figure 2. These images have not been beam corrected, given that WODEN does not include beam effects. We include an example WSClean command in Figure A1 to detail the settings used. We also include all software versions used during this research in Table B1.

Figure 2. MS CLEANed images of the WODEN simulation of Fornax A. Top: All Phase I configuration, imaged with briggs 0 weighing to balance the large number of short baselines present in Phase I data with a reasonable resolution. Bottom: A combination of both Phase I and Phase II configurations. These are imaged with uniform weighing to take advantage of the higher resolution Phase II data.

4.2. Modelling the Fornax A simulation

We fit the simulated Fornax A in the images in Figure 2 using SHAMFI as described in Section 2, and use Equation (5) to set $p_{\text{max}}$ . We set $\vartheta_{\text{max}} = 0.5^{\circ}$ , and set $\vartheta_{\text{min}}$ by oversampling the angular resolution of the array by 3, giving $\vartheta_{\text{min}} = \lambda / 3b_{\text{max}}$ , where $\lambda$ is the wavelength of the observation, and $b_{\text{max}}$ the maximum baseline of the array. We set the latter to 3 km for the Phase I array and 5.5 km for the Phase II array. This yields $p_{\text{max}} = 46$ for Phase I and $p_{\text{max}} = 86$ for Phase II.

Figure 3. Left: Example of initial Gaussian fit used to set $\phi_{\text{PA}}$ of the basis functions for one of the simulated Phase 1+2 lobe images. The red line shows the FWHM, with the white lines demonstrating the PA found. Right: Example of the grid-based approach for fitting $\beta_1$ and $\beta_2$ . The colour scale here represents the residuals in (Jy/pixel) $^2$ left after fitting the image in the left panel with all basis functions up to $p_{\text{max}} = 86$ .

Figure 4. Left: MS CLEANed Phase 1 simulated data Middle: Fitted shapelet model recreated in image space, using restoring beam convolved basis functions Right: The model subtracted from the data. The footprints of the two separate shapelet models (one for each lobe) are clearly visible in the residual plot.

During modelling, we found the choice of location of the zero pixel $x,y = 0,0$ of the basis functions severely affected the quality of fit. We found that for a very extended, double-lobed radio galaxy like Fornax A, the best solution was to split the galaxy into two lobes and fit each lobe separately. To avoid double-fitting flux, we divided the image in two by fitting each lobe with a normalised 2D Gaussian, which we label $N_1(x,y)$ , $N_2(x,y)$ , and then applying weights $w_1(x,y),w_2(x,y)$ to each pixel such that $w_1 = N_1 / (N_1 + N_2)$ , $w_2 = N_2 / (N_1 + N_2)$ . Figure 3 shows an example of fitting for $\phi$ and the grid-search optimising for $\beta_1,\beta_2$ . The fitting results for the simulated Phase I data are shown in Figure 4.

As is also prescribed in the LOFAR cookbook, we found it necessary to subtract any compact point source-like emission from the images including Phase II data, which reduced the number of higher-order shapelets required for a good fit. We do this manually by subtracting image-based Gaussian components, as we found traditional source finders struggle to deal with the complexity of Fornax A to successfully isolate the point-like details in the image. A full example of the image preparation needed to model Fornax A is given in Section 6.1. We note the tools to split a radio galaxy into Gaussian masked regions, and manually subtract Gaussians, are included in the SHAMFI package.

5.1. Peeling method on simulated data

To fit the gains, we follow the calibration scheme implemented in YANDAsoft.Footnote j As this scheme has yet to be published, we detail the formalism here. We label the ideal visibility model for baseline ij as $\mathcal{I}_{ij}$ . The visibility model $M_{ij}$ (over the range of time and frequency that antenna gains are assumed constant) is

(18) \begin{equation} M_{ij}(t,\nu) = g_ig_j^* \mathcal{I}_{ij}(t,\nu).\end{equation}

where $g_i$ is the gain on antenna i. The measured visibility $V_{ij}$ is

(19) \begin{equation} V_{ij}(t,\nu) = (g_i+e_i)(g_j+e_j)^*\mathcal{I}_{ij}(t,\nu),\end{equation}

where $e_i$ is an additive error on the gain $g_i$ . We label the residual visibilities as

(20) \begin{equation} R_{ij}(t,\nu) = V_{ij}(t,\nu) - M_{ij}(t,\nu).\end{equation}

The gains can be optimised by iteratively solving for and following local error gradients. These are found using linear least-squares with data $\Re(R)$ and $\Im(R)$ and free parameters $\Re(e)$ and $\Im(e)$ . The complex conjugation of gains in the equations above cannot be simply represented by a linear operator. So rather than solving for complex solutions, we instead solve for the real and imaginary parts separately and multiply the relevant coefficients by $-1$ . Assuming the error-squared term is small allows the problem to be linearised as $\mathbf{A}\mathbf{x} = \mathbf{b}$ , where $\mathbf{A}$ is the design matrix, $\mathbf{x}$ are the parameters, and $\mathbf{b}$ is the residual vector, e.g. the real parts of the residual vector $\mathbf{b}$ are

(21) \begin{equation} \Re[R_{ij}(t,\nu)] = \Re[ e_j^*g_iI_{ij}(t,\nu) + e_ig_j^*I_{ij}(t,\nu)].\end{equation}

This can be implemented to solve for all $e_i$ in the vector $\mathbf{x}$ as

(22) \begin{equation}\mathbf{x} = \left[ \mathbf{A}^\intercal \mathbf{A} \right]^{-1}\mathbf{A}^\intercal \mathbf{b}.\end{equation}

The assumption of linearity will improve as the system converges. We found that as the simulations had no instrumental or source side-lobe noise, the gains consistently converged after 10 iterations of the calibration loop, and so used 10 iterations for all the results shown in this and following sections. We fit a separate gain for each time and frequency step. Once we have the gains, we then fit a cubic spline across the entire bandwidth using the scipy.interpolate.interp1d Footnote k function, to ensure we apply spectrally smooth gains during peeling. This was found to be necessary for an EoR detection by Mouri Sardarabadi & Koopmans (Reference Mouri Sardarabadi and Koopmans2019), as without smooth gains, the EoR signal can be suppressed. The caveat to this being that we are assuming a linearised calibration case, where spectrally smoothing the gains does yield improvement. In the non-linear regime, this may not be true; however given that these simulations have zero noise and instrumental effects, the error terms in calibration are likely small.

Figure 5. Simulated peeling results. Plots (a) to (d) show MS CLEANed peel residuals of an integration over the full bandwidth and time span of a single simulated Phase 1 observation, each peeled with a different model for Fornax A. The models are: (a) a shapelet model from Phase I data; (b) a shapelet model from Phase 1 and 2 data; (c) an MS CLEAN model from Phase I data; (d) an MS CLEAN model from Phase I and II data. A contour plot of the unpeeled power is shown on each plot as a angular scale reference. Plot (f) shows the 1D PS as estimated with CHIPS for the models in (a) to (d). For reference, the power without peeling is shown for the unpeeled WODEN simulation in (e). The line style and marker for each 1D PS plotted in (f) is also plotted in each plot from (a) to (d) on the lower right for reference.

5.2. Peeling results

The results of peeling all four models from the full simulation are shown in Figures 5 and 6. We image the peeling residuals in Figure 5 using the same WSClean commands used to image the full simulation in Figure 2. We make 1D (Figure 5) and 2D (Figure 6) PS estimates using CHIPS (Trott et al. Reference Trott2016), and note that the values of the powers in Figures 5 and 6 are indicative, rather than absolute. This is because CHIPS grids using the primary beam shape of the MWA, and normalises the outputs based on the expected cosmological observing volume of the telescope. Neither are present in WODEN simulations. All PS estimates were run with identical settings however, and so all differences are solely due to the models used for peeling.

Figure 6. The 2D PS as estimated by CHIPS of the simulated peeling results. Plots (a) and (b) show the shapelet model peel results for Phase I and Phase I + II, with (c) showing the difference between the 2D PS in (a) and (b). The MS CLEAN peel results are shown in a similar manner in plots (d), (e), and (f). In the difference PS of (c) and (f), blue means the Phase I model subtracted less power than the Phase I + II model, and red means the Phase I model subtracted more power.

Figure 7. Results from compressing the MS CLEAN and shapelet models. Plots (a), (b), and (c) on the left show results from Phase I models, where (d), (e), and (f) on the right show results from Phase I + II models. The top row show residual 1D PS left after peeling MS CLEAN models (a) and (d) and shapelet models (b) and (e), at various levels of truncation. The percentage in the legend is the percentage of remaining components after truncation. We take slices through the truncation results at specific low k values, shown by the vertical dotted lines, and plot them as a function of GPU time in the bottom row (c) and (f). In both plots, dashed lines are MS CLEAN, and solid lines are shapelet results. Matching colours plot matching k-modes.

Figure 8. Stages of fitting for the Fornax A image using real MWA data as described in Section 6.1 and shown in (a). Four compact sources (outlined in red circles) were found through trial and error to cause fitting problems, and so were subtracted as Gaussians. The image was then split into two lobes, (b) and (c), and fitted separately. (d) shows the fitted shapelet model recreated in image space, using restoring beam convolved basis functions, as well as the subtracted Gaussians. (e) shows the fitting residuals.

The 1D PS is an average over the Fourier inverse (or k-space) of three dimensions; two angular dimensions across the sky, and one derived from the spectral dimension of the data. The 1D PS is currently used to generate upper limits on an EoR detection. In Figure 5, we include a fiducial EoR signal as a point of reference to the reader. This is taken from the public websiteFootnote l associated with the work of Mesinger, Greig, & Sobacchi (Reference Mesinger, Greig and Sobacchi2016). We select a redshift that matches the centre of the simulated frequency bandwidth. The greatest ratio of expected EoR signal to astrophysical foregrounds, after foreground removal/avoidance, is expected to lie at low k-modes in the 1D PS. In Figure 5, we concentrate on the smallest k-modes. We find little difference when peeling with the Phase I MS CLEAN and shapelet models (purple line with squares, and blue line with circles, respectively). We see a difference in behaviour in the Phase I + II models; however, including Phase II data in the model improves the subtraction at small k with the shapelet model, but actually worsens the subtraction with the MS CLEAN model.

To understand the difference, we turn to the 2D PS, shown in Figure 6. The 2D PS averages the k-modes obtained from angular scales upon the sky ( $k_\perp$ , horizontal axis), and plots them against the k-modes derived from the spectral response ( $k_\parallel$ , vertical axis). Broadly, astrophysical foregrounds are expected to be spectrally smooth, and so should inhabit the bottom of the 2D PS, where $k_\parallel$ is small. Instrumental chromaticity causes the foregrounds to be swept up at higher $k_\perp$ into the so-called ‘wedge’ on the lower right, with the upper left being known as the ‘window’, offering the cleanest area of k-space in which to attempt a detection. If the peeled residuals were perfectly spectrally smooth, we would expect this window to be devoid of emission; however, there is clearly power present. There are two processes that could have added in spectral structure: the calibration step during peeling, and the gridding kernel of CHIPS, which expects the data to have the imprint of an MWA beam pattern, which was not included in the simulation. Visibilities with quickly changing spectral structure have the effect of throwing power that would normally sit at low $k_\parallel$ up into high $k_\parallel$ , contaminating the window. This effect can happen in real data, and the net effect is more power in the wedge at a certain $k_\perp$ will also throw more power into the window at that $k_\perp$ .

This effect goes some way to explaining the difference in 1D power when adding higher angular resolution data into the modelling process between shapelets and MS CLEAN. When looking at the difference of the Phase I and Phase I + II shapelet model residuals [Figure 6(c)], we see more power has been subtracted out of the wedge at all $k_\perp$ beside the very smallest $k_\perp$ -bin, which has meant there is less power in the window also. Comparing the MS CLEAN models, the introduction of higher angular resolution data has yielded a large improvement in the power subtracted at large $k_\perp$ (small angular scales), but has worsened the subtraction at low $k_\perp$ , contaminating the window. It is difficult to isolate what might cause worsening subtraction at large angular scales with MS CLEAN. It is possible that as including the Phase II configuration reduces the angular size and shape of the synthesised PSF, an MS CLEAN is more biased to smaller angular scales. There are two WSClean parameters that can be manually set to change the MS CLEAN settings: -multiscale-scales, which forces the CLEAN to specific scales, and -multiscale-scale-bias, which biases the CLEAN to larger or smaller scales. We tested forcing the scales that were CLEAN-ed for the Phase I + II image to be the same as the scales CLEAN-ed for the Phase I alone (plus an extra smaller scale to account for the increase in resolution). We also changed -multiscale-scale-bias to favour larger scales. We propagated both these changes through to the PS, and found either no change in residuals, or worsened residuals at all scales. While not an exhaustive parameter search, we find the default settings give the best results, and so use them throughout the paper. What is clear is that the SHAMFI model was able to incorporate the increased angular resolution, while simultaneously still fitting the large-scale structure.

5.3. Model truncation

As discussed in Section 5, calculating a value for a single shapelet basis function requires more calculation than for a single point source or Gaussian. As a point of reference, for the simulated Phase I + II image, WSClean produced a total of 6 331 MS CLEAN components (4 540 point source and 1 791 Gaussian). Fitting a separate shapelet model for each lobe, each with $p_{\text{max}} = 86$ , gives a total of 7 656 basis functions to extrapolate. Therefore, for shapelets to be computationally competitive, some compression/truncation of the full model is necessary. We apply a basic truncation, by simply removing the MS CLEAN components or shapelet basis functions that contribute the smallest absolute flux density to the total number of components $N_{\text{tot}}$ , up until some fraction $f_{\textrm{comp}}$ of all components remain. Explicitly, we find the number of components $N_{\textrm{comp}}$ that satisfies

(23) \begin{align}\sum^{N_{\textrm{comp}}}_j |I_j| <= f_{\textrm{comp}} \sum^{N_\textrm{tot}}_j |I_j|, \end{align}

(24) \begin{align}\textrm{where}\quad\quad & |I_j| > |I_{j+1}|.\end{align}

For a point source or Gaussian, $I_j$ is simply the total integrated flux density, with the absolute value necessary as MS CLEAN components and shapelets can be positive or negative. To find $I_j \equiv I_{p_1,p_2}$ for the shapelet basis function $B_{p_1,p_2}$ , we image and sum the absolute of the restoring beam $G(\theta_{\text{bmaj}},\theta_{\text{bmin}})$ convolved basis function,

(25) \begin{equation} I_{p_1,p_2} = \sum_s \sum_t \left| C_{p_1,p_2} G(\theta_{\text{bmaj}},\theta_{\text{bmin}}) \ast B_{p_1,p_2}(x_s,y_t) \right|,\end{equation}

where $\ast$ is a convolution, for all $s \times t$ pixels used during the model fit.

We have implemented this type of truncation within SHAMFI, which has the necessary code base to allow the compressed (reduced) set of shapelet basis functions to be fit to the image again. We apply this extra step in the following results.

We apply three levels of truncation to both MS CLEAN and shapelets, for both the Phase I and Phase II models. We then simulate these models through WODEN for the same 2-min observation as used in Section 5, and apply the same peeling technique, fitting a single complex gain per antenna with 10 iterations of Equation (22). When running the simulations, we used nvprof,Footnote m the NVIDIA line profiler, to measure the time spent performing calculations and memory allocations on the GPU. For each level of truncation, we ran the peeled outputs through CHIPS, so we could compare the effect of peeling against time taken in generating the visibilities. The results are shown in Figure 7. As expected, the uncompressed shapelet model takes more than twice the time to generate than the MS CLEAN model. However, for both the Phase I and Phase I + II models, shapelets can be compressed to a level where they incur comparative computational expense, and still peel out more power at the largest spatial scales. By comparing the blue dashed lines in the bottom plots of Figure 7, we can see that the 100% MS CLEAN model subtracts more large-scale power when only including Phase I information, (bottom left), rather than including both Phases I and II. The amount of power subtracted by the shapelet model is comparable when incorporating just Phase I information or both Phases I and II.

Table 1. Observational parameters of the real data used to test peeling.

6.1. Imaging and modelling real MWA data

We use six 2-min MWA observations to create the image of Fornax A in Figure 8. Three observations were Phase I, taken in December 2014 as part of the GLEAM survey (Hurley-Walker et al. Reference Hurley-Walker2017; Wayth et al. Reference Wayth2015), and the other three were Phase II, taken in February 2018 as part of the GLEAM-X survey (Hurley-Walker Reference Hurley-Walker2020, in preparation). Calibration observations of the bright radio galaxy Pictor A were also taken at the same frequency on the same nights as the Fornax A observations. The Pictor A observations were used to solve for an initial set of calibration solutions, which were then applied to the Fornax A data. An image was made with the full 30.72 MHz bandwidth of the MWA, centred at 185 MHz using joint deconvolution with WSClean. This joint deconvolution was enabled by the implementation of Image Domain Gridding (IDG, van der Tol, Veenboer, & Offringa Reference van der Tol, Veenboer and Offringa2018) into WSClean, which uses the MWA primary beam (Sokolowski et al. Reference Sokolowski2017) to properly take into account the changing beam shapes for each snapshot observation in the gridding process. Several rounds of self-calibration were then performed with CALIBRATE (Offringa et al. Reference Offringa2016) to improve the antenna gain solutions, using the MS CLEAN component model produced by WSClean.

In Figure 8, we detail the image manipulation required when running SHAMFI to obtain a good model using real data, including the necessary subtraction of compact emission to reduce the need for higher-order (and therefore higher-resolution) shapelet basis functions. The process of fitting the real Fornax A image, along with all commands used, is available online as a tutorial in the documentation of SHAMFI.Footnote n

Figure 9. Real data peel results for Fornax A data detailed in Table 1: (a) the calibrated data before subtraction; (b) and (c) the residuals after peeling the shapelet model and MS CLEAN model, respectively; (d) and (e) the 2D PS of the residuals after peeling the shapelet model and MS CLEAN model, respectively; (f) the difference between the two 2D PS shown in (d) and (e), with the MS CLEAN PS subtracted from the shapelet PS. In (f), red means the shapelet model subtracted more power during the peel, and blue means the MS CLEAN model subtracted more.

Figure 10. Real data peel results for the EoR1 data detailed in Table 1: (a) and (b) the new Phase I + II shapelet model simulated using the RTS and the residuals after peeling with this model, respectively; (c) and (d) the old shapelet model and residuals after peeling, respectively; (e) an east-west polarisation difference 2D PS where the residuals of (d) (old shapelets) are subtracted from (b) (new shapelets); (f) 1D PS of the residuals left behind after peeling the new shapelet mode (blue line with squares) and the old model (orange with triangles), with an estimate of the thermal noise in the east-west polarisation; (g) the ratio of the 1D PS shown in (f), with the new model divide by the old model. In (e), red means the new model subtracted more power during the peel, and blue means the old model subtracted more. We only show east-west here, but the north-south display the same behaviour.

6.2. Calibrating and peeling real data

To carry out calibration and peeling on the real data sets, we use the RTS, software which is used by the Australian MWA EoR collaboration (see Jacobs et al. Reference Jacobs2016, for an overview of RTS processing). The RTS has the capacity to perform direction-dependent calibration, correct for first-order ionospheric effects, and peel sources (phase-rotate, calibrate, and subtract a source, see Noordam Reference Noordam and Oschmann2004). During the imaging of the real data, WSClean produced 2 537 MS CLEAN components. These components, as well as the shapelet model described in Figure 8, were used to peel with the RTS. We peel Fornax A from two data sets: the first set where Fornax A is near the primary beam centre, and therefore contributes a large fraction of the total power in the sky; the second set made of standard MWA EoR observations, where Fornax A is far from pointing centre, often at a total Stokes I beam power of $\sim 30\%$ . The latter data set is pointed towards the ‘EoR1’ field at $RA=4^h,\delta=-27^\circ$ . We use the ionospheric metrics detailed in Jordan et al. (Reference Jordan2017) to pick an ionospherically quiet night for the EoR1 data. The key points of both data sets are summarised in Table 1. For both data sets, we use the same sky model, and simply switch in the MS CLEAN or shapelet model of Fornax A. The sky model is mostly based on the GLEAM catalogue (Hurley-Walker et al. Reference Hurley-Walker2017), cross-matched to the following catalogues and frequencies: VLSSr, 74 MHz (Lane et al. Reference Lane, Cotton, Helmboldt and Kassim2012); TGSS ADR1, 150 MHz (Intema et al. Reference Intema, Jagannathan, Mooley and Frail2017); MRC, 408 MHz (Large et al. Reference Large, Mills, Little, Crawford and Sutton1981); SUMSS, 843 MHz (Mauch et al. Reference Mauch, Murphy, Buttery, Curran, Hunstead, Piestrzynski, Robertson and Sadler2003); NVSS, 1 400 MHz (Condon et al. Reference Condon, Cotton, Greisen, Yin, Perley, Taylor and Broderick1998). We use PUMA (Line et al. Reference Line, Webster, Pindor, Mitchell and Trott2017) for the cross-match, using higher-resolution data for position or source morphology where possible. We use the sky model developed in Procopio et al. (Reference Procopio2017) in the EoR1 field. This sky model has been used by the MWA EoR analysis for a number of years, most recently in the limits published in Trott et al. (Reference Trott2020).

We use the RTS to perform foreground subtraction in two steps: an initial direction-independent calibration to a sky model of 1 000 sources; a subtraction step where a number of sources are fitted for first-order ionospheric effects and peeled. Even though we are only comparing the Fornax A models, we need to subtract off a number of sources in the sky to be able to see the difference in residuals, lest it be swamped by the rest of the sky. Typically, we peel off 1 000 sources in the latter stage. As discussed in Section 3 however, the architecture of the RTS is currently such that creating a single calibrator source with thousands of components causes significant computational expense. We found as the Fornax A data set had a higher time resolution, when processing with the MS CLEAN model, we could only peel off 500 sources before running into GPU memory issues. The processing time using the MS CLEAN model was also around 4 times slower than when using the shapelet model. The results of peeling 500 sources from the Fornax A data are shown in Figure 9, which shows both image residuals, and 2D PS.

In image space, the MS CLEAN and shapelet models leave similar residuals. A limitation of the shapelet fitting routine developed in this paper is exposed in Figure 9(b). The residuals are surrounded by a negative over-peel, which we attribute to side-lobe noise from the original image that was fitted into the shapelet model. This adds in a ‘floor’ to the model, which the RTS then has to calibrate for, leading to a slight under-peel of the lobes, and a slight over-peel of the floor. Interestingly, a similar negative border surrounds the MS CLEAN residuals, which is possibly due to the largest scale fitted for during the CLEANing process, or a consequence of the peeling algorithms internal to the RTS. While the 2D difference PS in Figure 9 suggests that the shapelet model has subtracted more power at all scales, the difference in the window is negligible, with the 1D PS looking near identical.

Given the small difference seen between the two models, we decided to compare our shapelet model to an earlier shapelet model created from early Phase I data, for the EoR1 data. This model has historically been used by the Australian EoR team. Fornax A has long been suspected of causing poor limits in the 1D PS. The results of comparing these two models are shown in Figure 10, and again we see little difference in the resultant 1D PS, regardless of the stark contrast in residuals left behind in the image. We explore possible reasons for the small differences in the following section.

Figure 11. Residual 1D PS comparing peeling 1 000 simulated sources, including the new Fornax A MS CLEAN model, when peeling sources including the new shapelet model (blue with squares) and including the old shapelet model (orange with triangles).