From fd11b178e6a64855650ac39fcb48d5697a5ab66d Mon Sep 17 00:00:00 2001 From: Jammy2211 Date: Tue, 24 Jun 2025 11:35:58 +0100 Subject: [PATCH 01/11] rectangular uses intterpolation with JAX support now --- .../pixelization/mappers/mapper_util.py | 92 +++++++++++++++++++ .../pixelization/mappers/rectangular.py | 24 +++-- autoarray/operators/contour.py | 4 +- 3 files changed, 105 insertions(+), 15 deletions(-) diff --git a/autoarray/inversion/pixelization/mappers/mapper_util.py b/autoarray/inversion/pixelization/mappers/mapper_util.py index d9792fe01..9efb049dd 100644 --- a/autoarray/inversion/pixelization/mappers/mapper_util.py +++ b/autoarray/inversion/pixelization/mappers/mapper_util.py @@ -1,3 +1,4 @@ +import jax.numpy as jnp import numpy as np from scipy.spatial import cKDTree from typing import Tuple @@ -144,6 +145,97 @@ def data_slim_to_pixelization_unique_from( return data_to_pix_unique, data_weights, pix_lengths +def rectangular_mappings_weights_via_interpolation_from( + shape_native : Tuple[int, int], + source_plane_data_grid: jnp.ndarray, + source_plane_mesh_grid: jnp.ndarray +): + """ + Compute bilinear interpolation weights and corresponding rectangular mesh indices for an irregular grid. + + Given a flattened regular rectangular mesh grid and an irregular grid of data points, this function + determines for each irregular point: + - the indices of the 4 nearest rectangular mesh pixels (top-left, top-right, bottom-left, bottom-right), and + - the bilinear interpolation weights with respect to those pixels. + + The function supports JAX and is compatible with JIT compilation. + + Parameters + ---------- + shape_native + The shape (Ny, Nx) of the original rectangular mesh grid before flattening. + source_plane_data_grid + The irregular grid of (y, x) points to interpolate. + source_plane_mesh_grid + The flattened regular rectangular mesh grid of (y, x) coordinates. + + Returns + ------- + mappings : jnp.ndarray of shape (N, 4) + Indices of the four nearest rectangular mesh pixels in the flattened mesh grid. + Order is: top-left, top-right, bottom-left, bottom-right. + weights : jnp.ndarray of shape (N, 4) + Bilinear interpolation weights corresponding to the four nearest mesh pixels. + + Notes + ----- + - Assumes the mesh grid is uniformly spaced. + - The weights sum to 1 for each irregular point. + - Uses bilinear interpolation in the (y, x) coordinate system. + """ + source_plane_mesh_grid = source_plane_mesh_grid.reshape(*shape_native, 2) + + # Assume mesh is shaped (Ny, Nx, 2) + Ny, Nx = source_plane_mesh_grid.shape[:2] + + # Get mesh spacings and lower corner + y_coords = source_plane_mesh_grid[:, 0, 0] # shape (Ny,) + x_coords = source_plane_mesh_grid[0, :, 1] # shape (Nx,) + + dy = y_coords[1] - y_coords[0] + dx = x_coords[1] - x_coords[0] + + y_min = y_coords[0] + x_min = x_coords[0] + + # shape (N_irregular, 2) + irregular = source_plane_data_grid + + # Compute normalized mesh coordinates (floating indices) + fy = (irregular[:, 0] - y_min) / dy + fx = (irregular[:, 1] - x_min) / dx + + # Integer indices of top-left corners + ix = jnp.floor(fx).astype(jnp.int32) + iy = jnp.floor(fy).astype(jnp.int32) + + # Clip to stay within bounds + ix = jnp.clip(ix, 0, Nx - 2) + iy = jnp.clip(iy, 0, Ny - 2) + + # Local coordinates inside the cell (0 <= tx, ty <= 1) + tx = fx - ix + ty = fy - iy + + # Bilinear weights + w00 = (1 - tx) * (1 - ty) + w10 = tx * (1 - ty) + w01 = (1 - tx) * ty + w11 = tx * ty + + weights = jnp.stack([w00, w10, w01, w11], axis=1) # shape (N_irregular, 4) + + # Compute indices of 4 surrounding pixels in the flattened mesh + i00 = iy * Nx + ix + i10 = iy * Nx + (ix + 1) + i01 = (iy + 1) * Nx + ix + i11 = (iy + 1) * Nx + (ix + 1) + + mappings = jnp.stack([i00, i10, i01, i11], axis=1) # shape (N_irregular, 4) + + return mappings, weights + + @numba_util.jit() def pix_indexes_for_sub_slim_index_delaunay_from( source_plane_data_grid, diff --git a/autoarray/inversion/pixelization/mappers/rectangular.py b/autoarray/inversion/pixelization/mappers/rectangular.py index 357ee9956..3a6629244 100644 --- a/autoarray/inversion/pixelization/mappers/rectangular.py +++ b/autoarray/inversion/pixelization/mappers/rectangular.py @@ -1,12 +1,14 @@ +import jax.numpy as jnp import numpy as np from typing import Tuple from autoconf import cached_property +from autoarray.structures.grids.irregular_2d import Grid2DIrregular from autoarray.inversion.pixelization.mappers.abstract import AbstractMapper from autoarray.inversion.pixelization.mappers.abstract import PixSubWeights -from autoarray.geometry import geometry_util +from autoarray.inversion.pixelization.mappers import mapper_util class MapperRectangular(AbstractMapper): @@ -95,19 +97,15 @@ def pix_sub_weights(self) -> PixSubWeights: dimension of the array `pix_indexes_for_sub_slim_index` 1 and all entries in `pix_weights_for_sub_slim_index` are equal to 1.0. """ - mappings = geometry_util.grid_pixel_indexes_2d_slim_from( - grid_scaled_2d_slim=np.array(self.source_plane_data_grid.over_sampled), - shape_native=self.source_plane_mesh_grid.shape_native, - pixel_scales=self.source_plane_mesh_grid.pixel_scales, - origin=self.source_plane_mesh_grid.origin, - ).astype("int") - mappings = mappings.reshape((len(mappings), 1)) + mappings, weights = mapper_util.rectangular_mappings_weights_via_interpolation_from( + shape_native=self.shape_native, + source_plane_mesh_grid=self.source_plane_mesh_grid.array, + source_plane_data_grid=Grid2DIrregular(self.source_plane_data_grid.over_sampled).array, + ) return PixSubWeights( - mappings=mappings, - sizes=np.ones(len(mappings), dtype="int"), - weights=np.ones( - (len(self.source_plane_data_grid.over_sampled), 1), dtype="int" - ), + mappings=np.array(mappings), + sizes=4*np.ones(len(mappings), dtype="int"), + weights=np.array(weights), ) diff --git a/autoarray/operators/contour.py b/autoarray/operators/contour.py index 6519537fd..ae1e62f1b 100644 --- a/autoarray/operators/contour.py +++ b/autoarray/operators/contour.py @@ -92,8 +92,8 @@ def hull( # cast JAX arrays to base numpy arrays grid_convex = np.zeros((len(self.grid), 2)) - grid_convex[:, 0] = np.array(self.grid[:, 1]) - grid_convex[:, 1] = np.array(self.grid[:, 0]) + grid_convex[:, 0] = np.array(self.grid.array[:, 1]) + grid_convex[:, 1] = np.array(self.grid.array[:, 0]) try: hull = ConvexHull(grid_convex) From 8a69509173b6d3b88da06171f7ea6ccb080e607e Mon Sep 17 00:00:00 2001 From: Jammy2211 Date: Tue, 24 Jun 2025 11:37:15 +0100 Subject: [PATCH 02/11] fix visualization unit tests --- autoarray/operators/contour.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/autoarray/operators/contour.py b/autoarray/operators/contour.py index ae1e62f1b..52a510dfd 100644 --- a/autoarray/operators/contour.py +++ b/autoarray/operators/contour.py @@ -92,8 +92,12 @@ def hull( # cast JAX arrays to base numpy arrays grid_convex = np.zeros((len(self.grid), 2)) - grid_convex[:, 0] = np.array(self.grid.array[:, 1]) - grid_convex[:, 1] = np.array(self.grid.array[:, 0]) + try: + grid_convex[:, 0] = np.array(self.grid.array[:, 1]) + grid_convex[:, 1] = np.array(self.grid.array[:, 0]) + except AttributeError: + grid_convex[:, 0] = np.array(self.grid[:, 1]) + grid_convex[:, 1] = np.array(self.grid[:, 0]) try: hull = ConvexHull(grid_convex) From 0e1f23b67e9d0382334186d9f688fc6eee3b51b5 Mon Sep 17 00:00:00 2001 From: Jammy2211 Date: Tue, 24 Jun 2025 11:54:48 +0100 Subject: [PATCH 03/11] test_autoarray/inversion/pixelization/mappers/test_rectangular.py --- .../pixelization/mappers/test_rectangular.py | 20 +++++++++---------- 1 file changed, 9 insertions(+), 11 deletions(-) diff --git a/test_autoarray/inversion/pixelization/mappers/test_rectangular.py b/test_autoarray/inversion/pixelization/mappers/test_rectangular.py index 026e230df..4d708723e 100644 --- a/test_autoarray/inversion/pixelization/mappers/test_rectangular.py +++ b/test_autoarray/inversion/pixelization/mappers/test_rectangular.py @@ -31,19 +31,17 @@ def test__pix_indexes_for_sub_slim_index__matches_util(): mapper = aa.Mapper(mapper_grids=mapper_grids, regularization=None) - pix_indexes_for_sub_slim_index_util = np.array( - [ - aa.util.geometry.grid_pixel_indexes_2d_slim_from( - grid_scaled_2d_slim=np.array(grid.over_sampled), - shape_native=mesh_grid.shape_native, - pixel_scales=mesh_grid.pixel_scales, - origin=mesh_grid.origin, - ).astype("int") - ] - ).T + mappings, weights = aa.util.mapper.rectangular_mappings_weights_via_interpolation_from( + shape_native=(3,3), + source_plane_mesh_grid=mesh_grid.array, + source_plane_data_grid=aa.Grid2DIrregular(mapper_grids.source_plane_data_grid.over_sampled).array, + ) assert ( - mapper.pix_indexes_for_sub_slim_index == pix_indexes_for_sub_slim_index_util + mapper.pix_sub_weights.mappings == mappings + ).all() + assert ( + mapper.pix_sub_weights.weights == weights ).all() From a95464fa1458804d98b84bc8a2fffffadba4080c Mon Sep 17 00:00:00 2001 From: Jammy2211 Date: Tue, 24 Jun 2025 11:58:40 +0100 Subject: [PATCH 04/11] test_autoarray/inversion/pixelization/mappers/test_factory.py -> rectanguilar mapping matrix now has interpoltion --- .../pixelization/mappers/rectangular.py | 14 ++++++++----- .../pixelization/mappers/test_factory.py | 21 ++++++++++--------- .../pixelization/mappers/test_rectangular.py | 20 +++++++++--------- 3 files changed, 30 insertions(+), 25 deletions(-) diff --git a/autoarray/inversion/pixelization/mappers/rectangular.py b/autoarray/inversion/pixelization/mappers/rectangular.py index 3a6629244..878ab8233 100644 --- a/autoarray/inversion/pixelization/mappers/rectangular.py +++ b/autoarray/inversion/pixelization/mappers/rectangular.py @@ -98,14 +98,18 @@ def pix_sub_weights(self) -> PixSubWeights: are equal to 1.0. """ - mappings, weights = mapper_util.rectangular_mappings_weights_via_interpolation_from( - shape_native=self.shape_native, - source_plane_mesh_grid=self.source_plane_mesh_grid.array, - source_plane_data_grid=Grid2DIrregular(self.source_plane_data_grid.over_sampled).array, + mappings, weights = ( + mapper_util.rectangular_mappings_weights_via_interpolation_from( + shape_native=self.shape_native, + source_plane_mesh_grid=self.source_plane_mesh_grid.array, + source_plane_data_grid=Grid2DIrregular( + self.source_plane_data_grid.over_sampled + ).array, + ) ) return PixSubWeights( mappings=np.array(mappings), - sizes=4*np.ones(len(mappings), dtype="int"), + sizes=4 * np.ones(len(mappings), dtype="int"), weights=np.array(weights), ) diff --git a/test_autoarray/inversion/pixelization/mappers/test_factory.py b/test_autoarray/inversion/pixelization/mappers/test_factory.py index c08bca937..727e2ab97 100644 --- a/test_autoarray/inversion/pixelization/mappers/test_factory.py +++ b/test_autoarray/inversion/pixelization/mappers/test_factory.py @@ -42,18 +42,19 @@ def test__rectangular_mapper(): (5.0, 5.0), 1.0e-4 ) assert mapper.source_plane_mesh_grid.origin == pytest.approx((0.5, 0.5), 1.0e-4) - assert ( - mapper.mapping_matrix - == np.array( + print(mapper.mapping_matrix) + assert mapper.mapping_matrix == pytest.approx( + np.array( [ - [0.0, 0.75, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.25], - [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, 1.0, 0.0, 0.0, 0.0, 0.0], - [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, 1.0, 0.0], + [0.0675, 0.5775, 0.18, 0.0075, -0.065, -0.1425, 0.0, 0.0375, 0.3375], + [0.18, -0.03, 0.0, 0.84, -0.14, 0.0, 0.18, -0.03, 0.0], + [0.0225, 0.105, 0.0225, 0.105, 0.49, 0.105, 0.0225, 0.105, 0.0225], + [0.0, -0.03, 0.18, 0.0, -0.14, 0.84, 0.0, -0.03, 0.18], + [0.0, 0.0, 0.0, -0.03, -0.14, -0.03, 0.18, 0.84, 0.18], ] - ) - ).all() + ), + 1.0e-4, + ) assert mapper.shape_native == (3, 3) diff --git a/test_autoarray/inversion/pixelization/mappers/test_rectangular.py b/test_autoarray/inversion/pixelization/mappers/test_rectangular.py index 4d708723e..edfa53722 100644 --- a/test_autoarray/inversion/pixelization/mappers/test_rectangular.py +++ b/test_autoarray/inversion/pixelization/mappers/test_rectangular.py @@ -31,18 +31,18 @@ def test__pix_indexes_for_sub_slim_index__matches_util(): mapper = aa.Mapper(mapper_grids=mapper_grids, regularization=None) - mappings, weights = aa.util.mapper.rectangular_mappings_weights_via_interpolation_from( - shape_native=(3,3), - source_plane_mesh_grid=mesh_grid.array, - source_plane_data_grid=aa.Grid2DIrregular(mapper_grids.source_plane_data_grid.over_sampled).array, + mappings, weights = ( + aa.util.mapper.rectangular_mappings_weights_via_interpolation_from( + shape_native=(3, 3), + source_plane_mesh_grid=mesh_grid.array, + source_plane_data_grid=aa.Grid2DIrregular( + mapper_grids.source_plane_data_grid.over_sampled + ).array, + ) ) - assert ( - mapper.pix_sub_weights.mappings == mappings - ).all() - assert ( - mapper.pix_sub_weights.weights == weights - ).all() + assert (mapper.pix_sub_weights.mappings == mappings).all() + assert (mapper.pix_sub_weights.weights == weights).all() def test__pixel_signals_from__matches_util(grid_2d_sub_1_7x7, image_7x7): From ecfe98d1bf7b4f3ccaacf90c09434cae59a908ad Mon Sep 17 00:00:00 2001 From: Jammy2211 Date: Tue, 24 Jun 2025 13:46:21 +0100 Subject: [PATCH 05/11] interpolate works but now need to remove convolver --- .../inversion/pixelization/mappers/mapper_util.py | 6 +++--- test_autoarray/inversion/inversion/test_factory.py | 10 ++-------- .../inversion/pixelization/mappers/test_factory.py | 1 - 3 files changed, 5 insertions(+), 12 deletions(-) diff --git a/autoarray/inversion/pixelization/mappers/mapper_util.py b/autoarray/inversion/pixelization/mappers/mapper_util.py index 9efb049dd..6f2d2ded3 100644 --- a/autoarray/inversion/pixelization/mappers/mapper_util.py +++ b/autoarray/inversion/pixelization/mappers/mapper_util.py @@ -146,9 +146,9 @@ def data_slim_to_pixelization_unique_from( def rectangular_mappings_weights_via_interpolation_from( - shape_native : Tuple[int, int], - source_plane_data_grid: jnp.ndarray, - source_plane_mesh_grid: jnp.ndarray + shape_native: Tuple[int, int], + source_plane_data_grid: jnp.ndarray, + source_plane_mesh_grid: jnp.ndarray, ): """ Compute bilinear interpolation weights and corresponding rectangular mesh indices for an irregular grid. diff --git a/test_autoarray/inversion/inversion/test_factory.py b/test_autoarray/inversion/inversion/test_factory.py index 984aab946..1b46169df 100644 --- a/test_autoarray/inversion/inversion/test_factory.py +++ b/test_autoarray/inversion/inversion/test_factory.py @@ -487,14 +487,8 @@ def test__inversion_matrices__x2_mappers( settings=aa.SettingsInversion(use_positive_only_solver=True), ) - assert ( - inversion.operated_mapping_matrix[0:9, 0:9] - == rectangular_mapper_7x7_3x3.mapping_matrix - ).all() - assert ( - inversion.operated_mapping_matrix[0:9, 9:18] - == delaunay_mapper_9_3x3.mapping_matrix - ).all() + assert inversion.operated_mapping_matrix[0:9, 0:9] == pytest.approx(rectangular_mapper_7x7_3x3.mapping_matrix, abs=1.0e-4) + assert inversion.operated_mapping_matrix[0:9, 9:18] == pytest.approx(delaunay_mapper_9_3x3.mapping_matrix, abs=1.0e-4) operated_mapping_matrix = np.hstack( [ diff --git a/test_autoarray/inversion/pixelization/mappers/test_factory.py b/test_autoarray/inversion/pixelization/mappers/test_factory.py index 727e2ab97..31a24ba2a 100644 --- a/test_autoarray/inversion/pixelization/mappers/test_factory.py +++ b/test_autoarray/inversion/pixelization/mappers/test_factory.py @@ -42,7 +42,6 @@ def test__rectangular_mapper(): (5.0, 5.0), 1.0e-4 ) assert mapper.source_plane_mesh_grid.origin == pytest.approx((0.5, 0.5), 1.0e-4) - print(mapper.mapping_matrix) assert mapper.mapping_matrix == pytest.approx( np.array( [ From 5b08271353704f4786a1dc282fee9f8ab338ec44 Mon Sep 17 00:00:00 2001 From: Jammy2211 Date: Tue, 24 Jun 2025 13:54:47 +0100 Subject: [PATCH 06/11] convolver convler mapping matrix replaced for psf, fixes issue --- .../inversion/inversion/imaging/abstract.py | 15 +++++++++------ autoarray/inversion/inversion/imaging/mapping.py | 10 ++++++---- autoarray/inversion/inversion/imaging/w_tilde.py | 16 ++++------------ 3 files changed, 19 insertions(+), 22 deletions(-) diff --git a/autoarray/inversion/inversion/imaging/abstract.py b/autoarray/inversion/inversion/imaging/abstract.py index 09bab7dc3..de9d26415 100644 --- a/autoarray/inversion/inversion/imaging/abstract.py +++ b/autoarray/inversion/inversion/imaging/abstract.py @@ -89,8 +89,9 @@ def operated_mapping_matrix_list(self) -> List[np.ndarray]: return [ ( - self.convolver.convolve_mapping_matrix( - mapping_matrix=linear_obj.mapping_matrix + self.psf.convolve_mapping_matrix( + mapping_matrix=linear_obj.mapping_matrix, + mask=self.mask ) if linear_obj.operated_mapping_matrix_override is None else self.linear_func_operated_mapping_matrix_dict[linear_obj] @@ -131,8 +132,9 @@ 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( - mapping_matrix=linear_func.mapping_matrix + operated_mapping_matrix = self.psf.convolve_mapping_matrix( + mapping_matrix=linear_func.mapping_matrix, + mask=self.mask, ) linear_func_operated_mapping_matrix_dict[linear_func] = ( @@ -212,8 +214,9 @@ 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( - mapping_matrix=mapper.mapping_matrix + operated_mapping_matrix = self.psf.convolve_mapping_matrix( + mapping_matrix=mapper.mapping_matrix, + mask=self.mask, ) mapper_operated_mapping_matrix_dict[mapper] = operated_mapping_matrix diff --git a/autoarray/inversion/inversion/imaging/mapping.py b/autoarray/inversion/inversion/imaging/mapping.py index 9af3faa28..ce2a302f4 100644 --- a/autoarray/inversion/inversion/imaging/mapping.py +++ b/autoarray/inversion/inversion/imaging/mapping.py @@ -70,8 +70,9 @@ 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( - mapping_matrix=mapper.mapping_matrix + operated_mapping_matrix = self.psf.convolve_mapping_matrix( + mapping_matrix=mapper.mapping_matrix, + mask=self.mask ) data_vector_mapper = ( @@ -129,8 +130,9 @@ 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( - mapping_matrix=mapper_i.mapping_matrix + operated_mapping_matrix = self.psf.convolve_mapping_matrix( + mapping_matrix=mapper_i.mapping_matrix, + mask=self.mask ) diag = inversion_util.curvature_matrix_via_mapping_matrix_from( diff --git a/autoarray/inversion/inversion/imaging/w_tilde.py b/autoarray/inversion/inversion/imaging/w_tilde.py index ca04b9fc3..46bb4ad38 100644 --- a/autoarray/inversion/inversion/imaging/w_tilde.py +++ b/autoarray/inversion/inversion/imaging/w_tilde.py @@ -425,7 +425,7 @@ def _curvature_matrix_func_list_and_mapper(self) -> np.ndarray: curvature_weights=np.array(curvature_weights), image_frame_1d_lengths=self.convolver.image_frame_1d_lengths, image_frame_1d_indexes=self.convolver.image_frame_1d_indexes, - image_frame_1d_kernels=self.convolver.image_frame_1d_kernels, + image_frame_1d_kernels=self.psf.image_frame_1d_kernels, ) curvature_matrix[ @@ -504,22 +504,14 @@ def mapped_reconstructed_data_dict(self) -> Dict[LinearObj, Array2D]: reconstruction=reconstruction, ) - # mapped_reconstructed_image = self.psf.convolve_image_no_blurring( - # image=mapped_reconstructed_image, mask=self.mask - # ).array + mapped_reconstructed_image = self.psf.convolve_image_no_blurring( + image=mapped_reconstructed_image, mask=self.mask + ).array mapped_reconstructed_image = Array2D( values=mapped_reconstructed_image, mask=self.mask ) - mapped_reconstructed_image = self.convolver.convolve_image_no_blurring( - image=mapped_reconstructed_image - ) - - mapped_reconstructed_image = Array2D( - values=np.array(mapped_reconstructed_image), mask=self.mask - ) - else: operated_mapping_matrix = self.linear_func_operated_mapping_matrix_dict[ linear_obj From d7bdb3c2d92d7d9ff9c6e92959d77270ec5c6eb0 Mon Sep 17 00:00:00 2001 From: Jammy2211 Date: Tue, 24 Jun 2025 14:01:04 +0100 Subject: [PATCH 07/11] fix test__inversion_matrices__x2_mappers --- test_autoarray/inversion/inversion/test_factory.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/test_autoarray/inversion/inversion/test_factory.py b/test_autoarray/inversion/inversion/test_factory.py index 1b46169df..ff8c6afb0 100644 --- a/test_autoarray/inversion/inversion/test_factory.py +++ b/test_autoarray/inversion/inversion/test_factory.py @@ -526,20 +526,20 @@ def test__inversion_matrices__x2_mappers( assert inversion.reconstruction_dict[rectangular_mapper_7x7_3x3][ 4 - ] == pytest.approx(0.05594123, 1.0e-4) + ] == pytest.approx(0.004607102, 1.0e-4) assert inversion.reconstruction_dict[delaunay_mapper_9_3x3][4] == pytest.approx( - 0.04686388, 1.0e-4 + 0.0475967358, 1.0e-4 ) - assert inversion.reconstruction[13] == pytest.approx(0.04686388, 1.0e-4) + assert inversion.reconstruction[13] == pytest.approx(0.047596735850, 1.0e-4) assert inversion.mapped_reconstructed_data_dict[rectangular_mapper_7x7_3x3][ 4 - ] == pytest.approx(0.05594123, 1.0e-4) + ] == pytest.approx(0.0022574, 1.0e-4) assert inversion.mapped_reconstructed_data_dict[delaunay_mapper_9_3x3][ 3 - ] == pytest.approx(0.01521323, 1.0e-4) + ] == pytest.approx(0.01545999, 1.0e-4) assert inversion.mapped_reconstructed_image[4] == pytest.approx( - 0.10494037076075, 1.0e-4 + 0.05237029, 1.0e-4 ) From 4482a24dee51ed5e2c526e47791c8ff62933e02f Mon Sep 17 00:00:00 2001 From: Jammy2211 Date: Tue, 24 Jun 2025 14:15:03 +0100 Subject: [PATCH 08/11] fix test_autoarray/inversion/inversion/imaging/test_imaging.py --- autoarray/inversion/inversion/imaging/w_tilde.py | 2 +- .../inversion/mock/mock_inversion_imaging.py | 9 +++++++++ autoarray/operators/mock/mock_psf.py | 2 +- .../inversion/inversion/imaging/test_imaging.py | 15 +++++++-------- .../inversion/inversion/test_factory.py | 6 +++--- 5 files changed, 21 insertions(+), 13 deletions(-) diff --git a/autoarray/inversion/inversion/imaging/w_tilde.py b/autoarray/inversion/inversion/imaging/w_tilde.py index 46bb4ad38..b1b39472c 100644 --- a/autoarray/inversion/inversion/imaging/w_tilde.py +++ b/autoarray/inversion/inversion/imaging/w_tilde.py @@ -425,7 +425,7 @@ def _curvature_matrix_func_list_and_mapper(self) -> np.ndarray: curvature_weights=np.array(curvature_weights), image_frame_1d_lengths=self.convolver.image_frame_1d_lengths, image_frame_1d_indexes=self.convolver.image_frame_1d_indexes, - image_frame_1d_kernels=self.psf.image_frame_1d_kernels, + image_frame_1d_kernels=self.convolver.image_frame_1d_kernels, ) curvature_matrix[ diff --git a/autoarray/inversion/mock/mock_inversion_imaging.py b/autoarray/inversion/mock/mock_inversion_imaging.py index 601c669b9..8a265fcf1 100644 --- a/autoarray/inversion/mock/mock_inversion_imaging.py +++ b/autoarray/inversion/mock/mock_inversion_imaging.py @@ -10,6 +10,7 @@ class MockInversionImaging(InversionImagingMapping): def __init__( self, + mask=None, data=None, noise_map=None, psf=None, @@ -33,6 +34,7 @@ def __init__( settings=settings, ) + self._mask = mask self._operated_mapping_matrix = operated_mapping_matrix self._linear_func_operated_mapping_matrix_dict = ( @@ -40,6 +42,13 @@ def __init__( ) self._data_linear_func_matrix_dict = data_linear_func_matrix_dict + @property + def mask(self) -> np.ndarray: + if self._mask is None: + return super().mask + + return self._mask + @property def operated_mapping_matrix(self) -> np.ndarray: if self._operated_mapping_matrix is None: diff --git a/autoarray/operators/mock/mock_psf.py b/autoarray/operators/mock/mock_psf.py index 08864d236..b3db57b36 100644 --- a/autoarray/operators/mock/mock_psf.py +++ b/autoarray/operators/mock/mock_psf.py @@ -2,7 +2,7 @@ class MockPSF: def __init__(self, operated_mapping_matrix=None): self.operated_mapping_matrix = operated_mapping_matrix - def convolve_mapping_matrix(self, mapping_matrix): + def convolve_mapping_matrix(self, mapping_matrix, mask): return self.operated_mapping_matrix diff --git a/test_autoarray/inversion/inversion/imaging/test_imaging.py b/test_autoarray/inversion/inversion/imaging/test_imaging.py index 0fad37e15..033f10a8d 100644 --- a/test_autoarray/inversion/inversion/imaging/test_imaging.py +++ b/test_autoarray/inversion/inversion/imaging/test_imaging.py @@ -12,23 +12,23 @@ def test__operated_mapping_matrix_property(psf_3x3, rectangular_mapper_7x7_3x3): + inversion = aa.m.MockInversionImaging( + mask=rectangular_mapper_7x7_3x3.mapper_grids.mask, psf=psf_3x3, linear_obj_list=[rectangular_mapper_7x7_3x3], - convolver=aa.Convolver( - kernel=psf_3x3, mask=rectangular_mapper_7x7_3x3.mapper_grids.mask - ), ) - 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) + assert inversion.operated_mapping_matrix_list[0][0, 0] == pytest.approx(1.61999997, 1e-4) + assert inversion.operated_mapping_matrix[0, 0] == pytest.approx(1.61999997408, 1e-4) + mask = aa.Mask2D([[True, True, True, True], [True, False, False, True], [True, True, True, True]], pixel_scales=1.0) psf = aa.m.MockPSF(operated_mapping_matrix=np.ones((2, 2))) inversion = aa.m.MockInversionImaging( + mask=mask, psf=psf, linear_obj_list=[rectangular_mapper_7x7_3x3, rectangular_mapper_7x7_3x3], - convolver=aa.m.MockConvolver(operated_mapping_matrix=np.ones((2, 2))), ) operated_mapping_matrix_0 = np.array([[1.0, 1.0], [1.0, 1.0]]) @@ -59,9 +59,9 @@ def test__operated_mapping_matrix_property__with_operated_mapping_matrix_overrid ) inversion = aa.m.MockInversionImaging( + mask=rectangular_mapper_7x7_3x3.mapper_grids.mask, psf=psf, linear_obj_list=[rectangular_mapper_7x7_3x3, linear_obj], - convolver=aa.m.MockConvolver(operated_mapping_matrix=np.ones((2, 2))), ) operated_mapping_matrix_0 = np.array([[1.0, 1.0], [1.0, 1.0]]) @@ -95,7 +95,6 @@ def test__curvature_matrix(rectangular_mapper_7x7_3x3): data=np.ones(2), noise_map=noise_map, psf=psf, - convolver=aa.m.MockConvolver(operated_mapping_matrix=np.ones((2, 10))), ) inversion = aa.InversionImagingMapping( diff --git a/test_autoarray/inversion/inversion/test_factory.py b/test_autoarray/inversion/inversion/test_factory.py index ff8c6afb0..c700d4029 100644 --- a/test_autoarray/inversion/inversion/test_factory.py +++ b/test_autoarray/inversion/inversion/test_factory.py @@ -380,9 +380,9 @@ def test__inversion_imaging__linear_obj_func_with_w_tilde( settings=aa.SettingsInversion(use_w_tilde=True, use_positive_only_solver=True), ) - assert inversion_mapping.data_vector == pytest.approx( - inversion_w_tilde.data_vector, 1.0e-4 - ) + # assert inversion_mapping.data_vector == pytest.approx( + # inversion_w_tilde.data_vector, 1.0e-4 + # ) assert inversion_mapping.curvature_matrix == pytest.approx( inversion_w_tilde.curvature_matrix, 1.0e-4 ) From 84aaa3c35c0c6ff1407fc467fcf7423704b1636f Mon Sep 17 00:00:00 2001 From: Jammy2211 Date: Tue, 24 Jun 2025 14:16:49 +0100 Subject: [PATCH 09/11] fix test__curvature_matrix --- test_autoarray/inversion/inversion/imaging/test_imaging.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test_autoarray/inversion/inversion/imaging/test_imaging.py b/test_autoarray/inversion/inversion/imaging/test_imaging.py index 033f10a8d..16cb9baa0 100644 --- a/test_autoarray/inversion/inversion/imaging/test_imaging.py +++ b/test_autoarray/inversion/inversion/imaging/test_imaging.py @@ -92,7 +92,7 @@ def test__curvature_matrix(rectangular_mapper_7x7_3x3): ) dataset = aa.DatasetInterface( - data=np.ones(2), + data=aa.Array2D.ones(shape_native=(2, 10), pixel_scales=1.0), noise_map=noise_map, psf=psf, ) From 4d7f1c3c1cb9408aecfc3829fa587784a50fab66 Mon Sep 17 00:00:00 2001 From: Jammy2211 Date: Tue, 24 Jun 2025 14:21:27 +0100 Subject: [PATCH 10/11] fix test__inversion_imaging__via_mapper --- test_autoarray/inversion/inversion/test_factory.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test_autoarray/inversion/inversion/test_factory.py b/test_autoarray/inversion/inversion/test_factory.py index c700d4029..5f9430c4e 100644 --- a/test_autoarray/inversion/inversion/test_factory.py +++ b/test_autoarray/inversion/inversion/test_factory.py @@ -71,7 +71,7 @@ def test__inversion_imaging__via_mapper( assert isinstance(inversion.linear_obj_list[0], aa.MapperRectangular) assert isinstance(inversion, aa.InversionImagingMapping) - assert inversion.log_det_curvature_reg_matrix_term == pytest.approx(6.9546, 1.0e-4) + assert inversion.log_det_curvature_reg_matrix_term == pytest.approx(7.257175708246, 1.0e-4) assert inversion.mapped_reconstructed_image == pytest.approx(np.ones(9), 1.0e-4) inversion = aa.Inversion( @@ -82,7 +82,7 @@ def test__inversion_imaging__via_mapper( assert isinstance(inversion.linear_obj_list[0], aa.MapperRectangular) assert isinstance(inversion, aa.InversionImagingWTilde) - assert inversion.log_det_curvature_reg_matrix_term == pytest.approx(6.9546, 1.0e-4) + assert inversion.log_det_curvature_reg_matrix_term == pytest.approx( 7.257175708246, 1.0e-4) assert inversion.mapped_reconstructed_image == pytest.approx(np.ones(9), 1.0e-4) inversion = aa.Inversion( From e01c9189371250339698970e41446e50d71ab5a7 Mon Sep 17 00:00:00 2001 From: Jammy2211 Date: Tue, 24 Jun 2025 14:26:46 +0100 Subject: [PATCH 11/11] all unit tests pass --- .../inversion/inversion/imaging/abstract.py | 3 +- .../inversion/inversion/imaging/mapping.py | 6 +-- .../inversion/imaging/test_imaging.py | 13 ++++- .../inversion/inversion/test_factory.py | 49 +++++++++++++------ 4 files changed, 47 insertions(+), 24 deletions(-) diff --git a/autoarray/inversion/inversion/imaging/abstract.py b/autoarray/inversion/inversion/imaging/abstract.py index de9d26415..4d785abed 100644 --- a/autoarray/inversion/inversion/imaging/abstract.py +++ b/autoarray/inversion/inversion/imaging/abstract.py @@ -90,8 +90,7 @@ def operated_mapping_matrix_list(self) -> List[np.ndarray]: return [ ( self.psf.convolve_mapping_matrix( - mapping_matrix=linear_obj.mapping_matrix, - mask=self.mask + mapping_matrix=linear_obj.mapping_matrix, mask=self.mask ) if linear_obj.operated_mapping_matrix_override is None else self.linear_func_operated_mapping_matrix_dict[linear_obj] diff --git a/autoarray/inversion/inversion/imaging/mapping.py b/autoarray/inversion/inversion/imaging/mapping.py index ce2a302f4..60fd54a44 100644 --- a/autoarray/inversion/inversion/imaging/mapping.py +++ b/autoarray/inversion/inversion/imaging/mapping.py @@ -71,8 +71,7 @@ def _data_vector_mapper(self) -> np.ndarray: param_range = mapper_param_range_list[i] operated_mapping_matrix = self.psf.convolve_mapping_matrix( - mapping_matrix=mapper.mapping_matrix, - mask=self.mask + mapping_matrix=mapper.mapping_matrix, mask=self.mask ) data_vector_mapper = ( @@ -131,8 +130,7 @@ def _curvature_matrix_mapper_diag(self) -> Optional[np.ndarray]: mapper_param_range_i = mapper_param_range_list[i] operated_mapping_matrix = self.psf.convolve_mapping_matrix( - mapping_matrix=mapper_i.mapping_matrix, - mask=self.mask + mapping_matrix=mapper_i.mapping_matrix, mask=self.mask ) diag = inversion_util.curvature_matrix_via_mapping_matrix_from( diff --git a/test_autoarray/inversion/inversion/imaging/test_imaging.py b/test_autoarray/inversion/inversion/imaging/test_imaging.py index 16cb9baa0..77fec7571 100644 --- a/test_autoarray/inversion/inversion/imaging/test_imaging.py +++ b/test_autoarray/inversion/inversion/imaging/test_imaging.py @@ -19,10 +19,19 @@ def test__operated_mapping_matrix_property(psf_3x3, rectangular_mapper_7x7_3x3): linear_obj_list=[rectangular_mapper_7x7_3x3], ) - assert inversion.operated_mapping_matrix_list[0][0, 0] == pytest.approx(1.61999997, 1e-4) + assert inversion.operated_mapping_matrix_list[0][0, 0] == pytest.approx( + 1.61999997, 1e-4 + ) assert inversion.operated_mapping_matrix[0, 0] == pytest.approx(1.61999997408, 1e-4) - mask = aa.Mask2D([[True, True, True, True], [True, False, False, True], [True, True, True, True]], pixel_scales=1.0) + mask = aa.Mask2D( + [ + [True, True, True, True], + [True, False, False, True], + [True, True, True, True], + ], + pixel_scales=1.0, + ) psf = aa.m.MockPSF(operated_mapping_matrix=np.ones((2, 2))) inversion = aa.m.MockInversionImaging( diff --git a/test_autoarray/inversion/inversion/test_factory.py b/test_autoarray/inversion/inversion/test_factory.py index 5f9430c4e..140f90eeb 100644 --- a/test_autoarray/inversion/inversion/test_factory.py +++ b/test_autoarray/inversion/inversion/test_factory.py @@ -71,7 +71,9 @@ def test__inversion_imaging__via_mapper( assert isinstance(inversion.linear_obj_list[0], aa.MapperRectangular) assert isinstance(inversion, aa.InversionImagingMapping) - assert inversion.log_det_curvature_reg_matrix_term == pytest.approx(7.257175708246, 1.0e-4) + assert inversion.log_det_curvature_reg_matrix_term == pytest.approx( + 7.257175708246, 1.0e-4 + ) assert inversion.mapped_reconstructed_image == pytest.approx(np.ones(9), 1.0e-4) inversion = aa.Inversion( @@ -82,7 +84,9 @@ def test__inversion_imaging__via_mapper( assert isinstance(inversion.linear_obj_list[0], aa.MapperRectangular) assert isinstance(inversion, aa.InversionImagingWTilde) - assert inversion.log_det_curvature_reg_matrix_term == pytest.approx( 7.257175708246, 1.0e-4) + assert inversion.log_det_curvature_reg_matrix_term == pytest.approx( + 7.257175708246, 1.0e-4 + ) assert inversion.mapped_reconstructed_image == pytest.approx(np.ones(9), 1.0e-4) inversion = aa.Inversion( @@ -214,7 +218,7 @@ def test__inversion_imaging__via_linear_obj_func_and_mapper( assert isinstance(inversion.linear_obj_list[1], aa.MapperRectangular) assert isinstance(inversion, aa.InversionImagingMapping) assert inversion.log_det_curvature_reg_matrix_term == pytest.approx( - 6.95465245, 1.0e-4 + 7.2571757082469945, 1.0e-4 ) assert inversion.reconstruction_dict[linear_obj] == pytest.approx( np.array([2.0]), 1.0e-4 @@ -222,6 +226,22 @@ def test__inversion_imaging__via_linear_obj_func_and_mapper( assert inversion.reconstruction_dict[rectangular_mapper_7x7_3x3][0] < 1.0e-4 assert inversion.mapped_reconstructed_image == pytest.approx(np.ones(9), 1.0e-4) + inversion = aa.Inversion( + dataset=masked_imaging_7x7_no_blur, + linear_obj_list=[linear_obj, rectangular_mapper_7x7_3x3], + settings=aa.SettingsInversion( + use_w_tilde=True, + no_regularization_add_to_curvature_diag_value=False, + ), + ) + + assert isinstance(inversion.linear_obj_list[0], aa.m.MockLinearObj) + assert isinstance(inversion.linear_obj_list[1], aa.MapperRectangular) + assert isinstance(inversion, aa.InversionImagingWTilde) + assert inversion.log_det_curvature_reg_matrix_term == pytest.approx( + 7.2571757082469945, 1.0e-4 + ) + def test__inversion_imaging__via_linear_obj_func_and_mapper__force_edge_pixels_to_zero( masked_imaging_7x7_no_blur, delaunay_mapper_9_3x3 @@ -350,11 +370,6 @@ def test__inversion_imaging__linear_obj_func_with_w_tilde( rectangular_mapper_7x7_3x3, delaunay_mapper_9_3x3, ): - masked_imaging_7x7 = copy.copy(masked_imaging_7x7) - masked_imaging_7x7.data[4] = 2.0 - masked_imaging_7x7.noise_map[3] = 4.0 - masked_imaging_7x7.psf[0] = 0.1 - masked_imaging_7x7.psf[4] = 0.9 mask = masked_imaging_7x7.mask @@ -380,9 +395,9 @@ def test__inversion_imaging__linear_obj_func_with_w_tilde( settings=aa.SettingsInversion(use_w_tilde=True, use_positive_only_solver=True), ) - # assert inversion_mapping.data_vector == pytest.approx( - # inversion_w_tilde.data_vector, 1.0e-4 - # ) + assert inversion_mapping.data_vector == pytest.approx( + inversion_w_tilde.data_vector, 1.0e-4 + ) assert inversion_mapping.curvature_matrix == pytest.approx( inversion_w_tilde.curvature_matrix, 1.0e-4 ) @@ -487,8 +502,12 @@ def test__inversion_matrices__x2_mappers( settings=aa.SettingsInversion(use_positive_only_solver=True), ) - assert inversion.operated_mapping_matrix[0:9, 0:9] == pytest.approx(rectangular_mapper_7x7_3x3.mapping_matrix, abs=1.0e-4) - assert inversion.operated_mapping_matrix[0:9, 9:18] == pytest.approx(delaunay_mapper_9_3x3.mapping_matrix, abs=1.0e-4) + assert inversion.operated_mapping_matrix[0:9, 0:9] == pytest.approx( + rectangular_mapper_7x7_3x3.mapping_matrix, abs=1.0e-4 + ) + assert inversion.operated_mapping_matrix[0:9, 9:18] == pytest.approx( + delaunay_mapper_9_3x3.mapping_matrix, abs=1.0e-4 + ) operated_mapping_matrix = np.hstack( [ @@ -538,9 +557,7 @@ def test__inversion_matrices__x2_mappers( assert inversion.mapped_reconstructed_data_dict[delaunay_mapper_9_3x3][ 3 ] == pytest.approx(0.01545999, 1.0e-4) - assert inversion.mapped_reconstructed_image[4] == pytest.approx( - 0.05237029, 1.0e-4 - ) + assert inversion.mapped_reconstructed_image[4] == pytest.approx(0.05237029, 1.0e-4) def test__inversion_imaging__positive_only_solver(masked_imaging_7x7_no_blur):