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..e3ec74d3b 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 @@ -167,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( @@ -176,25 +178,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 +353,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 +446,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 +494,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..e546fb497 100644 --- a/autoarray/fixtures.py +++ b/autoarray/fixtures.py @@ -110,13 +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/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/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..199ba66b2 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, ) @@ -528,8 +525,8 @@ 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( - image=mapped_reconstructed_image + mapped_reconstructed_image = self.psf.convolve_image_no_blurring( + image=mapped_reconstructed_image, mask=self.mask ) else: 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/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/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/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/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/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 f80e3f6e3..bb62656a9 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 @@ -6,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 @@ -15,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): @@ -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. @@ -472,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), @@ -481,33 +485,90 @@ 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: + def convolve_image(self, image, blurring_image, jax_method="fft"): """ - Convolve an array with this Kernel2D + For a given 1D array and blurring array, convolve the two using this psf. Parameters ---------- image - An array representing the image the Kernel2D is convolved with. + 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 + 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`. + """ - Returns - ------- - convolved_image - An array representing the image after convolution. + slim_to_native = jnp.nonzero( + jnp.logical_not(image.mask.array), size=image.shape[0] + ) + slim_to_native_blurring = jnp.nonzero( + jnp.logical_not(blurring_image.mask.array), size=blurring_image.shape[0] + ) - Raises - ------ - KernelException if either Kernel2D psf dimension is odd + expanded_array_native = jnp.zeros(image.mask.shape) + + expanded_array_native = expanded_array_native.at[slim_to_native].set( + image.array + ) + expanded_array_native = expanded_array_native.at[slim_to_native_blurring].set( + blurring_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) + + def convolve_image_no_blurring(self, image, mask, jax_method="fft"): """ + For a given 1D array and blurring array, convolve the two using this psf. - if self.mask.shape[0] % 2 == 0 or self.mask.shape[1] % 2 == 0: - raise exc.KernelException("Kernel2D Kernel2D must be odd") + Parameters + ---------- + image + 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 + 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`. + """ - convolved_array_2d = scipy.signal.convolve2d(array, self.native, mode="same") + slim_to_native = jnp.nonzero(jnp.logical_not(mask.array), size=image.shape[0]) - convolved_array_1d = array_2d_util.array_2d_slim_from( - mask_2d=np.array(mask), - array_2d_native=np.array(convolved_array_2d), + expanded_array_native = jnp.zeros(mask.shape) + + expanded_array_native = expanded_array_native.at[slim_to_native].set(image) + + 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=mask) + + def convolve_mapping_matrix(self, mapping_matrix, mask): + """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 psf's PSF. + """ + 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/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..45ef0a80f 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", @@ -36,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_convolver=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_convolver=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) @@ -144,7 +142,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 +223,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() @@ -243,3 +239,9 @@ 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) 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_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..e05e13350 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,9 +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( - 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, @@ -303,7 +301,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 +354,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/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/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/files/array/output_test/array.fits b/test_autoarray/structures/arrays/files/array/output_test/array.fits index dab9da2fb..0cff0a8f7 100644 Binary files a/test_autoarray/structures/arrays/files/array/output_test/array.fits and b/test_autoarray/structures/arrays/files/array/output_test/array.fits differ diff --git a/test_autoarray/structures/arrays/test_kernel_2d.py b/test_autoarray/structures/arrays/test_kernel_2d.py index 9106b43a1..4cf8b92d7 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 @@ -248,108 +247,52 @@ def test__rescaled_with_odd_dimensions_from__different_scalings(): assert (kernel_2d.native == (1.0 / 15.0) * np.ones((5, 3))).all() -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() +def test__from_as_gaussian_via_alma_fits_header_parameters__identical_to_astropy_gaussian_model(): + pixel_scales = 0.1 - 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, + x_stddev = ( + 1.0e-5 + * (units.deg).to(units.arcsec) + / pixel_scales + / (2.0 * np.sqrt(2.0 * np.log(2.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 + y_stddev = ( + 2.0e-5 + * (units.deg).to(units.arcsec) + / pixel_scales + / (2.0 * np.sqrt(2.0 * np.log(2.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, - ) + theta_deg = 230.0 + theta = Angle(theta_deg, "deg").radian - 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 + 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, ) - 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() + 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) - 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.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, ) - 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 - ) + assert kernel_astropy == pytest.approx(kernel_2d.native._array, abs=1e-4) - 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() +def test__convolved_array_from(): array_2d = aa.Array2D.no_mask( [ @@ -408,53 +351,210 @@ 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_image(): - with pytest.raises(exc.KernelException): - kernel_2d.convolved_array_from(np.ones((5, 5))) + mask = aa.Mask2D.circular( + shape_native=(30, 30), pixel_scales=(1.0, 1.0), radius=4.0 + ) + import scipy.signal -def test__from_as_gaussian_via_alma_fits_header_parameters__identical_to_astropy_gaussian_model(): - pixel_scales = 0.1 + kernel = np.arange(49).reshape(7, 7) + image = np.arange(900).reshape(30, 30) - x_stddev = ( - 1.0e-5 - * (units.deg).to(units.arcsec) - / pixel_scales - / (2.0 * np.sqrt(2.0 * np.log(2.0))) + 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 ) - y_stddev = ( - 2.0e-5 - * (units.deg).to(units.arcsec) - / pixel_scales - / (2.0 * np.sqrt(2.0 * np.log(2.0))) + blurred_masked_image_via_scipy = aa.Array2D( + values=blurred_image_via_scipy.native, mask=mask ) - theta_deg = 230.0 - theta = Angle(theta_deg, "deg").radian + # Now reproduce this data using the convolve_image function - 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, + 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 ) - 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) + blurring_image = aa.Array2D(values=image.native, mask=blurring_mask) - 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, + blurred_masked_im_1 = kernel.convolve_image( + image=masked_image, blurring_image=blurring_image ) - assert kernel_astropy == pytest.approx(kernel_2d.native._array, abs=1e-4) + assert blurred_masked_image_via_scipy == pytest.approx( + blurred_masked_im_1.array, 1e-4 + ) + + +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. + + 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) + + 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" + ) + 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 + + 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, mask=mask + ) + + 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, + ) + + kernel = aa.Kernel2D.no_mask( + values=[[0, 0.0, 0], [0.4, 0.2, 0.3], [0, 0.1, 0]], pixel_scales=1.0 + ) + + 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 = kernel.convolve_mapping_matrix(mapping, mask) + + 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, + ) + + kernel = aa.Kernel2D.no_mask( + values=[[0, 0.0, 0], [0.4, 0.2, 0.3], [0, 0.1, 0]], pixel_scales=1.0 + ) + + 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 = kernel.convolve_mapping_matrix(mapping, mask) + + 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], + ] + ), + abs=1e-4, + ) 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]]],