diff --git a/autoarray/dataset/abstract/dataset.py b/autoarray/dataset/abstract/dataset.py index 25012376a..97801a4a2 100644 --- a/autoarray/dataset/abstract/dataset.py +++ b/autoarray/dataset/abstract/dataset.py @@ -134,34 +134,62 @@ def __init__( @property def grid(self): + """ + The primary coordinate grid of the dataset, equivalent to `grids.lp`. + + Returns the light-profile `Grid2D` aligned with the centres of all unmasked image pixels. + This is the grid used for the majority of model calculations (e.g. evaluating galaxy light + profiles). + """ return self.grids.lp @property def shape_native(self): + """ + The 2D shape of the dataset image in its native (unmasked) dimensions, e.g. (rows, columns). + """ return self.mask.shape_native @property def shape_slim(self): + """ + The 1D size of the dataset data array after masking, i.e. the number of unmasked pixels. + """ return self.data.shape_slim @property def pixel_scales(self): + """ + The (y, x) arcsecond-to-pixel conversion factor of the dataset, as a (float, float) tuple. + """ return self.mask.pixel_scales @property def mask(self) -> Union[Mask1D, Mask2D]: + """ + The mask of the dataset, derived from the mask of the `data` array. + """ return self.data.mask def apply_over_sampling(self): + """ + Apply new over-sampling sizes to the dataset grids. + + Subclasses must implement this method to rebuild the `GridsDataset` with updated + `over_sample_size_lp` and `over_sample_size_pixelization` values. + """ raise NotImplementedError @property def signal_to_noise_map(self) -> Structure: """ - The estimated signal-to-noise_maps mappers of the image. + The signal-to-noise map of the dataset, computed as `data / noise_map`. - Warnings arise when masked native noise-maps are used, whose masked entries are given values of 0.0. We - use the warnings module to suppress these RunTimeWarnings. + Values below zero are clamped to zero, as negative signal-to-noise is not physically + meaningful (it indicates the data is below zero due to noise, not a real negative signal). + + RuntimeWarnings from dividing by zero in masked pixels (where the noise map is 0.0) are + suppressed, as these masked values are never used in downstream calculations. """ warnings.filterwarnings("ignore") @@ -172,7 +200,7 @@ def signal_to_noise_map(self) -> Structure: @property def signal_to_noise_max(self) -> float: """ - The maximum value of signal-to-noise_maps in an image pixel in the image's signal-to-noise_maps mappers. + The maximum signal-to-noise value across all unmasked pixels in the dataset. """ return np.max(self.signal_to_noise_map) @@ -185,6 +213,27 @@ def noise_covariance_matrix_inv(self) -> np.ndarray: return np.linalg.inv(self.noise_covariance_matrix) def trimmed_after_convolution_from(self, kernel_shape) -> "AbstractDataset": + """ + Return a copy of the dataset with all arrays trimmed to remove the border pixels affected + by PSF convolution edge effects. + + When a model image is convolved with a PSF kernel, the pixels at the border of the image + cannot be correctly convolved because they lack sufficient neighbouring pixels. These border + pixels have unreliable values after convolution. This method trims the `data`, `noise_map`, + `over_sample_size_lp` and `over_sample_size_pixelization` arrays by the kernel half-width + on each side, so that only pixels with a complete convolution kernel neighbourhood remain. + + Parameters + ---------- + kernel_shape + The (rows, cols) shape of the PSF convolution kernel. The dataset arrays are trimmed + by `kernel_shape // 2` pixels on each side in each dimension. + + Returns + ------- + AbstractDataset + A shallow copy of the dataset with all arrays trimmed to the post-convolution shape. + """ dataset = copy.copy(self) dataset.data = dataset.data.trimmed_after_convolution_from( diff --git a/autoarray/dataset/grids.py b/autoarray/dataset/grids.py index 34c388c3f..6920145a4 100644 --- a/autoarray/dataset/grids.py +++ b/autoarray/dataset/grids.py @@ -44,6 +44,9 @@ def __init__( Parameters ---------- + mask + The 2D mask defining which pixels are included in the dataset. All grids are constructed + to align with the centres of the unmasked pixels in this mask. over_sample_size_lp The over sampling scheme size, which divides the grid into a sub grid of smaller pixels when computing values (e.g. images) from the grid to approximate the 2D line integral of the amount of light that falls @@ -53,7 +56,8 @@ def __init__( passed into the calculations performed in the `inversion` module. psf The Point Spread Function kernel of the image which accounts for diffraction due to the telescope optics - via 2D convolution. + via 2D convolution. Required to compute the blurring grid; if `None` the blurring grid + is not constructed. """ self.mask = mask self.over_sample_size_lp = over_sample_size_lp @@ -67,7 +71,17 @@ def __init__( @property def lp(self): + """ + The light-profile grid: a `Grid2D` of (y,x) Cartesian coordinates at the centre of every + unmasked image pixel, used for evaluating light profiles and other spatial calculations + during model fitting. + + The grid uses `over_sample_size_lp` to perform over-sampled sub-pixel integration, + approximating the 2D line integral of the light profile within each pixel. This grid is + what most model-fitting calculations use (e.g. computing galaxy images). + This property is lazily evaluated and cached on first access. + """ if self._lp is not None: return self._lp @@ -80,6 +94,16 @@ def lp(self): @property def pixelization(self): + """ + The pixelization grid: a `Grid2D` of (y,x) Cartesian coordinates at the centre of every + unmasked image pixel, dedicated to pixelized source reconstructions via the `inversion` module. + + This grid uses `over_sample_size_pixelization` which can differ from `over_sample_size_lp`, + allowing the pixelization to benefit from a different (e.g. lower) over-sampling resolution + than the light-profile grid. + + This property is lazily evaluated and cached on first access. + """ if self._pixelization is not None: return self._pixelization @@ -92,7 +116,18 @@ def pixelization(self): @property def blurring(self): + """ + The blurring grid: a `Grid2D` of (y,x) coordinates for pixels that lie just outside the + mask but whose light can be scattered into the unmasked region by the PSF. + + When convolving a model image with the PSF, pixels neighbouring the mask boundary can + contribute flux to unmasked pixels. The blurring grid provides the coordinates of these + border pixels so their light profile values can be evaluated and included in the convolution. + + Returns `None` if no PSF was supplied (i.e. no blurring is performed). + This property is lazily evaluated and cached on first access. + """ if self._blurring is not None: return self._blurring @@ -113,7 +148,16 @@ def blurring(self): @property def border_relocator(self) -> BorderRelocator: + """ + The border relocator for the pixelization grid. + + During pixelized source reconstruction, source-plane coordinates that map outside the + border of the pixelization mesh can cause numerical problems. The `BorderRelocator` + detects these coordinates and relocates them to the border of the mesh, preventing + ill-conditioned inversions. + This property is lazily evaluated and cached on first access. + """ if self._border_relocator is not None: return self._border_relocator @@ -133,6 +177,26 @@ def __init__( blurring=None, border_relocator=None, ): + """ + A lightweight plain-data container for pre-constructed dataset grids. + + Unlike `GridsDataset`, this class performs no computation — it simply holds grids that have + already been created elsewhere. It is used in test fixtures and mock datasets where a full + `GridsDataset` is not needed, but code that accesses `dataset.grids.lp` or + `dataset.grids.pixelization` still needs to work. + + Parameters + ---------- + lp + The light-profile `Grid2D` used for evaluating light profiles during model fitting. + pixelization + The pixelization `Grid2D` used for source reconstruction via the inversion module. + blurring + The blurring `Grid2D` for pixels outside the mask that contribute flux via PSF convolution. + border_relocator + The `BorderRelocator` used to remap out-of-bounds source-plane coordinates to the + pixelization mesh border. + """ self.lp = lp self.pixelization = pixelization self.blurring = blurring diff --git a/autoarray/dataset/imaging/dataset.py b/autoarray/dataset/imaging/dataset.py index c4ce96c83..82e60dfa8 100644 --- a/autoarray/dataset/imaging/dataset.py +++ b/autoarray/dataset/imaging/dataset.py @@ -68,6 +68,10 @@ def __init__( psf The Point Spread Function kernel of the image which accounts for diffraction due to the telescope optics via 2D convolution. + psf_setup_state + If `True`, a `ConvolverState` is precomputed from the PSF kernel and mask, storing the + convolution pair indices required for efficient 2D convolution. This is set automatically + to `True` when a mask is applied via `apply_mask()` and should not normally be set by hand. noise_covariance_matrix A noise-map covariance matrix representing the covariance between noise in every `data` value, which can be used via a bespoke fit to account for correlated noise in the data. @@ -236,13 +240,25 @@ def apply_mask(self, mask: Mask2D) -> "Imaging": Apply a mask to the imaging dataset, whereby the mask is applied to the image data, noise-map and other quantities one-by-one. - The `apply_mask` function cannot be called multiple times, if it is a mask may remove data, therefore - an exception is raised. If you wish to apply a new mask, reload the dataset from .fits files. + The mask is applied to the `data`, `noise_map`, `over_sample_size_lp` and + `over_sample_size_pixelization` arrays. If a `noise_covariance_matrix` is present, the rows + and columns corresponding to masked pixels are removed so it stays consistent with the + remaining unmasked pixels. The PSF `ConvolverState` is recomputed for the new mask. + + The `apply_mask` function cannot be called multiple times — a new mask cannot expand the + unmasked region beyond what was already unmasked, as the underlying data has already been + trimmed. An exception is raised if this is attempted. If you wish to apply a different mask, + reload the dataset from .fits files. Parameters ---------- mask The 2D mask that is applied to the image. + + Returns + ------- + Imaging + A new `Imaging` dataset with the mask applied to all arrays. """ invalid = np.logical_and(self.data.mask, np.logical_not(mask)) @@ -413,22 +429,27 @@ def apply_sparse_operator( batch_size: int = 128, ): """ + Precompute the PSF precision operator for efficient pixelized source reconstruction. + The sparse linear algebra formalism precomputes the convolution of every pair of masked - noise-map values given the PSF (see `inversion.inversion_util`). + noise-map values given the PSF (see `inversion.inversion_util`). This is the imaging + equivalent of the interferometer NUFFT precision matrix. - The `WTilde` object stores these precomputed values in the imaging dataset ensuring they are only computed once - per analysis. + The `ImagingSparseOperator` stores these precomputed values in the imaging dataset ensuring + they are only computed once per analysis, enabling fast repeated likelihood evaluations during + model fitting. - This uses lazy allocation such that the calculation is only performed when the wtilde matrices are used, - ensuring efficient set up of the `Imaging` class. + Parameters + ---------- + batch_size + The number of image pixels processed per batch when computing the sparse operator via + FFT-based convolution. Reducing this lowers peak memory usage at the cost of speed. Returns ------- - batch_size - The size of batches used to compute the w-tilde curvature matrix via FFT-based convolution, - which can be reduced to produce lower memory usage at the cost of speed - use_jax - Whether to use JAX to compute W-Tilde. This requires JAX to be installed. + Imaging + A new `Imaging` dataset with the precomputed `ImagingSparseOperator` attached, enabling + efficient pixelized source reconstruction via the sparse linear algebra formalism. """ logger.info( @@ -459,22 +480,20 @@ def apply_sparse_operator_cpu( self, ): """ - The sparse linear algebra formalism precomputes the convolution of every pair of masked - noise-map values given the PSF (see `inversion.inversion_util`). + Precompute the PSF precision operator using a CPU-only Numba implementation. - The `WTilde` object stores these precomputed values in the imaging dataset ensuring they are only computed once - per analysis. + This is the CPU alternative to `apply_sparse_operator()`, using Numba JIT compilation + for the convolution loop rather than JAX. It requires `numba` to be installed; an + `InversionException` is raised if it is not available. - This uses lazy allocation such that the calculation is only performed when the wtilde matrices are used, - ensuring efficient set up of the `Imaging` class. + The resulting `SparseLinAlgImagingNumba` operator is stored on the returned `Imaging` + dataset and used by `FitImaging` when performing pixelized source reconstructions. Returns ------- - batch_size - The size of batches used to compute the w-tilde curvature matrix via FFT-based convolution, - which can be reduced to produce lower memory usage at the cost of speed. - use_jax - Whether to use JAX to compute W-Tilde. This requires JAX to be installed. + Imaging + A new `Imaging` dataset with a precomputed Numba-based sparse operator attached, + enabling efficient pixelized source reconstruction on CPU hardware. """ try: import numba diff --git a/autoarray/dataset/imaging/simulator.py b/autoarray/dataset/imaging/simulator.py index 48b73b6b7..0f32b9add 100644 --- a/autoarray/dataset/imaging/simulator.py +++ b/autoarray/dataset/imaging/simulator.py @@ -33,8 +33,8 @@ def __init__( The simulation of an `Imaging` dataset uses the following steps: - 1) Receive as input a raw image of what the data looks like before any simulaiton process is applied. - 2) Include dirrection due to the telescope optics by convolve the image with an input Point Spread + 1) Receive as input a raw image of what the data looks like before any simulation process is applied. + 2) Include diffraction due to the telescope optics by convolving the image with an input Point Spread Function (PSF). 3) Use input values of the background sky level in every pixel of the image to add the background sky to the PSF convolved image. @@ -121,7 +121,18 @@ def via_image_from( Parameters ---------- image - The 2D image from which the Imaging dataset is simulated. + The 2D image from which the Imaging dataset is simulated (e.g. a model galaxy image + before any telescope effects are applied). Must be an `Array2D`. + over_sample_size + If provided, the returned dataset has its over-sampling updated via `apply_over_sampling`. + Should be an `Array2D` of integer sub-grid sizes with the same shape as the image. + xp + The array module to use for PSF convolution (default `np` for NumPy, or `jnp` for JAX). + + Returns + ------- + Imaging + The simulated imaging dataset with PSF convolution, noise and background sky applied. """ exposure_time_map = Array2D.full( diff --git a/autoarray/dataset/interferometer/dataset.py b/autoarray/dataset/interferometer/dataset.py index 7c52585e5..446f5b72d 100644 --- a/autoarray/dataset/interferometer/dataset.py +++ b/autoarray/dataset/interferometer/dataset.py @@ -79,6 +79,14 @@ def __init__( transformer_class The class of the Fourier Transform which maps images from real space to Fourier space visibilities and the uv-plane. + sparse_operator + A precomputed `InterferometerSparseOperator` containing the NUFFT precision matrix for efficient + pixelized source reconstruction. This is computed via `apply_sparse_operator()` and can be passed + here directly to avoid recomputing it (e.g. when loading a cached result from disk). + raise_error_dft_visibilities_limit + If `True`, an exception is raised if the dataset has more than 10,000 visibilities and + `transformer_class=TransformerDFT`. The DFT is too slow for large datasets and `TransformerNUFFT` + should be used instead. Set to `False` to suppress this check. """ self.real_space_mask = real_space_mask @@ -133,11 +141,46 @@ def from_fits( transformer_class=TransformerNUFFT, ): """ - Factory for loading the interferometer data_type from .fits files, as well as computing properties like the - noise-map, exposure-time map, etc. from the interferometer-data_type. + Load an interferometer dataset from multiple .fits files. - This factory also includes a number of routines for converting the interferometer-data_type from unit_label - not supported by PyAutoLens (e.g. adus, electrons) to electrons per second. + The visibilities (complex-valued Fourier-space data), noise map and uv_wavelengths (baseline + coordinates) are each loaded from separate .fits files. A real-space mask defining the sky + region used for Fourier transforms must be supplied separately. + + The visibilities are assumed to be stored as a 2D array of shape (total_visibilities, 2) where + column 0 is the real component and column 1 is the imaginary component. The noise map has the + same shape. The uv_wavelengths are a 2D array of shape (total_visibilities, 2) with columns + corresponding to the (u, v) baseline coordinates in units of wavelengths. + + Parameters + ---------- + data_path + The path to the .fits file containing the visibility data + (e.g. '/path/to/visibilities.fits'). + noise_map_path + The path to the .fits file containing the visibility noise map + (e.g. '/path/to/noise_map.fits'). + uv_wavelengths_path + The path to the .fits file containing the (u, v) baseline coordinates in units of + wavelengths (e.g. '/path/to/uv_wavelengths.fits'). + real_space_mask + A `Mask2D` defining the real-space region of the sky that contains signal. This mask + determines the pixel grid used by the Fourier transformer and the coordinate grids + associated with the dataset. + visibilities_hdu + The HDU index in the visibilities .fits file from which data is loaded. + noise_map_hdu + The HDU index in the noise map .fits file from which data is loaded. + uv_wavelengths_hdu + The HDU index in the uv_wavelengths .fits file from which data is loaded. + transformer_class + The class of the Fourier Transform which maps images from real space to Fourier space + visibilities. Defaults to `TransformerNUFFT` for efficiency with large datasets. + + Returns + ------- + Interferometer + The interferometer dataset loaded from the .fits files. """ visibilities = Visibilities.from_fits(file_path=data_path, hdu=visibilities_hdu) @@ -168,28 +211,44 @@ def apply_sparse_operator( use_jax: bool = False, ): """ - The sparse linear algebra equations precomputes the Fourier Transform of all the visibilities - given the `uv_wavelengths` (see `inversion.inversion_util`). + Precompute the NUFFT precision operator for efficient pixelized source reconstruction. + + The sparse linear algebra formalism precomputes the Fourier Transform response matrix for all + visibility baselines, enabling fast repeated evaluation during model fitting. This avoids + recomputing the full NUFFT on every likelihood call. - The `WTilde` object stores these precomputed values in the interferometer dataset ensuring they are only - computed once per analysis. + The resulting `InterferometerSparseOperator` is stored on the returned `Interferometer` dataset + and is used automatically by `FitInterferometer` when performing pixelized reconstructions via + the inversion module. - This uses lazy allocation such that the calculation is only performed when the wtilde matrices are used, - ensuring efficient set up of the `Interferometer` class. + Computing the NUFFT precision matrix from scratch can be very slow (runtime scales with both + the number of visibilities and the real-space mask resolution — potentially hours for large + datasets). The result can be cached to disk and reloaded to avoid recomputation. Parameters ---------- nufft_precision_operator - An already computed curvature preload matrix for this dataset (e.g. loaded from hard-disk), to prevent - long recalculations of this matrix for large datasets. + An already computed NUFFT precision matrix for this dataset (e.g. loaded from disk via + `np.load`) to avoid an expensive recomputation. If `None` it is computed from scratch + by calling `psf_precision_operator_from()`. batch_size - The size of batches used to compute the w-tilde curvature matrix via FFT-based convolution, - which can be reduced to produce lower memory usage at the cost of speed. + The number of real-space pixels processed per batch when building the sparse operator. + Reducing this lowers peak memory usage at the cost of speed. + chunk_k + The number of visibilities processed per chunk when computing the NUFFT precision matrix + inside `psf_precision_operator_from()`. Reducing this lowers peak memory usage. + show_progress + If `True`, a progress bar is displayed while computing the NUFFT precision matrix. + show_memory + If `True`, memory usage statistics are printed while computing the NUFFT precision matrix. + use_jax + If `True`, JAX is used to accelerate the NUFFT precision matrix computation. Returns ------- - WTildeInterferometer - Precomputed values used for the w tilde formalism of linear algebra calculations. + Interferometer + A new `Interferometer` dataset with the precomputed `InterferometerSparseOperator` attached, + enabling efficient pixelized source reconstruction via the sparse linear algebra formalism. """ if nufft_precision_operator is None: @@ -233,6 +292,37 @@ def psf_precision_operator_from( show_memory: bool = False, use_jax: bool = False, ): + """ + Compute the NUFFT precision matrix for this interferometer dataset. + + The precision matrix encodes the response of every real-space pixel to every visibility + baseline, weighted by the noise map. It is the core precomputed quantity required for + efficient pixelized source reconstruction via the sparse linear algebra formalism. + + This computation can be very slow for large datasets (runtime scales with the number of + visibilities multiplied by the number of unmasked real-space pixels). For datasets with + tens of thousands of visibilities and high-resolution masks, computation can take hours + on a CPU. The result should be saved to disk and reloaded rather than recomputed on each + run. Use `apply_sparse_operator(nufft_precision_operator=...)` to attach a cached result. + + Parameters + ---------- + chunk_k + The number of visibilities processed per chunk. Reducing this lowers peak memory + usage during computation at the cost of speed. + show_progress + If `True`, a progress bar is shown during computation. + show_memory + If `True`, memory usage statistics are printed during computation. + use_jax + If `True`, JAX is used to accelerate the computation. + + Returns + ------- + np.ndarray + The NUFFT precision matrix of shape (total_pixels, total_pixels) where total_pixels + is the number of unmasked real-space pixels. + """ return inversion_interferometer_util.nufft_precision_operator_from( noise_map_real=self.noise_map.array.real, uv_wavelengths=self.uv_wavelengths, @@ -246,36 +336,81 @@ def psf_precision_operator_from( @property def mask(self): + """ + The real-space mask of the interferometer dataset. + + For an interferometer, this is the `real_space_mask` which defines the region of sky that + contains signal. It is used as the spatial domain for the Fourier transform, determining + the pixel grid size and coordinate grids. + """ return self.real_space_mask @property def amplitudes(self): + """ + The amplitudes of the complex visibilities, defined as the absolute value of each visibility: + amplitude = sqrt(real^2 + imag^2). + """ return self.data.amplitudes @property def phases(self): + """ + The phases of the complex visibilities in radians, defined as arctan(imag / real) for + each visibility. + """ return self.data.phases @property def uv_distances(self): + """ + The radial distance of each visibility baseline from the origin of the UV-plane, in units + of wavelengths. Computed as sqrt(u^2 + v^2) for each (u, v) baseline pair. + """ return np.sqrt( np.square(self.uv_wavelengths[:, 0]) + np.square(self.uv_wavelengths[:, 1]) ) @property def dirty_image(self): + """ + The dirty image, computed as the inverse Fourier transform of the observed visibilities. + + This is the raw image obtained by back-projecting the visibilities without any deconvolution. + It provides a quick visual representation of the data but is convolved with the synthesized + beam (the Fourier transform of the UV-plane sampling function). + """ return self.transformer.image_from(visibilities=self.data) @property def dirty_noise_map(self): + """ + The dirty noise map, computed as the inverse Fourier transform of the noise map visibilities. + + Provides a real-space representation of the noise levels in the dirty image. + """ return self.transformer.image_from(visibilities=self.noise_map) @property def dirty_signal_to_noise_map(self): + """ + The dirty signal-to-noise map, computed as the inverse Fourier transform of the + complex signal-to-noise visibility map. + """ return self.transformer.image_from(visibilities=self.signal_to_noise_map) @property def signal_to_noise_map(self): + """ + The complex signal-to-noise map of the visibilities. + + Computed separately for the real and imaginary components as data / noise_map. Values + below zero are clamped to zero, as negative signal-to-noise is not physically meaningful. + + Unlike the base class implementation (which operates on real-valued data), this override + handles the complex nature of interferometric visibilities by treating the real and + imaginary parts independently. + """ signal_to_noise_map_real = np.divide( np.real(self.data.array), np.real(self.noise_map.array) ) @@ -296,6 +431,27 @@ def output_to_fits( uv_wavelengths_path=None, overwrite=False, ): + """ + Output the interferometer dataset to multiple .fits files. + + Each component (visibilities, noise map, uv_wavelengths) is saved to its own .fits file. + Any path set to `None` means that component is not saved. + + Parameters + ---------- + data_path + The path where the visibility data is saved (e.g. '/path/to/visibilities.fits'). + If `None`, the visibilities are not saved. + noise_map_path + The path where the noise map is saved (e.g. '/path/to/noise_map.fits'). + If `None`, the noise map is not saved. + uv_wavelengths_path + The path where the uv_wavelengths array is saved (e.g. '/path/to/uv_wavelengths.fits'). + If `None`, the uv_wavelengths are not saved. + overwrite + If `True`, existing .fits files are overwritten. If `False`, an exception is raised + if a file already exists at the given path. + """ if data_path is not None: self.data.output_to_fits(file_path=data_path, overwrite=overwrite) @@ -311,4 +467,11 @@ def output_to_fits( @property def psf(self): + """ + Returns `None` for interferometer datasets. + + Interferometers do not have a Point Spread Function in the same sense as imaging datasets. + The equivalent quantity is the synthesized beam, which is determined by the UV-plane coverage + and is not stored explicitly. This property exists to satisfy the `AbstractDataset` interface. + """ return None diff --git a/autoarray/dataset/interferometer/simulator.py b/autoarray/dataset/interferometer/simulator.py index 1f8ea2210..52455cddb 100644 --- a/autoarray/dataset/interferometer/simulator.py +++ b/autoarray/dataset/interferometer/simulator.py @@ -18,21 +18,47 @@ def __init__( noise_if_add_noise_false=0.1, noise_seed=-1, ): - """A class representing a Imaging observation, using the shape of the image, the pixel scale, - psf, exposure time, etc. + """ + Simulates observations of `Interferometer` data, including transforming a real-space image to + complex-valued visibilities in Fourier space and optionally adding complex Gaussian noise. + + The simulation of an `Interferometer` dataset uses the following steps: + + 1) Receive as input a real-space image (e.g. a model galaxy or lens system) and a set of UV-plane + baselines (uv_wavelengths) describing the interferometer configuration. + 2) Fourier transform the real-space image to the UV-plane using the configured transformer class + (DFT or NUFFT) to produce model visibilities. + 3) Optionally add complex Gaussian noise to the visibilities, with the noise level controlled by + `noise_sigma`. The noise is added independently to the real and imaginary parts. + 4) Create a constant noise map with value `noise_sigma` for every visibility. If noise is not + added (`noise_sigma=None`), the noise map is filled with `noise_if_add_noise_false` instead. + + The returned `Interferometer` dataset contains the simulated visibilities, noise map, + uv_wavelengths and real-space mask, and can be used directly with `FitInterferometer`. Parameters ---------- - real_space_shape_native - The shape of the observation. Note that we do not simulator a full Imaging array (e.g. 2000 x 2000 pixels for \ - Hubble imaging), but instead just a cut-out around the strong lens. - real_space_pixel_scales - The (y,x) arcsecond-to-pixel units conversion factor of every pixel. If this is input as a `float`, - it is converted to a (float, float). - psf : PSF - An arrays describing the PSF kernel of the image. - exposure_time_map - The exposure time of an observation using this data_type. + uv_wavelengths + The (u, v) baseline coordinates of the interferometer in units of wavelengths. This is a + 2D array of shape (total_visibilities, 2) where each row is a (u, v) baseline pair. These + define the Fourier-space sampling of the observation. + exposure_time + The exposure time of the simulated interferometer observation in seconds. + transformer_class + The class used to perform the Fourier transform between real space and the UV-plane. The + default `TransformerDFT` is suitable for small datasets (fewer than ~10,000 visibilities). + For larger datasets use `TransformerNUFFT` for efficiency. + noise_sigma + The standard deviation of the complex Gaussian noise added to each visibility. Noise is + added independently to the real and imaginary components. If `None`, no noise is added to + the visibilities but a noise map is still created using `noise_if_add_noise_false`. + noise_if_add_noise_false + The noise value assigned to every visibility in the noise map when `noise_sigma=None` + (i.e. when no noise is added to the data). This gives the noise map a non-zero value + so that downstream fits remain well-defined. + noise_seed + The random seed used for noise generation. A value of -1 uses a different random seed + on every run, producing different noise realisations each time. """ self.uv_wavelengths = uv_wavelengths @@ -44,23 +70,25 @@ def __init__( def via_image_from(self, image): """ - Returns a realistic simulated image by applying effects to a plain simulated image. + Simulate an `Interferometer` dataset from an input real-space image. + + The steps of the simulation process are described in the `SimulatorInterferometer` `__init__` + docstring. In brief: the image is Fourier transformed to visibilities using the configured + transformer class and the uv_wavelengths baselines, then complex Gaussian noise is optionally + added and a constant noise map is created. Parameters ---------- - real_space_image - The image before simulating (e.g. the lens and source galaxies before optics blurring and UVPlane read-out). - real_space_pixel_scales - The scale of each pixel in scaled units - exposure_time_map - An array representing the effective exposure time of each pixel. - psf: PSF - An array describing the PSF the simulated image is blurred with. - add_poisson_noise_to_data: Bool - If `True` poisson noise_maps is simulated and added to the image, based on the total counts in each image - pixel - noise_seed: int - A seed for random noise_maps generation + image + The 2D real-space image from which the interferometer dataset is simulated (e.g. the + surface brightness of a galaxy or lens system). Must be an `Array2D` with an associated + mask that defines the real-space region used for the Fourier transform. + + Returns + ------- + Interferometer + The simulated interferometer dataset containing visibilities, noise map, uv_wavelengths + and the real-space mask derived from the input image. """ transformer = self.transformer_class( diff --git a/autoarray/dataset/mock/mock_dataset.py b/autoarray/dataset/mock/mock_dataset.py index c425c6a27..ab9ee927b 100644 --- a/autoarray/dataset/mock/mock_dataset.py +++ b/autoarray/dataset/mock/mock_dataset.py @@ -3,6 +3,30 @@ class MockDataset: + """ + A lightweight mock dataset for use in tests and unit testing fixtures. + + Provides a minimal dataset-like interface with `grids`, `psf` and `mask` attributes, + without any real data arrays or numerical computation. Code that only needs to access + `dataset.grids.pixelization` or similar can use this instead of constructing a full + `Imaging` or `Interferometer` dataset. + + If no `grids` are provided, a default `GridsInterface` is constructed with a single + pixelization pixel at coordinate (1.0, 1.0) and pixel_scales=1.0. This acts as a + minimal placeholder sufficient for tests that exercise inversion or pixelization code + paths without real imaging data. + + Parameters + ---------- + grids + A `GridsInterface` (or compatible object) providing `lp`, `pixelization`, `blurring` + and `border_relocator` attributes. If `None`, a default single-pixel grid is used. + psf + An optional PSF kernel object. If `None`, no PSF is available on this mock dataset. + mask + An optional mask object. If `None`, no mask is available on this mock dataset. + """ + def __init__(self, grids=None, psf=None, mask=None): if grids is None: self.grids = GridsInterface( diff --git a/autoarray/dataset/preprocess.py b/autoarray/dataset/preprocess.py index d3c0de84a..52a07512c 100644 --- a/autoarray/dataset/preprocess.py +++ b/autoarray/dataset/preprocess.py @@ -7,7 +7,7 @@ def array_with_new_shape(array, new_shape): """Resize an input array around its centre to a new input shape. If a new_shape dimension is smaller than the array's current dimension, the data at the edges is trimmed and - removedd. If it is larger, the data is padded with zeros. + removed. If it is larger, the data is padded with zeros. If the array has even sized dimensions, the central pixel around which data is trimmed / padded is chosen as the top-left pixel of the central quadrant of pixels. @@ -44,7 +44,7 @@ def array_eps_to_counts(array_eps, exposure_time_map): def array_counts_to_eps(array_counts, exposure_time_map): """ - Convert an array in units of electrons per second to counts, using an exposure time map containing the exposure + Convert an array in units of counts to electrons per second, using an exposure time map containing the exposure time at every point in the array. The conversion from counts to electrons per second is: