From e2ddf0f2e7ad565e6eae90cbfee7a4634c868ae7 Mon Sep 17 00:00:00 2001 From: James Nightingale Date: Wed, 2 Apr 2025 19:14:15 +0100 Subject: [PATCH 01/11] removal of convolver and switch to psf where required --- autoarray/__init__.py | 2 - autoarray/dataset/imaging/dataset.py | 34 +- autoarray/dataset/interferometer/dataset.py | 2 +- autoarray/dataset/preprocess.py | 2 +- autoarray/fixtures.py | 6 - .../inversion/inversion/dataset_interface.py | 6 +- autoarray/inversion/inversion/factory.py | 3 +- .../inversion/inversion/imaging/abstract.py | 15 +- .../inversion/inversion/imaging/mapping.py | 10 +- .../inversion/inversion/imaging/w_tilde.py | 9 +- autoarray/inversion/linear_obj/linear_obj.py | 2 +- .../inversion/mock/mock_inversion_imaging.py | 8 +- autoarray/mask/derive/mask_2d.py | 4 +- autoarray/mock.py | 1 - autoarray/operators/convolver.py | 592 --------- autoarray/operators/mock/mock_convolver.py | 2 +- autoarray/structures/arrays/kernel_2d.py | 25 +- test_autoarray/conftest.py | 5 - .../dataset/imaging/test_dataset.py | 6 +- .../inversion/imaging/test_imaging.py | 20 +- .../imaging/test_inversion_imaging_util.py | 12 +- test_autoarray/operators/test_convolver.py | 1164 ----------------- .../structures/arrays/test_kernel_2d.py | 334 ++++- 23 files changed, 374 insertions(+), 1890 deletions(-) delete mode 100644 autoarray/operators/convolver.py delete mode 100644 test_autoarray/operators/test_convolver.py diff --git a/autoarray/__init__.py b/autoarray/__init__.py index a212002f5..1a6d79c53 100644 --- a/autoarray/__init__.py +++ b/autoarray/__init__.py @@ -54,8 +54,6 @@ from .mask.derive.grid_2d import DeriveGrid2D from .mask.mask_1d import Mask1D from .mask.mask_2d import Mask2D -from .operators.convolver import Convolver -from .operators.convolver import Convolver from .operators.transformer import TransformerDFT from .operators.transformer import TransformerNUFFT from .operators.over_sampling.decorator import over_sample diff --git a/autoarray/dataset/imaging/dataset.py b/autoarray/dataset/imaging/dataset.py index 20bb8b0be..7fa3c515f 100644 --- a/autoarray/dataset/imaging/dataset.py +++ b/autoarray/dataset/imaging/dataset.py @@ -9,7 +9,6 @@ from autoarray.dataset.grids import GridsDataset from autoarray.dataset.imaging.w_tilde import WTildeImaging from autoarray.structures.arrays.uniform_2d import Array2D -from autoarray.operators.convolver import Convolver from autoarray.structures.arrays.kernel_2d import Kernel2D from autoarray.mask.mask_2d import Mask2D from autoarray import type as ty @@ -30,7 +29,7 @@ def __init__( noise_covariance_matrix: Optional[np.ndarray] = None, over_sample_size_lp: Union[int, Array2D] = 4, over_sample_size_pixelization: Union[int, Array2D] = 4, - pad_for_convolver: bool = False, + pad_for_psf: bool = False, use_normalized_psf: Optional[bool] = True, check_noise_map: bool = True, ): @@ -77,7 +76,7 @@ def __init__( over_sample_size_pixelization How over sampling is performed for the grid which is associated with a pixelization, which is therefore passed into the calculations performed in the `inversion` module. - pad_for_convolver + pad_for_psf The PSF convolution may extend beyond the edges of the image mask, which can lead to edge effects in the convolved image. If `True`, the image and noise-map are padded to ensure the PSF convolution does not extend beyond the edge of the image. @@ -90,9 +89,9 @@ def __init__( self.unmasked = None - self.pad_for_convolver = pad_for_convolver + self.pad_for_psf = pad_for_psf - if pad_for_convolver and psf is not None: + if pad_for_psf and psf is not None: try: data.mask.derive_mask.blurring_from( kernel_shape_native=psf.shape_native @@ -176,25 +175,6 @@ def grids(self): psf=self.psf, ) - @cached_property - def convolver(self): - """ - Returns a `Convolver` from a mask and 2D PSF kernel. - - The `Convolver` stores in memory the array indexing between the mask and PSF, enabling efficient 2D PSF - convolution of images and matrices used for linear algebra calculations (see `operators.convolver`). - - This uses lazy allocation such that the calculation is only performed when the convolver is used, ensuring - efficient set up of the `Imaging` class. - - Returns - ------- - Convolver - The convolver given the masked imaging data's mask and PSF. - """ - - return Convolver(mask=self.mask, kernel=Kernel2D(values=self.psf._array, mask=self.psf.mask, header=self.psf.header)) - @cached_property def w_tilde(self): """ @@ -370,7 +350,7 @@ def apply_mask(self, mask: Mask2D) -> "Imaging": noise_covariance_matrix=noise_covariance_matrix, over_sample_size_lp=over_sample_size_lp, over_sample_size_pixelization=over_sample_size_pixelization, - pad_for_convolver=True, + pad_for_psf=True, ) dataset.unmasked = unmasked_dataset @@ -463,7 +443,7 @@ def apply_noise_scaling( noise_covariance_matrix=self.noise_covariance_matrix, over_sample_size_lp=self.over_sample_size_lp, over_sample_size_pixelization=self.over_sample_size_pixelization, - pad_for_convolver=False, + pad_for_psf=False, check_noise_map=False, ) @@ -511,7 +491,7 @@ def apply_over_sampling( over_sample_size_lp=over_sample_size_lp or self.over_sample_size_lp, over_sample_size_pixelization=over_sample_size_pixelization or self.over_sample_size_pixelization, - pad_for_convolver=False, + pad_for_psf=False, check_noise_map=False, ) diff --git a/autoarray/dataset/interferometer/dataset.py b/autoarray/dataset/interferometer/dataset.py index 21dd600a6..06892ea68 100644 --- a/autoarray/dataset/interferometer/dataset.py +++ b/autoarray/dataset/interferometer/dataset.py @@ -276,5 +276,5 @@ def output_to_fits( ) @property - def convolver(self): + def psf(self): return None diff --git a/autoarray/dataset/preprocess.py b/autoarray/dataset/preprocess.py index 046ab3ea8..5c7338204 100644 --- a/autoarray/dataset/preprocess.py +++ b/autoarray/dataset/preprocess.py @@ -328,7 +328,7 @@ def background_noise_map_via_edges_from(image, no_edges): def psf_with_odd_dimensions_from(psf): """ If the PSF kernel has one or two even-sized dimensions, return a PSF object where the kernel has odd-sized - dimensions (odd-sized dimensions are required by a *Convolver*). + dimensions (odd-sized dimensions are required for 2D convolution). Kernels are rescaled using the scikit-image routine rescale, which performs rescaling via an interpolation routine. This may lead to loss of accuracy in the PSF kernel and it is advised that users, where possible, diff --git a/autoarray/fixtures.py b/autoarray/fixtures.py index acfd277df..68a8a73f8 100644 --- a/autoarray/fixtures.py +++ b/autoarray/fixtures.py @@ -110,12 +110,6 @@ def make_blurring_grid_2d_7x7(): return aa.Grid2D.from_mask(mask=make_blurring_mask_2d_7x7()) -# CONVOLVERS # - - -def make_convolver_7x7(): - return aa.Convolver(mask=make_mask_2d_7x7(), kernel=make_psf_3x3()) - def make_image_7x7(): return aa.Array2D.ones(shape_native=(7, 7), pixel_scales=(1.0, 1.0)) diff --git a/autoarray/inversion/inversion/dataset_interface.py b/autoarray/inversion/inversion/dataset_interface.py index abc780411..0e417e71a 100644 --- a/autoarray/inversion/inversion/dataset_interface.py +++ b/autoarray/inversion/inversion/dataset_interface.py @@ -4,7 +4,7 @@ def __init__( data, noise_map, grids=None, - convolver=None, + psf=None, transformer=None, w_tilde=None, noise_covariance_matrix=None, @@ -41,7 +41,7 @@ def __init__( border_relocator The border relocator, which relocates coordinates outside the border of the source-plane data grid to its edge. - convolver + psf Perform 2D convolution of the imaigng data's PSF when computing the operated mapping matrix. transformer Performs a Fourier transform of the image-data from real-space to visibilities when computing the @@ -59,7 +59,7 @@ def __init__( self.data = data self.noise_map = noise_map self.grids = grids - self.convolver = convolver + self.psf = psf self.transformer = transformer self.w_tilde = w_tilde self.noise_covariance_matrix = noise_covariance_matrix diff --git a/autoarray/inversion/inversion/factory.py b/autoarray/inversion/inversion/factory.py index 31c69d57b..327262786 100644 --- a/autoarray/inversion/inversion/factory.py +++ b/autoarray/inversion/inversion/factory.py @@ -172,8 +172,7 @@ def inversion_interferometer_from( dataset The dataset (e.g. `Interferometer`) whose data is reconstructed via the `Inversion`. w_tilde - Object which uses the `Imaging` dataset's PSF / `Convolver` operateor to perform the `Inversion` using the - w-tilde formalism. + Object which uses the `Imaging` dataset's PSF to perform the `Inversion` using the w-tilde formalism. linear_obj_list The list of linear objects (e.g. analytic functions, a mapper with a pixelized grid) which reconstruct the input dataset's data and whose values are solved for via the inversion. diff --git a/autoarray/inversion/inversion/imaging/abstract.py b/autoarray/inversion/inversion/imaging/abstract.py index 1ea826f88..5efc4d0a9 100644 --- a/autoarray/inversion/inversion/imaging/abstract.py +++ b/autoarray/inversion/inversion/imaging/abstract.py @@ -32,8 +32,7 @@ def __init__( of the linear object parameters that best reconstruct the dataset to be solved, via linear matrix algebra. This object contains matrices and vectors which perform an inversion for fits to an `Imaging` dataset. This - includes operations which use a PSF / `Convolver` in order to incorporate blurring into the solved for - linear object pixels. + includes operations which use a PSF in order to incorporate blurring into the solved for linear object pixels. The inversion may be regularized, whereby the parameters of the linear objects used to reconstruct the data are smoothed with one another such that their solved for values conform to certain properties (e.g. smoothness @@ -76,8 +75,8 @@ def __init__( ) @property - def convolver(self): - return self.dataset.convolver + def psf(self): + return self.dataset.psf @property def operated_mapping_matrix_list(self) -> List[np.ndarray]: @@ -88,7 +87,7 @@ def operated_mapping_matrix_list(self) -> List[np.ndarray]: This is used to construct the simultaneous linear equations which reconstruct the data. This property returns the a list of each linear object's blurred mapping matrix, which is computed by - blurring each linear object's `mapping_matrix` property with the `Convolver` operator. + blurring each linear object's `mapping_matrix` property with the `psf` operator. A linear object may have a `operated_mapping_matrix_override` property, which bypasses the `mapping_matrix` computation and convolution operator and is directly placed in the `operated_mapping_matrix_list`. @@ -96,7 +95,7 @@ def operated_mapping_matrix_list(self) -> List[np.ndarray]: return [ ( - self.convolver.convolve_mapping_matrix( + self.psf.convolve_mapping_matrix( mapping_matrix=linear_obj.mapping_matrix ) if linear_obj.operated_mapping_matrix_override is None @@ -139,7 +138,7 @@ def linear_func_operated_mapping_matrix_dict(self) -> Dict: if linear_func.operated_mapping_matrix_override is not None: operated_mapping_matrix = linear_func.operated_mapping_matrix_override else: - operated_mapping_matrix = self.convolver.convolve_mapping_matrix( + operated_mapping_matrix = self.psf.convolve_mapping_matrix( mapping_matrix=linear_func.mapping_matrix ) @@ -221,7 +220,7 @@ def mapper_operated_mapping_matrix_dict(self) -> Dict: mapper_operated_mapping_matrix_dict = {} for mapper in self.cls_list_from(cls=AbstractMapper): - operated_mapping_matrix = self.convolver.convolve_mapping_matrix( + operated_mapping_matrix = self.psf.convolve_mapping_matrix( mapping_matrix=mapper.mapping_matrix ) diff --git a/autoarray/inversion/inversion/imaging/mapping.py b/autoarray/inversion/inversion/imaging/mapping.py index 5b3e40966..03d73ff63 100644 --- a/autoarray/inversion/inversion/imaging/mapping.py +++ b/autoarray/inversion/inversion/imaging/mapping.py @@ -39,10 +39,8 @@ def __init__( Parameters ---------- - noise_map - The noise-map of the observed imaging data which values are solved for. - convolver - The convolver which performs a 2D convolution on the mapping matrix with the imaging data's PSF. + dataset + The dataset containing the image data, noise-map and psf which is fitted by the inversion. linear_obj_list The linear objects used to reconstruct the data's observed values. If multiple linear objects are passed the simultaneous linear equations are combined and solved simultaneously. @@ -80,7 +78,7 @@ def _data_vector_mapper(self) -> np.ndarray: mapper = mapper_list[i] param_range = mapper_param_range_list[i] - operated_mapping_matrix = self.convolver.convolve_mapping_matrix( + operated_mapping_matrix = self.psf.convolve_mapping_matrix( mapping_matrix=mapper.mapping_matrix ) @@ -142,7 +140,7 @@ def _curvature_matrix_mapper_diag(self) -> Optional[np.ndarray]: mapper_i = mapper_list[i] mapper_param_range_i = mapper_param_range_list[i] - operated_mapping_matrix = self.convolver.convolve_mapping_matrix( + operated_mapping_matrix = self.psf.convolve_mapping_matrix( mapping_matrix=mapper_i.mapping_matrix ) diff --git a/autoarray/inversion/inversion/imaging/w_tilde.py b/autoarray/inversion/inversion/imaging/w_tilde.py index 5f791873c..c249f4ce9 100644 --- a/autoarray/inversion/inversion/imaging/w_tilde.py +++ b/autoarray/inversion/inversion/imaging/w_tilde.py @@ -42,11 +42,8 @@ def __init__( Parameters ---------- - noise_map - The noise-map of the observed imaging data which values are solved for. - convolver - The convolver used to perform 2D convolution of the imaigng data's PSF when computing the operated - mapping matrix. + dataset + The dataset containing the image data, noise-map and psf which is fitted by the inversion. w_tilde An object containing matrices that construct the linear equations via the w-tilde formalism which bypasses the mapping matrix. @@ -76,7 +73,7 @@ def w_tilde_data(self): return inversion_imaging_util.w_tilde_data_imaging_from( image_native=np.array(self.data.native), noise_map_native=np.array(self.noise_map.native), - kernel_native=np.array(self.convolver.kernel.native), + kernel_native=np.array(self.psf.kernel.native), native_index_for_slim_index=self.data.mask.derive_indexes.native_for_slim, ) diff --git a/autoarray/inversion/linear_obj/linear_obj.py b/autoarray/inversion/linear_obj/linear_obj.py index 3bcc1ae57..1e0bacc85 100644 --- a/autoarray/inversion/linear_obj/linear_obj.py +++ b/autoarray/inversion/linear_obj/linear_obj.py @@ -124,7 +124,7 @@ def pixel_signals_from(self, signal_scale) -> np.ndarray: def operated_mapping_matrix_override(self) -> Optional[np.ndarray]: """ An `Inversion` takes the `mapping_matrix` of each linear object and combines it with the data's operators - (e.g. a `Convolver` for `Imaging` data) to compute the `operated_mapping_matrix`. + (e.g. a PSF for `Imaging` data) to compute the `operated_mapping_matrix`. If this property is overwritten this operation is not performed, with the `operated_mapping_matrix` output by this property automatically used instead. diff --git a/autoarray/inversion/mock/mock_inversion_imaging.py b/autoarray/inversion/mock/mock_inversion_imaging.py index 64076db56..de7f2baa1 100644 --- a/autoarray/inversion/mock/mock_inversion_imaging.py +++ b/autoarray/inversion/mock/mock_inversion_imaging.py @@ -12,7 +12,7 @@ def __init__( self, data=None, noise_map=None, - convolver=None, + psf=None, linear_obj_list=None, operated_mapping_matrix=None, linear_func_operated_mapping_matrix_dict=None, @@ -22,7 +22,7 @@ def __init__( dataset = DatasetInterface( data=data, noise_map=noise_map, - convolver=convolver, + psf=psf, ) super().__init__( @@ -70,7 +70,7 @@ def __init__( self, data=None, noise_map=None, - convolver=None, + psf=None, w_tilde=None, linear_obj_list=None, curvature_matrix_mapper_diag=None, @@ -79,7 +79,7 @@ def __init__( dataset = DatasetInterface( data=data, noise_map=noise_map, - convolver=convolver, + psf=psf, ) super().__init__( diff --git a/autoarray/mask/derive/mask_2d.py b/autoarray/mask/derive/mask_2d.py index 4d3090a52..9332e273d 100644 --- a/autoarray/mask/derive/mask_2d.py +++ b/autoarray/mask/derive/mask_2d.py @@ -146,8 +146,8 @@ def blurring_from(self, kernel_shape_native: Tuple[int, int]) -> Mask2D: Returns a blurring ``Mask2D``, representing all masked pixels (given by ``True``) whose values are blurred into unmasked pixels (given by ``False``) when a 2D convolution is performed. - This mask is used by the ``Convolver2D`` object to ensure that 2D convolution can be performed on masked - data structures without missing values. + This mask is used by the PSF to ensure that 2D convolution can be performed on masked data structures without + missing values. For example, for the following ``Mask2D``: diff --git a/autoarray/mock.py b/autoarray/mock.py index 78a857e7e..2261ba5e2 100644 --- a/autoarray/mock.py +++ b/autoarray/mock.py @@ -15,7 +15,6 @@ from autoarray.fit.mock.mock_fit_imaging import MockFitImaging from autoarray.fit.mock.mock_fit_interferometer import MockFitInterferometer from autoarray.mask.mock.mock_mask import MockMask -from autoarray.operators.mock.mock_convolver import MockConvolver from autoarray.structures.mock.mock_grid import MockGrid2DMesh from autoarray.structures.mock.mock_grid import MockMeshGrid from autoarray.structures.mock.mock_decorators import MockGridRadialMinimum diff --git a/autoarray/operators/convolver.py b/autoarray/operators/convolver.py deleted file mode 100644 index 73d3767d5..000000000 --- a/autoarray/operators/convolver.py +++ /dev/null @@ -1,592 +0,0 @@ -from autoarray import numba_util -import numpy as np - -from autoarray.structures.arrays.uniform_2d import Array2D - -from autoarray import exc -from autoarray.mask import mask_2d_util - -from os import environ - -use_jax = environ.get("USE_JAX", "0") == "1" - -if use_jax: - import jax - import jax.numpy as jnp - - -class Convolver: - def __init__(self, mask, kernel): - """ - Class to setup the 1D convolution of an / mapping matrix. - - Take a simple 3x3 and masks: - - [[2, 8, 2], - [5, 7, 5], - [3, 1, 4]] - - [[True, False, True], (True means that the value is masked) - [False, False, False], - [True, False, True]] - - A set of values in a corresponding 1d array of this might be represented as: - - [2, 8, 2, 5, 7, 5, 3, 1, 4] - - and after masking as: - - [8, 5, 7, 5, 1] - - Setup is required to perform 2D real-space convolution on the masked array. This module finds the \ - relationship between the unmasked 2D data, masked data and kernel, so that 2D real-space convolutions \ - can be efficiently applied to reduced 1D masked structures. - - This calculation also accounts for the blurring of light outside of the masked regions which blurs into \ - the masked region. - - - **IMAGE FRAMES:** - - For a masked in 2D, one can compute for every pixel all of the unmasked pixels it will blur light into for \ - a given PSF kernel size, e.g.: - - IxIxIxIxIxIxIxIxIxIxI - IxIxIxIxIxIxIxIxIxIxI This is an imaging.Mask2D, where: - IxIxIxIxIxIxIxIxIxIxI - IxIxIxIxIxIxIxIxIxIxI x = `True` (Pixel is masked and excluded from lens) - IxIxIxIoIoIoIxIxIxIxI o = `False` (Pixel is not masked and included in lens) - IxIxIxIoIoIoIxIxIxIxI - IxIxIxIoIoIoIxIxIxIxI - IxIxIxIxIxIxIxIxIxIxI - IxIxIxIxIxIxIxIxIxIxI - IxIxIxIxIxIxIxIxIxIxI - - Here, there are 9 unmasked pixels. Indexing of each unmasked pixel goes from the top-left corner right and \ - downwards, therefore: - - IxIxIxIxIxIxIxIxIxIxI - IxIxIxIxIxIxIxIxIxIxI - IxIxIxIxIxIxIxIxIxIxI - IxIxIxIxIxIxIxIxIxIxI - IxIxIxI0I1I2IxIxIxIxI - IxIxIxI3I4I5IxIxIxIxI - IxIxIxI6I7I8IxIxIxIxI - IxIxIxIxIxIxIxIxIxIxI - IxIxIxIxIxIxIxIxIxIxI - IxIxIxIxIxIxIxIxIxIxI - - For every unmasked pixel, the Convolver over-lays the PSF and computes three quantities; - - image_frame_indexes - The indexes of all masked pixels it will blur light into. - image_frame_kernels - The kernel values that overlap each masked pixel it will blur light into. - image_frame_length - The number of masked pixels it will blur light into (unmasked pixels are excluded) - - For example, if we had the following 3x3 kernel: - - I0.1I0.2I0.3I - I0.4I0.5I0.6I - I0.7I0.8I0.9I - - For pixel 0 above, when we overlap the kernel 4 unmasked pixels overlap this kernel, such that: - - image_frame_indexes = [0, 1, 3, 4] - image_frame_kernels = [0.5, 0.6, 0.8, 0.9] - image_frame_length = 4 - - Noting that the other 5 kernel values (0.1, 0.2, 0.3, 0.4, 0.7) overlap masked pixels and are thus discarded. - - For pixel 1, we get the following results: - - image_frame_indexes = [0, 1, 2, 3, 4, 5] - image_frame_kernels = [0.4, 0.5, 0.6, 0.7, 0.8, 0.9] - image_frame_lengths = 6 - - In the majority of cases, the kernel will overlap only unmasked pixels. This is the case above for \ - central pixel 4, where: - - image_frame_indexes = [0, 1, 2, 3, 4, 5, 6, 7, 8] - image_frame_kernels = [0,1, 0.2, 0,3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9] - image_frame_lengths = 9 - - Once we have set up all these quantities, the convolution routine simply uses them to convolve a 1D array of a - masked or the masked of a util in the inversion module. - - - **BLURRING FRAMES:** - - Whilst the scheme above accounts for all blurred light within the masks, it does not account for the fact that \ - pixels outside of the masks will also blur light into it. This effect is accounted for using blurring frames. - - It is omitted for mapping matrix blurring, as an inversion does not fit data outside of the masks. - - First, a blurring masks is computed from a masks, which describes all pixels which are close enough to the masks \ - to blur light into it for a given kernel size. Following the example above, the following blurring masks is \ - computed: - - IxIxIxIxIxIxIxIxIxIxI - IxIxIxIxIxIxIxIxIxIxI This is an example grid.Mask2D, where: - IxIxIxIxIxIxIxIxIxIxI - IxIxIoIoIoIoIoIxIxIxI x = `True` (Pixel is masked and excluded from lens) - IxIxIoIxIxIxIoIxIxIxI o = `False` (Pixel is not masked and included in lens) - IxIxIoIxIxIxIoIxIxIxI - IxIxIoIxIxIxIoIxIxIxI - IxIxIoIoIoIoIoIxIxIxI - IxIxIxIxIxIxIxIxIxIxI - IxIxIxIxIxIxIxIxIxIxI - - Indexing again goes from the top-left corner right and downwards: - - IxIxI xI xI xI xI xIxIxIxI - IxIxI xI xI xI xI xIxIxIxI - IxIxI xI xI xI xI xIxIxIxI - IxIxI 0I 1I 2I 3I 4IxIxIxI - IxIxI 5I xI xI xI 6IxIxIxI - IxIxI 7I xI xI xI 8IxIxIxI - IxIxI 9I xI xI xI10IxIxIxI - IxIxI11I12I13I14I15IxIxIxI - IxIxI xI xI xI xI xIxIxIxI - IxIxI xI xI xI xI xIxIxIxI - - For every unmasked blurring-pixel, the Convolver over-lays the PSF kernel and computes three quantities; - - blurring_frame_indexes - The indexes of all unmasked pixels (not unmasked blurring pixels) it will \ - blur light into. - bluring_frame_kernels - The kernel values that overlap each pixel it will blur light into. - blurring_frame_length - The number of pixels it will blur light into. - - The blurring frame therefore does not perform any blurring which blurs light into other blurring pixels. \ - It only performs computations which add light inside of the masks. - - For pixel 0 above, when we overlap the 3x3 kernel above only 1 unmasked pixels overlaps the kernel, such that: - - blurring_frame_indexes = [0] (This 0 refers to pixel 0 within the masks, not blurring_frame_pixel 0) - blurring_frame_kernels = [0.9] - blurring_frame_length = 1 - - For pixel 1 above, when we overlap the 3x3 kernel above 2 unmasked pixels overlap the kernel, such that: - - blurring_frame_indexes = [0, 1] (This 0 and 1 refer to pixels 0 and 1 within the masks) - blurring_frame_kernels = [0.8, 0.9] - blurring_frame_length = 2 - - For pixel 3 above, when we overlap the 3x3 kernel above 3 unmasked pixels overlap the kernel, such that: - - blurring_frame_indexes = [0, 1, 2] (Again, these are pixels 0, 1 and 2) - blurring_frame_kernels = [0.7, 0.8, 0.9] - blurring_frame_length = 3 - - Parameters - ---------- - mask - The mask within which the convolved signal is calculated. - blurring_mask - A masks of pixels outside the masks but whose light blurs into it after PSF convolution. - kernel : grid.PSF or ndarray - An array representing a PSF. - """ - if kernel.shape_native[0] % 2 == 0 or kernel.shape_native[1] % 2 == 0: - raise exc.KernelException("PSF kernel must be odd") - - self.mask = mask - - self.mask_index_array = np.full(mask.shape, -1) - self.pixels_in_mask = int(np.size(mask) - np.sum(mask)) - - count = 0 - for x in range(mask.shape[0]): - for y in range(mask.shape[1]): - if not mask[x, y]: - self.mask_index_array[x, y] = count - count += 1 - - self.kernel = kernel - self.kernel_max_size = self.kernel.shape_native[0] * self.kernel.shape_native[1] - - mask_1d_index = 0 - self.image_frame_1d_indexes = np.zeros( - (self.pixels_in_mask, self.kernel_max_size), dtype="int" - ) - self.image_frame_1d_kernels = np.zeros( - (self.pixels_in_mask, self.kernel_max_size) - ) - self.image_frame_1d_lengths = np.zeros((self.pixels_in_mask), dtype="int") - for x in range(self.mask_index_array.shape[0]): - for y in range(self.mask_index_array.shape[1]): - if not mask[x][y]: - ( - image_frame_1d_indexes, - image_frame_1d_kernels, - ) = self.frame_at_coordinates_jit( - coordinates=(x, y), - mask=np.array(mask), - mask_index_array=self.mask_index_array, - kernel_2d=self.kernel.native, - ) - self.image_frame_1d_indexes[mask_1d_index, :] = ( - image_frame_1d_indexes - ) - self.image_frame_1d_kernels[mask_1d_index, :] = ( - image_frame_1d_kernels - ) - self.image_frame_1d_lengths[mask_1d_index] = image_frame_1d_indexes[ - image_frame_1d_indexes >= 0 - ].shape[0] - mask_1d_index += 1 - - self.blurring_mask = mask_2d_util.blurring_mask_2d_from( - mask_2d=np.array(mask), - kernel_shape_native=kernel.shape_native, - ) - - self.pixels_in_blurring_mask = int( - np.size(self.blurring_mask) - np.sum(self.blurring_mask) - ) - - mask_1d_index = 0 - self.blurring_frame_1d_indexes = np.zeros( - (self.pixels_in_blurring_mask, self.kernel_max_size), dtype="int" - ) - self.blurring_frame_1d_kernels = np.zeros( - (self.pixels_in_blurring_mask, self.kernel_max_size) - ) - self.blurring_frame_1d_lengths = np.zeros( - (self.pixels_in_blurring_mask), dtype="int" - ) - for x in range(mask.shape[0]): - for y in range(mask.shape[1]): - if mask[x][y] and not self.blurring_mask[x, y]: - ( - image_frame_1d_indexes, - image_frame_1d_kernels, - ) = self.frame_at_coordinates_jit( - coordinates=(x, y), - mask=np.array(mask), - mask_index_array=np.array(self.mask_index_array), - kernel_2d=np.array(self.kernel.native), - ) - self.blurring_frame_1d_indexes[mask_1d_index, :] = ( - image_frame_1d_indexes - ) - self.blurring_frame_1d_kernels[mask_1d_index, :] = ( - image_frame_1d_kernels - ) - self.blurring_frame_1d_lengths[mask_1d_index] = ( - image_frame_1d_indexes[image_frame_1d_indexes >= 0].shape[0] - ) - mask_1d_index += 1 - - @staticmethod - @numba_util.jit() - def frame_at_coordinates_jit(coordinates, mask, mask_index_array, kernel_2d): - """ - Returns the frame (indexes of pixels light is blurred into) and kernel_frame (kernel kernel values of those \ - pixels) for a given coordinate in a masks and its PSF. - - Parameters - ---------- - coordinates: Tuple[int, int] - The coordinates of mask_index_array on which the frame should be centred - kernel_shape_native: Tuple[int, int] - The shape of the kernel for which this frame will be used - """ - - kernel_shape_native = kernel_2d.shape - kernel_max_size = kernel_shape_native[0] * kernel_shape_native[1] - - half_x = int(kernel_shape_native[0] / 2) - half_y = int(kernel_shape_native[1] / 2) - - frame = -1 * np.ones((kernel_max_size)) - kernel_frame = -1.0 * np.ones((kernel_max_size)) - - count = 0 - for i in range(kernel_shape_native[0]): - for j in range(kernel_shape_native[1]): - x = coordinates[0] - half_x + i - y = coordinates[1] - half_y + j - if ( - 0 <= x < mask_index_array.shape[0] - and 0 <= y < mask_index_array.shape[1] - ): - value = mask_index_array[x, y] - if value >= 0 and not mask[x, y]: - frame[count] = value - kernel_frame[count] = kernel_2d[i, j] - count += 1 - - return frame, kernel_frame - - def jax_convolve(self, image, blurring_image, method="auto"): - slim_to_2D_index_image = jnp.nonzero( - jnp.logical_not(self.mask.array), size=image.shape[0] - ) - slim_to_2D_index_blurring = jnp.nonzero( - jnp.logical_not(self.blurring_mask), size=blurring_image.shape[0] - ) - expanded_image_native = jnp.zeros(self.mask.shape) - expanded_image_native = expanded_image_native.at[slim_to_2D_index_image].set( - image.array - ) - expanded_image_native = expanded_image_native.at[slim_to_2D_index_blurring].set( - blurring_image.array - ) - kernel = np.array(self.kernel.native.array) - convolve_native = jax.scipy.signal.convolve( - expanded_image_native, kernel, mode="same", method=method - ) - convolve_slim = convolve_native[slim_to_2D_index_image] - return convolve_slim - - def convolve_image(self, image, blurring_image, jax_method="fft"): - """ - For a given 1D array and blurring array, convolve the two using this convolver. - - Parameters - ---------- - image - 1D array of the values which are to be blurred with the convolver's PSF. - blurring_image - 1D array of the blurring values which blur into the array after PSF convolution. - jax_method - If JAX is enabled this keyword will indicate what method is used for the PSF - convolution. Can be either `direct` to calculate it in real space or `fft` - to calculated it via a fast Fourier transform. `fft` is typically faster for - kernels that are more than about 5x5. Default is `fft`. - """ - - def exception_message(): - raise exc.KernelException( - "You cannot use the convolve_image function of a Convolver if the Convolver was" - "not created with a blurring_mask." - ) - - if use_jax: - jax.lax.cond( - self.blurring_mask is None, - lambda _: jax.debug.callback(exception_message), - lambda _: None, - None, - ) - - return self.jax_convolve(image, blurring_image, method=jax_method) - - else: - if self.blurring_mask is None: - exception_message() - - convolved_image = self.convolve_jit( - image_1d_array=np.array(image.slim), - image_frame_1d_indexes=self.image_frame_1d_indexes, - image_frame_1d_kernels=self.image_frame_1d_kernels, - image_frame_1d_lengths=self.image_frame_1d_lengths, - blurring_1d_array=np.array(blurring_image.slim), - blurring_frame_1d_indexes=self.blurring_frame_1d_indexes, - blurring_frame_1d_kernels=self.blurring_frame_1d_kernels, - blurring_frame_1d_lengths=self.blurring_frame_1d_lengths, - ) - - return Array2D(values=convolved_image, mask=self.mask) - - @staticmethod - @numba_util.jit() - def convolve_jit( - image_1d_array, - image_frame_1d_indexes, - image_frame_1d_kernels, - image_frame_1d_lengths, - blurring_1d_array, - blurring_frame_1d_indexes, - blurring_frame_1d_kernels, - blurring_frame_1d_lengths, - ): - blurred_image_1d = np.zeros(image_1d_array.shape) - - for image_1d_index in range(len(image_1d_array)): - frame_1d_indexes = image_frame_1d_indexes[image_1d_index] - frame_1d_kernel = image_frame_1d_kernels[image_1d_index] - frame_1d_length = image_frame_1d_lengths[image_1d_index] - image_value = image_1d_array[image_1d_index] - - for kernel_1d_index in range(frame_1d_length): - vector_index = frame_1d_indexes[kernel_1d_index] - kernel_value = frame_1d_kernel[kernel_1d_index] - blurred_image_1d[vector_index] += image_value * kernel_value - - for blurring_1d_index in range(len(blurring_1d_array)): - frame_1d_indexes = blurring_frame_1d_indexes[blurring_1d_index] - frame_1d_kernel = blurring_frame_1d_kernels[blurring_1d_index] - frame_1d_length = blurring_frame_1d_lengths[blurring_1d_index] - image_value = blurring_1d_array[blurring_1d_index] - - for kernel_1d_index in range(frame_1d_length): - vector_index = frame_1d_indexes[kernel_1d_index] - kernel_value = frame_1d_kernel[kernel_1d_index] - blurred_image_1d[vector_index] += image_value * kernel_value - - return blurred_image_1d - - def convolve_image_no_blurring(self, image): - """For a given 1D array and blurring array, convolve the two using this convolver. - - Parameters - ---------- - image - 1D array of the values which are to be blurred with the convolver's PSF. - """ - - convolved_image = self.convolve_no_blurring_jit( - image_1d_array=np.array(image.slim), - image_frame_1d_indexes=self.image_frame_1d_indexes, - image_frame_1d_kernels=self.image_frame_1d_kernels, - image_frame_1d_lengths=self.image_frame_1d_lengths, - ) - - return Array2D(values=convolved_image, mask=self.mask) - - def convolve_image_no_blurring_interpolation(self, image): - """For a given 1D array and blurring array, convolve the two using this convolver. - - Parameters - ---------- - image - 1D array of the values which are to be blurred with the convolver's PSF. - """ - - convolved_image = self.convolve_no_blurring_jit( - image_1d_array=image, - image_frame_1d_indexes=self.image_frame_1d_indexes, - image_frame_1d_kernels=self.image_frame_1d_kernels, - image_frame_1d_lengths=self.image_frame_1d_lengths, - ) - - return Array2D(values=convolved_image, mask=self.mask) - - @staticmethod - @numba_util.jit() - def convolve_no_blurring_jit( - image_1d_array, - image_frame_1d_indexes, - image_frame_1d_kernels, - image_frame_1d_lengths, - ): - blurred_image_1d = np.zeros(image_1d_array.shape) - - for image_1d_index in range(len(image_1d_array)): - frame_1d_indexes = image_frame_1d_indexes[image_1d_index] - frame_1d_kernel = image_frame_1d_kernels[image_1d_index] - frame_1d_length = image_frame_1d_lengths[image_1d_index] - image_value = image_1d_array[image_1d_index] - - for kernel_1d_index in range(frame_1d_length): - vector_index = frame_1d_indexes[kernel_1d_index] - kernel_value = frame_1d_kernel[kernel_1d_index] - blurred_image_1d[vector_index] += image_value * kernel_value - - return blurred_image_1d - - def convolve_mapping_matrix(self, mapping_matrix): - """For a given inversion mapping matrix, convolve every pixel's mapped with the PSF kernel. - - A mapping matrix provides non-zero entries in all elements which map two pixels to one another - (see *inversions.mappers*). - - For example, lets take an which is masked using a 'cross' of 5 pixels: - - [[ True, False, True]], - [[False, False, False]], - [[ True, False, True]] - - As example mapping matrix of this cross is as follows (5 pixels x 3 source pixels): - - [1, 0, 0] [0->0] - [1, 0, 0] [1->0] - [0, 1, 0] [2->1] - [0, 1, 0] [3->1] - [0, 0, 1] [4->2] - - For each source-pixel, we can create an of its unit-surface brightnesses by util the non-zero - entries back to masks. For example, doing this for source pixel 1 gives: - - [[0.0, 1.0, 0.0]], - [[1.0, 0.0, 0.0]] - [[0.0, 0.0, 0.0]] - - And source pixel 2: - - [[0.0, 0.0, 0.0]], - [[0.0, 1.0, 1.0]] - [[0.0, 0.0, 0.0]] - - We then convolve each of these with our PSF kernel, in 2 dimensions, like we would a grid. For - example, using the kernel below: - - kernel: - - [[0.0, 0.1, 0.0]] - [[0.1, 0.6, 0.1]] - [[0.0, 0.1, 0.0]] - - Blurred Source Pixel 1 (we don't need to perform the convolution into masked pixels): - - [[0.0, 0.6, 0.0]], - [[0.6, 0.0, 0.0]], - [[0.0, 0.0, 0.0]] - - Blurred Source pixel 2: - - [[0.0, 0.0, 0.0]], - [[0.0, 0.7, 0.7]], - [[0.0, 0.0, 0.0]] - - Finally, we map each of these blurred back to a blurred mapping matrix, which is analogous to the - mapping matrix. - - [0.6, 0.0, 0.0] [0->0] - [0.6, 0.0, 0.0] [1->0] - [0.0, 0.7, 0.0] [2->1] - [0.0, 0.7, 0.0] [3->1] - [0.0, 0.0, 0.6] [4->2] - - If the mapping matrix is sub-gridded, we perform the convolution on the fractional surface brightnesses in an - identical fashion to above. - - Parameters - ---------- - mapping_matrix - The 2D mapping matrix describing how every inversion pixel maps to a pixel on the data pixel. - """ - return self.convolve_matrix_jit( - mapping_matrix=mapping_matrix, - image_frame_1d_indexes=self.image_frame_1d_indexes, - image_frame_1d_kernels=self.image_frame_1d_kernels, - image_frame_1d_lengths=self.image_frame_1d_lengths, - ) - - @staticmethod - @numba_util.jit() - def convolve_matrix_jit( - mapping_matrix, - image_frame_1d_indexes, - image_frame_1d_kernels, - image_frame_1d_lengths, - ): - blurred_mapping_matrix = np.zeros(mapping_matrix.shape) - - for pixel_1d_index in range(mapping_matrix.shape[1]): - for image_1d_index in range(mapping_matrix.shape[0]): - value = mapping_matrix[image_1d_index, pixel_1d_index] - - if value > 0: - frame_1d_indexes = image_frame_1d_indexes[image_1d_index] - frame_1d_kernel = image_frame_1d_kernels[image_1d_index] - frame_1d_length = image_frame_1d_lengths[image_1d_index] - - for kernel_1d_index in range(frame_1d_length): - vector_index = frame_1d_indexes[kernel_1d_index] - kernel_value = frame_1d_kernel[kernel_1d_index] - blurred_mapping_matrix[vector_index, pixel_1d_index] += ( - value * kernel_value - ) - - return blurred_mapping_matrix diff --git a/autoarray/operators/mock/mock_convolver.py b/autoarray/operators/mock/mock_convolver.py index 290be228f..7dc5dfcbb 100644 --- a/autoarray/operators/mock/mock_convolver.py +++ b/autoarray/operators/mock/mock_convolver.py @@ -1,4 +1,4 @@ -class MockConvolver: +class MockPSF: def __init__(self, operated_mapping_matrix=None): self.operated_mapping_matrix = operated_mapping_matrix diff --git a/autoarray/structures/arrays/kernel_2d.py b/autoarray/structures/arrays/kernel_2d.py index f80e3f6e3..09ff331ef 100644 --- a/autoarray/structures/arrays/kernel_2d.py +++ b/autoarray/structures/arrays/kernel_2d.py @@ -1,4 +1,6 @@ from astropy import units +import jax +import jax.numpy as jnp import numpy as np import scipy.signal from pathlib import Path @@ -361,7 +363,7 @@ def rescaled_with_odd_dimensions_from( ) -> "Kernel2D": """ If the PSF kernel has one or two even-sized dimensions, return a PSF object where the kernel has odd-sized - dimensions (odd-sized dimensions are required by a *Convolver*). + dimensions (odd-sized dimensions are required for 2D convolution). The PSF can be scaled to larger / smaller sizes than the input size, if the rescale factor uses values that deviate furher from 1.0. @@ -511,3 +513,24 @@ def convolved_array_with_mask_from(self, array: Array2D, mask: Mask2D) -> Array2 ) return Array2D(values=convolved_array_1d, mask=mask) + + def jax_convolve(self, image, blurring_image, method="auto"): + slim_to_2D_index_image = jnp.nonzero( + jnp.logical_not(self.mask.array), size=image.shape[0] + ) + slim_to_2D_index_blurring = jnp.nonzero( + jnp.logical_not(self.blurring_mask), size=blurring_image.shape[0] + ) + expanded_image_native = jnp.zeros(self.mask.shape) + expanded_image_native = expanded_image_native.at[slim_to_2D_index_image].set( + image.array + ) + expanded_image_native = expanded_image_native.at[slim_to_2D_index_blurring].set( + blurring_image.array + ) + kernel = np.array(self.kernel.native.array) + convolve_native = jax.scipy.signal.convolve( + expanded_image_native, kernel, mode="same", method=method + ) + convolve_slim = convolve_native[slim_to_2D_index_image] + return convolve_slim \ No newline at end of file diff --git a/test_autoarray/conftest.py b/test_autoarray/conftest.py index 1dbe19e19..657ac4b0d 100644 --- a/test_autoarray/conftest.py +++ b/test_autoarray/conftest.py @@ -119,11 +119,6 @@ def make_psf_3x3(): return fixtures.make_psf_3x3() -@pytest.fixture(name="convolver_7x7") -def make_convolver_7x7(): - return fixtures.make_convolver_7x7() - - @pytest.fixture(name="imaging_7x7") def make_imaging_7x7(): return fixtures.make_imaging_7x7() diff --git a/test_autoarray/dataset/imaging/test_dataset.py b/test_autoarray/dataset/imaging/test_dataset.py index 6025733f5..9c079b94d 100644 --- a/test_autoarray/dataset/imaging/test_dataset.py +++ b/test_autoarray/dataset/imaging/test_dataset.py @@ -37,14 +37,14 @@ def test__psf_and_mask_hit_edge__automatically_pads_image_and_noise_map(): psf = aa.Kernel2D.ones(shape_native=(3, 3), pixel_scales=1.0) dataset = aa.Imaging( - data=image, noise_map=noise_map, psf=psf, pad_for_convolver=False + data=image, noise_map=noise_map, psf=psf, pad_for_psf=False ) assert dataset.data.shape_native == (3, 3) assert dataset.noise_map.shape_native == (3, 3) dataset = aa.Imaging( - data=image, noise_map=noise_map, psf=psf, pad_for_convolver=True + data=image, noise_map=noise_map, psf=psf, pad_for_psf=True ) assert dataset.data.shape_native == (5, 5) @@ -144,7 +144,6 @@ def test__apply_mask(imaging_7x7, mask_2d_7x7, psf_3x3): assert (masked_imaging_7x7.psf.slim == (1.0 / 3.0) * psf_3x3.slim).all() assert type(masked_imaging_7x7.psf) == aa.Kernel2D - assert type(masked_imaging_7x7.convolver) == aa.Convolver assert masked_imaging_7x7.w_tilde.curvature_preload.shape == (35,) assert masked_imaging_7x7.w_tilde.indexes.shape == (35,) assert masked_imaging_7x7.w_tilde.lengths.shape == (9,) @@ -226,7 +225,6 @@ def test__different_imaging_without_mock_objects__customize_constructor_inputs() assert masked_dataset.psf.native == pytest.approx( (1.0 / 49.0) * np.ones((7, 7)), 1.0e-4 ) - assert masked_dataset.convolver.kernel.shape_native == (7, 7) assert (masked_dataset.data == np.array([1.0])).all() assert (masked_dataset.noise_map == np.array([2.0])).all() diff --git a/test_autoarray/inversion/inversion/imaging/test_imaging.py b/test_autoarray/inversion/inversion/imaging/test_imaging.py index 46d726bad..df615ad85 100644 --- a/test_autoarray/inversion/inversion/imaging/test_imaging.py +++ b/test_autoarray/inversion/inversion/imaging/test_imaging.py @@ -11,18 +11,18 @@ directory = path.dirname(path.realpath(__file__)) -def test__operated_mapping_matrix_property(convolver_7x7, rectangular_mapper_7x7_3x3): +def test__operated_mapping_matrix_property(psf_7x7, rectangular_mapper_7x7_3x3): inversion = aa.m.MockInversionImaging( - convolver=convolver_7x7, linear_obj_list=[rectangular_mapper_7x7_3x3] + psf=psf_7x7, linear_obj_list=[rectangular_mapper_7x7_3x3] ) assert inversion.operated_mapping_matrix_list[0][0, 0] == pytest.approx(1.0, 1e-4) assert inversion.operated_mapping_matrix[0, 0] == pytest.approx(1.0, 1e-4) - convolver = aa.m.MockConvolver(operated_mapping_matrix=np.ones((2, 2))) + psf = aa.m.MockPSF(operated_mapping_matrix=np.ones((2, 2))) inversion = aa.m.MockInversionImaging( - convolver=convolver, + psf=psf, linear_obj_list=[rectangular_mapper_7x7_3x3, rectangular_mapper_7x7_3x3], ) @@ -42,9 +42,9 @@ def test__operated_mapping_matrix_property(convolver_7x7, rectangular_mapper_7x7 def test__operated_mapping_matrix_property__with_operated_mapping_matrix_override( - convolver_7x7, rectangular_mapper_7x7_3x3 + psf_7x7, rectangular_mapper_7x7_3x3 ): - convolver = aa.m.MockConvolver(operated_mapping_matrix=np.ones((2, 2))) + psf = aa.m.MockPSF(operated_mapping_matrix=np.ones((2, 2))) operated_mapping_matrix_override = np.array([[1.0, 2.0], [3.0, 4.0]]) @@ -54,7 +54,7 @@ def test__operated_mapping_matrix_property__with_operated_mapping_matrix_overrid ) inversion = aa.m.MockInversionImaging( - convolver=convolver, linear_obj_list=[rectangular_mapper_7x7_3x3, linear_obj] + psf=psf, linear_obj_list=[rectangular_mapper_7x7_3x3, linear_obj] ) operated_mapping_matrix_0 = np.array([[1.0, 1.0], [1.0, 1.0]]) @@ -73,7 +73,7 @@ def test__operated_mapping_matrix_property__with_operated_mapping_matrix_overrid def test__curvature_matrix(rectangular_mapper_7x7_3x3): noise_map = np.ones(2) - convolver = aa.m.MockConvolver(operated_mapping_matrix=np.ones((2, 10))) + psf = aa.m.MockPSF(operated_mapping_matrix=np.ones((2, 10))) operated_mapping_matrix_override = np.array([[1.0, 2.0], [3.0, 4.0]]) @@ -87,7 +87,7 @@ def test__curvature_matrix(rectangular_mapper_7x7_3x3): dataset = aa.DatasetInterface( data=np.ones(2), noise_map=noise_map, - convolver=convolver, + psf=psf, ) inversion = aa.InversionImagingMapping( @@ -135,7 +135,7 @@ def test__w_tilde_checks_noise_map_and_raises_exception_if_preloads_dont_match_n dataset = aa.DatasetInterface( data=np.ones(9), noise_map=np.ones(9), - convolver=aa.m.MockConvolver(matrix_shape), + psf=aa.m.MockPSF(matrix_shape), ) # noinspection PyTypeChecker diff --git a/test_autoarray/inversion/inversion/imaging/test_inversion_imaging_util.py b/test_autoarray/inversion/inversion/imaging/test_inversion_imaging_util.py index 711135ebc..b616a95d4 100644 --- a/test_autoarray/inversion/inversion/imaging/test_inversion_imaging_util.py +++ b/test_autoarray/inversion/inversion/imaging/test_inversion_imaging_util.py @@ -181,7 +181,7 @@ def test__data_vector_via_w_tilde_data_two_methods_agree(): shape_native=(7, 7), pixel_scales=mask.pixel_scales, sigma=1.0, normalize=True ) - convolver = aa.Convolver(mask=mask, kernel=kernel) + psf = kernel pixelization = aa.mesh.Rectangular(shape=(20, 20)) @@ -203,7 +203,7 @@ def test__data_vector_via_w_tilde_data_two_methods_agree(): mapping_matrix = mapper.mapping_matrix - blurred_mapping_matrix = convolver.convolve_mapping_matrix( + blurred_mapping_matrix = psf.convolve_mapping_matrix( mapping_matrix=mapping_matrix ) @@ -258,7 +258,7 @@ def test__curvature_matrix_via_w_tilde_two_methods_agree(): shape_native=(7, 7), pixel_scales=mask.pixel_scales, sigma=1.0, normalize=True ) - convolver = aa.Convolver(mask=mask, kernel=kernel) + psf = kernel pixelization = aa.mesh.Rectangular(shape=(20, 20)) @@ -282,7 +282,7 @@ def test__curvature_matrix_via_w_tilde_two_methods_agree(): w_tilde=w_tilde, mapping_matrix=mapping_matrix ) - blurred_mapping_matrix = convolver.convolve_mapping_matrix( + blurred_mapping_matrix = psf.convolve_mapping_matrix( mapping_matrix=mapping_matrix ) @@ -303,7 +303,7 @@ def test__curvature_matrix_via_w_tilde_preload_two_methods_agree(): shape_native=(7, 7), pixel_scales=mask.pixel_scales, sigma=1.0, normalize=True ) - convolver = aa.Convolver(mask=mask, kernel=kernel) + psf = kernel pixelization = aa.mesh.Rectangular(shape=(20, 20)) @@ -356,7 +356,7 @@ def test__curvature_matrix_via_w_tilde_preload_two_methods_agree(): pix_pixels=pixelization.pixels, ) - blurred_mapping_matrix = convolver.convolve_mapping_matrix( + blurred_mapping_matrix = psf.convolve_mapping_matrix( mapping_matrix=mapping_matrix ) diff --git a/test_autoarray/operators/test_convolver.py b/test_autoarray/operators/test_convolver.py deleted file mode 100644 index c298dda31..000000000 --- a/test_autoarray/operators/test_convolver.py +++ /dev/null @@ -1,1164 +0,0 @@ -import numpy as np -import pytest - -import autoarray as aa -from autoarray import exc - - -@pytest.fixture(name="simple_mask_2d_7x7") -def make_simple_mask_2d_7x7(): - mask = [ - [True, True, True, True, True, True, True], - [True, True, True, True, True, True, True], - [True, True, False, False, False, True, True], - [True, True, False, False, False, True, True], - [True, True, False, False, False, True, True], - [True, True, True, True, True, True, True], - [True, True, True, True, True, True, True], - ] - - return aa.Mask2D(mask=mask, pixel_scales=1.0) - - -@pytest.fixture(name="simple_mask_5x5") -def make_simple_mask_5x5(): - mask = [ - [True, True, True, True, True], - [True, False, False, False, True], - [True, False, False, False, True], - [True, False, False, False, True], - [True, True, True, True, True], - ] - - return aa.Mask2D(mask=mask, pixel_scales=1.0) - - -@pytest.fixture(name="cross_mask") -def make_cross_mask(): - mask = np.full((5, 5), True) - - mask[2, 2] = False - mask[1, 2] = False - mask[3, 2] = False - mask[2, 1] = False - mask[2, 3] = False - - return aa.Mask2D(mask=mask, pixel_scales=1.0) - - -def test__numbering__uses_mask_correctly(simple_mask_5x5, cross_mask): - convolver = aa.Convolver( - mask=simple_mask_5x5, - kernel=aa.Kernel2D.ones(shape_native=(1, 1), pixel_scales=1.0), - ) - - mask_index_array = convolver.mask_index_array - - assert mask_index_array.shape == (5, 5) - # noinspection PyUnresolvedReferences - assert ( - mask_index_array - == np.array( - [ - [-1, -1, -1, -1, -1], - [-1, 0, 1, 2, -1], - [-1, 3, 4, 5, -1], - [-1, 6, 7, 8, -1], - [-1, -1, -1, -1, -1], - ] - ) - ).all() - - convolver = aa.Convolver( - mask=cross_mask, kernel=aa.Kernel2D.ones(shape_native=(1, 1), pixel_scales=1.0) - ) - - assert ( - convolver.mask_index_array - == np.array( - [ - [-1, -1, -1, -1, -1], - [-1, -1, 0, -1, -1], - [-1, 1, 2, 3, -1], - [-1, -1, 4, -1, -1], - [-1, -1, -1, -1, -1], - ] - ) - ).all() - - -def test__even_kernel_failure(): - with pytest.raises(exc.KernelException): - aa.Convolver( - mask=np.full((3, 3), False), - kernel=aa.Kernel2D.ones(shape_native=(2, 2), pixel_scales=1.0), - ) - - -def test__frame_extraction__frame_and_kernel_frame_at_coords(simple_mask_5x5): - convolver = aa.Convolver( - mask=simple_mask_5x5, - kernel=aa.Kernel2D.no_mask( - [[1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [7.0, 8.0, 9.0]], pixel_scales=1.0 - ), - ) - - frame, kernel_frame = convolver.frame_at_coordinates_jit( - coordinates=(2, 2), - mask=np.array(simple_mask_5x5), - mask_index_array=convolver.mask_index_array, - kernel_2d=np.array(convolver.kernel.native), - ) - - assert (frame == np.array([i for i in range(9)])).all() - - assert ( - kernel_frame == np.array([[1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0]]) - ).all() - - corner_frame = np.array([0, 1, 3, 4, -1, -1, -1, -1, -1]) - - frame, kernel_frame = convolver.frame_at_coordinates_jit( - coordinates=(1, 1), - mask=np.array(simple_mask_5x5), - mask_index_array=convolver.mask_index_array, - kernel_2d=np.array(convolver.kernel.native), - ) - - assert (frame == corner_frame).all() - - frame, kernel_frame = convolver.frame_at_coordinates_jit( - coordinates=(1, 1), - mask=np.array(simple_mask_5x5), - mask_index_array=convolver.mask_index_array, - kernel_2d=np.array(convolver.kernel.native), - ) - - assert (kernel_frame == np.array([5.0, 6.0, 8.0, 9.0, -1, -1, -1, -1, -1])).all() - - assert 9 == len(convolver.image_frame_1d_indexes) - - assert ( - convolver.image_frame_1d_indexes[4] == np.array([i for i in range(9)]) - ).all() - - -def test__frame_extraction__more_complicated_frames(simple_mask_2d_7x7): - convolver = aa.Convolver( - mask=simple_mask_2d_7x7, - kernel=aa.Kernel2D.no_mask( - [ - [1.0, 2.0, 3.0, 4.0, 5.0], - [6.0, 7.0, 8.0, 9.0, 10.0], - [11.0, 12.0, 13.0, 14.0, 15.0], - [16.0, 17.0, 18.0, 19.0, 20.0], - [21.0, 22.0, 23.0, 24.0, 25.0], - ], - pixel_scales=1.0, - ), - ) - - frame, kernel_frame = convolver.frame_at_coordinates_jit( - coordinates=(2, 2), - mask=np.array(simple_mask_2d_7x7), - mask_index_array=convolver.mask_index_array, - kernel_2d=np.array(convolver.kernel.native), - ) - - assert ( - kernel_frame - == np.array( - [ - 13.0, - 14.0, - 15.0, - 18.0, - 19.0, - 20.0, - 23.0, - 24.0, - 25.0, - -1.0, - -1.0, - -1.0, - -1.0, - -1.0, - -1.0, - -1.0, - -1.0, - -1.0, - -1.0, - -1.0, - -1.0, - -1.0, - -1.0, - -1.0, - -1.0, - ] - ) - ).all() - - frame, kernel_frame = convolver.frame_at_coordinates_jit( - coordinates=(3, 2), - mask=np.array(simple_mask_2d_7x7), - mask_index_array=convolver.mask_index_array, - kernel_2d=np.array(convolver.kernel.native), - ) - - assert ( - kernel_frame - == np.array( - [ - 8.0, - 9.0, - 10.0, - 13.0, - 14.0, - 15.0, - 18.0, - 19.0, - 20.0, - -1.0, - -1.0, - -1.0, - -1.0, - -1.0, - -1.0, - -1.0, - -1.0, - -1.0, - -1.0, - -1.0, - -1.0, - -1.0, - -1.0, - -1.0, - -1, - ] - ) - ).all() - - frame, kernel_frame = convolver.frame_at_coordinates_jit( - coordinates=(3, 3), - mask=np.array(simple_mask_2d_7x7), - mask_index_array=convolver.mask_index_array, - kernel_2d=np.array(convolver.kernel.native), - ) - - assert ( - kernel_frame - == np.array( - [ - 7.0, - 8.0, - 9.0, - 12.0, - 13.0, - 14.0, - 17.0, - 18.0, - 19.0, - -1.0, - -1.0, - -1.0, - -1.0, - -1.0, - -1.0, - -1.0, - -1.0, - -1.0, - -1.0, - -1.0, - -1.0, - -1.0, - -1.0, - -1.0, - -1, - ] - ) - ).all() - - -def test__image_frame_indexes__for_different_masks(cross_mask, simple_mask_2d_7x7): - convolver = aa.Convolver( - mask=cross_mask, - kernel=aa.Kernel2D.no_mask( - [[1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [7.0, 8.0, 9.0]], pixel_scales=1.0 - ), - ) - - assert 5 == len(convolver.image_frame_1d_indexes) - - assert ( - convolver.image_frame_1d_indexes[0] - == np.array([0, 1, 2, 3, -1, -1, -1, -1, -1]) - ).all() - assert ( - convolver.image_frame_1d_indexes[1] - == np.array([0, 1, 2, 4, -1, -1, -1, -1, -1]) - ).all() - assert ( - convolver.image_frame_1d_indexes[2] == np.array([0, 1, 2, 3, 4, -1, -1, -1, -1]) - ).all() - assert ( - convolver.image_frame_1d_indexes[3] - == np.array([0, 2, 3, 4, -1, -1, -1, -1, -1]) - ).all() - assert ( - convolver.image_frame_1d_indexes[4] - == np.array([1, 2, 3, 4, -1, -1, -1, -1, -1]) - ).all() - - convolver = aa.Convolver( - mask=simple_mask_2d_7x7, - kernel=aa.Kernel2D.ones(shape_native=(3, 5), pixel_scales=1.0), - ) - - assert ( - convolver.image_frame_1d_indexes[0] - == np.array([0, 1, 2, 3, 4, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1]) - ).all() - assert ( - convolver.image_frame_1d_indexes[1] - == np.array([0, 1, 2, 3, 4, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1]) - ).all() - assert ( - convolver.image_frame_1d_indexes[2] - == np.array([0, 1, 2, 3, 4, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1]) - ).all() - assert ( - convolver.image_frame_1d_indexes[3] - == np.array([0, 1, 2, 3, 4, 5, 6, 7, 8, -1, -1, -1, -1, -1, -1]) - ).all() - assert ( - convolver.image_frame_1d_indexes[4] - == np.array([0, 1, 2, 3, 4, 5, 6, 7, 8, -1, -1, -1, -1, -1, -1]) - ).all() - assert ( - convolver.image_frame_1d_indexes[5] - == np.array([0, 1, 2, 3, 4, 5, 6, 7, 8, -1, -1, -1, -1, -1, -1]) - ).all() - assert ( - convolver.image_frame_1d_indexes[6] - == np.array([3, 4, 5, 6, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1]) - ).all() - assert ( - convolver.image_frame_1d_indexes[7] - == np.array([3, 4, 5, 6, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1]) - ).all() - assert ( - convolver.image_frame_1d_indexes[8] - == np.array([3, 4, 5, 6, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1]) - ).all() - - convolver = aa.Convolver( - mask=simple_mask_2d_7x7, - kernel=aa.Kernel2D.ones(shape_native=(5, 3), pixel_scales=1.0), - ) - - assert ( - convolver.image_frame_1d_indexes[0] - == np.array([0, 1, 3, 4, 6, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1]) - ).all() - assert ( - convolver.image_frame_1d_indexes[1] - == np.array([0, 1, 2, 3, 4, 5, 6, 7, 8, -1, -1, -1, -1, -1, -1]) - ).all() - assert ( - convolver.image_frame_1d_indexes[2] - == np.array([1, 2, 4, 5, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1]) - ).all() - assert ( - convolver.image_frame_1d_indexes[3] - == np.array([0, 1, 3, 4, 6, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1]) - ).all() - assert ( - convolver.image_frame_1d_indexes[4] - == np.array([0, 1, 2, 3, 4, 5, 6, 7, 8, -1, -1, -1, -1, -1, -1]) - ).all() - assert ( - convolver.image_frame_1d_indexes[5] - == np.array([1, 2, 4, 5, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1]) - ).all() - assert ( - convolver.image_frame_1d_indexes[6] - == np.array([0, 1, 3, 4, 6, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1]) - ).all() - assert ( - convolver.image_frame_1d_indexes[7] - == np.array([0, 1, 2, 3, 4, 5, 6, 7, 8, -1, -1, -1, -1, -1, -1]) - ).all() - assert ( - convolver.image_frame_1d_indexes[8] - == np.array([1, 2, 4, 5, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1]) - ).all() - - convolver = aa.Convolver( - mask=simple_mask_2d_7x7, - kernel=aa.Kernel2D.ones(shape_native=(5, 5), pixel_scales=1.0), - ) - - assert ( - convolver.image_frame_1d_indexes[0] - == np.array( - [ - 0, - 1, - 2, - 3, - 4, - 5, - 6, - 7, - 8, - -1, - -1, - -1, - -1, - -1, - -1, - -1, - -1, - -1, - -1, - -1, - -1, - -1, - -1, - -1, - -1, - ] - ) - ).all() - assert ( - convolver.image_frame_1d_indexes[1] - == np.array( - [ - 0, - 1, - 2, - 3, - 4, - 5, - 6, - 7, - 8, - -1, - -1, - -1, - -1, - -1, - -1, - -1, - -1, - -1, - -1, - -1, - -1, - -1, - -1, - -1, - -1, - ] - ) - ).all() - assert ( - convolver.image_frame_1d_indexes[2] - == np.array( - [ - 0, - 1, - 2, - 3, - 4, - 5, - 6, - 7, - 8, - -1, - -1, - -1, - -1, - -1, - -1, - -1, - -1, - -1, - -1, - -1, - -1, - -1, - -1, - -1, - -1, - ] - ) - ).all() - assert ( - convolver.image_frame_1d_indexes[3] - == np.array( - [ - 0, - 1, - 2, - 3, - 4, - 5, - 6, - 7, - 8, - -1, - -1, - -1, - -1, - -1, - -1, - -1, - -1, - -1, - -1, - -1, - -1, - -1, - -1, - -1, - -1, - ] - ) - ).all() - assert ( - convolver.image_frame_1d_indexes[4] - == np.array( - [ - 0, - 1, - 2, - 3, - 4, - 5, - 6, - 7, - 8, - -1, - -1, - -1, - -1, - -1, - -1, - -1, - -1, - -1, - -1, - -1, - -1, - -1, - -1, - -1, - -1, - ] - ) - ).all() - assert ( - convolver.image_frame_1d_indexes[5] - == np.array( - [ - 0, - 1, - 2, - 3, - 4, - 5, - 6, - 7, - 8, - -1, - -1, - -1, - -1, - -1, - -1, - -1, - -1, - -1, - -1, - -1, - -1, - -1, - -1, - -1, - -1, - ] - ) - ).all() - assert ( - convolver.image_frame_1d_indexes[6] - == np.array( - [ - 0, - 1, - 2, - 3, - 4, - 5, - 6, - 7, - 8, - -1, - -1, - -1, - -1, - -1, - -1, - -1, - -1, - -1, - -1, - -1, - -1, - -1, - -1, - -1, - -1, - ] - ) - ).all() - assert ( - convolver.image_frame_1d_indexes[7] - == np.array( - [ - 0, - 1, - 2, - 3, - 4, - 5, - 6, - 7, - 8, - -1, - -1, - -1, - -1, - -1, - -1, - -1, - -1, - -1, - -1, - -1, - -1, - -1, - -1, - -1, - -1, - ] - ) - ).all() - assert ( - convolver.image_frame_1d_indexes[8] - == np.array( - [ - 0, - 1, - 2, - 3, - 4, - 5, - 6, - 7, - 8, - -1, - -1, - -1, - -1, - -1, - -1, - -1, - -1, - -1, - -1, - -1, - -1, - -1, - -1, - -1, - -1, - ] - ) - ).all() - - -def test_image_frame_kernels__different_shape_masks( - simple_mask_5x5, simple_mask_2d_7x7 -): - convolver = aa.Convolver( - mask=simple_mask_5x5, - kernel=aa.Kernel2D.no_mask( - [[1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [7.0, 8.0, 9.0]], pixel_scales=1.0 - ), - ) - - assert 9 == len(convolver.image_frame_1d_indexes) - - assert ( - convolver.image_frame_1d_kernels[4] - == np.array([[1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0]]) - ).all() - - convolver = aa.Convolver( - mask=simple_mask_2d_7x7, - kernel=aa.Kernel2D.no_mask( - [ - [1.0, 2.0, 3.0, 4.0, 5.0], - [6.0, 7.0, 8.0, 9.0, 10.0], - [11.0, 12.0, 13.0, 14.0, 15.0], - ], - pixel_scales=1.0, - ), - ) - - assert ( - convolver.image_frame_1d_kernels[0] - == np.array( - [8.0, 9.0, 10.0, 13.0, 14.0, 15.0, -1, -1, -1, -1, -1, -1, -1, -1, -1] - ) - ).all() - assert ( - convolver.image_frame_1d_kernels[1] - == np.array( - [7.0, 8.0, 9.0, 12.0, 13.0, 14.0, -1, -1, -1, -1, -1, -1, -1, -1, -1] - ) - ).all() - assert ( - convolver.image_frame_1d_kernels[2] - == np.array( - [6.0, 7.0, 8.0, 11.0, 12.0, 13.0, -1, -1, -1, -1, -1, -1, -1, -1, -1] - ) - ).all() - assert ( - convolver.image_frame_1d_kernels[3] - == np.array( - [3.0, 4.0, 5.0, 8.0, 9.0, 10.0, 13.0, 14.0, 15.0, -1, -1, -1, -1, -1, -1] - ) - ).all() - assert ( - convolver.image_frame_1d_kernels[4] - == np.array( - [2.0, 3.0, 4.0, 7.0, 8.0, 9.0, 12.0, 13.0, 14.0, -1, -1, -1, -1, -1, -1] - ) - ).all() - assert ( - convolver.image_frame_1d_kernels[5] - == np.array( - [1.0, 2.0, 3.0, 6.0, 7.0, 8.0, 11.0, 12.0, 13.0, -1, -1, -1, -1, -1, -1] - ) - ).all() - assert ( - convolver.image_frame_1d_kernels[6] - == np.array([3.0, 4.0, 5.0, 8.0, 9.0, 10.0, -1, -1, -1, -1, -1, -1, -1, -1, -1]) - ).all() - assert ( - convolver.image_frame_1d_kernels[7] - == np.array([2.0, 3.0, 4.0, 7.0, 8.0, 9.0, -1, -1, -1, -1, -1, -1, -1, -1, -1]) - ).all() - assert ( - convolver.image_frame_1d_kernels[8] - == np.array([1.0, 2.0, 3.0, 6.0, 7.0, 8.0, -1, -1, -1, -1, -1, -1, -1, -1, -1]) - ).all() - - convolver = aa.Convolver( - mask=simple_mask_2d_7x7, - kernel=aa.Kernel2D.no_mask( - [ - [1.0, 2.0, 3.0], - [4.0, 5.0, 6.0], - [7.0, 8.0, 9.0], - [10.0, 11.0, 12.0], - [13.0, 14.0, 15.0], - ], - pixel_scales=1.0, - ), - ) - - assert ( - convolver.image_frame_1d_kernels[0] - == np.array( - [8.0, 9.0, 11.0, 12.0, 14.0, 15.0, -1, -1, -1, -1, -1, -1, -1, -1, -1] - ) - ).all() - assert ( - convolver.image_frame_1d_kernels[1] - == np.array( - [7.0, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0, -1, -1, -1, -1, -1, -1] - ) - ).all() - assert ( - convolver.image_frame_1d_kernels[2] - == np.array( - [7.0, 8.0, 10.0, 11.0, 13.0, 14.0, -1, -1, -1, -1, -1, -1, -1, -1, -1] - ) - ).all() - assert ( - convolver.image_frame_1d_kernels[3] - == np.array( - [5.0, 6.0, 8.0, 9.0, 11.0, 12.0, -1, -1, -1, -1, -1, -1, -1, -1, -1] - ) - ).all() - assert ( - convolver.image_frame_1d_kernels[4] - == np.array( - [4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0, -1, -1, -1, -1, -1, -1] - ) - ).all() - assert ( - convolver.image_frame_1d_kernels[5] - == np.array( - [4.0, 5.0, 7.0, 8.0, 10.0, 11.0, -1, -1, -1, -1, -1, -1, -1, -1, -1] - ) - ).all() - assert ( - convolver.image_frame_1d_kernels[6] - == np.array([2.0, 3.0, 5.0, 6.0, 8.0, 9.0, -1, -1, -1, -1, -1, -1, -1, -1, -1]) - ).all() - assert ( - convolver.image_frame_1d_kernels[7] - == np.array( - [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, -1, -1, -1, -1, -1, -1] - ) - ).all() - assert ( - convolver.image_frame_1d_kernels[8] - == np.array([1.0, 2.0, 4.0, 5.0, 7.0, 8.0, -1, -1, -1, -1, -1, -1, -1, -1, -1]) - ).all() - - -def test__blurring_frame_indexes__blurring_region_3x3_kernel(cross_mask): - convolver = aa.Convolver( - mask=cross_mask, kernel=aa.Kernel2D.ones(shape_native=(3, 3), pixel_scales=1.0) - ) - - assert ( - convolver.blurring_frame_1d_indexes[4] - == np.array([0, 1, 2, -1, -1, -1, -1, -1, -1]) - ).all() - assert ( - convolver.blurring_frame_1d_indexes[5] - == np.array([0, 2, 3, -1, -1, -1, -1, -1, -1]) - ).all() - assert ( - convolver.blurring_frame_1d_indexes[10] - == np.array([1, 2, 4, -1, -1, -1, -1, -1, -1]) - ).all() - assert ( - convolver.blurring_frame_1d_indexes[11] - == np.array([2, 3, 4, -1, -1, -1, -1, -1, -1]) - ).all() - - -def test__blurring_frame_kernels__blurring_region_3x3_kernel(cross_mask): - convolver = aa.Convolver( - mask=cross_mask, - kernel=aa.Kernel2D.no_mask( - [[1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [7.0, 8.0, 9.0]], pixel_scales=1.0 - ), - ) - - assert ( - convolver.blurring_frame_1d_kernels[4] - == np.array([6.0, 8.0, 9.0, -1, -1, -1, -1, -1, -1]) - ).all() - assert ( - convolver.blurring_frame_1d_kernels[5] - == np.array([4.0, 7.0, 8.0, -1, -1, -1, -1, -1, -1]) - ).all() - assert ( - convolver.blurring_frame_1d_kernels[10] - == np.array([2.0, 3.0, 6.0, -1, -1, -1, -1, -1, -1]) - ).all() - assert ( - convolver.blurring_frame_1d_kernels[11] - == np.array([1.0, 2.0, 4.0, -1, -1, -1, -1, -1, -1]) - ).all() - - -def test__frame_lengths__frames_are_from_examples_above__lengths_are_right( - simple_mask_2d_7x7, -): - convolver = aa.Convolver( - mask=simple_mask_2d_7x7, - kernel=aa.Kernel2D.ones(shape_native=(3, 5), pixel_scales=1.0), - ) - - # convolver_image.image_frame_indexes[0] == np.array([0, 1, 2, 3, 4, 5]) - # convolver_image.image_frame_indexes[1] == np.array([0, 1, 2, 3, 4, 5]) - # convolver_image.image_frame_indexes[2] == np.array([0, 1, 2, 3, 4, 5]) - # convolver_image.image_frame_indexes[3] == np.array([0, 1, 2, 3, 4, 5, 6, 7, 8]) - # (convolver_image.image_frame_indexes[4] == np.array([0, 1, 2, 3, 4, 5, 6, 7, 8]) - # convolver_image.image_frame_indexes[5] == np.array([0, 1, 2, 3, 4, 5, 6, 7, 8]) - # convolver_image.image_frame_indexes[6] == np.array([3, 4, 5, 6, 7, 8]) - # convolver_image.image_frame_indexes[7] == np.array([3, 4, 5, 6, 7, 8]) - # convolver_image.image_frame_indexes[8] == np.array([3, 4, 5, 6, 7, 8]) - - assert ( - convolver.image_frame_1d_lengths == np.array([6, 6, 6, 9, 9, 9, 6, 6, 6]) - ).all() - - -def test__convolve_mapping_matrix__asymetric_convolver__matrix_blurred_correctly(): - mask = np.array( - [ - [True, True, True, True, True, True], - [True, False, False, False, False, True], - [True, False, False, False, False, True], - [True, False, False, False, False, True], - [True, False, False, False, False, True], - [True, True, True, True, True, True], - ] - ) - - asymmetric_kernel = aa.Kernel2D.no_mask( - values=[[0, 0.0, 0], [0.4, 0.2, 0.3], [0, 0.1, 0]], pixel_scales=1.0 - ) - - convolver = aa.Convolver(mask=mask, kernel=asymmetric_kernel) - - mapping = np.array( - [ - [0, 0, 0], - [0, 0, 0], - [0, 0, 0], - [0, 0, 0], - [0, 0, 0], - [0, 0, 0], - [0, 0, 0], - [ - 0, - 1, - 0, - ], # The 0.3 should be 'chopped' from this pixel as it is on the right-most edge - [0, 0, 0], - [1, 0, 0], - [0, 0, 1], - [0, 0, 0], - [0, 0, 0], - [0, 0, 0], - [0, 0, 0], - [0, 0, 0], - ] - ) - - blurred_mapping = convolver.convolve_mapping_matrix(mapping) - - assert ( - blurred_mapping - == np.array( - [ - [0, 0, 0], - [0, 0, 0], - [0, 0, 0], - [0, 0, 0], - [0, 0, 0], - [0, 0, 0], - [0, 0.4, 0], - [0, 0.2, 0], - [0.4, 0, 0], - [0.2, 0, 0.4], - [0.3, 0, 0.2], - [0, 0.1, 0.3], - [0, 0, 0], - [0.1, 0, 0], - [0, 0, 0.1], - [0, 0, 0], - ] - ) - ).all() - - asymmetric_kernel = aa.Kernel2D.no_mask( - values=[[0, 0.0, 0], [0.4, 0.2, 0.3], [0, 0.1, 0]], pixel_scales=1.0 - ) - - convolver = aa.Convolver(mask=mask, kernel=asymmetric_kernel) - - mapping = np.array( - [ - [0, 1, 0], - [0, 1, 0], - [0, 1, 0], - [0, 0, 0], - [0, 0, 0], - [0, 0, 0], - [0, 0, 0], - [ - 0, - 1, - 0, - ], # The 0.3 should be 'chopped' from this pixel as it is on the right-most edge - [1, 0, 0], - [1, 0, 0], - [0, 0, 1], - [0, 0, 0], - [0, 0, 0], - [0, 0, 0], - [0, 0, 0], - [0, 0, 0], - ] - ) - - blurred_mapping = convolver.convolve_mapping_matrix(mapping) - - assert blurred_mapping == pytest.approx( - np.array( - [ - [0, 0.6, 0], - [0, 0.9, 0], - [0, 0.5, 0], - [0, 0.3, 0], - [0, 0.1, 0], - [0, 0.1, 0], - [0, 0.5, 0], - [0, 0.2, 0], - [0.6, 0, 0], - [0.5, 0, 0.4], - [0.3, 0, 0.2], - [0, 0.1, 0.3], - [0.1, 0, 0], - [0.1, 0, 0], - [0, 0, 0.1], - [0, 0, 0], - ] - ), - 1e-4, - ) - - -def test__convolution__cross_mask_with_blurring_entries__returns_array(): - cross_mask = aa.Mask2D( - mask=[ - [True, True, True, True, True], - [True, True, False, True, True], - [True, False, False, False, True], - [True, True, False, True, True], - [True, True, True, True, True], - ], - pixel_scales=0.1, - ) - - kernel = aa.Kernel2D.no_mask( - values=[[0, 0.2, 0], [0.2, 0.4, 0.2], [0, 0.2, 0]], pixel_scales=0.1 - ) - - convolver = aa.Convolver(mask=cross_mask, kernel=kernel) - - image_array = aa.Array2D(values=[1, 0, 0, 0, 0], mask=cross_mask) - - blurring_mask = cross_mask.derive_mask.blurring_from( - kernel_shape_native=kernel.shape_native - ) - - blurring_array = aa.Array2D( - values=[0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], mask=blurring_mask - ) - - result = convolver.convolve_image(image=image_array, blurring_image=blurring_array) - - assert (np.round(result, 1) == np.array([0.6, 0.2, 0.2, 0.0, 0.0])).all() - - -def test__compare_to_full_2d_convolution(): - # Setup a blurred data, using the PSF to perform the convolution in 2D, then masks it to make a 1d array. - - import scipy.signal - - mask = aa.Mask2D.circular( - shape_native=(30, 30), pixel_scales=(1.0, 1.0), radius=4.0 - ) - kernel = aa.Kernel2D.no_mask(values=np.arange(49).reshape(7, 7), pixel_scales=1.0) - image = aa.Array2D.no_mask(values=np.arange(900).reshape(30, 30), pixel_scales=1.0) - - blurred_image_via_scipy = scipy.signal.convolve2d( - image.native, kernel.native, mode="same" - ) - blurred_image_via_scipy = aa.Array2D.no_mask( - values=blurred_image_via_scipy, pixel_scales=1.0 - ) - blurred_masked_image_via_scipy = aa.Array2D( - values=blurred_image_via_scipy.native, mask=mask - ) - - # Now reproduce this data using the frame convolver_image - - masked_image = aa.Array2D(values=image.native, mask=mask) - - blurring_mask = mask.derive_mask.blurring_from( - kernel_shape_native=kernel.shape_native - ) - - convolver = aa.Convolver(mask=mask, kernel=kernel) - - blurring_image = aa.Array2D(values=image.native, mask=blurring_mask) - - blurred_masked_im_1 = convolver.convolve_image( - image=masked_image, blurring_image=blurring_image - ) - - assert blurred_masked_image_via_scipy == pytest.approx(blurred_masked_im_1, 1e-4) - - -def test__compare_to_full_2d_convolution__no_blurring_image(): - # Setup a blurred data, using the PSF to perform the convolution in 2D, then masks it to make a 1d array. - - import scipy.signal - - mask = aa.Mask2D.circular( - shape_native=(30, 30), pixel_scales=(1.0, 1.0), radius=4.0 - ) - kernel = aa.Kernel2D.no_mask(values=np.arange(49).reshape(7, 7), pixel_scales=1.0) - image = aa.Array2D.no_mask(values=np.arange(900).reshape(30, 30), pixel_scales=1.0) - - blurring_mask = mask.derive_mask.blurring_from( - kernel_shape_native=kernel.shape_native - ) - blurred_image_via_scipy = scipy.signal.convolve2d( - image.native * blurring_mask, kernel.native, mode="same" - ) - blurred_image_via_scipy = aa.Array2D.no_mask( - values=blurred_image_via_scipy, pixel_scales=1.0 - ) - blurred_masked_image_via_scipy = aa.Array2D( - values=blurred_image_via_scipy.native, mask=mask - ) - - # Now reproduce this data using the frame convolver_image - - masked_image = aa.Array2D(values=image.native, mask=mask) - - convolver = aa.Convolver(mask=mask, kernel=kernel) - - blurred_masked_im_1 = convolver.convolve_image_no_blurring(image=masked_image) - - assert blurred_masked_image_via_scipy == pytest.approx(blurred_masked_im_1, 1e-4) - - -def test__summed_convolved_array_from(): - mask = aa.Mask2D( - mask=[ - [True, True, True, True, True], - [True, True, True, True, True], - [True, False, False, False, True], - [True, True, True, True, True], - [True, True, True, True, True], - ], - pixel_scales=0.1, - ) - - kernel = aa.Kernel2D.no_mask( - values=[[0, 0.0, 0], [0.5, 1.0, 0.5], [0, 0.0, 0]], pixel_scales=0.1 - ) - - convolver = aa.Convolver(mask=mask, kernel=kernel) - - image_array = aa.Array2D(values=[1.0, 2.0, 3.0], mask=mask) - - summed_convolved_array = convolver.convolve_image_no_blurring(image=image_array) - - assert summed_convolved_array == pytest.approx(np.array([2.0, 4.0, 4.0]), 1.0e-4) diff --git a/test_autoarray/structures/arrays/test_kernel_2d.py b/test_autoarray/structures/arrays/test_kernel_2d.py index 9106b43a1..23f66de19 100644 --- a/test_autoarray/structures/arrays/test_kernel_2d.py +++ b/test_autoarray/structures/arrays/test_kernel_2d.py @@ -248,6 +248,58 @@ def test__rescaled_with_odd_dimensions_from__different_scalings(): assert (kernel_2d.native == (1.0 / 15.0) * np.ones((5, 3))).all() +def test__from_as_gaussian_via_alma_fits_header_parameters__identical_to_astropy_gaussian_model(): + pixel_scales = 0.1 + + x_stddev = ( + 1.0e-5 + * (units.deg).to(units.arcsec) + / pixel_scales + / (2.0 * np.sqrt(2.0 * np.log(2.0))) + ) + y_stddev = ( + 2.0e-5 + * (units.deg).to(units.arcsec) + / pixel_scales + / (2.0 * np.sqrt(2.0 * np.log(2.0))) + ) + + theta_deg = 230.0 + theta = Angle(theta_deg, "deg").radian + + gaussian_astropy = functional_models.Gaussian2D( + amplitude=1.0, + x_mean=1.0, + y_mean=1.0, + x_stddev=x_stddev, + y_stddev=y_stddev, + theta=theta, + ) + + shape = (3, 3) + y, x = np.mgrid[0 : shape[1], 0 : shape[0]] + kernel_astropy = gaussian_astropy(x, y) + kernel_astropy /= np.sum(kernel_astropy) + + kernel_2d = aa.Kernel2D.from_as_gaussian_via_alma_fits_header_parameters( + shape_native=shape, + pixel_scales=pixel_scales, + y_stddev=2.0e-5, + x_stddev=1.0e-5, + theta=theta_deg, + normalize=True, + ) + + assert kernel_astropy == pytest.approx(kernel_2d.native._array, abs=1e-4) + + +def test__convolved_array_from__not_odd_x_odd_kernel__raises_error(): + kernel_2d = aa.Kernel2D.no_mask(values=[[0.0, 1.0], [1.0, 2.0]], pixel_scales=1.0) + + with pytest.raises(exc.KernelException): + kernel_2d.convolved_array_from(np.ones((5, 5))) + + def test__convolved_array_from(): array_2d = aa.Array2D.no_mask( values=[[0.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 0.0]], pixel_scales=1.0 @@ -408,53 +460,261 @@ def test__convolved_array_from(): ).all() -def test__convolved_array_from__not_odd_x_odd_kernel__raises_error(): - kernel_2d = aa.Kernel2D.no_mask(values=[[0.0, 1.0], [1.0, 2.0]], pixel_scales=1.0) +def test__convolve_mapping_matrix__asymetric_convolver__matrix_blurred_correctly(): + mask = np.array( + [ + [True, True, True, True, True, True], + [True, False, False, False, False, True], + [True, False, False, False, False, True], + [True, False, False, False, False, True], + [True, False, False, False, False, True], + [True, True, True, True, True, True], + ] + ) - with pytest.raises(exc.KernelException): - kernel_2d.convolved_array_from(np.ones((5, 5))) + asymmetric_kernel = aa.Kernel2D.no_mask( + values=[[0, 0.0, 0], [0.4, 0.2, 0.3], [0, 0.1, 0]], pixel_scales=1.0 + ) + convolver = aa.Convolver(mask=mask, kernel=asymmetric_kernel) -def test__from_as_gaussian_via_alma_fits_header_parameters__identical_to_astropy_gaussian_model(): - pixel_scales = 0.1 + mapping = np.array( + [ + [0, 0, 0], + [0, 0, 0], + [0, 0, 0], + [0, 0, 0], + [0, 0, 0], + [0, 0, 0], + [0, 0, 0], + [ + 0, + 1, + 0, + ], # The 0.3 should be 'chopped' from this pixel as it is on the right-most edge + [0, 0, 0], + [1, 0, 0], + [0, 0, 1], + [0, 0, 0], + [0, 0, 0], + [0, 0, 0], + [0, 0, 0], + [0, 0, 0], + ] + ) + + blurred_mapping = convolver.convolve_mapping_matrix(mapping) - x_stddev = ( - 1.0e-5 - * (units.deg).to(units.arcsec) - / pixel_scales - / (2.0 * np.sqrt(2.0 * np.log(2.0))) + assert ( + blurred_mapping + == np.array( + [ + [0, 0, 0], + [0, 0, 0], + [0, 0, 0], + [0, 0, 0], + [0, 0, 0], + [0, 0, 0], + [0, 0.4, 0], + [0, 0.2, 0], + [0.4, 0, 0], + [0.2, 0, 0.4], + [0.3, 0, 0.2], + [0, 0.1, 0.3], + [0, 0, 0], + [0.1, 0, 0], + [0, 0, 0.1], + [0, 0, 0], + ] + ) + ).all() + + asymmetric_kernel = aa.Kernel2D.no_mask( + values=[[0, 0.0, 0], [0.4, 0.2, 0.3], [0, 0.1, 0]], pixel_scales=1.0 ) - y_stddev = ( - 2.0e-5 - * (units.deg).to(units.arcsec) - / pixel_scales - / (2.0 * np.sqrt(2.0 * np.log(2.0))) + + convolver = aa.Convolver(mask=mask, kernel=asymmetric_kernel) + + mapping = np.array( + [ + [0, 1, 0], + [0, 1, 0], + [0, 1, 0], + [0, 0, 0], + [0, 0, 0], + [0, 0, 0], + [0, 0, 0], + [ + 0, + 1, + 0, + ], # The 0.3 should be 'chopped' from this pixel as it is on the right-most edge + [1, 0, 0], + [1, 0, 0], + [0, 0, 1], + [0, 0, 0], + [0, 0, 0], + [0, 0, 0], + [0, 0, 0], + [0, 0, 0], + ] + ) + + blurred_mapping = convolver.convolve_mapping_matrix(mapping) + + assert blurred_mapping == pytest.approx( + np.array( + [ + [0, 0.6, 0], + [0, 0.9, 0], + [0, 0.5, 0], + [0, 0.3, 0], + [0, 0.1, 0], + [0, 0.1, 0], + [0, 0.5, 0], + [0, 0.2, 0], + [0.6, 0, 0], + [0.5, 0, 0.4], + [0.3, 0, 0.2], + [0, 0.1, 0.3], + [0.1, 0, 0], + [0.1, 0, 0], + [0, 0, 0.1], + [0, 0, 0], + ] + ), + 1e-4, ) - theta_deg = 230.0 - theta = Angle(theta_deg, "deg").radian - gaussian_astropy = functional_models.Gaussian2D( - amplitude=1.0, - x_mean=1.0, - y_mean=1.0, - x_stddev=x_stddev, - y_stddev=y_stddev, - theta=theta, +def test__convolution__cross_mask_with_blurring_entries__returns_array(): + cross_mask = aa.Mask2D( + mask=[ + [True, True, True, True, True], + [True, True, False, True, True], + [True, False, False, False, True], + [True, True, False, True, True], + [True, True, True, True, True], + ], + pixel_scales=0.1, ) - shape = (3, 3) - y, x = np.mgrid[0 : shape[1], 0 : shape[0]] - kernel_astropy = gaussian_astropy(x, y) - kernel_astropy /= np.sum(kernel_astropy) + kernel = aa.Kernel2D.no_mask( + values=[[0, 0.2, 0], [0.2, 0.4, 0.2], [0, 0.2, 0]], pixel_scales=0.1 + ) - kernel_2d = aa.Kernel2D.from_as_gaussian_via_alma_fits_header_parameters( - shape_native=shape, - pixel_scales=pixel_scales, - y_stddev=2.0e-5, - x_stddev=1.0e-5, - theta=theta_deg, - normalize=True, + convolver = aa.Convolver(mask=cross_mask, kernel=kernel) + + image_array = aa.Array2D(values=[1, 0, 0, 0, 0], mask=cross_mask) + + blurring_mask = cross_mask.derive_mask.blurring_from( + kernel_shape_native=kernel.shape_native ) - assert kernel_astropy == pytest.approx(kernel_2d.native._array, abs=1e-4) + blurring_array = aa.Array2D( + values=[0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], mask=blurring_mask + ) + + result = convolver.convolve_image(image=image_array, blurring_image=blurring_array) + + assert (np.round(result, 1) == np.array([0.6, 0.2, 0.2, 0.0, 0.0])).all() + + +def test__compare_to_full_2d_convolution(): + # Setup a blurred data, using the PSF to perform the convolution in 2D, then masks it to make a 1d array. + + import scipy.signal + + mask = aa.Mask2D.circular( + shape_native=(30, 30), pixel_scales=(1.0, 1.0), radius=4.0 + ) + kernel = aa.Kernel2D.no_mask(values=np.arange(49).reshape(7, 7), pixel_scales=1.0) + image = aa.Array2D.no_mask(values=np.arange(900).reshape(30, 30), pixel_scales=1.0) + + blurred_image_via_scipy = scipy.signal.convolve2d( + image.native, kernel.native, mode="same" + ) + blurred_image_via_scipy = aa.Array2D.no_mask( + values=blurred_image_via_scipy, pixel_scales=1.0 + ) + blurred_masked_image_via_scipy = aa.Array2D( + values=blurred_image_via_scipy.native, mask=mask + ) + + # Now reproduce this data using the frame convolver_image + + masked_image = aa.Array2D(values=image.native, mask=mask) + + blurring_mask = mask.derive_mask.blurring_from( + kernel_shape_native=kernel.shape_native + ) + + convolver = aa.Convolver(mask=mask, kernel=kernel) + + blurring_image = aa.Array2D(values=image.native, mask=blurring_mask) + + blurred_masked_im_1 = convolver.convolve_image( + image=masked_image, blurring_image=blurring_image + ) + + assert blurred_masked_image_via_scipy == pytest.approx(blurred_masked_im_1, 1e-4) + + +def test__compare_to_full_2d_convolution__no_blurring_image(): + # Setup a blurred data, using the PSF to perform the convolution in 2D, then masks it to make a 1d array. + + import scipy.signal + + mask = aa.Mask2D.circular( + shape_native=(30, 30), pixel_scales=(1.0, 1.0), radius=4.0 + ) + kernel = aa.Kernel2D.no_mask(values=np.arange(49).reshape(7, 7), pixel_scales=1.0) + image = aa.Array2D.no_mask(values=np.arange(900).reshape(30, 30), pixel_scales=1.0) + + blurring_mask = mask.derive_mask.blurring_from( + kernel_shape_native=kernel.shape_native + ) + blurred_image_via_scipy = scipy.signal.convolve2d( + image.native * blurring_mask, kernel.native, mode="same" + ) + blurred_image_via_scipy = aa.Array2D.no_mask( + values=blurred_image_via_scipy, pixel_scales=1.0 + ) + blurred_masked_image_via_scipy = aa.Array2D( + values=blurred_image_via_scipy.native, mask=mask + ) + + # Now reproduce this data using the frame convolver_image + + masked_image = aa.Array2D(values=image.native, mask=mask) + + convolver = aa.Convolver(mask=mask, kernel=kernel) + + blurred_masked_im_1 = convolver.convolve_image_no_blurring(image=masked_image) + + assert blurred_masked_image_via_scipy == pytest.approx(blurred_masked_im_1, 1e-4) + + +def test__summed_convolved_array_from(): + mask = aa.Mask2D( + mask=[ + [True, True, True, True, True], + [True, True, True, True, True], + [True, False, False, False, True], + [True, True, True, True, True], + [True, True, True, True, True], + ], + pixel_scales=0.1, + ) + + kernel = aa.Kernel2D.no_mask( + values=[[0, 0.0, 0], [0.5, 1.0, 0.5], [0, 0.0, 0]], pixel_scales=0.1 + ) + + convolver = aa.Convolver(mask=mask, kernel=kernel) + + image_array = aa.Array2D(values=[1.0, 2.0, 3.0], mask=mask) + + summed_convolved_array = convolver.convolve_image_no_blurring(image=image_array) + + assert summed_convolved_array == pytest.approx(np.array([2.0, 4.0, 4.0]), 1.0e-4) From 0b4763ce2d1838b97479e46eb800e319a3ae841d Mon Sep 17 00:00:00 2001 From: James Nightingale Date: Wed, 2 Apr 2025 19:17:04 +0100 Subject: [PATCH 02/11] cleaned up test_kernel_2d --- .../structures/arrays/test_kernel_2d.py | 181 +----------------- 1 file changed, 7 insertions(+), 174 deletions(-) diff --git a/test_autoarray/structures/arrays/test_kernel_2d.py b/test_autoarray/structures/arrays/test_kernel_2d.py index 23f66de19..bd52796a0 100644 --- a/test_autoarray/structures/arrays/test_kernel_2d.py +++ b/test_autoarray/structures/arrays/test_kernel_2d.py @@ -301,107 +301,6 @@ def test__convolved_array_from__not_odd_x_odd_kernel__raises_error(): def test__convolved_array_from(): - array_2d = aa.Array2D.no_mask( - values=[[0.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 0.0]], pixel_scales=1.0 - ) - - kernel_2d = aa.Kernel2D.no_mask( - values=[[0.0, 1.0, 0.0], [1.0, 2.0, 1.0], [0.0, 1.0, 0.0]], pixel_scales=1.0 - ) - - blurred_array_2d = kernel_2d.convolved_array_from(array_2d) - - assert (blurred_array_2d == kernel_2d).all() - - array_2d = aa.Array2D.no_mask( - values=[ - [0.0, 0.0, 0.0, 0.0], - [0.0, 1.0, 0.0, 0.0], - [0.0, 0.0, 0.0, 0.0], - [0.0, 0.0, 0.0, 0.0], - ], - pixel_scales=1.0, - ) - - kernel_2d = aa.Kernel2D.no_mask( - values=[[0.0, 1.0, 0.0], [1.0, 2.0, 1.0], [0.0, 1.0, 0.0]], pixel_scales=1.0 - ) - - blurred_array_2d = kernel_2d.convolved_array_from(array=array_2d) - - assert ( - blurred_array_2d.native - == np.array( - [ - [0.0, 1.0, 0.0, 0.0], - [1.0, 2.0, 1.0, 0.0], - [0.0, 1.0, 0.0, 0.0], - [0.0, 0.0, 0.0, 0.0], - ] - ) - ).all() - - array_2d = aa.Array2D.no_mask( - values=[[0.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]], - pixel_scales=1.0, - ) - - kernel_2d = aa.Kernel2D.no_mask( - values=[[0.0, 1.0, 0.0], [1.0, 2.0, 1.0], [0.0, 1.0, 0.0]], pixel_scales=1.0 - ) - - blurred_array_2d = kernel_2d.convolved_array_from(array_2d) - - assert ( - blurred_array_2d.native - == np.array( - [[0.0, 1.0, 0.0], [1.0, 2.0, 1.0], [0.0, 1.0, 0.0], [0.0, 0.0, 0.0]] - ) - ).all() - - array_2d = aa.Array2D.no_mask( - values=[[0.0, 0.0, 0.0, 0.0], [0.0, 1.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0]], - pixel_scales=1.0, - ) - - kernel_2d = aa.Kernel2D.no_mask( - values=[[0.0, 1.0, 0.0], [1.0, 2.0, 1.0], [0.0, 1.0, 0.0]], pixel_scales=1.0 - ) - - blurred_array_2d = kernel_2d.convolved_array_from(array_2d) - - assert ( - blurred_array_2d.native - == np.array([[0.0, 1.0, 0.0, 0.0], [1.0, 2.0, 1.0, 0.0], [0.0, 1.0, 0.0, 0.0]]) - ).all() - - array_2d = aa.Array2D.no_mask( - values=[ - [0.0, 0.0, 0.0, 0.0], - [0.0, 1.0, 0.0, 0.0], - [0.0, 0.0, 1.0, 0.0], - [0.0, 0.0, 0.0, 0.0], - ], - pixel_scales=1.0, - ) - - kernel_2d = aa.Kernel2D.no_mask( - values=[[1.0, 1.0, 1.0], [2.0, 2.0, 1.0], [1.0, 3.0, 3.0]], pixel_scales=1.0 - ) - - blurred_array_2d = kernel_2d.convolved_array_from(array_2d) - - assert ( - blurred_array_2d.native - == np.array( - [ - [1.0, 1.0, 1.0, 0.0], - [2.0, 3.0, 2.0, 1.0], - [1.0, 5.0, 5.0, 1.0], - [0.0, 1.0, 3.0, 3.0], - ] - ) - ).all() array_2d = aa.Array2D.no_mask( [ @@ -460,7 +359,7 @@ def test__convolved_array_from(): ).all() -def test__convolve_mapping_matrix__asymetric_convolver__matrix_blurred_correctly(): +def test__convolve_mapping_matrix(): mask = np.array( [ [True, True, True, True, True, True], @@ -472,12 +371,10 @@ def test__convolve_mapping_matrix__asymetric_convolver__matrix_blurred_correctly ] ) - asymmetric_kernel = aa.Kernel2D.no_mask( + kernel = aa.Kernel2D.no_mask( values=[[0, 0.0, 0], [0.4, 0.2, 0.3], [0, 0.1, 0]], pixel_scales=1.0 ) - convolver = aa.Convolver(mask=mask, kernel=asymmetric_kernel) - mapping = np.array( [ [0, 0, 0], @@ -503,7 +400,7 @@ def test__convolve_mapping_matrix__asymetric_convolver__matrix_blurred_correctly ] ) - blurred_mapping = convolver.convolve_mapping_matrix(mapping) + blurred_mapping = kernel.convolve_mapping_matrix(mapping) assert ( blurred_mapping @@ -529,12 +426,10 @@ def test__convolve_mapping_matrix__asymetric_convolver__matrix_blurred_correctly ) ).all() - asymmetric_kernel = aa.Kernel2D.no_mask( + kernel = aa.Kernel2D.no_mask( values=[[0, 0.0, 0], [0.4, 0.2, 0.3], [0, 0.1, 0]], pixel_scales=1.0 ) - convolver = aa.Convolver(mask=mask, kernel=asymmetric_kernel) - mapping = np.array( [ [0, 1, 0], @@ -560,7 +455,7 @@ def test__convolve_mapping_matrix__asymetric_convolver__matrix_blurred_correctly ] ) - blurred_mapping = convolver.convolve_mapping_matrix(mapping) + blurred_mapping = kernel.convolve_mapping_matrix(mapping) assert blurred_mapping == pytest.approx( np.array( @@ -587,39 +482,6 @@ def test__convolve_mapping_matrix__asymetric_convolver__matrix_blurred_correctly ) -def test__convolution__cross_mask_with_blurring_entries__returns_array(): - cross_mask = aa.Mask2D( - mask=[ - [True, True, True, True, True], - [True, True, False, True, True], - [True, False, False, False, True], - [True, True, False, True, True], - [True, True, True, True, True], - ], - pixel_scales=0.1, - ) - - kernel = aa.Kernel2D.no_mask( - values=[[0, 0.2, 0], [0.2, 0.4, 0.2], [0, 0.2, 0]], pixel_scales=0.1 - ) - - convolver = aa.Convolver(mask=cross_mask, kernel=kernel) - - image_array = aa.Array2D(values=[1, 0, 0, 0, 0], mask=cross_mask) - - blurring_mask = cross_mask.derive_mask.blurring_from( - kernel_shape_native=kernel.shape_native - ) - - blurring_array = aa.Array2D( - values=[0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], mask=blurring_mask - ) - - result = convolver.convolve_image(image=image_array, blurring_image=blurring_array) - - assert (np.round(result, 1) == np.array([0.6, 0.2, 0.2, 0.0, 0.0])).all() - - def test__compare_to_full_2d_convolution(): # Setup a blurred data, using the PSF to perform the convolution in 2D, then masks it to make a 1d array. @@ -649,11 +511,9 @@ def test__compare_to_full_2d_convolution(): kernel_shape_native=kernel.shape_native ) - convolver = aa.Convolver(mask=mask, kernel=kernel) - blurring_image = aa.Array2D(values=image.native, mask=blurring_mask) - blurred_masked_im_1 = convolver.convolve_image( + blurred_masked_im_1 = kernel.convolve_image( image=masked_image, blurring_image=blurring_image ) @@ -688,33 +548,6 @@ def test__compare_to_full_2d_convolution__no_blurring_image(): masked_image = aa.Array2D(values=image.native, mask=mask) - convolver = aa.Convolver(mask=mask, kernel=kernel) - - blurred_masked_im_1 = convolver.convolve_image_no_blurring(image=masked_image) + blurred_masked_im_1 = kernel.convolve_image_no_blurring(image=masked_image) assert blurred_masked_image_via_scipy == pytest.approx(blurred_masked_im_1, 1e-4) - - -def test__summed_convolved_array_from(): - mask = aa.Mask2D( - mask=[ - [True, True, True, True, True], - [True, True, True, True, True], - [True, False, False, False, True], - [True, True, True, True, True], - [True, True, True, True, True], - ], - pixel_scales=0.1, - ) - - kernel = aa.Kernel2D.no_mask( - values=[[0, 0.0, 0], [0.5, 1.0, 0.5], [0, 0.0, 0]], pixel_scales=0.1 - ) - - convolver = aa.Convolver(mask=mask, kernel=kernel) - - image_array = aa.Array2D(values=[1.0, 2.0, 3.0], mask=mask) - - summed_convolved_array = convolver.convolve_image_no_blurring(image=image_array) - - assert summed_convolved_array == pytest.approx(np.array([2.0, 4.0, 4.0]), 1.0e-4) From 5762fcf320794b12f66359d5d29a56ddc3eb130c Mon Sep 17 00:00:00 2001 From: James Nightingale Date: Wed, 2 Apr 2025 19:19:43 +0100 Subject: [PATCH 03/11] remopve convolved_array_with_mask_From --- autoarray/structures/arrays/kernel_2d.py | 31 ------------------------ 1 file changed, 31 deletions(-) diff --git a/autoarray/structures/arrays/kernel_2d.py b/autoarray/structures/arrays/kernel_2d.py index 09ff331ef..d384f8a4b 100644 --- a/autoarray/structures/arrays/kernel_2d.py +++ b/autoarray/structures/arrays/kernel_2d.py @@ -483,37 +483,6 @@ def convolved_array_from(self, array: Array2D) -> Array2D: return Array2D(values=convolved_array_1d, mask=array_2d.mask) - def convolved_array_with_mask_from(self, array: Array2D, mask: Mask2D) -> Array2D: - """ - Convolve an array with this Kernel2D - - Parameters - ---------- - image - An array representing the image the Kernel2D is convolved with. - - Returns - ------- - convolved_image - An array representing the image after convolution. - - Raises - ------ - KernelException if either Kernel2D psf dimension is odd - """ - - if self.mask.shape[0] % 2 == 0 or self.mask.shape[1] % 2 == 0: - raise exc.KernelException("Kernel2D Kernel2D must be odd") - - convolved_array_2d = scipy.signal.convolve2d(array, self.native, mode="same") - - convolved_array_1d = array_2d_util.array_2d_slim_from( - mask_2d=np.array(mask), - array_2d_native=np.array(convolved_array_2d), - ) - - return Array2D(values=convolved_array_1d, mask=mask) - def jax_convolve(self, image, blurring_image, method="auto"): slim_to_2D_index_image = jnp.nonzero( jnp.logical_not(self.mask.array), size=image.shape[0] From 6fa35a6f838fe89e4fb0209bc25c9f960e99e152 Mon Sep 17 00:00:00 2001 From: James Nightingale Date: Wed, 2 Apr 2025 19:21:20 +0100 Subject: [PATCH 04/11] simplify jax_convolve --- autoarray/structures/arrays/kernel_2d.py | 17 ++++++++++++----- 1 file changed, 12 insertions(+), 5 deletions(-) diff --git a/autoarray/structures/arrays/kernel_2d.py b/autoarray/structures/arrays/kernel_2d.py index d384f8a4b..def73f067 100644 --- a/autoarray/structures/arrays/kernel_2d.py +++ b/autoarray/structures/arrays/kernel_2d.py @@ -484,22 +484,29 @@ def convolved_array_from(self, array: Array2D) -> Array2D: return Array2D(values=convolved_array_1d, mask=array_2d.mask) def jax_convolve(self, image, blurring_image, method="auto"): - slim_to_2D_index_image = jnp.nonzero( + + slim_to_native = jnp.nonzero( jnp.logical_not(self.mask.array), size=image.shape[0] ) - slim_to_2D_index_blurring = jnp.nonzero( + slim_to_native_blurring = jnp.nonzero( jnp.logical_not(self.blurring_mask), size=blurring_image.shape[0] ) + expanded_image_native = jnp.zeros(self.mask.shape) - expanded_image_native = expanded_image_native.at[slim_to_2D_index_image].set( + + expanded_image_native = expanded_image_native.at[slim_to_native].set( image.array ) - expanded_image_native = expanded_image_native.at[slim_to_2D_index_blurring].set( + expanded_image_native = expanded_image_native.at[slim_to_native_blurring].set( blurring_image.array ) + kernel = np.array(self.kernel.native.array) + convolve_native = jax.scipy.signal.convolve( expanded_image_native, kernel, mode="same", method=method ) - convolve_slim = convolve_native[slim_to_2D_index_image] + + convolve_slim = convolve_native[slim_to_native] + return convolve_slim \ No newline at end of file From 0386bddb1860869c91766b3ead0672f3b87418e7 Mon Sep 17 00:00:00 2001 From: James Nightingale Date: Wed, 2 Apr 2025 19:47:45 +0100 Subject: [PATCH 05/11] test__convolve_image --- autoarray/structures/arrays/kernel_2d.py | 66 ++++++++++++++++--- .../structures/arrays/test_kernel_2d.py | 52 ++++++++++++--- 2 files changed, 98 insertions(+), 20 deletions(-) diff --git a/autoarray/structures/arrays/kernel_2d.py b/autoarray/structures/arrays/kernel_2d.py index def73f067..28ee92725 100644 --- a/autoarray/structures/arrays/kernel_2d.py +++ b/autoarray/structures/arrays/kernel_2d.py @@ -8,7 +8,6 @@ from autoconf.fitsable import header_obj_from -from autoarray.mask.mask_2d import Mask2D from autoarray.structures.arrays.uniform_2d import AbstractArray2D from autoarray.structures.arrays.uniform_2d import Array2D from autoarray.structures.grids.uniform_2d import Grid2D @@ -17,6 +16,7 @@ from autoarray import exc from autoarray import type as ty from autoarray.structures.arrays import array_2d_util +from autoarray.mask.mask_2d import mask_2d_util class Kernel2D(AbstractArray2D): @@ -483,28 +483,74 @@ def convolved_array_from(self, array: Array2D) -> Array2D: return Array2D(values=convolved_array_1d, mask=array_2d.mask) - def jax_convolve(self, image, blurring_image, method="auto"): + def convolve_image(self, image, blurring_image, jax_method="fft"): + """ + For a given 1D array and blurring array, convolve the two using this convolver. + + Parameters + ---------- + image + 1D array of the values which are to be blurred with the convolver's PSF. + blurring_image + 1D array of the blurring values which blur into the array after PSF convolution. + jax_method + If JAX is enabled this keyword will indicate what method is used for the PSF + convolution. Can be either `direct` to calculate it in real space or `fft` + to calculated it via a fast Fourier transform. `fft` is typically faster for + kernels that are more than about 5x5. Default is `fft`. + """ + + if self.mask.shape[0] % 2 == 0 or self.mask.shape[1] % 2 == 0: + raise exc.KernelException("Kernel2D Kernel2D must be odd") + + print(type(image.native + blurring_image.native)) + print(type(self.native)) + + convolved_array_2d = scipy.signal.convolve2d((image.native + blurring_image.native)._array, self.native._array, mode="same") + convolved_array_1d = array_2d_util.array_2d_slim_from( + mask_2d=np.array(image.mask), + array_2d_native=convolved_array_2d, + ) + + return Array2D(values=convolved_array_1d, mask=image.mask) + + def convolve_image_jax_from(self, array, blurring_array, method="auto"): + """ + For a given 1D array and blurring array, convolve the two using this convolver. + + Parameters + ---------- + array + 1D array of the values which are to be blurred with the convolver's PSF. + blurring_array + 1D array of the blurring values which blur into the array after PSF convolution. + jax_method + If JAX is enabled this keyword will indicate what method is used for the PSF + convolution. Can be either `direct` to calculate it in real space or `fft` + to calculated it via a fast Fourier transform. `fft` is typically faster for + kernels that are more than about 5x5. Default is `fft`. + """ slim_to_native = jnp.nonzero( - jnp.logical_not(self.mask.array), size=image.shape[0] + jnp.logical_not(self.mask.array), size=array.shape[0] ) slim_to_native_blurring = jnp.nonzero( - jnp.logical_not(self.blurring_mask), size=blurring_image.shape[0] + jnp.logical_not(self.blurring_mask), size=blurring_array.shape[0] ) - expanded_image_native = jnp.zeros(self.mask.shape) + expanded_array_native = jnp.zeros(self.mask.shape) - expanded_image_native = expanded_image_native.at[slim_to_native].set( - image.array + expanded_array_native = expanded_array_native.at[slim_to_native].set( + array.array ) - expanded_image_native = expanded_image_native.at[slim_to_native_blurring].set( - blurring_image.array + expanded_array_native = expanded_array_native.at[slim_to_native_blurring].set( + blurring_array.array ) kernel = np.array(self.kernel.native.array) convolve_native = jax.scipy.signal.convolve( - expanded_image_native, kernel, mode="same", method=method + expanded_array_native, kernel, mode="same", method=method ) convolve_slim = convolve_native[slim_to_native] diff --git a/test_autoarray/structures/arrays/test_kernel_2d.py b/test_autoarray/structures/arrays/test_kernel_2d.py index bd52796a0..50c931ae5 100644 --- a/test_autoarray/structures/arrays/test_kernel_2d.py +++ b/test_autoarray/structures/arrays/test_kernel_2d.py @@ -1,11 +1,10 @@ -from astropy.io import fits from astropy import units from astropy.modeling import functional_models from astropy.coordinates import Angle +import jax.numpy as jnp import numpy as np import pytest from os import path -import os import autoarray as aa from autoarray import exc @@ -359,6 +358,36 @@ def test__convolved_array_from(): ).all() +def test__convolved_array_from__input_jax_array(): + + array_2d = jnp.array( + [ + [0.0, 0.0, 0.0, 0.0], + [1.0, 0.0, 0.0, 0.0], + [0.0, 0.0, 0.0, 1.0], + [0.0, 0.0, 0.0, 0.0], + ]) + + kernel_2d = aa.Kernel2D.no_mask( + values=[[1.0, 1.0, 1.0], [2.0, 2.0, 1.0], [1.0, 3.0, 3.0]], pixel_scales=1.0 + ) + + blurred_array_2d = kernel_2d.convolved_array_from(array_2d) + + assert ( + blurred_array_2d.native + == np.array( + [ + [1.0, 1.0, 0.0, 0.0], + [2.0, 1.0, 1.0, 1.0], + [3.0, 3.0, 2.0, 2.0], + [0.0, 0.0, 1.0, 3.0], + ] + ) + ).all() + + + def test__convolve_mapping_matrix(): mask = np.array( [ @@ -482,19 +511,19 @@ def test__convolve_mapping_matrix(): ) -def test__compare_to_full_2d_convolution(): - # Setup a blurred data, using the PSF to perform the convolution in 2D, then masks it to make a 1d array. - - import scipy.signal +def test__convolve_image(): mask = aa.Mask2D.circular( shape_native=(30, 30), pixel_scales=(1.0, 1.0), radius=4.0 ) - kernel = aa.Kernel2D.no_mask(values=np.arange(49).reshape(7, 7), pixel_scales=1.0) - image = aa.Array2D.no_mask(values=np.arange(900).reshape(30, 30), pixel_scales=1.0) + + import scipy.signal + + kernel = np.arange(49).reshape(7, 7) + image = np.arange(900).reshape(30, 30) blurred_image_via_scipy = scipy.signal.convolve2d( - image.native, kernel.native, mode="same" + image, kernel, mode="same" ) blurred_image_via_scipy = aa.Array2D.no_mask( values=blurred_image_via_scipy, pixel_scales=1.0 @@ -503,7 +532,10 @@ def test__compare_to_full_2d_convolution(): values=blurred_image_via_scipy.native, mask=mask ) - # Now reproduce this data using the frame convolver_image + # Now reproduce this data using the convolve_image function + + image = aa.Array2D.no_mask(values=np.arange(900).reshape(30, 30), pixel_scales=1.0) + kernel = aa.Kernel2D.no_mask(values=np.arange(49).reshape(7, 7), pixel_scales=1.0) masked_image = aa.Array2D(values=image.native, mask=mask) From 537b5ef8683d4a4ef4168fc0e9583dc94604c787 Mon Sep 17 00:00:00 2001 From: James Nightingale Date: Wed, 2 Apr 2025 19:54:04 +0100 Subject: [PATCH 06/11] convolve_image now only uses JAX --- autoarray/structures/arrays/kernel_2d.py | 52 +++++-------------- .../structures/arrays/test_kernel_2d.py | 5 +- 2 files changed, 16 insertions(+), 41 deletions(-) diff --git a/autoarray/structures/arrays/kernel_2d.py b/autoarray/structures/arrays/kernel_2d.py index 28ee92725..3f34c49bc 100644 --- a/autoarray/structures/arrays/kernel_2d.py +++ b/autoarray/structures/arrays/kernel_2d.py @@ -54,6 +54,9 @@ def __init__( store_native=store_native, ) + if self.mask.shape[0] % 2 == 0 or self.mask.shape[1] % 2 == 0: + raise exc.KernelException("Kernel2D Kernel2D must be odd") + if normalize: self._array = np.divide(self._array, np.sum(self._array)) @@ -500,59 +503,28 @@ def convolve_image(self, image, blurring_image, jax_method="fft"): kernels that are more than about 5x5. Default is `fft`. """ - if self.mask.shape[0] % 2 == 0 or self.mask.shape[1] % 2 == 0: - raise exc.KernelException("Kernel2D Kernel2D must be odd") - - print(type(image.native + blurring_image.native)) - print(type(self.native)) - - convolved_array_2d = scipy.signal.convolve2d((image.native + blurring_image.native)._array, self.native._array, mode="same") - - convolved_array_1d = array_2d_util.array_2d_slim_from( - mask_2d=np.array(image.mask), - array_2d_native=convolved_array_2d, - ) - - return Array2D(values=convolved_array_1d, mask=image.mask) - - def convolve_image_jax_from(self, array, blurring_array, method="auto"): - """ - For a given 1D array and blurring array, convolve the two using this convolver. - - Parameters - ---------- - array - 1D array of the values which are to be blurred with the convolver's PSF. - blurring_array - 1D array of the blurring values which blur into the array after PSF convolution. - jax_method - If JAX is enabled this keyword will indicate what method is used for the PSF - convolution. Can be either `direct` to calculate it in real space or `fft` - to calculated it via a fast Fourier transform. `fft` is typically faster for - kernels that are more than about 5x5. Default is `fft`. - """ slim_to_native = jnp.nonzero( - jnp.logical_not(self.mask.array), size=array.shape[0] + jnp.logical_not(image.mask.array), size=image.shape[0] ) slim_to_native_blurring = jnp.nonzero( - jnp.logical_not(self.blurring_mask), size=blurring_array.shape[0] + jnp.logical_not(blurring_image.mask.array), size=blurring_image.shape[0] ) - expanded_array_native = jnp.zeros(self.mask.shape) + expanded_array_native = jnp.zeros(image.mask.shape) expanded_array_native = expanded_array_native.at[slim_to_native].set( - array.array + image.array ) expanded_array_native = expanded_array_native.at[slim_to_native_blurring].set( - blurring_array.array + blurring_image.array ) - kernel = np.array(self.kernel.native.array) + kernel = np.array(self.native.array) convolve_native = jax.scipy.signal.convolve( - expanded_array_native, kernel, mode="same", method=method + expanded_array_native, kernel, mode="same", method=jax_method ) - convolve_slim = convolve_native[slim_to_native] + convolved_array_1d = convolve_native[slim_to_native] - return convolve_slim \ No newline at end of file + return Array2D(values=convolved_array_1d, mask=image.mask) \ No newline at end of file diff --git a/test_autoarray/structures/arrays/test_kernel_2d.py b/test_autoarray/structures/arrays/test_kernel_2d.py index 50c931ae5..11dab836d 100644 --- a/test_autoarray/structures/arrays/test_kernel_2d.py +++ b/test_autoarray/structures/arrays/test_kernel_2d.py @@ -549,7 +549,10 @@ def test__convolve_image(): image=masked_image, blurring_image=blurring_image ) - assert blurred_masked_image_via_scipy == pytest.approx(blurred_masked_im_1, 1e-4) + print(blurred_masked_image_via_scipy) + print(blurred_masked_im_1) + + assert blurred_masked_image_via_scipy == pytest.approx(blurred_masked_im_1.array, 1e-4) def test__compare_to_full_2d_convolution__no_blurring_image(): From 27bf06a9b15b6a20755ffef749a0a05a2dcef70a Mon Sep 17 00:00:00 2001 From: James Nightingale Date: Wed, 2 Apr 2025 19:56:18 +0100 Subject: [PATCH 07/11] fix test on array shape --- .../structures/arrays/test_kernel_2d.py | 170 +++++++----------- 1 file changed, 68 insertions(+), 102 deletions(-) diff --git a/test_autoarray/structures/arrays/test_kernel_2d.py b/test_autoarray/structures/arrays/test_kernel_2d.py index 11dab836d..d70854486 100644 --- a/test_autoarray/structures/arrays/test_kernel_2d.py +++ b/test_autoarray/structures/arrays/test_kernel_2d.py @@ -292,11 +292,10 @@ def test__from_as_gaussian_via_alma_fits_header_parameters__identical_to_astropy assert kernel_astropy == pytest.approx(kernel_2d.native._array, abs=1e-4) -def test__convolved_array_from__not_odd_x_odd_kernel__raises_error(): - kernel_2d = aa.Kernel2D.no_mask(values=[[0.0, 1.0], [1.0, 2.0]], pixel_scales=1.0) +def test__not_odd_x_odd_kernel__raises_error(): with pytest.raises(exc.KernelException): - kernel_2d.convolved_array_from(np.ones((5, 5))) + aa.Kernel2D.no_mask(values=[[0.0, 1.0], [1.0, 2.0]], pixel_scales=1.0) def test__convolved_array_from(): @@ -358,34 +357,78 @@ def test__convolved_array_from(): ).all() -def test__convolved_array_from__input_jax_array(): +def test__convolve_image(): - array_2d = jnp.array( - [ - [0.0, 0.0, 0.0, 0.0], - [1.0, 0.0, 0.0, 0.0], - [0.0, 0.0, 0.0, 1.0], - [0.0, 0.0, 0.0, 0.0], - ]) + mask = aa.Mask2D.circular( + shape_native=(30, 30), pixel_scales=(1.0, 1.0), radius=4.0 + ) - kernel_2d = aa.Kernel2D.no_mask( - values=[[1.0, 1.0, 1.0], [2.0, 2.0, 1.0], [1.0, 3.0, 3.0]], pixel_scales=1.0 + import scipy.signal + + kernel = np.arange(49).reshape(7, 7) + image = np.arange(900).reshape(30, 30) + + blurred_image_via_scipy = scipy.signal.convolve2d( + image, kernel, mode="same" + ) + blurred_image_via_scipy = aa.Array2D.no_mask( + values=blurred_image_via_scipy, pixel_scales=1.0 + ) + blurred_masked_image_via_scipy = aa.Array2D( + values=blurred_image_via_scipy.native, mask=mask ) - blurred_array_2d = kernel_2d.convolved_array_from(array_2d) + # Now reproduce this data using the convolve_image function - assert ( - blurred_array_2d.native - == np.array( - [ - [1.0, 1.0, 0.0, 0.0], - [2.0, 1.0, 1.0, 1.0], - [3.0, 3.0, 2.0, 2.0], - [0.0, 0.0, 1.0, 3.0], - ] - ) - ).all() + image = aa.Array2D.no_mask(values=np.arange(900).reshape(30, 30), pixel_scales=1.0) + kernel = aa.Kernel2D.no_mask(values=np.arange(49).reshape(7, 7), pixel_scales=1.0) + + masked_image = aa.Array2D(values=image.native, mask=mask) + + blurring_mask = mask.derive_mask.blurring_from( + kernel_shape_native=kernel.shape_native + ) + + blurring_image = aa.Array2D(values=image.native, mask=blurring_mask) + + blurred_masked_im_1 = kernel.convolve_image( + image=masked_image, blurring_image=blurring_image + ) + + assert blurred_masked_image_via_scipy == pytest.approx(blurred_masked_im_1.array, 1e-4) + + +def test__compare_to_full_2d_convolution__no_blurring_image(): + # Setup a blurred data, using the PSF to perform the convolution in 2D, then masks it to make a 1d array. + + import scipy.signal + + mask = aa.Mask2D.circular( + shape_native=(30, 30), pixel_scales=(1.0, 1.0), radius=4.0 + ) + kernel = aa.Kernel2D.no_mask(values=np.arange(49).reshape(7, 7), pixel_scales=1.0) + image = aa.Array2D.no_mask(values=np.arange(900).reshape(30, 30), pixel_scales=1.0) + + blurring_mask = mask.derive_mask.blurring_from( + kernel_shape_native=kernel.shape_native + ) + blurred_image_via_scipy = scipy.signal.convolve2d( + image.native * blurring_mask, kernel.native, mode="same" + ) + blurred_image_via_scipy = aa.Array2D.no_mask( + values=blurred_image_via_scipy, pixel_scales=1.0 + ) + blurred_masked_image_via_scipy = aa.Array2D( + values=blurred_image_via_scipy.native, mask=mask + ) + + # Now reproduce this data using the frame convolver_image + masked_image = aa.Array2D(values=image.native, mask=mask) + + blurred_masked_im_1 = kernel.convolve_image_no_blurring(image=masked_image) + + assert blurred_masked_image_via_scipy == pytest.approx(blurred_masked_im_1, 1e-4) def test__convolve_mapping_matrix(): @@ -509,80 +552,3 @@ def test__convolve_mapping_matrix(): ), 1e-4, ) - - -def test__convolve_image(): - - mask = aa.Mask2D.circular( - shape_native=(30, 30), pixel_scales=(1.0, 1.0), radius=4.0 - ) - - import scipy.signal - - kernel = np.arange(49).reshape(7, 7) - image = np.arange(900).reshape(30, 30) - - blurred_image_via_scipy = scipy.signal.convolve2d( - image, kernel, mode="same" - ) - blurred_image_via_scipy = aa.Array2D.no_mask( - values=blurred_image_via_scipy, pixel_scales=1.0 - ) - blurred_masked_image_via_scipy = aa.Array2D( - values=blurred_image_via_scipy.native, mask=mask - ) - - # Now reproduce this data using the convolve_image function - - image = aa.Array2D.no_mask(values=np.arange(900).reshape(30, 30), pixel_scales=1.0) - kernel = aa.Kernel2D.no_mask(values=np.arange(49).reshape(7, 7), pixel_scales=1.0) - - masked_image = aa.Array2D(values=image.native, mask=mask) - - blurring_mask = mask.derive_mask.blurring_from( - kernel_shape_native=kernel.shape_native - ) - - blurring_image = aa.Array2D(values=image.native, mask=blurring_mask) - - blurred_masked_im_1 = kernel.convolve_image( - image=masked_image, blurring_image=blurring_image - ) - - print(blurred_masked_image_via_scipy) - print(blurred_masked_im_1) - - assert blurred_masked_image_via_scipy == pytest.approx(blurred_masked_im_1.array, 1e-4) - - -def test__compare_to_full_2d_convolution__no_blurring_image(): - # Setup a blurred data, using the PSF to perform the convolution in 2D, then masks it to make a 1d array. - - import scipy.signal - - mask = aa.Mask2D.circular( - shape_native=(30, 30), pixel_scales=(1.0, 1.0), radius=4.0 - ) - kernel = aa.Kernel2D.no_mask(values=np.arange(49).reshape(7, 7), pixel_scales=1.0) - image = aa.Array2D.no_mask(values=np.arange(900).reshape(30, 30), pixel_scales=1.0) - - blurring_mask = mask.derive_mask.blurring_from( - kernel_shape_native=kernel.shape_native - ) - blurred_image_via_scipy = scipy.signal.convolve2d( - image.native * blurring_mask, kernel.native, mode="same" - ) - blurred_image_via_scipy = aa.Array2D.no_mask( - values=blurred_image_via_scipy, pixel_scales=1.0 - ) - blurred_masked_image_via_scipy = aa.Array2D( - values=blurred_image_via_scipy.native, mask=mask - ) - - # Now reproduce this data using the frame convolver_image - - masked_image = aa.Array2D(values=image.native, mask=mask) - - blurred_masked_im_1 = kernel.convolve_image_no_blurring(image=masked_image) - - assert blurred_masked_image_via_scipy == pytest.approx(blurred_masked_im_1, 1e-4) From 44cd4157b070a56b26077189352eb7aa105d967b Mon Sep 17 00:00:00 2001 From: James Nightingale Date: Wed, 2 Apr 2025 20:01:05 +0100 Subject: [PATCH 08/11] convolve_image_no_blurring --- autoarray/structures/arrays/kernel_2d.py | 37 +++++++++++++++++++ .../structures/arrays/test_kernel_2d.py | 20 ++++++---- 2 files changed, 49 insertions(+), 8 deletions(-) diff --git a/autoarray/structures/arrays/kernel_2d.py b/autoarray/structures/arrays/kernel_2d.py index 3f34c49bc..cf2f451f9 100644 --- a/autoarray/structures/arrays/kernel_2d.py +++ b/autoarray/structures/arrays/kernel_2d.py @@ -527,4 +527,41 @@ def convolve_image(self, image, blurring_image, jax_method="fft"): convolved_array_1d = convolve_native[slim_to_native] + return Array2D(values=convolved_array_1d, mask=image.mask) + + def convolve_image_no_blurring(self, image, jax_method="fft"): + """ + For a given 1D array and blurring array, convolve the two using this convolver. + + Parameters + ---------- + image + 1D array of the values which are to be blurred with the convolver's PSF. + blurring_image + 1D array of the blurring values which blur into the array after PSF convolution. + jax_method + If JAX is enabled this keyword will indicate what method is used for the PSF + convolution. Can be either `direct` to calculate it in real space or `fft` + to calculated it via a fast Fourier transform. `fft` is typically faster for + kernels that are more than about 5x5. Default is `fft`. + """ + + slim_to_native = jnp.nonzero( + jnp.logical_not(image.mask.array), size=image.shape[0] + ) + + expanded_array_native = jnp.zeros(image.mask.shape) + + expanded_array_native = expanded_array_native.at[slim_to_native].set( + image.array + ) + + kernel = np.array(self.native.array) + + convolve_native = jax.scipy.signal.convolve( + expanded_array_native, kernel, mode="same", method=jax_method + ) + + convolved_array_1d = convolve_native[slim_to_native] + return Array2D(values=convolved_array_1d, mask=image.mask) \ No newline at end of file diff --git a/test_autoarray/structures/arrays/test_kernel_2d.py b/test_autoarray/structures/arrays/test_kernel_2d.py index d70854486..c24092ff1 100644 --- a/test_autoarray/structures/arrays/test_kernel_2d.py +++ b/test_autoarray/structures/arrays/test_kernel_2d.py @@ -398,22 +398,23 @@ def test__convolve_image(): assert blurred_masked_image_via_scipy == pytest.approx(blurred_masked_im_1.array, 1e-4) -def test__compare_to_full_2d_convolution__no_blurring_image(): +def test__convolve_image_no_blurring(): # Setup a blurred data, using the PSF to perform the convolution in 2D, then masks it to make a 1d array. - import scipy.signal - mask = aa.Mask2D.circular( shape_native=(30, 30), pixel_scales=(1.0, 1.0), radius=4.0 ) - kernel = aa.Kernel2D.no_mask(values=np.arange(49).reshape(7, 7), pixel_scales=1.0) - image = aa.Array2D.no_mask(values=np.arange(900).reshape(30, 30), pixel_scales=1.0) + + import scipy.signal + + kernel = np.arange(49).reshape(7, 7) + image = np.arange(900).reshape(30, 30) blurring_mask = mask.derive_mask.blurring_from( - kernel_shape_native=kernel.shape_native + kernel_shape_native=kernel.shape ) blurred_image_via_scipy = scipy.signal.convolve2d( - image.native * blurring_mask, kernel.native, mode="same" + image * blurring_mask, kernel, mode="same" ) blurred_image_via_scipy = aa.Array2D.no_mask( values=blurred_image_via_scipy, pixel_scales=1.0 @@ -424,11 +425,14 @@ def test__compare_to_full_2d_convolution__no_blurring_image(): # Now reproduce this data using the frame convolver_image + kernel = aa.Kernel2D.no_mask(values=np.arange(49).reshape(7, 7), pixel_scales=1.0) + image = aa.Array2D.no_mask(values=np.arange(900).reshape(30, 30), pixel_scales=1.0) + masked_image = aa.Array2D(values=image.native, mask=mask) blurred_masked_im_1 = kernel.convolve_image_no_blurring(image=masked_image) - assert blurred_masked_image_via_scipy == pytest.approx(blurred_masked_im_1, 1e-4) + assert blurred_masked_image_via_scipy == pytest.approx(blurred_masked_im_1.array, 1e-4) def test__convolve_mapping_matrix(): From 4e0b92596ad40168b9908bec6d356a7a4bcce421 Mon Sep 17 00:00:00 2001 From: James Nightingale Date: Wed, 2 Apr 2025 20:55:00 +0100 Subject: [PATCH 09/11] maapping matrix convolve works --- autoarray/dataset/imaging/dataset.py | 3 +++ .../inversion/inversion/imaging/w_tilde.py | 3 ++- autoarray/structures/arrays/kernel_2d.py | 23 +++++++++++------- .../dataset/imaging/test_dataset.py | 7 ++++++ .../structures/arrays/test_kernel_2d.py | 24 ++++++++----------- 5 files changed, 37 insertions(+), 23 deletions(-) diff --git a/autoarray/dataset/imaging/dataset.py b/autoarray/dataset/imaging/dataset.py index 7fa3c515f..e3ec74d3b 100644 --- a/autoarray/dataset/imaging/dataset.py +++ b/autoarray/dataset/imaging/dataset.py @@ -166,6 +166,9 @@ def __init__( self.psf = psf + if psf.mask.shape[0] % 2 == 0 or psf.mask.shape[1] % 2 == 0: + raise exc.KernelException("Kernel2D Kernel2D must be odd") + @cached_property def grids(self): return GridsDataset( diff --git a/autoarray/inversion/inversion/imaging/w_tilde.py b/autoarray/inversion/inversion/imaging/w_tilde.py index c249f4ce9..888be987e 100644 --- a/autoarray/inversion/inversion/imaging/w_tilde.py +++ b/autoarray/inversion/inversion/imaging/w_tilde.py @@ -526,7 +526,8 @@ def mapped_reconstructed_data_dict(self) -> Dict[LinearObj, Array2D]: ) mapped_reconstructed_image = self.convolver.convolve_image_no_blurring( - image=mapped_reconstructed_image + image=mapped_reconstructed_image, + mask=self.mask ) else: diff --git a/autoarray/structures/arrays/kernel_2d.py b/autoarray/structures/arrays/kernel_2d.py index cf2f451f9..a9c148709 100644 --- a/autoarray/structures/arrays/kernel_2d.py +++ b/autoarray/structures/arrays/kernel_2d.py @@ -54,9 +54,6 @@ def __init__( store_native=store_native, ) - if self.mask.shape[0] % 2 == 0 or self.mask.shape[1] % 2 == 0: - raise exc.KernelException("Kernel2D Kernel2D must be odd") - if normalize: self._array = np.divide(self._array, np.sum(self._array)) @@ -529,7 +526,7 @@ def convolve_image(self, image, blurring_image, jax_method="fft"): return Array2D(values=convolved_array_1d, mask=image.mask) - def convolve_image_no_blurring(self, image, jax_method="fft"): + def convolve_image_no_blurring(self, image, mask, jax_method="fft"): """ For a given 1D array and blurring array, convolve the two using this convolver. @@ -547,13 +544,13 @@ def convolve_image_no_blurring(self, image, jax_method="fft"): """ slim_to_native = jnp.nonzero( - jnp.logical_not(image.mask.array), size=image.shape[0] + jnp.logical_not(mask.array), size=image.shape[0] ) - expanded_array_native = jnp.zeros(image.mask.shape) + expanded_array_native = jnp.zeros(mask.shape) expanded_array_native = expanded_array_native.at[slim_to_native].set( - image.array + image ) kernel = np.array(self.native.array) @@ -564,4 +561,14 @@ def convolve_image_no_blurring(self, image, jax_method="fft"): convolved_array_1d = convolve_native[slim_to_native] - return Array2D(values=convolved_array_1d, mask=image.mask) \ No newline at end of file + return Array2D(values=convolved_array_1d, mask=mask) + + def convolve_mapping_matrix(self, mapping_matrix, mask): + """For a given 1D array and blurring array, convolve the two using this convolver. + + Parameters + ---------- + image + 1D array of the values which are to be blurred with the convolver's PSF. + """ + return jax.vmap(self.convolve_image_no_blurring, in_axes=(1, None))(mapping_matrix, mask).T \ No newline at end of file diff --git a/test_autoarray/dataset/imaging/test_dataset.py b/test_autoarray/dataset/imaging/test_dataset.py index 9c079b94d..9758f7cdb 100644 --- a/test_autoarray/dataset/imaging/test_dataset.py +++ b/test_autoarray/dataset/imaging/test_dataset.py @@ -8,6 +8,8 @@ import autoarray as aa +from autoarray import exc + test_data_path = path.join( "{}".format(path.dirname(path.realpath(__file__))), "files", @@ -241,3 +243,8 @@ def test__noise_map_unmasked_has_zeros_or_negative__raises_exception(): with pytest.raises(aa.exc.DatasetException): aa.Imaging(data=array, noise_map=noise_map) + +def test__psf_not_odd_x_odd_kernel__raises_error(): + + with pytest.raises(exc.KernelException): + aa.Kernel2D.no_mask(values=[[0.0, 1.0], [1.0, 2.0]], pixel_scales=1.0) \ No newline at end of file diff --git a/test_autoarray/structures/arrays/test_kernel_2d.py b/test_autoarray/structures/arrays/test_kernel_2d.py index c24092ff1..2941801ba 100644 --- a/test_autoarray/structures/arrays/test_kernel_2d.py +++ b/test_autoarray/structures/arrays/test_kernel_2d.py @@ -292,12 +292,6 @@ def test__from_as_gaussian_via_alma_fits_header_parameters__identical_to_astropy assert kernel_astropy == pytest.approx(kernel_2d.native._array, abs=1e-4) -def test__not_odd_x_odd_kernel__raises_error(): - - with pytest.raises(exc.KernelException): - aa.Kernel2D.no_mask(values=[[0.0, 1.0], [1.0, 2.0]], pixel_scales=1.0) - - def test__convolved_array_from(): array_2d = aa.Array2D.no_mask( @@ -430,13 +424,13 @@ def test__convolve_image_no_blurring(): masked_image = aa.Array2D(values=image.native, mask=mask) - blurred_masked_im_1 = kernel.convolve_image_no_blurring(image=masked_image) + blurred_masked_im_1 = kernel.convolve_image_no_blurring(image=masked_image, mask=mask) assert blurred_masked_image_via_scipy == pytest.approx(blurred_masked_im_1.array, 1e-4) def test__convolve_mapping_matrix(): - mask = np.array( + mask = aa.Mask2D(mask=np.array( [ [True, True, True, True, True, True], [True, False, False, False, False, True], @@ -445,7 +439,7 @@ def test__convolve_mapping_matrix(): [True, False, False, False, False, True], [True, True, True, True, True, True], ] - ) + ), pixel_scales=1.0) kernel = aa.Kernel2D.no_mask( values=[[0, 0.0, 0], [0.4, 0.2, 0.3], [0, 0.1, 0]], pixel_scales=1.0 @@ -476,11 +470,11 @@ def test__convolve_mapping_matrix(): ] ) - blurred_mapping = kernel.convolve_mapping_matrix(mapping) + blurred_mapping = kernel.convolve_mapping_matrix(mapping, mask) assert ( blurred_mapping - == np.array( + == pytest.approx(np.array( [ [0, 0, 0], [0, 0, 0], @@ -500,7 +494,7 @@ def test__convolve_mapping_matrix(): [0, 0, 0], ] ) - ).all() + ), 1.0e-4) kernel = aa.Kernel2D.no_mask( values=[[0, 0.0, 0], [0.4, 0.2, 0.3], [0, 0.1, 0]], pixel_scales=1.0 @@ -531,7 +525,9 @@ def test__convolve_mapping_matrix(): ] ) - blurred_mapping = kernel.convolve_mapping_matrix(mapping) + blurred_mapping = kernel.convolve_mapping_matrix(mapping, mask) + + print(blurred_mapping) assert blurred_mapping == pytest.approx( np.array( @@ -554,5 +550,5 @@ def test__convolve_mapping_matrix(): [0, 0, 0], ] ), - 1e-4, + abs=1e-4, ) From 731fcb5914e2b11abaf9b6bb520e277c140ac9e0 Mon Sep 17 00:00:00 2001 From: James Nightingale Date: Wed, 2 Apr 2025 20:58:07 +0100 Subject: [PATCH 10/11] black --- autoarray/fixtures.py | 1 - autoarray/geometry/geometry_util.py | 21 ++--- .../inversion/inversion/imaging/w_tilde.py | 3 +- autoarray/mask/derive/indexes_2d.py | 6 +- autoarray/mask/mask_1d.py | 1 + autoarray/mask/mask_2d_util.py | 33 ++++--- .../over_sampling/over_sample_util.py | 10 +-- .../operators/over_sampling/over_sampler.py | 4 +- autoarray/structures/arrays/array_2d_util.py | 13 ++- autoarray/structures/arrays/kernel_2d.py | 16 ++-- autoarray/structures/grids/uniform_2d.py | 6 +- .../dataset/imaging/test_dataset.py | 11 +-- test_autoarray/geometry/test_geometry_util.py | 4 +- .../imaging/test_inversion_imaging_util.py | 4 +- test_autoarray/mask/test_mask_2d_util.py | 1 - .../structures/arrays/test_kernel_2d.py | 90 ++++++++++--------- .../structures/grids/test_uniform_2d.py | 2 - 17 files changed, 110 insertions(+), 116 deletions(-) diff --git a/autoarray/fixtures.py b/autoarray/fixtures.py index 68a8a73f8..e546fb497 100644 --- a/autoarray/fixtures.py +++ b/autoarray/fixtures.py @@ -110,7 +110,6 @@ def make_blurring_grid_2d_7x7(): return aa.Grid2D.from_mask(mask=make_blurring_mask_2d_7x7()) - def make_image_7x7(): return aa.Array2D.ones(shape_native=(7, 7), pixel_scales=(1.0, 1.0)) diff --git a/autoarray/geometry/geometry_util.py b/autoarray/geometry/geometry_util.py index b646c7d08..e2c6c8898 100644 --- a/autoarray/geometry/geometry_util.py +++ b/autoarray/geometry/geometry_util.py @@ -181,6 +181,7 @@ def convert_pixel_scales_2d(pixel_scales: ty.PixelScales) -> Tuple[float, float] return pixel_scales + @numba_util.jit() def central_pixel_coordinates_2d_numba_from( shape_native: Tuple[int, int], @@ -205,6 +206,7 @@ def central_pixel_coordinates_2d_numba_from( """ return (float(shape_native[0] - 1) / 2, float(shape_native[1] - 1) / 2) + @numba_util.jit() def central_scaled_coordinate_2d_numba_from( shape_native: Tuple[int, int], @@ -305,6 +307,7 @@ def central_scaled_coordinate_2d_from( return (y_pixel, x_pixel) + def pixel_coordinates_2d_from( scaled_coordinates_2d: Tuple[float, float], shape_native: Tuple[int, int], @@ -589,9 +592,9 @@ def grid_pixel_centres_2d_slim_from( centres_scaled = np.array(centres_scaled) pixel_scales = np.array(pixel_scales) sign = np.array([-1.0, 1.0]) - return ( - (sign * grid_scaled_2d_slim / pixel_scales) + centres_scaled + 0.5 - ).astype(int) + return ((sign * grid_scaled_2d_slim / pixel_scales) + centres_scaled + 0.5).astype( + int + ) def grid_pixel_indexes_2d_slim_from( @@ -647,9 +650,7 @@ def grid_pixel_indexes_2d_slim_from( ) return ( - (grid_pixels_2d_slim * np.array([shape_native[1], 1])) - .sum(axis=1) - .astype(int) + (grid_pixels_2d_slim * np.array([shape_native[1], 1])).sum(axis=1).astype(int) ) @@ -698,9 +699,7 @@ def grid_scaled_2d_slim_from( centres_scaled = np.array(centres_scaled) pixel_scales = np.array(pixel_scales) sign = np.array([-1, 1]) - return ( - (grid_pixels_2d_slim - centres_scaled - 0.5) * pixel_scales * sign - ) + return (grid_pixels_2d_slim - centres_scaled - 0.5) * pixel_scales * sign def grid_pixel_centres_2d_from( @@ -750,9 +749,7 @@ def grid_pixel_centres_2d_from( centres_scaled = np.array(centres_scaled) pixel_scales = np.array(pixel_scales) sign = np.array([-1.0, 1.0]) - return ( - (sign * grid_scaled_2d / pixel_scales) + centres_scaled + 0.5 - ).astype(int) + return ((sign * grid_scaled_2d / pixel_scales) + centres_scaled + 0.5).astype(int) def extent_symmetric_from( diff --git a/autoarray/inversion/inversion/imaging/w_tilde.py b/autoarray/inversion/inversion/imaging/w_tilde.py index 888be987e..aef314586 100644 --- a/autoarray/inversion/inversion/imaging/w_tilde.py +++ b/autoarray/inversion/inversion/imaging/w_tilde.py @@ -526,8 +526,7 @@ def mapped_reconstructed_data_dict(self) -> Dict[LinearObj, Array2D]: ) mapped_reconstructed_image = self.convolver.convolve_image_no_blurring( - image=mapped_reconstructed_image, - mask=self.mask + image=mapped_reconstructed_image, mask=self.mask ) else: diff --git a/autoarray/mask/derive/indexes_2d.py b/autoarray/mask/derive/indexes_2d.py index 4807799d9..062c8e664 100644 --- a/autoarray/mask/derive/indexes_2d.py +++ b/autoarray/mask/derive/indexes_2d.py @@ -198,9 +198,9 @@ def edge_slim(self) -> np.ndarray: print(derive_indexes_2d.edge_slim) """ - return mask_2d_util.edge_1d_indexes_from(mask_2d=np.array(self.mask).astype("bool")).astype( - "int" - ) + return mask_2d_util.edge_1d_indexes_from( + mask_2d=np.array(self.mask).astype("bool") + ).astype("int") @property def edge_native(self) -> np.ndarray: diff --git a/autoarray/mask/mask_1d.py b/autoarray/mask/mask_1d.py index 9330e62ee..8c36d8866 100644 --- a/autoarray/mask/mask_1d.py +++ b/autoarray/mask/mask_1d.py @@ -25,6 +25,7 @@ class Mask1DKeys(Enum): PIXSCA = "PIXSCA" ORIGIN = "ORIGIN" + class Mask1D(Mask): def __init__( self, diff --git a/autoarray/mask/mask_2d_util.py b/autoarray/mask/mask_2d_util.py index 9ea4e35c0..462448073 100644 --- a/autoarray/mask/mask_2d_util.py +++ b/autoarray/mask/mask_2d_util.py @@ -6,6 +6,7 @@ from autoarray import exc from autoarray.numpy_wrapper import np as jnp + def native_index_for_slim_index_2d_from( mask_2d: np.ndarray, ) -> np.ndarray: @@ -400,6 +401,7 @@ def mask_2d_via_pixel_coordinates_from( return mask_2d return buffed_mask_2d_from(mask_2d=mask_2d, buffer=buffer) # Apply buf + def min_false_distance_to_edge(mask: np.ndarray) -> Tuple[int, int]: """ Compute the minimum 1D distance in the y and x directions from any `False` value at the mask's extreme positions @@ -618,14 +620,18 @@ def edge_1d_indexes_from(mask_2d: np.ndarray) -> np.ndarray: array([0, 1, 2, 3, 5, 6, 7, 8]) """ # Pad the mask to handle edge cases without index errors - padded_mask = np.pad(mask_2d, pad_width=1, mode='constant', constant_values=True) + padded_mask = np.pad(mask_2d, pad_width=1, mode="constant", constant_values=True) # Identify neighbors in 3x3 regions around each pixel neighbors = ( - padded_mask[:-2, 1:-1] | padded_mask[2:, 1:-1] | # Up, Down - padded_mask[1:-1, :-2] | padded_mask[1:-1, 2:] | # Left, Right - padded_mask[:-2, :-2] | padded_mask[:-2, 2:] | # Top-left, Top-right - padded_mask[2:, :-2] | padded_mask[2:, 2:] # Bottom-left, Bottom-right + padded_mask[:-2, 1:-1] + | padded_mask[2:, 1:-1] # Up, Down + | padded_mask[1:-1, :-2] + | padded_mask[1:-1, 2:] # Left, Right + | padded_mask[:-2, :-2] + | padded_mask[:-2, 2:] # Top-left, Top-right + | padded_mask[2:, :-2] + | padded_mask[2:, 2:] # Bottom-left, Bottom-right ) # Identify edge pixels: False values with at least one True neighbor @@ -708,10 +714,10 @@ def border_slim_indexes_from(mask_2d: np.ndarray) -> np.ndarray: # Identify border pixels: where the full length in any direction is True border_mask = ( - (up_sums == np.arange(height)[:, None]) | - (down_sums == np.arange(height - 1, -1, -1)[:, None]) | - (left_sums == np.arange(width)[None, :]) | - (right_sums == np.arange(width - 1, -1, -1)[None, :]) + (up_sums == np.arange(height)[:, None]) + | (down_sums == np.arange(height - 1, -1, -1)[:, None]) + | (left_sums == np.arange(width)[None, :]) + | (right_sums == np.arange(width - 1, -1, -1)[None, :]) ) & ~mask_2d # Create an index array where False entries get sequential 1D indices @@ -767,14 +773,16 @@ def buffed_mask_2d_from(mask_2d: np.ndarray, buffer: int = 1) -> np.ndarray: buffer_range = np.arange(-buffer, buffer + 1) # Generate all possible neighbors for each False entry - dy, dx = np.meshgrid(buffer_range, buffer_range, indexing='ij') + dy, dx = np.meshgrid(buffer_range, buffer_range, indexing="ij") neighbors = np.stack([dy.ravel(), dx.ravel()], axis=-1) # Calculate all neighboring positions for all False coordinates all_neighbors = np.add(np.array(false_coords).T[:, np.newaxis], neighbors) # Clip the neighbors to stay within the bounds of the mask - valid_neighbors = np.clip(all_neighbors, [0, 0], [mask_2d.shape[0] - 1, mask_2d.shape[1] - 1]) + valid_neighbors = np.clip( + all_neighbors, [0, 0], [mask_2d.shape[0] - 1, mask_2d.shape[1] - 1] + ) # Update the buffed mask: set all the neighbors to False buffed_mask_2d[valid_neighbors[:, :, 0], valid_neighbors[:, :, 1]] = False @@ -833,6 +841,3 @@ def rescaled_mask_2d_from(mask_2d: np.ndarray, rescale_factor: float) -> np.ndar rescaled_mask_2d[:, 0] = 1 rescaled_mask_2d[:, rescaled_mask_2d.shape[1] - 1] = 1 return np.isclose(rescaled_mask_2d, 1) - - - diff --git a/autoarray/operators/over_sampling/over_sample_util.py b/autoarray/operators/over_sampling/over_sample_util.py index 8966e785e..084d0bd0d 100644 --- a/autoarray/operators/over_sampling/over_sample_util.py +++ b/autoarray/operators/over_sampling/over_sample_util.py @@ -463,16 +463,10 @@ def grid_2d_slim_over_sampled_via_mask_from( # ) # else: grid_slim[sub_index, 0] = -( - y_scaled - - y_sub_half - + y1 * y_sub_step - + (y_sub_step / 2.0) + y_scaled - y_sub_half + y1 * y_sub_step + (y_sub_step / 2.0) ) grid_slim[sub_index, 1] = ( - x_scaled - - x_sub_half - + x1 * x_sub_step - + (x_sub_step / 2.0) + x_scaled - x_sub_half + x1 * x_sub_step + (x_sub_step / 2.0) ) sub_index += 1 diff --git a/autoarray/operators/over_sampling/over_sampler.py b/autoarray/operators/over_sampling/over_sampler.py index 0aafbe008..6f12f4b9f 100644 --- a/autoarray/operators/over_sampling/over_sampler.py +++ b/autoarray/operators/over_sampling/over_sampler.py @@ -226,7 +226,9 @@ def binned_array_2d_from(self, array: Array2D) -> "Array2D": # sub_size=np.array(self.sub_size).astype("int"), # ) - binned_array_2d = array.reshape(self.mask.shape_slim, self.sub_size[0]**2).mean(axis=1) + binned_array_2d = array.reshape( + self.mask.shape_slim, self.sub_size[0] ** 2 + ).mean(axis=1) return Array2D( values=binned_array_2d, diff --git a/autoarray/structures/arrays/array_2d_util.py b/autoarray/structures/arrays/array_2d_util.py index 0e694846d..d8a0f7e10 100644 --- a/autoarray/structures/arrays/array_2d_util.py +++ b/autoarray/structures/arrays/array_2d_util.py @@ -27,10 +27,7 @@ def convert_array(array: Union[np.ndarray, List]) -> np.ndarray: array = np.asarray(array) elif isinstance(array, jnp.ndarray): array = jax.lax.cond( - type(array) is list, - lambda _: jnp.asarray(array), - lambda _: array, - None + type(array) is list, lambda _: jnp.asarray(array), lambda _: array, None ) return array @@ -41,6 +38,7 @@ def check_array_2d(array_2d: np.ndarray): "An array input into the Array2D.__new__ method is not of shape 1." ) + def check_array_2d_and_mask_2d(array_2d: np.ndarray, mask_2d: Mask2D): """ The `manual` classmethods in the `Array2D` object take as input a list or ndarray which is returned as an @@ -90,6 +88,7 @@ def check_array_2d_and_mask_2d(array_2d: np.ndarray, mask_2d: Mask2D): """ ) + def convert_array_2d( array_2d: Union[np.ndarray, List], mask_2d: Mask2D, @@ -489,6 +488,7 @@ def index_slim_for_index_2d_from(indexes_2d: np.ndarray, shape_native) -> np.nda return index_slim_for_index_native_2d + def array_2d_slim_from( array_2d_native: np.ndarray, mask_2d: np.ndarray, @@ -534,6 +534,7 @@ def array_2d_slim_from( """ return array_2d_native[~mask_2d.astype(bool)] + def array_2d_native_from( array_2d_slim: np.ndarray, mask_2d: np.ndarray, @@ -620,9 +621,7 @@ def array_2d_via_indexes_from( The native 2D array of values mapped from the slimmed array with dimensions (total_values, total_values). """ return ( - jnp.zeros(shape) - .at[tuple(native_index_for_slim_index_2d.T)] - .set(array_2d_slim) + jnp.zeros(shape).at[tuple(native_index_for_slim_index_2d.T)].set(array_2d_slim) ) diff --git a/autoarray/structures/arrays/kernel_2d.py b/autoarray/structures/arrays/kernel_2d.py index a9c148709..5938198cf 100644 --- a/autoarray/structures/arrays/kernel_2d.py +++ b/autoarray/structures/arrays/kernel_2d.py @@ -474,7 +474,9 @@ def convolved_array_from(self, array: Array2D) -> Array2D: array_2d = array.native - convolved_array_2d = scipy.signal.convolve2d(array_2d._array, np.array(self.native._array), mode="same") + convolved_array_2d = scipy.signal.convolve2d( + array_2d._array, np.array(self.native._array), mode="same" + ) convolved_array_1d = array_2d_util.array_2d_slim_from( mask_2d=np.array(array_2d.mask), @@ -543,15 +545,11 @@ def convolve_image_no_blurring(self, image, mask, jax_method="fft"): kernels that are more than about 5x5. Default is `fft`. """ - slim_to_native = jnp.nonzero( - jnp.logical_not(mask.array), size=image.shape[0] - ) + slim_to_native = jnp.nonzero(jnp.logical_not(mask.array), size=image.shape[0]) expanded_array_native = jnp.zeros(mask.shape) - expanded_array_native = expanded_array_native.at[slim_to_native].set( - image - ) + expanded_array_native = expanded_array_native.at[slim_to_native].set(image) kernel = np.array(self.native.array) @@ -571,4 +569,6 @@ def convolve_mapping_matrix(self, mapping_matrix, mask): image 1D array of the values which are to be blurred with the convolver's PSF. """ - return jax.vmap(self.convolve_image_no_blurring, in_axes=(1, None))(mapping_matrix, mask).T \ No newline at end of file + return jax.vmap(self.convolve_image_no_blurring, in_axes=(1, None))( + mapping_matrix, mask + ).T diff --git a/autoarray/structures/grids/uniform_2d.py b/autoarray/structures/grids/uniform_2d.py index 4d565dd27..670cbfcbe 100644 --- a/autoarray/structures/grids/uniform_2d.py +++ b/autoarray/structures/grids/uniform_2d.py @@ -840,9 +840,9 @@ def squared_distances_to_coordinate_from( The (y,x) coordinate from which the squared distance of every grid (y,x) coordinate is computed. """ if isinstance(self, jnp.ndarray): - squared_distances = jnp.square(self.array[:, 0] - coordinate[0]) + jnp.square( - self.array[:, 1] - coordinate[1] - ) + squared_distances = jnp.square( + self.array[:, 0] - coordinate[0] + ) + jnp.square(self.array[:, 1] - coordinate[1]) else: squared_distances = np.square(self[:, 0] - coordinate[0]) + np.square( self[:, 1] - coordinate[1] diff --git a/test_autoarray/dataset/imaging/test_dataset.py b/test_autoarray/dataset/imaging/test_dataset.py index 9758f7cdb..45ef0a80f 100644 --- a/test_autoarray/dataset/imaging/test_dataset.py +++ b/test_autoarray/dataset/imaging/test_dataset.py @@ -38,16 +38,12 @@ def test__psf_and_mask_hit_edge__automatically_pads_image_and_noise_map(): noise_map = aa.Array2D.ones(shape_native=(3, 3), pixel_scales=1.0) psf = aa.Kernel2D.ones(shape_native=(3, 3), pixel_scales=1.0) - dataset = aa.Imaging( - data=image, noise_map=noise_map, psf=psf, pad_for_psf=False - ) + dataset = aa.Imaging(data=image, noise_map=noise_map, psf=psf, pad_for_psf=False) assert dataset.data.shape_native == (3, 3) assert dataset.noise_map.shape_native == (3, 3) - dataset = aa.Imaging( - data=image, noise_map=noise_map, psf=psf, pad_for_psf=True - ) + dataset = aa.Imaging(data=image, noise_map=noise_map, psf=psf, pad_for_psf=True) assert dataset.data.shape_native == (5, 5) assert dataset.noise_map.shape_native == (5, 5) @@ -244,7 +240,8 @@ def test__noise_map_unmasked_has_zeros_or_negative__raises_exception(): with pytest.raises(aa.exc.DatasetException): aa.Imaging(data=array, noise_map=noise_map) + def test__psf_not_odd_x_odd_kernel__raises_error(): with pytest.raises(exc.KernelException): - aa.Kernel2D.no_mask(values=[[0.0, 1.0], [1.0, 2.0]], pixel_scales=1.0) \ No newline at end of file + aa.Kernel2D.no_mask(values=[[0.0, 1.0], [1.0, 2.0]], pixel_scales=1.0) diff --git a/test_autoarray/geometry/test_geometry_util.py b/test_autoarray/geometry/test_geometry_util.py index 4bc2706dc..e65608bf1 100644 --- a/test_autoarray/geometry/test_geometry_util.py +++ b/test_autoarray/geometry/test_geometry_util.py @@ -994,7 +994,8 @@ def test__transform_2d_grid_to_reference_frame(): [0.0, np.sqrt(2)], [np.sqrt(2) / 2.0, np.sqrt(2) / 2.0], ] - ), abs=1.0e-4 + ), + abs=1.0e-4, ) transformed_grid_2d = aa.util.geometry.transform_grid_2d_to_reference_frame( @@ -1071,7 +1072,6 @@ def test__transform_2d_grid_from_reference_frame(): grid_2d=transformed_grid_2d, centre=(8.0, 5.0), angle=137.0 ) - assert grid_2d == pytest.approx(original_grid_2d, abs=1.0e-4) diff --git a/test_autoarray/inversion/inversion/imaging/test_inversion_imaging_util.py b/test_autoarray/inversion/inversion/imaging/test_inversion_imaging_util.py index b616a95d4..e05e13350 100644 --- a/test_autoarray/inversion/inversion/imaging/test_inversion_imaging_util.py +++ b/test_autoarray/inversion/inversion/imaging/test_inversion_imaging_util.py @@ -282,9 +282,7 @@ def test__curvature_matrix_via_w_tilde_two_methods_agree(): w_tilde=w_tilde, mapping_matrix=mapping_matrix ) - blurred_mapping_matrix = psf.convolve_mapping_matrix( - mapping_matrix=mapping_matrix - ) + blurred_mapping_matrix = psf.convolve_mapping_matrix(mapping_matrix=mapping_matrix) curvature_matrix = aa.util.inversion.curvature_matrix_via_mapping_matrix_from( mapping_matrix=blurred_mapping_matrix, diff --git a/test_autoarray/mask/test_mask_2d_util.py b/test_autoarray/mask/test_mask_2d_util.py index bd2ebd84a..f3db2938e 100644 --- a/test_autoarray/mask/test_mask_2d_util.py +++ b/test_autoarray/mask/test_mask_2d_util.py @@ -723,7 +723,6 @@ def test__mask_1d_indexes_from(): assert masked_slim[-1] == 48 - def test__edge_1d_indexes_from(): mask = np.array( [ diff --git a/test_autoarray/structures/arrays/test_kernel_2d.py b/test_autoarray/structures/arrays/test_kernel_2d.py index 2941801ba..4cf8b92d7 100644 --- a/test_autoarray/structures/arrays/test_kernel_2d.py +++ b/test_autoarray/structures/arrays/test_kernel_2d.py @@ -362,9 +362,7 @@ def test__convolve_image(): kernel = np.arange(49).reshape(7, 7) image = np.arange(900).reshape(30, 30) - blurred_image_via_scipy = scipy.signal.convolve2d( - image, kernel, mode="same" - ) + blurred_image_via_scipy = scipy.signal.convolve2d(image, kernel, mode="same") blurred_image_via_scipy = aa.Array2D.no_mask( values=blurred_image_via_scipy, pixel_scales=1.0 ) @@ -389,7 +387,9 @@ def test__convolve_image(): image=masked_image, blurring_image=blurring_image ) - assert blurred_masked_image_via_scipy == pytest.approx(blurred_masked_im_1.array, 1e-4) + assert blurred_masked_image_via_scipy == pytest.approx( + blurred_masked_im_1.array, 1e-4 + ) def test__convolve_image_no_blurring(): @@ -404,9 +404,7 @@ def test__convolve_image_no_blurring(): kernel = np.arange(49).reshape(7, 7) image = np.arange(900).reshape(30, 30) - blurring_mask = mask.derive_mask.blurring_from( - kernel_shape_native=kernel.shape - ) + blurring_mask = mask.derive_mask.blurring_from(kernel_shape_native=kernel.shape) blurred_image_via_scipy = scipy.signal.convolve2d( image * blurring_mask, kernel, mode="same" ) @@ -424,22 +422,29 @@ def test__convolve_image_no_blurring(): masked_image = aa.Array2D(values=image.native, mask=mask) - blurred_masked_im_1 = kernel.convolve_image_no_blurring(image=masked_image, mask=mask) + blurred_masked_im_1 = kernel.convolve_image_no_blurring( + image=masked_image, mask=mask + ) - assert blurred_masked_image_via_scipy == pytest.approx(blurred_masked_im_1.array, 1e-4) + assert blurred_masked_image_via_scipy == pytest.approx( + blurred_masked_im_1.array, 1e-4 + ) def test__convolve_mapping_matrix(): - mask = aa.Mask2D(mask=np.array( - [ - [True, True, True, True, True, True], - [True, False, False, False, False, True], - [True, False, False, False, False, True], - [True, False, False, False, False, True], - [True, False, False, False, False, True], - [True, True, True, True, True, True], - ] - ), pixel_scales=1.0) + mask = aa.Mask2D( + mask=np.array( + [ + [True, True, True, True, True, True], + [True, False, False, False, False, True], + [True, False, False, False, False, True], + [True, False, False, False, False, True], + [True, False, False, False, False, True], + [True, True, True, True, True, True], + ] + ), + pixel_scales=1.0, + ) kernel = aa.Kernel2D.no_mask( values=[[0, 0.0, 0], [0.4, 0.2, 0.3], [0, 0.1, 0]], pixel_scales=1.0 @@ -474,27 +479,30 @@ def test__convolve_mapping_matrix(): assert ( blurred_mapping - == pytest.approx(np.array( - [ - [0, 0, 0], - [0, 0, 0], - [0, 0, 0], - [0, 0, 0], - [0, 0, 0], - [0, 0, 0], - [0, 0.4, 0], - [0, 0.2, 0], - [0.4, 0, 0], - [0.2, 0, 0.4], - [0.3, 0, 0.2], - [0, 0.1, 0.3], - [0, 0, 0], - [0.1, 0, 0], - [0, 0, 0.1], - [0, 0, 0], - ] - ) - ), 1.0e-4) + == pytest.approx( + np.array( + [ + [0, 0, 0], + [0, 0, 0], + [0, 0, 0], + [0, 0, 0], + [0, 0, 0], + [0, 0, 0], + [0, 0.4, 0], + [0, 0.2, 0], + [0.4, 0, 0], + [0.2, 0, 0.4], + [0.3, 0, 0.2], + [0, 0.1, 0.3], + [0, 0, 0], + [0.1, 0, 0], + [0, 0, 0.1], + [0, 0, 0], + ] + ) + ), + 1.0e-4, + ) kernel = aa.Kernel2D.no_mask( values=[[0, 0.0, 0], [0.4, 0.2, 0.3], [0, 0.1, 0]], pixel_scales=1.0 @@ -527,8 +535,6 @@ def test__convolve_mapping_matrix(): blurred_mapping = kernel.convolve_mapping_matrix(mapping, mask) - print(blurred_mapping) - assert blurred_mapping == pytest.approx( np.array( [ diff --git a/test_autoarray/structures/grids/test_uniform_2d.py b/test_autoarray/structures/grids/test_uniform_2d.py index 686721c70..78813329e 100644 --- a/test_autoarray/structures/grids/test_uniform_2d.py +++ b/test_autoarray/structures/grids/test_uniform_2d.py @@ -575,7 +575,6 @@ def test__grid_2d_radial_projected_shape_slim_from(): pixel_scales=grid_2d.pixel_scales, ) - assert grid_radii == pytest.approx(grid_radii_util, 1.0e-4) assert grid_radial_shape_slim == grid_radii_util.shape[0] @@ -782,7 +781,6 @@ def test__grid_with_coordinates_within_distance_removed_from(): ).all() - def test__recursive_shape_storage(): grid_2d = aa.Grid2D.no_mask( values=[[[1.0, 2.0], [3.0, 4.0]], [[5.0, 6.0], [7.0, 8.0]]], From ef5ba9efef1ea3c6234fcc5557c380e9b8e20db6 Mon Sep 17 00:00:00 2001 From: James Nightingale Date: Wed, 2 Apr 2025 21:08:38 +0100 Subject: [PATCH 11/11] finish --- .../inversion/inversion/imaging/w_tilde.py | 2 +- autoarray/structures/arrays/kernel_2d.py | 12 ++++++------ .../arrays/files/array/output_test/array.fits | Bin 5760 -> 5760 bytes 3 files changed, 7 insertions(+), 7 deletions(-) diff --git a/autoarray/inversion/inversion/imaging/w_tilde.py b/autoarray/inversion/inversion/imaging/w_tilde.py index aef314586..199ba66b2 100644 --- a/autoarray/inversion/inversion/imaging/w_tilde.py +++ b/autoarray/inversion/inversion/imaging/w_tilde.py @@ -525,7 +525,7 @@ def mapped_reconstructed_data_dict(self) -> Dict[LinearObj, Array2D]: values=mapped_reconstructed_image, mask=self.mask ) - mapped_reconstructed_image = self.convolver.convolve_image_no_blurring( + mapped_reconstructed_image = self.psf.convolve_image_no_blurring( image=mapped_reconstructed_image, mask=self.mask ) diff --git a/autoarray/structures/arrays/kernel_2d.py b/autoarray/structures/arrays/kernel_2d.py index 5938198cf..bb62656a9 100644 --- a/autoarray/structures/arrays/kernel_2d.py +++ b/autoarray/structures/arrays/kernel_2d.py @@ -487,12 +487,12 @@ def convolved_array_from(self, array: Array2D) -> Array2D: def convolve_image(self, image, blurring_image, jax_method="fft"): """ - For a given 1D array and blurring array, convolve the two using this convolver. + For a given 1D array and blurring array, convolve the two using this psf. Parameters ---------- image - 1D array of the values which are to be blurred with the convolver's PSF. + 1D array of the values which are to be blurred with the psf's PSF. blurring_image 1D array of the blurring values which blur into the array after PSF convolution. jax_method @@ -530,12 +530,12 @@ def convolve_image(self, image, blurring_image, jax_method="fft"): def convolve_image_no_blurring(self, image, mask, jax_method="fft"): """ - For a given 1D array and blurring array, convolve the two using this convolver. + For a given 1D array and blurring array, convolve the two using this psf. Parameters ---------- image - 1D array of the values which are to be blurred with the convolver's PSF. + 1D array of the values which are to be blurred with the psf's PSF. blurring_image 1D array of the blurring values which blur into the array after PSF convolution. jax_method @@ -562,12 +562,12 @@ def convolve_image_no_blurring(self, image, mask, jax_method="fft"): return Array2D(values=convolved_array_1d, mask=mask) def convolve_mapping_matrix(self, mapping_matrix, mask): - """For a given 1D array and blurring array, convolve the two using this convolver. + """For a given 1D array and blurring array, convolve the two using this psf. Parameters ---------- image - 1D array of the values which are to be blurred with the convolver's PSF. + 1D array of the values which are to be blurred with the psf's PSF. """ return jax.vmap(self.convolve_image_no_blurring, in_axes=(1, None))( mapping_matrix, mask diff --git a/test_autoarray/structures/arrays/files/array/output_test/array.fits b/test_autoarray/structures/arrays/files/array/output_test/array.fits index dab9da2fb676efc023e53b73faa54fe78fc4000c..0cff0a8f7db90d3ded28c644cd66e3791627feeb 100644 GIT binary patch delta 84 hcmZqBZP1;N!(?W%G4C>$;|B&XuqT_|T*&>83jh#~5uX46 delta 56 kcmZqBZP1;N!(?o$Vgmz%Jzl(dBKLc)$qTsk0j$6e2mk;8