From 8f0079570a17220b5fcefc5c3996342ba5bcd49d Mon Sep 17 00:00:00 2001 From: Jammy2211 Date: Mon, 15 Sep 2025 12:34:25 +0100 Subject: [PATCH 01/17] adpative stuff implemented --- .../pixelization/mappers/mapper_util.py | 141 ++++++++++++++++++ .../pixelization/mappers/rectangular.py | 32 ++-- .../inversion/inversion/test_factory.py | 2 +- 3 files changed, 153 insertions(+), 22 deletions(-) diff --git a/autoarray/inversion/pixelization/mappers/mapper_util.py b/autoarray/inversion/pixelization/mappers/mapper_util.py index c3e8d470a..3ef01e479 100644 --- a/autoarray/inversion/pixelization/mappers/mapper_util.py +++ b/autoarray/inversion/pixelization/mappers/mapper_util.py @@ -105,6 +105,147 @@ def data_slim_to_pixelization_unique_from( return data_to_pix_unique, data_weights, pix_lengths +import jax +import jax.numpy as jnp + +from functools import partial + + +def forward_interp(xp, yp, x): + return jax.vmap(jnp.interp, in_axes=(1, 1, None, None, None))(x, xp, yp, 0, 1).T + + +def reverse_interp(xp, yp, x): + return jax.vmap(jnp.interp, in_axes=(1, None, 1))(x, xp, yp).T + + +def create_transforms(traced_points): + # make functions that takes a set of traced points + # stored in a (N, 2) array and return functions that + # take in (N, 2) arrays and transform the values into + # the range (0, 1) and the inverse transform + N = traced_points.shape[0] + t = jnp.arange(1, N + 1) / (N + 1) + + + sort_points = jnp.sort(traced_points, axis=0) + + transform = partial(forward_interp, sort_points, t) + inv_transform = partial(reverse_interp, t, sort_points) + return transform, inv_transform + + +def adaptive_rectangular_mappings_weights_via_interpolation_from( + source_grid_size: int, + source_plane_data_grid, + source_plane_data_grid_over_sampled, +): + """ + Compute bilinear interpolation indices and weights for mapping an oversampled + source-plane grid onto a regular rectangular pixelization. + + This function takes a set of irregularly-sampled source-plane coordinates and + builds an adaptive mapping onto a `source_grid_size x source_grid_size` rectangular + pixelization using bilinear interpolation. The interpolation is expressed as: + + f(x, y) ≈ w_bl * f(ix_down, iy_down) + + w_br * f(ix_up, iy_down) + + w_tl * f(ix_down, iy_up) + + w_tr * f(ix_up, iy_up) + + where `(ix_down, ix_up, iy_down, iy_up)` are the integer grid coordinates + surrounding the continuous position `(x, y)`. + + Steps performed: + 1. Normalize the source-plane grid by subtracting its mean and dividing by + the minimum axis standard deviation (to balance scaling). + 2. Construct forward/inverse transforms which map the grid into the unit square [0,1]^2. + 3. Transform the oversampled source-plane grid into [0,1]^2, then scale it + to index space `[0, source_grid_size)`. + 4. Compute floor/ceil along x and y axes to find the enclosing rectangular cell. + 5. Build the four corner indices: bottom-left (bl), bottom-right (br), + top-left (tl), and top-right (tr). + 6. Flatten the 2D indices into 1D indices suitable for scatter operations, + with a flipped row-major convention: row = source_grid_size - i, col = j. + 7. Compute bilinear interpolation weights (`w_bl, w_br, w_tl, w_tr`). + 8. Return arrays of flattened indices and weights of shape `(N, 4)`, where + `N` is the number of oversampled coordinates. + + Parameters + ---------- + source_grid_size : int + The number of pixels along one dimension of the rectangular pixelization. + The grid is square: (source_grid_size x source_grid_size). + source_plane_data_grid : (M, 2) ndarray + The base source-plane coordinates, used to define normalization and transforms. + source_plane_data_grid_over_sampled : (N, 2) ndarray + Oversampled source-plane coordinates to be interpolated onto the rectangular grid. + + Returns + ------- + flat_indices : (N, 4) int ndarray + The flattened indices of the four neighboring pixel corners for each oversampled point. + Order: [bl, br, tl, tr]. + weights : (N, 4) float ndarray + The bilinear interpolation weights for each of the four neighboring pixels. + Order: [w_bl, w_br, w_tl, w_tr]. + """ + + # --- Step 1. Normalize grid --- + mu = source_plane_data_grid.mean(axis=0) + scale = source_plane_data_grid.std(axis=0).min() + grid_scaled = (source_plane_data_grid - mu) / scale + grid_over_scaled = (source_plane_data_grid_over_sampled - mu) / scale + + # --- Step 2. Build transforms --- + transform, inv_transform = create_transforms(grid_scaled) + + # --- Step 3. Transform oversampled grid into index space --- + grid_over_transformed = transform(grid_over_scaled) + grid_over_index = source_grid_size * grid_over_transformed + + # --- Step 4. Floor/ceil indices --- + ix_down = jnp.floor(grid_over_index[:, 0]) + ix_up = jnp.ceil(grid_over_index[:, 0]) + iy_down = jnp.floor(grid_over_index[:, 1]) + iy_up = jnp.ceil(grid_over_index[:, 1]) + + # --- Step 5. Four corners --- + idx_bl = jnp.stack([ix_down, iy_down], axis=1) # bottom-left + idx_br = jnp.stack([ix_up, iy_down], axis=1) # bottom-right + idx_tl = jnp.stack([ix_down, iy_up], axis=1) # top-left + idx_tr = jnp.stack([ix_up, iy_up], axis=1) # top-right + + # --- Step 6. Flatten indices --- + def flatten(idx, n): + row = n - idx[:, 0] + col = idx[:, 1] + return row * n + col + + flat_bl = flatten(idx_bl, source_grid_size) + flat_br = flatten(idx_br, source_grid_size) + flat_tl = flatten(idx_tl, source_grid_size) + flat_tr = flatten(idx_tr, source_grid_size) + + # --- Step 7. Bilinear interpolation weights --- + idx_down = jnp.stack([ix_down, iy_down], axis=1) + idx_up = jnp.stack([ix_up, iy_up], axis=1) + + delta_up = idx_up - grid_over_index + delta_down = grid_over_index - idx_down + + w_tl = delta_up[:, 0] * delta_up[:, 1] + w_bl = delta_up[:, 0] * delta_down[:, 1] + w_tr = delta_down[:, 0] * delta_up[:, 1] + w_br = delta_down[:, 0] * delta_down[:, 1] + + # --- Step 8. Stack outputs --- + flat_indices = jnp.stack([flat_bl, flat_br, flat_tl, flat_tr], axis=1).astype("int64") + weights = jnp.stack([w_bl, w_br, w_tl, w_tr], axis=1) + + return flat_indices, weights + + def rectangular_mappings_weights_via_interpolation_from( shape_native: Tuple[int, int], source_plane_data_grid: jnp.ndarray, diff --git a/autoarray/inversion/pixelization/mappers/rectangular.py b/autoarray/inversion/pixelization/mappers/rectangular.py index 79d2ad774..9dc0d6688 100644 --- a/autoarray/inversion/pixelization/mappers/rectangular.py +++ b/autoarray/inversion/pixelization/mappers/rectangular.py @@ -97,30 +97,20 @@ 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. """ - # from autoarray.geometry import geometry_util - # - # 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)) - # - # 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, 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=self.source_plane_data_grid.over_sampled, + # ) # ) 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=self.source_plane_data_grid.over_sampled, + mapper_util.adaptive_rectangular_mappings_weights_via_interpolation_from( + source_grid_size=self.shape_native[0], + source_plane_data_grid=self.source_plane_mesh_grid.array, + source_plane_data_grid_over_sampled=jnp.array(self.source_plane_data_grid.over_sampled), ) ) diff --git a/test_autoarray/inversion/inversion/test_factory.py b/test_autoarray/inversion/inversion/test_factory.py index ed3e6fa53..5905f570f 100644 --- a/test_autoarray/inversion/inversion/test_factory.py +++ b/test_autoarray/inversion/inversion/test_factory.py @@ -72,7 +72,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( - 7.257175708246, 1.0e-4 + 4.609440907938719, 1.0e-4 ) assert inversion.mapped_reconstructed_image == pytest.approx(np.ones(9), 1.0e-4) From e3332b68f9a2d42476e624de8a482929fcea2db6 Mon Sep 17 00:00:00 2001 From: Jammy2211 Date: Mon, 15 Sep 2025 17:24:15 +0100 Subject: [PATCH 02/17] edges unit test passes --- .../inversion/pixelization/mesh/mesh_util.py | 41 +++++++++++++++++ autoarray/structures/mesh/abstract_2d.py | 9 ++++ autoarray/structures/mesh/rectangular_2d.py | 21 ++++++--- autoarray/structures/mesh/triangulation_2d.py | 8 ---- .../structures/mesh/test_rectangular.py | 44 +++++++++++++++++++ 5 files changed, 108 insertions(+), 15 deletions(-) diff --git a/autoarray/inversion/pixelization/mesh/mesh_util.py b/autoarray/inversion/pixelization/mesh/mesh_util.py index 420d9c2ab..d8e78f18f 100644 --- a/autoarray/inversion/pixelization/mesh/mesh_util.py +++ b/autoarray/inversion/pixelization/mesh/mesh_util.py @@ -1,3 +1,4 @@ +import jax.numpy as jnp import numpy as np from typing import List, Tuple, Union @@ -300,6 +301,46 @@ def rectangular_central_neighbors( return neighbors, neighbors_sizes +def rectangular_edges_from(shape_native, pixel_scales): + """ + Returns all pixel edges for a rectangular grid as a JAX array of shape (N, 4, 2, 2), + where N = Ny * Nx. Edge order per pixel matches the user's convention: + + 0: (x1, y0) -> (x1, y1) + 1: (x1, y1) -> (x0, y1) + 2: (x0, y1) -> (x0, y0) + 3: (x0, y0) -> (x1, y0) + + Notes + ----- + - x is flipped so that the leftmost column has the largest +x (e.g. centres start at x=+1.0). + - y increases upward (top row has the most negative y when dy>0). + """ + Ny, Nx = shape_native + dy, dx = pixel_scales + + # Grid edge coordinates. Flip x so leftmost column has largest +x, matching your convention. + x_edges = ((jnp.arange(Nx + 1) - Nx / 2) * dx)[::-1] + y_edges = (jnp.arange(Ny + 1) - Ny / 2) * dy + + edges_list = [] + + # Pixel order: row-major (y outer, x inner). If you want column-major, swap the loop nesting. + for j in range(Ny): + for i in range(Nx): + y0, y1 = y_edges[i], y_edges[i + 1] + xa, xb = x_edges[j], x_edges[j + 1] # xa is the "right" boundary in your convention + + # Edge order to match your pytest: [(xa,y0)->(xa,y1), (xa,y1)->(xb,y1), (xb,y1)->(xb,y0), (xb,y0)->(xa,y0)] + e0 = jnp.array([[xa, y0], [xa, y1]]) # "top" in your test (vertical at x=xa) + e1 = jnp.array([[xa, y1], [xb, y1]]) # "right" in your test (horizontal at y=y1) + e2 = jnp.array([[xb, y1], [xb, y0]]) # "bottom" in your test (vertical at x=xb) + e3 = jnp.array([[xb, y0], [xa, y0]]) # "left" in your test (horizontal at y=y0) + + edges_list.append(jnp.stack([e0, e1, e2, e3], axis=0)) + + return jnp.stack(edges_list, axis=0) + def rectangular_edge_pixel_list_from(shape_native: Tuple[int, int]) -> List[int]: """ diff --git a/autoarray/structures/mesh/abstract_2d.py b/autoarray/structures/mesh/abstract_2d.py index 55fb7b9b8..cf630443e 100644 --- a/autoarray/structures/mesh/abstract_2d.py +++ b/autoarray/structures/mesh/abstract_2d.py @@ -5,6 +5,15 @@ class Abstract2DMesh(Structure): + + @property + def slim(self) -> "Structure": + raise NotImplementedError() + + @property + def native(self) -> Structure: + raise NotImplementedError() + @property def parameters(self) -> int: return self.pixels diff --git a/autoarray/structures/mesh/rectangular_2d.py b/autoarray/structures/mesh/rectangular_2d.py index e8c7a8a82..845ae5739 100644 --- a/autoarray/structures/mesh/rectangular_2d.py +++ b/autoarray/structures/mesh/rectangular_2d.py @@ -15,13 +15,6 @@ class Mesh2DRectangular(Abstract2DMesh): - @property - def slim(self) -> "Structure": - raise NotImplementedError() - - @property - def native(self) -> Structure: - raise NotImplementedError() def __init__( self, @@ -134,6 +127,20 @@ def neighbors(self) -> Neighbors: return Neighbors(arr=neighbors.astype("int"), sizes=sizes.astype("int")) + @cached_property + def edges(self) -> Neighbors: + """ + A class packing the ndarrays describing the neighbors of every pixel in the rectangular pixelization (see + `Neighbors` for a complete description of the neighboring scheme). + + The neighbors of a rectangular pixelization are computed by exploiting the uniform and symmetric nature of the + rectangular grid, as described in the method `mesh_util.rectangular_neighbors_from`. + """ + return mesh_util.rectangular_edges_from( + shape_native=self.shape_native, + pixel_scales=self.pixel_scales, + ) + @cached_property def edge_pixel_list(self) -> List: return mesh_util.rectangular_edge_pixel_list_from( diff --git a/autoarray/structures/mesh/triangulation_2d.py b/autoarray/structures/mesh/triangulation_2d.py index 3ddb95c01..3eeaa249d 100644 --- a/autoarray/structures/mesh/triangulation_2d.py +++ b/autoarray/structures/mesh/triangulation_2d.py @@ -14,14 +14,6 @@ class Abstract2DMeshTriangulation(Abstract2DMesh): - @property - def slim(self) -> "Structure": - raise NotImplementedError() - - @property - def native(self) -> Structure: - raise NotImplementedError() - def __init__( self, values: Union[np.ndarray, List], diff --git a/test_autoarray/structures/mesh/test_rectangular.py b/test_autoarray/structures/mesh/test_rectangular.py index a63489733..fd5f6fe1e 100644 --- a/test_autoarray/structures/mesh/test_rectangular.py +++ b/test_autoarray/structures/mesh/test_rectangular.py @@ -24,6 +24,50 @@ def test__neighbors__compare_to_mesh_util(): assert (mesh.neighbors.sizes == neighbors_sizes_util).all() +def test__edges(): + + grid = np.array( + [ + [-1.5, -1.5], + [-1.5, 0.0], + [-1.5, 1.5], + [0.0, -1.5], + [0.0, 0.0], + [0.0, 1.5], + [1.5, -1.5], + [1.5, 0.0], + [1.5, 1.5], + ] + ) + + mesh = aa.Mesh2DRectangular.overlay_grid( + shape_native=(3, 3), grid=grid, buffer=1e-8 + ) + + print(mesh[1]) + print(mesh.edges) + + assert mesh.edges[0] == pytest.approx( + np.array([[[1.5, -1.5], [1.5, -0.5]], # top + [[1.5, -0.5], [0.5, -0.5]], # right + [[0.5, -0.5], [0.5, -1.5]], # bottom + [[0.5, -1.5], [1.5, -1.5]]]), # left + ) + assert mesh.edges[1] == pytest.approx( + np.array( [[[1.5, -0.5], [1.5, 0.5]], # top + [[1.5, 0.5], [0.5, 0.5]], # right + [[0.5, 0.5], [0.5, -0.5]], # bottom + [[0.5, -0.5], [1.5, -0.5]]]), # left + ) + assert mesh.edges[5] == pytest.approx( + np.array( [[[0.5, 0.5], [0.5, 1.5]], # top + [[0.5, 1.5], [-0.5, 1.5]], # right + [[-0.5, 1.5], [-0.5, 0.5]], # bottom + [[-0.5, 0.5], [0.5, 0.5]]]), # left + ) + + + def test__edge_pixel_list(): grid = np.array( [ From 972690bdd7a6f717be0c6aefe75d938e3e26496e Mon Sep 17 00:00:00 2001 From: Jammy2211 Date: Tue, 16 Sep 2025 14:22:54 +0100 Subject: [PATCH 03/17] rectangular edges --- .../pixelization/mappers/mapper_util.py | 11 ++--- .../pixelization/mappers/rectangular.py | 4 +- .../inversion/pixelization/mesh/mesh_util.py | 24 ++++++++--- .../structures/mesh/test_rectangular.py | 40 +++++++++++-------- 4 files changed, 51 insertions(+), 28 deletions(-) diff --git a/autoarray/inversion/pixelization/mappers/mapper_util.py b/autoarray/inversion/pixelization/mappers/mapper_util.py index 3ef01e479..559c748f7 100644 --- a/autoarray/inversion/pixelization/mappers/mapper_util.py +++ b/autoarray/inversion/pixelization/mappers/mapper_util.py @@ -127,7 +127,6 @@ def create_transforms(traced_points): N = traced_points.shape[0] t = jnp.arange(1, N + 1) / (N + 1) - sort_points = jnp.sort(traced_points, axis=0) transform = partial(forward_interp, sort_points, t) @@ -136,9 +135,9 @@ def create_transforms(traced_points): def adaptive_rectangular_mappings_weights_via_interpolation_from( - source_grid_size: int, - source_plane_data_grid, - source_plane_data_grid_over_sampled, + source_grid_size: int, + source_plane_data_grid, + source_plane_data_grid_over_sampled, ): """ Compute bilinear interpolation indices and weights for mapping an oversampled @@ -240,7 +239,9 @@ def flatten(idx, n): w_br = delta_down[:, 0] * delta_down[:, 1] # --- Step 8. Stack outputs --- - flat_indices = jnp.stack([flat_bl, flat_br, flat_tl, flat_tr], axis=1).astype("int64") + flat_indices = jnp.stack([flat_bl, flat_br, flat_tl, flat_tr], axis=1).astype( + "int64" + ) weights = jnp.stack([w_bl, w_br, w_tl, w_tr], axis=1) return flat_indices, weights diff --git a/autoarray/inversion/pixelization/mappers/rectangular.py b/autoarray/inversion/pixelization/mappers/rectangular.py index 9dc0d6688..5c6d96e19 100644 --- a/autoarray/inversion/pixelization/mappers/rectangular.py +++ b/autoarray/inversion/pixelization/mappers/rectangular.py @@ -110,7 +110,9 @@ def pix_sub_weights(self) -> PixSubWeights: mapper_util.adaptive_rectangular_mappings_weights_via_interpolation_from( source_grid_size=self.shape_native[0], source_plane_data_grid=self.source_plane_mesh_grid.array, - source_plane_data_grid_over_sampled=jnp.array(self.source_plane_data_grid.over_sampled), + source_plane_data_grid_over_sampled=jnp.array( + self.source_plane_data_grid.over_sampled + ), ) ) diff --git a/autoarray/inversion/pixelization/mesh/mesh_util.py b/autoarray/inversion/pixelization/mesh/mesh_util.py index d8e78f18f..514fb5e1d 100644 --- a/autoarray/inversion/pixelization/mesh/mesh_util.py +++ b/autoarray/inversion/pixelization/mesh/mesh_util.py @@ -301,6 +301,7 @@ def rectangular_central_neighbors( return neighbors, neighbors_sizes + def rectangular_edges_from(shape_native, pixel_scales): """ Returns all pixel edges for a rectangular grid as a JAX array of shape (N, 4, 2, 2), @@ -321,7 +322,7 @@ def rectangular_edges_from(shape_native, pixel_scales): # Grid edge coordinates. Flip x so leftmost column has largest +x, matching your convention. x_edges = ((jnp.arange(Nx + 1) - Nx / 2) * dx)[::-1] - y_edges = (jnp.arange(Ny + 1) - Ny / 2) * dy + y_edges = (jnp.arange(Ny + 1) - Ny / 2) * dy edges_list = [] @@ -329,13 +330,24 @@ def rectangular_edges_from(shape_native, pixel_scales): for j in range(Ny): for i in range(Nx): y0, y1 = y_edges[i], y_edges[i + 1] - xa, xb = x_edges[j], x_edges[j + 1] # xa is the "right" boundary in your convention + xa, xb = ( + x_edges[j], + x_edges[j + 1], + ) # xa is the "right" boundary in your convention # Edge order to match your pytest: [(xa,y0)->(xa,y1), (xa,y1)->(xb,y1), (xb,y1)->(xb,y0), (xb,y0)->(xa,y0)] - e0 = jnp.array([[xa, y0], [xa, y1]]) # "top" in your test (vertical at x=xa) - e1 = jnp.array([[xa, y1], [xb, y1]]) # "right" in your test (horizontal at y=y1) - e2 = jnp.array([[xb, y1], [xb, y0]]) # "bottom" in your test (vertical at x=xb) - e3 = jnp.array([[xb, y0], [xa, y0]]) # "left" in your test (horizontal at y=y0) + e0 = jnp.array( + [[xa, y0], [xa, y1]] + ) # "top" in your test (vertical at x=xa) + e1 = jnp.array( + [[xa, y1], [xb, y1]] + ) # "right" in your test (horizontal at y=y1) + e2 = jnp.array( + [[xb, y1], [xb, y0]] + ) # "bottom" in your test (vertical at x=xb) + e3 = jnp.array( + [[xb, y0], [xa, y0]] + ) # "left" in your test (horizontal at y=y0) edges_list.append(jnp.stack([e0, e1, e2, e3], axis=0)) diff --git a/test_autoarray/structures/mesh/test_rectangular.py b/test_autoarray/structures/mesh/test_rectangular.py index fd5f6fe1e..6cdfef609 100644 --- a/test_autoarray/structures/mesh/test_rectangular.py +++ b/test_autoarray/structures/mesh/test_rectangular.py @@ -44,30 +44,38 @@ def test__edges(): shape_native=(3, 3), grid=grid, buffer=1e-8 ) - print(mesh[1]) - print(mesh.edges) - assert mesh.edges[0] == pytest.approx( - np.array([[[1.5, -1.5], [1.5, -0.5]], # top - [[1.5, -0.5], [0.5, -0.5]], # right - [[0.5, -0.5], [0.5, -1.5]], # bottom - [[0.5, -1.5], [1.5, -1.5]]]), # left + np.array( + [ + [[1.5, -1.5], [1.5, -0.5]], # top + [[1.5, -0.5], [0.5, -0.5]], # right + [[0.5, -0.5], [0.5, -1.5]], # bottom + [[0.5, -1.5], [1.5, -1.5]], + ] + ), # left ) assert mesh.edges[1] == pytest.approx( - np.array( [[[1.5, -0.5], [1.5, 0.5]], # top - [[1.5, 0.5], [0.5, 0.5]], # right - [[0.5, 0.5], [0.5, -0.5]], # bottom - [[0.5, -0.5], [1.5, -0.5]]]), # left + np.array( + [ + [[1.5, -0.5], [1.5, 0.5]], # top + [[1.5, 0.5], [0.5, 0.5]], # right + [[0.5, 0.5], [0.5, -0.5]], # bottom + [[0.5, -0.5], [1.5, -0.5]], + ] + ), # left ) assert mesh.edges[5] == pytest.approx( - np.array( [[[0.5, 0.5], [0.5, 1.5]], # top - [[0.5, 1.5], [-0.5, 1.5]], # right - [[-0.5, 1.5], [-0.5, 0.5]], # bottom - [[-0.5, 0.5], [0.5, 0.5]]]), # left + np.array( + [ + [[0.5, 0.5], [0.5, 1.5]], # top + [[0.5, 1.5], [-0.5, 1.5]], # right + [[-0.5, 1.5], [-0.5, 0.5]], # bottom + [[-0.5, 0.5], [0.5, 0.5]], + ] + ), # left ) - def test__edge_pixel_list(): grid = np.array( [ From c939bd03bbfd770015fc6c21f81f0af78174ca7c Mon Sep 17 00:00:00 2001 From: Jammy2211 Date: Tue, 16 Sep 2025 18:35:19 +0100 Subject: [PATCH 04/17] edges_transformed --- .../pixelization/mappers/mapper_util.py | 26 ++++++++++-- autoarray/structures/mesh/rectangular_2d.py | 29 ++++++++++++-- .../structures/mesh/test_rectangular.py | 40 +++++++++++++++++++ 3 files changed, 88 insertions(+), 7 deletions(-) diff --git a/autoarray/inversion/pixelization/mappers/mapper_util.py b/autoarray/inversion/pixelization/mappers/mapper_util.py index 559c748f7..428c7b3a0 100644 --- a/autoarray/inversion/pixelization/mappers/mapper_util.py +++ b/autoarray/inversion/pixelization/mappers/mapper_util.py @@ -134,6 +134,26 @@ def create_transforms(traced_points): return transform, inv_transform + +def adaptive_rectangular_transformed_grid_from( + source_plane_data_grid, + grid +): + + mu = source_plane_data_grid.mean(axis=0) + scale = source_plane_data_grid.std(axis=0).min() + + grid_scaled = (grid - mu) / scale + + transform, inv_transform = create_transforms(grid_scaled) + + def inv_full(U): + return inv_transform(U) * scale + mu + + return inv_full(grid) + + + def adaptive_rectangular_mappings_weights_via_interpolation_from( source_grid_size: int, source_plane_data_grid, @@ -194,14 +214,14 @@ def adaptive_rectangular_mappings_weights_via_interpolation_from( mu = source_plane_data_grid.mean(axis=0) scale = source_plane_data_grid.std(axis=0).min() grid_scaled = (source_plane_data_grid - mu) / scale - grid_over_scaled = (source_plane_data_grid_over_sampled - mu) / scale # --- Step 2. Build transforms --- transform, inv_transform = create_transforms(grid_scaled) # --- Step 3. Transform oversampled grid into index space --- - grid_over_transformed = transform(grid_over_scaled) - grid_over_index = source_grid_size * grid_over_transformed + grid_over_sampled_scaled = (source_plane_data_grid_over_sampled - mu) / scale + grid_over_sampled_transformed = transform(grid_over_sampled_scaled) + grid_over_index = source_grid_size * grid_over_sampled_transformed # --- Step 4. Floor/ceil indices --- ix_down = jnp.floor(grid_over_index[:, 0]) diff --git a/autoarray/structures/mesh/rectangular_2d.py b/autoarray/structures/mesh/rectangular_2d.py index 845ae5739..df1fdbd57 100644 --- a/autoarray/structures/mesh/rectangular_2d.py +++ b/autoarray/structures/mesh/rectangular_2d.py @@ -3,16 +3,19 @@ from typing import List, Optional, Tuple +from autoconf import cached_property + from autoarray import type as ty from autoarray.inversion.linear_obj.neighbors import Neighbors -from autoarray.inversion.pixelization.mesh import mesh_util from autoarray.mask.mask_2d import Mask2D from autoarray.structures.abstract_structure import Structure from autoarray.structures.arrays.uniform_2d import Array2D -from autoarray.structures.grids import grid_2d_util + from autoarray.structures.mesh.abstract_2d import Abstract2DMesh -from autoconf import cached_property +from autoarray.inversion.pixelization.mappers import mapper_util +from autoarray.inversion.pixelization.mesh import mesh_util +from autoarray.structures.grids import grid_2d_util class Mesh2DRectangular(Abstract2DMesh): @@ -128,7 +131,7 @@ def neighbors(self) -> Neighbors: return Neighbors(arr=neighbors.astype("int"), sizes=sizes.astype("int")) @cached_property - def edges(self) -> Neighbors: + def edges(self): """ A class packing the ndarrays describing the neighbors of every pixel in the rectangular pixelization (see `Neighbors` for a complete description of the neighboring scheme). @@ -141,6 +144,24 @@ def edges(self) -> Neighbors: pixel_scales=self.pixel_scales, ) + @cached_property + def edges_transformed(self): + """ + A class packing the ndarrays describing the neighbors of every pixel in the rectangular pixelization (see + `Neighbors` for a complete description of the neighboring scheme). + + The neighbors of a rectangular pixelization are computed by exploiting the uniform and symmetric nature of the + rectangular grid, as described in the method `mesh_util.rectangular_neighbors_from`. + """ + + edges_transformed = mapper_util.adaptive_rectangular_transformed_grid_from( + source_plane_data_grid=self.array, + grid=self.edges.reshape(-1, 2), + ) + + + return edges_transformed.reshape(self.edges.shape) + @cached_property def edge_pixel_list(self) -> List: return mesh_util.rectangular_edge_pixel_list_from( diff --git a/test_autoarray/structures/mesh/test_rectangular.py b/test_autoarray/structures/mesh/test_rectangular.py index 6cdfef609..80bcce919 100644 --- a/test_autoarray/structures/mesh/test_rectangular.py +++ b/test_autoarray/structures/mesh/test_rectangular.py @@ -75,6 +75,46 @@ def test__edges(): ), # left ) +def test__edges_transformed(): + + grid = np.array( + [ + [-1.5, -1.5], + [-1.5, 0.0], + [-1.5, 1.5], + [0.0, -1.5], + [0.0, 0.0], + [0.0, 1.5], + [1.5, -1.5], + [1.5, 0.0], + [1.5, 1.5], + ] + ) + + mesh = aa.Mesh2DRectangular.overlay_grid( + shape_native=(3, 3), grid=grid, buffer=1e-8 + ) + + assert mesh.edges_transformed[0] == pytest.approx( + np.array( + [ + [[1.5, -1.5], [1.5, -0.5]], # top + [[1.5, -0.5], [0.5, -0.5]], # right + [[0.5, -0.5], [0.5, -1.5]], # bottom + [[0.5, -1.5], [1.5, -1.5]], + ] + ), # left + ) + assert mesh.edges[1] == pytest.approx( + np.array( + [ + [[1.5, -0.5], [1.5, 0.5]], # top + [[1.5, 0.5], [0.5, 0.5]], # right + [[0.5, 0.5], [0.5, -0.5]], # bottom + [[0.5, -0.5], [1.5, -0.5]], + ] + ), # left + ) def test__edge_pixel_list(): grid = np.array( From 5c352bba2501f1978322d8dbe352fa374eb4e53f Mon Sep 17 00:00:00 2001 From: Jammy2211 Date: Thu, 18 Sep 2025 08:48:10 +0100 Subject: [PATCH 05/17] remove InversionException tests --- .../pixelization/mappers/rectangular.py | 28 +++++++++---------- autoarray/structures/mesh/rectangular_2d.py | 1 - .../inversion/inversion/test_abstract.py | 11 -------- .../inversion/inversion/test_factory.py | 5 +++- .../inversion/test_inversion_util.py | 17 ----------- 5 files changed, 18 insertions(+), 44 deletions(-) diff --git a/autoarray/inversion/pixelization/mappers/rectangular.py b/autoarray/inversion/pixelization/mappers/rectangular.py index 5c6d96e19..c54f5f8e0 100644 --- a/autoarray/inversion/pixelization/mappers/rectangular.py +++ b/autoarray/inversion/pixelization/mappers/rectangular.py @@ -98,24 +98,24 @@ 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=self.source_plane_data_grid.over_sampled, - # ) - # ) - mappings, weights = ( - mapper_util.adaptive_rectangular_mappings_weights_via_interpolation_from( - source_grid_size=self.shape_native[0], - source_plane_data_grid=self.source_plane_mesh_grid.array, - source_plane_data_grid_over_sampled=jnp.array( - self.source_plane_data_grid.over_sampled - ), + 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=self.source_plane_data_grid.over_sampled, ) ) + # mappings, weights = ( + # mapper_util.adaptive_rectangular_mappings_weights_via_interpolation_from( + # source_grid_size=self.shape_native[0], + # source_plane_data_grid=self.source_plane_mesh_grid.array, + # source_plane_data_grid_over_sampled=jnp.array( + # self.source_plane_data_grid.over_sampled + # ), + # ) + # ) + return PixSubWeights( mappings=mappings, sizes=4 * jnp.ones(len(mappings), dtype="int"), diff --git a/autoarray/structures/mesh/rectangular_2d.py b/autoarray/structures/mesh/rectangular_2d.py index df1fdbd57..744e06b16 100644 --- a/autoarray/structures/mesh/rectangular_2d.py +++ b/autoarray/structures/mesh/rectangular_2d.py @@ -159,7 +159,6 @@ def edges_transformed(self): grid=self.edges.reshape(-1, 2), ) - return edges_transformed.reshape(self.edges.shape) @cached_property diff --git a/test_autoarray/inversion/inversion/test_abstract.py b/test_autoarray/inversion/inversion/test_abstract.py index 8880b544c..84338c4ae 100644 --- a/test_autoarray/inversion/inversion/test_abstract.py +++ b/test_autoarray/inversion/inversion/test_abstract.py @@ -426,17 +426,6 @@ def test__data_subtracted_dict(): assert (inversion.data_subtracted_dict[linear_obj_1] == 2.0 * np.ones(3)).all() -def test__reconstruction_raises_exception_for_linalg_error(): - # noinspection PyTypeChecker - inversion = aa.m.MockInversion( - data_vector=np.ones(3), curvature_reg_matrix=np.ones((3, 3)) - ) - - with pytest.raises(exc.InversionException): - # noinspection PyStatementEffect - inversion.reconstruction - - def test__regularization_term(): reconstruction = np.array([1.0, 1.0, 1.0]) diff --git a/test_autoarray/inversion/inversion/test_factory.py b/test_autoarray/inversion/inversion/test_factory.py index 5905f570f..602001429 100644 --- a/test_autoarray/inversion/inversion/test_factory.py +++ b/test_autoarray/inversion/inversion/test_factory.py @@ -72,8 +72,11 @@ 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( - 4.609440907938719, 1.0e-4 + 7.2571757082, 1.0e-4 ) + # assert inversion.log_det_curvature_reg_matrix_term == pytest.approx( + # 4.609440907938719, 1.0e-4 + # ) assert inversion.mapped_reconstructed_image == pytest.approx(np.ones(9), 1.0e-4) inversion = aa.Inversion( diff --git a/test_autoarray/inversion/inversion/test_inversion_util.py b/test_autoarray/inversion/inversion/test_inversion_util.py index 86b722812..e2eb96a3d 100644 --- a/test_autoarray/inversion/inversion/test_inversion_util.py +++ b/test_autoarray/inversion/inversion/test_inversion_util.py @@ -86,23 +86,6 @@ def test__reconstruction_positive_negative_from(): assert reconstruction == pytest.approx(np.array([1.0, -1.0, 3.0]), 1.0e-4) - -def test__reconstruction_positive_negative_from__check_solution_raises_error_cause_all_values_identical(): - data_vector = np.array([1.0, 1.0, 1.0]) - - curvature_reg_matrix = np.array([[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]]) - - # reconstruction = np.array([1.0, 1.0, 1.0]) - - with pytest.raises(aa.exc.InversionException): - aa.util.inversion.reconstruction_positive_negative_from( - data_vector=data_vector, - curvature_reg_matrix=curvature_reg_matrix, - mapper_param_range_list=[[0, 3]], - force_check_reconstruction=True, - ) - - def test__mapped_reconstructed_data_via_mapping_matrix_from(): mapping_matrix = np.array([[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]]) From c1fff0e17d6c12c7da50866548b330c6e5e4749d Mon Sep 17 00:00:00 2001 From: Jammy2211 Date: Thu, 18 Sep 2025 08:51:15 +0100 Subject: [PATCH 06/17] comment out edges test --- .../structures/mesh/test_rectangular.py | 80 +++++++++---------- 1 file changed, 40 insertions(+), 40 deletions(-) diff --git a/test_autoarray/structures/mesh/test_rectangular.py b/test_autoarray/structures/mesh/test_rectangular.py index 80bcce919..8d22c5dae 100644 --- a/test_autoarray/structures/mesh/test_rectangular.py +++ b/test_autoarray/structures/mesh/test_rectangular.py @@ -75,46 +75,46 @@ def test__edges(): ), # left ) -def test__edges_transformed(): - - grid = np.array( - [ - [-1.5, -1.5], - [-1.5, 0.0], - [-1.5, 1.5], - [0.0, -1.5], - [0.0, 0.0], - [0.0, 1.5], - [1.5, -1.5], - [1.5, 0.0], - [1.5, 1.5], - ] - ) - - mesh = aa.Mesh2DRectangular.overlay_grid( - shape_native=(3, 3), grid=grid, buffer=1e-8 - ) - - assert mesh.edges_transformed[0] == pytest.approx( - np.array( - [ - [[1.5, -1.5], [1.5, -0.5]], # top - [[1.5, -0.5], [0.5, -0.5]], # right - [[0.5, -0.5], [0.5, -1.5]], # bottom - [[0.5, -1.5], [1.5, -1.5]], - ] - ), # left - ) - assert mesh.edges[1] == pytest.approx( - np.array( - [ - [[1.5, -0.5], [1.5, 0.5]], # top - [[1.5, 0.5], [0.5, 0.5]], # right - [[0.5, 0.5], [0.5, -0.5]], # bottom - [[0.5, -0.5], [1.5, -0.5]], - ] - ), # left - ) +# def test__edges_transformed(): +# +# grid = np.array( +# [ +# [-1.5, -1.5], +# [-1.5, 0.0], +# [-1.5, 1.5], +# [0.0, -1.5], +# [0.0, 0.0], +# [0.0, 1.5], +# [1.5, -1.5], +# [1.5, 0.0], +# [1.5, 1.5], +# ] +# ) +# +# mesh = aa.Mesh2DRectangular.overlay_grid( +# shape_native=(3, 3), grid=grid, buffer=1e-8 +# ) +# +# assert mesh.edges_transformed[0] == pytest.approx( +# np.array( +# [ +# [[1.5, -1.5], [1.5, -0.5]], # top +# [[1.5, -0.5], [0.5, -0.5]], # right +# [[0.5, -0.5], [0.5, -1.5]], # bottom +# [[0.5, -1.5], [1.5, -1.5]], +# ] +# ), # left +# ) +# assert mesh.edges[1] == pytest.approx( +# np.array( +# [ +# [[1.5, -0.5], [1.5, 0.5]], # top +# [[1.5, 0.5], [0.5, 0.5]], # right +# [[0.5, 0.5], [0.5, -0.5]], # bottom +# [[0.5, -0.5], [1.5, -0.5]], +# ] +# ), # left +# ) def test__edge_pixel_list(): grid = np.array( From dde10c4d3ef86fc158f0617fa21ac1f395b38329 Mon Sep 17 00:00:00 2001 From: Jammy2211 Date: Thu, 18 Sep 2025 09:56:55 +0100 Subject: [PATCH 07/17] data_vector_via_w_tilde_data_imaging_from converted to JAX --- .../imaging/inversion_imaging_util.py | 46 ++++++++++--------- .../inversion/inversion/imaging/w_tilde.py | 14 +++--- 2 files changed, 30 insertions(+), 30 deletions(-) diff --git a/autoarray/inversion/inversion/imaging/inversion_imaging_util.py b/autoarray/inversion/inversion/imaging/inversion_imaging_util.py index 699361638..6b8bd756c 100644 --- a/autoarray/inversion/inversion/imaging/inversion_imaging_util.py +++ b/autoarray/inversion/inversion/imaging/inversion_imaging_util.py @@ -1,10 +1,11 @@ +from scipy.signal import convolve2d +import jax.numpy as jnp import numpy as np from typing import Tuple from autoarray import numba_util +from scipy.signal import correlate2d - -@numba_util.jit() def w_tilde_data_imaging_from( image_native: np.ndarray, noise_map_native: np.ndarray, @@ -44,32 +45,33 @@ def w_tilde_data_imaging_from( efficient calculation of the data vector. """ - kernel_shift_y = -(kernel_native.shape[1] // 2) - kernel_shift_x = -(kernel_native.shape[0] // 2) - - image_pixels = len(native_index_for_slim_index) - - w_tilde_data = np.zeros((image_pixels,)) - - weight_map_native = image_native / noise_map_native**2.0 + # 1) weight map = image / noise^2 (safe where noise==0) + weight_map = jnp.where(noise_map_native > 0.0, + image_native / (noise_map_native ** 2), + 0.0) - for ip0 in range(image_pixels): - ip0_y, ip0_x = native_index_for_slim_index[ip0] + Ky, Kx = kernel_native.shape + ph, pw = Ky // 2, Kx // 2 - value = 0.0 + # 2) pad so neighbourhood gathers never go OOB + padded = jnp.pad(weight_map, ((ph, ph), (pw, pw)), mode="constant", constant_values=0.0) - for k0_y in range(kernel_native.shape[0]): - for k0_x in range(kernel_native.shape[1]): - weight_value = weight_map_native[ - ip0_y + k0_y + kernel_shift_y, ip0_x + k0_x + kernel_shift_x - ] + # 3) build broadcasted neighbourhood indices for all requested pixels + # shift pixel coords into the padded frame + ys = native_index_for_slim_index[:, 0] + ph # (N,) + xs = native_index_for_slim_index[:, 1] + pw # (N,) - if not np.isnan(weight_value): - value += kernel_native[k0_y, k0_x] * weight_value + # kernel-relative offsets + dy = jnp.arange(Ky) - ph # (Ky,) + dx = jnp.arange(Kx) - pw # (Kx,) - w_tilde_data[ip0] = value + # broadcast to (N, Ky, Kx) + Y = ys[:, None, None] + dy[None, :, None] + X = xs[:, None, None] + dx[None, None, :] - return w_tilde_data + # 4) gather patches and correlate (no kernel flip) + patches = padded[Y, X] # (N, Ky, Kx) + return jnp.sum(patches * kernel_native[None, :, :], axis=(1, 2)) # (N,) @numba_util.jit() diff --git a/autoarray/inversion/inversion/imaging/w_tilde.py b/autoarray/inversion/inversion/imaging/w_tilde.py index b725adb6c..e4a18d00b 100644 --- a/autoarray/inversion/inversion/imaging/w_tilde.py +++ b/autoarray/inversion/inversion/imaging/w_tilde.py @@ -75,12 +75,10 @@ def __init__( @cached_property def w_tilde_data(self): return inversion_imaging_util.w_tilde_data_imaging_from( - image_native=np.array(self.data.native.array).astype("float"), - noise_map_native=np.array(self.noise_map.native.array).astype("float"), - kernel_native=np.array(self.psf.native.array).astype("float"), - native_index_for_slim_index=np.array( - self.data.mask.derive_indexes.native_for_slim - ).astype("int"), + image_native=self.data.native.array, + noise_map_native=self.noise_map.native.array, + kernel_native=self.psf.native.array, + native_index_for_slim_index=self.data.mask.derive_indexes.native_for_slim ) @property @@ -104,7 +102,7 @@ def _data_vector_mapper(self) -> np.ndarray: for mapper_index, mapper in enumerate(mapper_list): data_vector_mapper = ( inversion_imaging_util.data_vector_via_w_tilde_data_imaging_from( - w_tilde_data=self.w_tilde_data, + w_tilde_data=np.array(self.w_tilde_data), data_to_pix_unique=np.array( mapper.unique_mappings.data_to_pix_unique ), @@ -152,7 +150,7 @@ def _data_vector_x1_mapper(self) -> np.ndarray: linear_obj = self.linear_obj_list[0] return inversion_imaging_util.data_vector_via_w_tilde_data_imaging_from( - w_tilde_data=self.w_tilde_data, + w_tilde_data=np.array(self.w_tilde_data), data_to_pix_unique=linear_obj.unique_mappings.data_to_pix_unique, data_weights=linear_obj.unique_mappings.data_weights, pix_lengths=linear_obj.unique_mappings.pix_lengths, From 6fd813bfcec589bb46e373f7f67684478fab619f Mon Sep 17 00:00:00 2001 From: Jammy2211 Date: Thu, 18 Sep 2025 10:15:15 +0100 Subject: [PATCH 08/17] _data_vector_mapper in JAX --- .../inversion/inversion/imaging/w_tilde.py | 57 ++++++++++--------- .../imaging/test_inversion_imaging_util.py | 2 +- 2 files changed, 31 insertions(+), 28 deletions(-) diff --git a/autoarray/inversion/inversion/imaging/w_tilde.py b/autoarray/inversion/inversion/imaging/w_tilde.py index e4a18d00b..d95ed07ae 100644 --- a/autoarray/inversion/inversion/imaging/w_tilde.py +++ b/autoarray/inversion/inversion/imaging/w_tilde.py @@ -1,3 +1,4 @@ +import jax.numpy as jnp import numpy as np from typing import Dict, List, Optional, Union @@ -91,31 +92,33 @@ def _data_vector_mapper(self) -> np.ndarray: in the inversion, and is separated into a separate method to enable preloading of the mapper `data_vector`. """ - if not self.has(cls=AbstractMapper): - return None - - data_vector = np.zeros(self.total_params) - - mapper_list = self.cls_list_from(cls=AbstractMapper) - mapper_param_range = self.param_range_list_from(cls=AbstractMapper) - - for mapper_index, mapper in enumerate(mapper_list): - data_vector_mapper = ( - inversion_imaging_util.data_vector_via_w_tilde_data_imaging_from( - w_tilde_data=np.array(self.w_tilde_data), - data_to_pix_unique=np.array( - mapper.unique_mappings.data_to_pix_unique - ), - data_weights=np.array(mapper.unique_mappings.data_weights), - pix_lengths=np.array(mapper.unique_mappings.pix_lengths), - pix_pixels=mapper.params, - ) - ) - param_range = mapper_param_range[mapper_index] - - data_vector[param_range[0] : param_range[1],] = data_vector_mapper - - return data_vector + return jnp.dot(self.mapping_matrix.T, self.w_tilde_data) + + # if not self.has(cls=AbstractMapper): + # return None + # + # data_vector = np.zeros(self.total_params) + # + # mapper_list = self.cls_list_from(cls=AbstractMapper) + # mapper_param_range = self.param_range_list_from(cls=AbstractMapper) + # + # for mapper_index, mapper in enumerate(mapper_list): + # data_vector_mapper = ( + # inversion_imaging_util.data_vector_via_w_tilde_data_imaging_from( + # w_tilde_data=np.array(self.w_tilde_data), + # data_to_pix_unique=np.array( + # mapper.unique_mappings.data_to_pix_unique + # ), + # data_weights=np.array(mapper.unique_mappings.data_weights), + # pix_lengths=np.array(mapper.unique_mappings.pix_lengths), + # pix_pixels=mapper.params, + # ) + # ) + # param_range = mapper_param_range[mapper_index] + # + # data_vector[param_range[0] : param_range[1],] = data_vector_mapper + # + # return data_vector @cached_property def data_vector(self) -> np.ndarray: @@ -170,7 +173,7 @@ def _data_vector_multi_mapper(self) -> np.ndarray: return np.concatenate( [ inversion_imaging_util.data_vector_via_w_tilde_data_imaging_from( - w_tilde_data=self.w_tilde_data, + w_tilde_data=np.array(self.w_tilde_data), data_to_pix_unique=linear_obj.unique_mappings.data_to_pix_unique, data_weights=linear_obj.unique_mappings.data_weights, pix_lengths=linear_obj.unique_mappings.pix_lengths, @@ -193,7 +196,7 @@ def _data_vector_func_list_and_mapper(self) -> np.ndarray: separation of functions enables the `data_vector` to be preloaded in certain circumstances. """ - data_vector = self._data_vector_mapper + data_vector = np.array(self._data_vector_mapper) linear_func_param_range = self.param_range_list_from( cls=AbstractLinearObjFuncList 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 f96731aa9..e5abd0b4b 100644 --- a/test_autoarray/inversion/inversion/imaging/test_inversion_imaging_util.py +++ b/test_autoarray/inversion/inversion/imaging/test_inversion_imaging_util.py @@ -243,7 +243,7 @@ def test__data_vector_via_w_tilde_data_two_methods_agree(): data_vector_via_w_tilde = ( aa.util.inversion_imaging.data_vector_via_w_tilde_data_imaging_from( - w_tilde_data=w_tilde_data, + w_tilde_data=np.array(w_tilde_data), data_to_pix_unique=data_to_pix_unique.astype("int"), data_weights=data_weights, pix_lengths=pix_lengths.astype("int"), From 9eea7c2a930b83fa79269ffc5b7a68b74931f6f9 Mon Sep 17 00:00:00 2001 From: Jammy2211 Date: Thu, 18 Sep 2025 10:53:33 +0100 Subject: [PATCH 09/17] _curvature_matrix_x1_mapper --- autoarray/dataset/imaging/dataset.py | 3 ++ autoarray/dataset/imaging/w_tilde.py | 48 +++++++++++++++++++ .../inversion/inversion/imaging/w_tilde.py | 9 +++- .../inversion/inversion/inversion_util.py | 5 +- .../inversion/imaging/test_imaging.py | 2 +- 5 files changed, 60 insertions(+), 7 deletions(-) diff --git a/autoarray/dataset/imaging/dataset.py b/autoarray/dataset/imaging/dataset.py index b851089bb..d5d146a46 100644 --- a/autoarray/dataset/imaging/dataset.py +++ b/autoarray/dataset/imaging/dataset.py @@ -252,6 +252,9 @@ def w_tilde(self): indexes=indexes.astype("int"), lengths=lengths.astype("int"), noise_map_value=self.noise_map[0], + noise_map=self.noise_map, + psf=self.psf, + mask=self.mask, ) @classmethod diff --git a/autoarray/dataset/imaging/w_tilde.py b/autoarray/dataset/imaging/w_tilde.py index 985caeeb8..818941162 100644 --- a/autoarray/dataset/imaging/w_tilde.py +++ b/autoarray/dataset/imaging/w_tilde.py @@ -2,8 +2,12 @@ import logging import numpy as np +from autoconf import cached_property + from autoarray.dataset.abstract.w_tilde import AbstractWTilde +from autoarray.inversion.inversion.imaging import inversion_imaging_util + logger = logging.getLogger(__name__) @@ -13,6 +17,9 @@ def __init__( curvature_preload: np.ndarray, indexes: np.ndim, lengths: np.ndarray, + noise_map : np.ndarray, + psf : np.ndarray, + mask : np.ndarray, noise_map_value: float, ): """ @@ -44,3 +51,44 @@ def __init__( self.indexes = indexes self.lengths = lengths + self.noise_map = noise_map + self.psf = psf + self.mask = mask + + @cached_property + def w_matrix(self): + """ + The matrix `w_tilde_curvature` is a matrix of dimensions [image_pixels, image_pixels] that encodes the PSF + convolution of every pair of image pixels given the noise map. This can be used to efficiently compute the + curvature matrix via the mappings between image and source pixels, in a way that omits having to perform the + PSF convolution on every individual source pixel. This provides a significant speed up for inversions of imaging + datasets. + + The limitation of this matrix is that the dimensions of [image_pixels, image_pixels] can exceed many 10s of GB's, + making it impossible to store in memory and its use in linear algebra calculations extremely. The method + `w_tilde_curvature_preload_imaging_from` describes a compressed representation that overcomes this hurdles. It is + advised `w_tilde` and this method are only used for testing. + + Parameters + ---------- + noise_map_native + The two dimensional masked noise-map of values which w_tilde is computed from. + kernel_native + The two dimensional PSF kernel that w_tilde encodes the convolution of. + native_index_for_slim_index + An array of shape [total_x_pixels*sub_size] that maps pixels from the slimmed array to the native array. + + Returns + ------- + ndarray + A matrix that encodes the PSF convolution values between the noise map that enables efficient calculation of + the curvature matrix. + """ + + return inversion_imaging_util.w_tilde_curvature_imaging_from( + noise_map_native=np.array(self.noise_map.native.array).astype("float64"), + kernel_native=np.array(self.psf.native.array).astype("float64"), + native_index_for_slim_index=np.array( + self.mask.derive_indexes.native_for_slim + ).astype("int"), + ) \ No newline at end of file diff --git a/autoarray/inversion/inversion/imaging/w_tilde.py b/autoarray/inversion/inversion/imaging/w_tilde.py index d95ed07ae..6f97b9285 100644 --- a/autoarray/inversion/inversion/imaging/w_tilde.py +++ b/autoarray/inversion/inversion/imaging/w_tilde.py @@ -91,7 +91,6 @@ def _data_vector_mapper(self) -> np.ndarray: This method is used to compute part of the `data_vector` if there are also linear function list objects in the inversion, and is separated into a separate method to enable preloading of the mapper `data_vector`. """ - return jnp.dot(self.mapping_matrix.T, self.w_tilde_data) # if not self.has(cls=AbstractMapper): @@ -360,7 +359,13 @@ def _curvature_matrix_x1_mapper(self) -> np.ndarray: This method computes the `curvature_matrix` when there is a single mapper object in the `Inversion`, which circumvents `block_diag` for speed up. """ - return self._curvature_matrix_mapper_diag + + return inversion_util.curvature_matrix_via_w_tilde_from( + w_tilde=self.w_tilde.w_matrix, + mapping_matrix=self.mapping_matrix + ) + +# return self._curvature_matrix_mapper_diag @property def _curvature_matrix_multi_mapper(self) -> np.ndarray: diff --git a/autoarray/inversion/inversion/inversion_util.py b/autoarray/inversion/inversion/inversion_util.py index b2eb4f26c..be0b314dc 100644 --- a/autoarray/inversion/inversion/inversion_util.py +++ b/autoarray/inversion/inversion/inversion_util.py @@ -4,8 +4,6 @@ from typing import List, Optional, Type -from autoconf import conf - from autoarray.inversion.inversion.settings import SettingsInversion from autoarray import numba_util @@ -36,8 +34,7 @@ def curvature_matrix_via_w_tilde_from( ndarray The curvature matrix `F` (see Warren & Dye 2003). """ - - return np.dot(mapping_matrix.T, np.dot(w_tilde, mapping_matrix)) + return jnp.dot(mapping_matrix.T, jnp.dot(w_tilde, mapping_matrix)) def curvature_matrix_with_added_to_diag_from( diff --git a/test_autoarray/inversion/inversion/imaging/test_imaging.py b/test_autoarray/inversion/inversion/imaging/test_imaging.py index bd54c35f1..c77bd0862 100644 --- a/test_autoarray/inversion/inversion/imaging/test_imaging.py +++ b/test_autoarray/inversion/inversion/imaging/test_imaging.py @@ -144,7 +144,7 @@ def test__w_tilde_checks_noise_map_and_raises_exception_if_preloads_dont_match_n grid = aa.Grid2D.from_mask(mask=mask) w_tilde = WTildeImaging( - curvature_preload=None, indexes=None, lengths=None, noise_map_value=2.0 + curvature_preload=None, indexes=None, lengths=None, noise_map_value=2.0, noise_map=None, psf=None, mask=mask ) with pytest.raises(exc.InversionException): From 74c8d3903190d5c437b8a4ce9a9e20d0b37892eb Mon Sep 17 00:00:00 2001 From: Jammy2211 Date: Thu, 18 Sep 2025 11:42:35 +0100 Subject: [PATCH 10/17] mapped_reconstructed_data_via_w_tilde_from --- autoarray/dataset/imaging/w_tilde.py | 12 ++++ .../imaging/inversion_imaging_util.py | 57 +++++++++++++++++++ .../inversion/inversion/imaging/w_tilde.py | 16 +++++- .../inversion/inversion/inversion_util.py | 28 +++++++++ 4 files changed, 112 insertions(+), 1 deletion(-) diff --git a/autoarray/dataset/imaging/w_tilde.py b/autoarray/dataset/imaging/w_tilde.py index 818941162..7fc1fdf51 100644 --- a/autoarray/dataset/imaging/w_tilde.py +++ b/autoarray/dataset/imaging/w_tilde.py @@ -91,4 +91,16 @@ def w_matrix(self): native_index_for_slim_index=np.array( self.mask.derive_indexes.native_for_slim ).astype("int"), + ) + + @cached_property + def psf_operator_matrix_dense(self): + + return inversion_imaging_util.psf_operator_matrix_dense_from( + kernel_native=np.array(self.psf.native.array).astype("float64"), + native_index_for_slim_index=np.array( + self.mask.derive_indexes.native_for_slim + ).astype("int"), + native_shape=self.noise_map.shape_native, + correlate=False ) \ No newline at end of file diff --git a/autoarray/inversion/inversion/imaging/inversion_imaging_util.py b/autoarray/inversion/inversion/imaging/inversion_imaging_util.py index 6b8bd756c..009db6526 100644 --- a/autoarray/inversion/inversion/imaging/inversion_imaging_util.py +++ b/autoarray/inversion/inversion/imaging/inversion_imaging_util.py @@ -6,6 +6,63 @@ from autoarray import numba_util from scipy.signal import correlate2d +import numpy as np + +def psf_operator_matrix_dense_from( + kernel_native: np.ndarray, + native_index_for_slim_index: np.ndarray, # shape (N_pix, 2), native (y,x) coords of masked pixels + native_shape: tuple[int, int], + correlate: bool = True, +) -> np.ndarray: + """ + Construct a dense PSF operator W (N_pix x N_pix) that maps masked image pixels to masked image pixels. + + Parameters + ---------- + kernel_native : (Ky, Kx) PSF kernel. + native_index_for_slim_index : (N_pix, 2) array of int + Native (y, x) coords for each masked pixel. + native_shape : (Ny, Nx) + Native 2D image shape. + correlate : bool, default True + If True, use correlation convention (no kernel flip). + If False, use convolution convention (flip kernel). + + Returns + ------- + W : ndarray, shape (N_pix, N_pix) + Dense PSF operator. + """ + Ky, Kx = kernel_native.shape + ph, pw = Ky // 2, Kx // 2 + Ny, Nx = native_shape + N_pix = native_index_for_slim_index.shape[0] + + ker = kernel_native if correlate else kernel_native[::-1, ::-1] + + # Padded index grid: -1 everywhere, slim index where masked + index_padded = -np.ones((Ny + 2*ph, Nx + 2*pw), dtype=np.int64) + for p, (y, x) in enumerate(native_index_for_slim_index): + index_padded[y + ph, x + pw] = p + + # Neighborhood offsets + dy = np.arange(Ky) - ph + dx = np.arange(Kx) - pw + + W = np.zeros((N_pix, N_pix), dtype=float) + + for i, (y, x) in enumerate(native_index_for_slim_index): + yp = y + ph + xp = x + pw + for j, dy_ in enumerate(dy): + for k, dx_ in enumerate(dx): + neigh = index_padded[yp + dy_, xp + dx_] + if neigh >= 0: + W[i, neigh] += ker[j, k] + + return W + + def w_tilde_data_imaging_from( image_native: np.ndarray, noise_map_native: np.ndarray, diff --git a/autoarray/inversion/inversion/imaging/w_tilde.py b/autoarray/inversion/inversion/imaging/w_tilde.py index 6f97b9285..e3c78cd8f 100644 --- a/autoarray/inversion/inversion/imaging/w_tilde.py +++ b/autoarray/inversion/inversion/imaging/w_tilde.py @@ -517,7 +517,8 @@ def mapped_reconstructed_data_dict(self) -> Dict[LinearObj, Array2D]: for linear_obj in self.linear_obj_list: reconstruction = reconstruction_dict[linear_obj] - if isinstance(linear_obj, AbstractMapper): + if isinstance(linear_obj, AbstractMapper) and self.has(cls=AbstractLinearObjFuncList): + mapped_reconstructed_image = inversion_util.mapped_reconstructed_data_via_image_to_pix_unique_from( data_to_pix_unique=linear_obj.unique_mappings.data_to_pix_unique, data_weights=linear_obj.unique_mappings.data_weights, @@ -533,7 +534,20 @@ def mapped_reconstructed_data_dict(self) -> Dict[LinearObj, Array2D]: values=mapped_reconstructed_image, mask=self.mask ) + elif isinstance(linear_obj, AbstractMapper) and not self.has(cls=AbstractLinearObjFuncList): + + mapped_reconstructed_image = inversion_util.mapped_reconstructed_data_via_w_tilde_from( + w_tilde=self.w_tilde.psf_operator_matrix_dense, + mapping_matrix=self.mapping_matrix, + reconstruction=reconstruction, + ) + + mapped_reconstructed_image = Array2D( + values=mapped_reconstructed_image, mask=self.mask + ) + else: + operated_mapping_matrix = self.linear_func_operated_mapping_matrix_dict[ linear_obj ] diff --git a/autoarray/inversion/inversion/inversion_util.py b/autoarray/inversion/inversion/inversion_util.py index be0b314dc..6f9d7d580 100644 --- a/autoarray/inversion/inversion/inversion_util.py +++ b/autoarray/inversion/inversion/inversion_util.py @@ -187,6 +187,34 @@ def mapped_reconstructed_data_via_mapping_matrix_from( return jnp.dot(mapping_matrix, reconstruction) +def mapped_reconstructed_data_via_w_tilde_from( + w_tilde: np.ndarray, mapping_matrix: np.ndarray, reconstruction: np.ndarray +) -> np.ndarray: + """ + Returns the reconstructed data vector from the unblurred mapping matrix `M`, + the reconstruction vector `s`, and the PSF convolution operator `w_tilde`. + + Equivalent to: + reconstructed = (W @ M) @ s + = W @ (M @ s) + + Parameters + ---------- + w_tilde + Array of shape [image_pixels, image_pixels], the PSF convolution operator. + mapping_matrix + Array of shape [image_pixels, source_pixels], unblurred mapping matrix. + reconstruction + Array of shape [source_pixels], solution vector. + + Returns + ------- + ndarray + The reconstructed data vector of shape [image_pixels]. + """ + return w_tilde @ (mapping_matrix @ reconstruction) + + def reconstruction_positive_negative_from( data_vector: np.ndarray, curvature_reg_matrix: np.ndarray, From c38ed0744159af6d908be7f973b88c731555d71c Mon Sep 17 00:00:00 2001 From: Jammy2211 Date: Sun, 21 Sep 2025 13:38:49 +0100 Subject: [PATCH 11/17] lots of hacky stuff --- .../inversion/inversion/imaging/w_tilde.py | 19 +-- autoarray/structures/arrays/kernel_2d.py | 112 ++++++++++++++++++ 2 files changed, 122 insertions(+), 9 deletions(-) diff --git a/autoarray/inversion/inversion/imaging/w_tilde.py b/autoarray/inversion/inversion/imaging/w_tilde.py index e3c78cd8f..a61b311d1 100644 --- a/autoarray/inversion/inversion/imaging/w_tilde.py +++ b/autoarray/inversion/inversion/imaging/w_tilde.py @@ -148,16 +148,17 @@ def _data_vector_x1_mapper(self) -> np.ndarray: This method computes the `data_vector` whenthere is a single mapper object in the `Inversion`, which circumvents `np.concatenate` for speed up. """ + return self._data_vector_mapper - linear_obj = self.linear_obj_list[0] - - return inversion_imaging_util.data_vector_via_w_tilde_data_imaging_from( - w_tilde_data=np.array(self.w_tilde_data), - data_to_pix_unique=linear_obj.unique_mappings.data_to_pix_unique, - data_weights=linear_obj.unique_mappings.data_weights, - pix_lengths=linear_obj.unique_mappings.pix_lengths, - pix_pixels=linear_obj.params, - ) + # linear_obj = self.linear_obj_list[0] + # + # return inversion_imaging_util.data_vector_via_w_tilde_data_imaging_from( + # w_tilde_data=self.w_tilde_data, + # data_to_pix_unique=linear_obj.unique_mappings.data_to_pix_unique, + # data_weights=linear_obj.unique_mappings.data_weights, + # pix_lengths=linear_obj.unique_mappings.pix_lengths, + # pix_pixels=linear_obj.params, + # ) @property def _data_vector_multi_mapper(self) -> np.ndarray: diff --git a/autoarray/structures/arrays/kernel_2d.py b/autoarray/structures/arrays/kernel_2d.py index 3ae0bb1cc..44c7625ff 100644 --- a/autoarray/structures/arrays/kernel_2d.py +++ b/autoarray/structures/arrays/kernel_2d.py @@ -1,5 +1,6 @@ import jax import jax.numpy as jnp +from jax import lax import numpy as np from pathlib import Path from typing import List, Tuple, Union @@ -16,6 +17,7 @@ from autoarray.structures.arrays import array_2d_util + class Kernel2D(AbstractArray2D): def __init__( self, @@ -698,6 +700,115 @@ def convolve_image_no_blurring_for_mapping(self, image, mask, jax_method="direct return Array2D(values=convolved_array_1d, mask=mask) + # def convolve_mapping_matrix(self, mapping_matrix, mask, jax_method="direct", chan_chunk=128): + # """ + # mapping_matrix : (N_masked, S) rows correspond to mask == False in row-major order + # mask : (H, W) boolean (False = kept, True = masked) + # chan_chunk : number of source columns processed at once (avoid large workspaces) + # to_f32 : cast inputs to float32 for smaller workspace & faster conv + # returns : (N_masked, S) convolved, re-masked mapping matrix + # """ + # H, W = mask.shape + # S = mapping_matrix.shape[1] + # N = H * W + # + # psf_kernel = self.stored_native.array + # + # # Indices of kept pixels (False in the mask), STATIC size for jit + # flat_mask = mask.reshape(-1) + # K = mapping_matrix.shape[0] # number of unmasked pixels + # keep_idx = jnp.where(~flat_mask, size=K, fill_value=0)[0] # (K,) + # + # # Expand to full grid (N, S): zeros everywhere, fill only kept rows + # full = jnp.zeros((N, S), dtype=mapping_matrix.dtype) + # full = full.at[keep_idx, :].set(mapping_matrix) + # + # out_cols = [] + # kH, kW = psf_kernel.shape + # + # for c0 in range(0, S, chan_chunk): + # c1 = min(c0 + chan_chunk, S) + # G = c1 - c0 # groups in this chunk + # + # # Treat the G columns as G input channels: images shape [1, G, H, W] + # imgs = full[:, c0:c1].T.reshape(1, G, H, W) + # + # # Depthwise/grouped conv: kernel shape [G, 1, kH, kW] + # # Use broadcast_to to avoid actual data replication. + # kernel = jnp.broadcast_to(psf_kernel[None, None, :, :], (G, 1, kH, kW)) + # + # y = lax.conv_general_dilated( + # imgs, + # kernel, + # window_strides=(1, 1), + # padding="SAME", + # dimension_numbers=("NCHW", "OIHW", "NCHW"), + # feature_group_count=G, # depthwise: each channel convolved independently + # ) # [1, G, H, W] + # + # # Back to (N, G) then re-mask + # y_full = y.reshape(G, N).T + # y_masked = y_full[keep_idx, :] # (K, G) + # out_cols.append(y_masked) + # + # return jnp.concatenate(out_cols, axis=1) # (K, S) + # + # + # def convolve_mapping_matrix(self, mapping_matrix, mask, jax_method="direct"): + # """ + # mapping_matrix: [N_masked, S] (rows correspond to mask==False in row-major order) + # mask : [H, W] boolean (False=kept, True=masked) + # returns : [N_masked, S] convolved, re-masked mapping matrix + # """ + # psf_kernel = self.stored_native.array + # H, W = mask.shape + # S = mapping_matrix.shape[1] + # N = H * W + # + # # Indices of kept pixels (False in the mask), with STATIC size for jit. + # flat_mask = mask.reshape(-1) + # K = mapping_matrix.shape[0] # number of unmasked pixels + # keep_idx = jnp.where(~flat_mask, size=K, fill_value=0)[0] # shape (K,) + # + # # Expand to full grid: zeros everywhere, fill only kept rows. + # full = jnp.zeros((N, S), dtype=mapping_matrix.dtype) + # full = full.at[keep_idx, :].set(mapping_matrix) + # + # # Batch conv: [S, 1, H, W] * [1, 1, kH, kW] -> [S, 1, H, W] + # images = full.T.reshape(S, 1, H, W) + # kernel = psf_kernel[jnp.newaxis, jnp.newaxis, :, :] + # + # convolved = lax.conv_general_dilated( + # images, + # kernel, + # window_strides=(1, 1), + # padding="SAME", + # dimension_numbers=("NCHW", "OIHW", "NCHW"), + # ) + # + # # Back to masked slim grid + # convolved_full = convolved.reshape(S, N).T + # convolved_masked = convolved_full[keep_idx, :] # shape (K, S) + # + # return convolved_masked + + + # def convolve_mapping_matrix(self, mapping_matrix, mask, jax_method="direct"): + + # mapping_matrix_list = [] + # + # for i in range(mapping_matrix.shape[1]): + # + # blurred_image_2d = self.convolve_image_no_blurring_for_mapping( + # image=mapping_matrix[:,i], + # mask=mask + # ) + # + # mapping_matrix_list.append(blurred_image_2d.array) + # + # return jnp.stack(mapping_matrix_list, axis=1) + + def convolve_mapping_matrix(self, mapping_matrix, mask, jax_method="direct"): """For a given 1D array and blurring array, convolve the two using this psf. @@ -709,3 +820,4 @@ def convolve_mapping_matrix(self, mapping_matrix, mask, jax_method="direct"): return jax.vmap( self.convolve_image_no_blurring_for_mapping, in_axes=(1, None, None) )(mapping_matrix, mask, jax_method).T + From e3a8341539f1c1f1da3ac3efa2a947763ad30ec4 Mon Sep 17 00:00:00 2001 From: Jammy2211 Date: Mon, 22 Sep 2025 19:06:42 +0100 Subject: [PATCH 12/17] remove commented code --- autoarray/structures/arrays/kernel_2d.py | 109 ----------------------- 1 file changed, 109 deletions(-) diff --git a/autoarray/structures/arrays/kernel_2d.py b/autoarray/structures/arrays/kernel_2d.py index 44c7625ff..64f680c89 100644 --- a/autoarray/structures/arrays/kernel_2d.py +++ b/autoarray/structures/arrays/kernel_2d.py @@ -700,115 +700,6 @@ def convolve_image_no_blurring_for_mapping(self, image, mask, jax_method="direct return Array2D(values=convolved_array_1d, mask=mask) - # def convolve_mapping_matrix(self, mapping_matrix, mask, jax_method="direct", chan_chunk=128): - # """ - # mapping_matrix : (N_masked, S) rows correspond to mask == False in row-major order - # mask : (H, W) boolean (False = kept, True = masked) - # chan_chunk : number of source columns processed at once (avoid large workspaces) - # to_f32 : cast inputs to float32 for smaller workspace & faster conv - # returns : (N_masked, S) convolved, re-masked mapping matrix - # """ - # H, W = mask.shape - # S = mapping_matrix.shape[1] - # N = H * W - # - # psf_kernel = self.stored_native.array - # - # # Indices of kept pixels (False in the mask), STATIC size for jit - # flat_mask = mask.reshape(-1) - # K = mapping_matrix.shape[0] # number of unmasked pixels - # keep_idx = jnp.where(~flat_mask, size=K, fill_value=0)[0] # (K,) - # - # # Expand to full grid (N, S): zeros everywhere, fill only kept rows - # full = jnp.zeros((N, S), dtype=mapping_matrix.dtype) - # full = full.at[keep_idx, :].set(mapping_matrix) - # - # out_cols = [] - # kH, kW = psf_kernel.shape - # - # for c0 in range(0, S, chan_chunk): - # c1 = min(c0 + chan_chunk, S) - # G = c1 - c0 # groups in this chunk - # - # # Treat the G columns as G input channels: images shape [1, G, H, W] - # imgs = full[:, c0:c1].T.reshape(1, G, H, W) - # - # # Depthwise/grouped conv: kernel shape [G, 1, kH, kW] - # # Use broadcast_to to avoid actual data replication. - # kernel = jnp.broadcast_to(psf_kernel[None, None, :, :], (G, 1, kH, kW)) - # - # y = lax.conv_general_dilated( - # imgs, - # kernel, - # window_strides=(1, 1), - # padding="SAME", - # dimension_numbers=("NCHW", "OIHW", "NCHW"), - # feature_group_count=G, # depthwise: each channel convolved independently - # ) # [1, G, H, W] - # - # # Back to (N, G) then re-mask - # y_full = y.reshape(G, N).T - # y_masked = y_full[keep_idx, :] # (K, G) - # out_cols.append(y_masked) - # - # return jnp.concatenate(out_cols, axis=1) # (K, S) - # - # - # def convolve_mapping_matrix(self, mapping_matrix, mask, jax_method="direct"): - # """ - # mapping_matrix: [N_masked, S] (rows correspond to mask==False in row-major order) - # mask : [H, W] boolean (False=kept, True=masked) - # returns : [N_masked, S] convolved, re-masked mapping matrix - # """ - # psf_kernel = self.stored_native.array - # H, W = mask.shape - # S = mapping_matrix.shape[1] - # N = H * W - # - # # Indices of kept pixels (False in the mask), with STATIC size for jit. - # flat_mask = mask.reshape(-1) - # K = mapping_matrix.shape[0] # number of unmasked pixels - # keep_idx = jnp.where(~flat_mask, size=K, fill_value=0)[0] # shape (K,) - # - # # Expand to full grid: zeros everywhere, fill only kept rows. - # full = jnp.zeros((N, S), dtype=mapping_matrix.dtype) - # full = full.at[keep_idx, :].set(mapping_matrix) - # - # # Batch conv: [S, 1, H, W] * [1, 1, kH, kW] -> [S, 1, H, W] - # images = full.T.reshape(S, 1, H, W) - # kernel = psf_kernel[jnp.newaxis, jnp.newaxis, :, :] - # - # convolved = lax.conv_general_dilated( - # images, - # kernel, - # window_strides=(1, 1), - # padding="SAME", - # dimension_numbers=("NCHW", "OIHW", "NCHW"), - # ) - # - # # Back to masked slim grid - # convolved_full = convolved.reshape(S, N).T - # convolved_masked = convolved_full[keep_idx, :] # shape (K, S) - # - # return convolved_masked - - - # def convolve_mapping_matrix(self, mapping_matrix, mask, jax_method="direct"): - - # mapping_matrix_list = [] - # - # for i in range(mapping_matrix.shape[1]): - # - # blurred_image_2d = self.convolve_image_no_blurring_for_mapping( - # image=mapping_matrix[:,i], - # mask=mask - # ) - # - # mapping_matrix_list.append(blurred_image_2d.array) - # - # return jnp.stack(mapping_matrix_list, axis=1) - - def convolve_mapping_matrix(self, mapping_matrix, mask, jax_method="direct"): """For a given 1D array and blurring array, convolve the two using this psf. From 9b1c91cf2c8615aac476f688c51454c07fc0d349 Mon Sep 17 00:00:00 2001 From: Jammy2211 Date: Tue, 23 Sep 2025 13:10:46 +0100 Subject: [PATCH 13/17] fixed mappings and weights --- .../pixelization/mappers/mapper_util.py | 59 +++++++-------- .../pixelization/mappers/rectangular.py | 47 ++++++++---- autoarray/plot/mat_plot/two_d.py | 75 +++++++++++++++---- autoarray/structures/mesh/rectangular_2d.py | 31 -------- .../structures/mesh/test_rectangular.py | 70 ++--------------- 5 files changed, 129 insertions(+), 153 deletions(-) diff --git a/autoarray/inversion/pixelization/mappers/mapper_util.py b/autoarray/inversion/pixelization/mappers/mapper_util.py index 428c7b3a0..f41aa07cf 100644 --- a/autoarray/inversion/pixelization/mappers/mapper_util.py +++ b/autoarray/inversion/pixelization/mappers/mapper_util.py @@ -124,28 +124,25 @@ def create_transforms(traced_points): # stored in a (N, 2) array and return functions that # take in (N, 2) arrays and transform the values into # the range (0, 1) and the inverse transform - N = traced_points.shape[0] + N = traced_points.shape[0] #// 2 t = jnp.arange(1, N + 1) / (N + 1) - sort_points = jnp.sort(traced_points, axis=0) + sort_points = jnp.sort(traced_points, axis=0)#[::2] transform = partial(forward_interp, sort_points, t) inv_transform = partial(reverse_interp, t, sort_points) return transform, inv_transform - def adaptive_rectangular_transformed_grid_from( - source_plane_data_grid, - grid + source_plane_data_grid, + grid ): - mu = source_plane_data_grid.mean(axis=0) scale = source_plane_data_grid.std(axis=0).min() - - grid_scaled = (grid - mu) / scale + source_grid_scaled = (source_plane_data_grid - mu) / scale - transform, inv_transform = create_transforms(grid_scaled) + transform, inv_transform = create_transforms(source_grid_scaled) def inv_full(U): return inv_transform(U) * scale + mu @@ -154,6 +151,7 @@ def inv_full(U): + def adaptive_rectangular_mappings_weights_via_interpolation_from( source_grid_size: int, source_plane_data_grid, @@ -213,10 +211,10 @@ def adaptive_rectangular_mappings_weights_via_interpolation_from( # --- Step 1. Normalize grid --- mu = source_plane_data_grid.mean(axis=0) scale = source_plane_data_grid.std(axis=0).min() - grid_scaled = (source_plane_data_grid - mu) / scale + source_grid_scaled = (source_plane_data_grid - mu) / scale # --- Step 2. Build transforms --- - transform, inv_transform = create_transforms(grid_scaled) + transform, inv_transform = create_transforms(source_grid_scaled) # --- Step 3. Transform oversampled grid into index space --- grid_over_sampled_scaled = (source_plane_data_grid_over_sampled - mu) / scale @@ -230,10 +228,10 @@ def adaptive_rectangular_mappings_weights_via_interpolation_from( iy_up = jnp.ceil(grid_over_index[:, 1]) # --- Step 5. Four corners --- - idx_bl = jnp.stack([ix_down, iy_down], axis=1) # bottom-left - idx_br = jnp.stack([ix_up, iy_down], axis=1) # bottom-right - idx_tl = jnp.stack([ix_down, iy_up], axis=1) # top-left - idx_tr = jnp.stack([ix_up, iy_up], axis=1) # top-right + idx_tl = jnp.stack([ix_up, iy_down], axis=1) + idx_tr = jnp.stack([ix_up, iy_up], axis=1) + idx_br = jnp.stack([ix_down, iy_up], axis=1) + idx_bl = jnp.stack([ix_down, iy_down], axis=1) # --- Step 6. Flatten indices --- def flatten(idx, n): @@ -241,28 +239,25 @@ def flatten(idx, n): col = idx[:, 1] return row * n + col - flat_bl = flatten(idx_bl, source_grid_size) - flat_br = flatten(idx_br, source_grid_size) flat_tl = flatten(idx_tl, source_grid_size) flat_tr = flatten(idx_tr, source_grid_size) + flat_bl = flatten(idx_bl, source_grid_size) + flat_br = flatten(idx_br, source_grid_size) - # --- Step 7. Bilinear interpolation weights --- - idx_down = jnp.stack([ix_down, iy_down], axis=1) - idx_up = jnp.stack([ix_up, iy_up], axis=1) - - delta_up = idx_up - grid_over_index - delta_down = grid_over_index - idx_down - - w_tl = delta_up[:, 0] * delta_up[:, 1] - w_bl = delta_up[:, 0] * delta_down[:, 1] - w_tr = delta_down[:, 0] * delta_up[:, 1] - w_br = delta_down[:, 0] * delta_down[:, 1] - - # --- Step 8. Stack outputs --- - flat_indices = jnp.stack([flat_bl, flat_br, flat_tl, flat_tr], axis=1).astype( + flat_indices = jnp.stack([flat_tl, flat_tr, flat_bl, flat_br], axis=1).astype( "int64" ) - weights = jnp.stack([w_bl, w_br, w_tl, w_tr], axis=1) + + # --- Step 7. Bilinear interpolation weights --- + t_row = (grid_over_index[:, 0] - ix_down) / (ix_up - ix_down + 1e-12) + t_col = (grid_over_index[:, 1] - iy_down) / (iy_up - iy_down + 1e-12) + + # Weights + w_tl = (1 - t_row) * (1 - t_col) + w_tr = (1 - t_row) * t_col + w_bl = t_row * (1 - t_col) + w_br = t_row * t_col + weights = jnp.stack([w_tl, w_tr, w_bl, w_br], axis=1) return flat_indices, weights diff --git a/autoarray/inversion/pixelization/mappers/rectangular.py b/autoarray/inversion/pixelization/mappers/rectangular.py index c54f5f8e0..2ace32ba9 100644 --- a/autoarray/inversion/pixelization/mappers/rectangular.py +++ b/autoarray/inversion/pixelization/mappers/rectangular.py @@ -98,26 +98,45 @@ 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=self.source_plane_data_grid.over_sampled, - ) - ) - # mappings, weights = ( - # mapper_util.adaptive_rectangular_mappings_weights_via_interpolation_from( - # source_grid_size=self.shape_native[0], - # source_plane_data_grid=self.source_plane_mesh_grid.array, - # source_plane_data_grid_over_sampled=jnp.array( - # self.source_plane_data_grid.over_sampled - # ), + # 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=self.source_plane_data_grid.over_sampled, # ) # ) + mappings, weights = ( + mapper_util.adaptive_rectangular_mappings_weights_via_interpolation_from( + source_grid_size=self.shape_native[0], + source_plane_data_grid=self.source_plane_data_grid.array, + source_plane_data_grid_over_sampled=jnp.array( + self.source_plane_data_grid.over_sampled + ), + ) + ) + return PixSubWeights( mappings=mappings, sizes=4 * jnp.ones(len(mappings), dtype="int"), weights=weights, ) + + @cached_property + def edges_transformed(self): + """ + A class packing the ndarrays describing the neighbors of every pixel in the rectangular pixelization (see + `Neighbors` for a complete description of the neighboring scheme). + + The neighbors of a rectangular pixelization are computed by exploiting the uniform and symmetric nature of the + rectangular grid, as described in the method `mesh_util.rectangular_neighbors_from`. + """ + + # edges defined in 0 -> 1 space, there is one more edge than pixel centers on each side + edges = jnp.linspace(0, 1, self.shape_native[0] + 1) + edges_reshaped = jnp.stack([edges, edges]).T + + return mapper_util.adaptive_rectangular_transformed_grid_from( + source_plane_data_grid=self.source_plane_data_grid.array, + grid=edges_reshaped + ) \ No newline at end of file diff --git a/autoarray/plot/mat_plot/two_d.py b/autoarray/plot/mat_plot/two_d.py index 5e38bd896..1aa8e84d1 100644 --- a/autoarray/plot/mat_plot/two_d.py +++ b/autoarray/plot/mat_plot/two_d.py @@ -564,18 +564,68 @@ def _plot_rectangular_mapper( else: ax = self.setup_subplot(aspect=aspect_inv) + shape_native = mapper.source_plane_mesh_grid.shape_native + if pixel_values is not None: - self.plot_array( - array=pixel_values, - visuals_2d=visuals_2d, - auto_labels=auto_labels, - bypass=True, + + # self.plot_array( + # array=pixel_values, + # visuals_2d=visuals_2d, + # auto_labels=auto_labels, + # bypass=True, + # ) + + norm = self.cmap.norm_from(array=pixel_values.array, use_log10=self.use_log10) + + edges_transformed = mapper.edges_transformed + + edges_transformed_dense = np.moveaxis( + np.stack(np.meshgrid(*edges_transformed.T)), + 0, + 2 ) - self.axis.set(extent=extent, grid=mapper.source_plane_mesh_grid) + plt.pcolormesh( + edges_transformed_dense[..., 0], + edges_transformed_dense[..., 1], + pixel_values.array.reshape(shape_native), + shading='flat', + norm=norm, + cmap=self.cmap.cmap, + ) - self.yticks.set(min_value=extent[2], max_value=extent[3], units=self.units) - self.xticks.set(min_value=extent[0], max_value=extent[1], units=self.units) + if self.colorbar is not False: + + cb = self.colorbar.set( + units=self.units, + ax=ax, + norm=norm, + cb_unit=auto_labels.cb_unit, + use_log10=self.use_log10, + ) + self.colorbar_tickparams.set(cb=cb) + + extent_axis = self.axis.config_dict.get("extent") + + if extent_axis is None: + extent_axis = extent + + self.axis.set(extent=extent_axis) + + self.tickparams.set() + self.yticks.set( + min_value=extent_axis[2], + max_value=extent_axis[3], + units=self.units, + pixels=shape_native[0], + ) + + self.xticks.set( + min_value=extent_axis[0], + max_value=extent_axis[1], + units=self.units, + pixels=shape_native[1], + ) if not isinstance(self.text, list): self.text.set() @@ -587,13 +637,12 @@ def _plot_rectangular_mapper( else: [annotate.set() for annotate in self.annotate] - self.grid_plot.plot_rectangular_grid_lines( - extent=mapper.source_plane_mesh_grid.geometry.extent, - shape_native=mapper.shape_native, - ) + # self.grid_plot.plot_rectangular_grid_lines( + # extent=mapper.source_plane_mesh_grid.geometry.extent, + # shape_native=mapper.shape_native, + # ) self.title.set(auto_title=auto_labels.title) - self.tickparams.set() self.ylabel.set() self.xlabel.set() diff --git a/autoarray/structures/mesh/rectangular_2d.py b/autoarray/structures/mesh/rectangular_2d.py index 744e06b16..59736c1e5 100644 --- a/autoarray/structures/mesh/rectangular_2d.py +++ b/autoarray/structures/mesh/rectangular_2d.py @@ -130,37 +130,6 @@ def neighbors(self) -> Neighbors: return Neighbors(arr=neighbors.astype("int"), sizes=sizes.astype("int")) - @cached_property - def edges(self): - """ - A class packing the ndarrays describing the neighbors of every pixel in the rectangular pixelization (see - `Neighbors` for a complete description of the neighboring scheme). - - The neighbors of a rectangular pixelization are computed by exploiting the uniform and symmetric nature of the - rectangular grid, as described in the method `mesh_util.rectangular_neighbors_from`. - """ - return mesh_util.rectangular_edges_from( - shape_native=self.shape_native, - pixel_scales=self.pixel_scales, - ) - - @cached_property - def edges_transformed(self): - """ - A class packing the ndarrays describing the neighbors of every pixel in the rectangular pixelization (see - `Neighbors` for a complete description of the neighboring scheme). - - The neighbors of a rectangular pixelization are computed by exploiting the uniform and symmetric nature of the - rectangular grid, as described in the method `mesh_util.rectangular_neighbors_from`. - """ - - edges_transformed = mapper_util.adaptive_rectangular_transformed_grid_from( - source_plane_data_grid=self.array, - grid=self.edges.reshape(-1, 2), - ) - - return edges_transformed.reshape(self.edges.shape) - @cached_property def edge_pixel_list(self) -> List: return mesh_util.rectangular_edge_pixel_list_from( diff --git a/test_autoarray/structures/mesh/test_rectangular.py b/test_autoarray/structures/mesh/test_rectangular.py index 8d22c5dae..42ab9d9ce 100644 --- a/test_autoarray/structures/mesh/test_rectangular.py +++ b/test_autoarray/structures/mesh/test_rectangular.py @@ -23,8 +23,7 @@ def test__neighbors__compare_to_mesh_util(): assert (mesh.neighbors == neighbors_util).all() assert (mesh.neighbors.sizes == neighbors_sizes_util).all() - -def test__edges(): +def test__edges_transformed(): grid = np.array( [ @@ -44,17 +43,13 @@ def test__edges(): shape_native=(3, 3), grid=grid, buffer=1e-8 ) - assert mesh.edges[0] == pytest.approx( + print(mesh.edges_transformed[0]) + + assert mesh.edges_transformed[0] == pytest.approx( np.array( - [ - [[1.5, -1.5], [1.5, -0.5]], # top - [[1.5, -0.5], [0.5, -0.5]], # right - [[0.5, -0.5], [0.5, -1.5]], # bottom - [[0.5, -1.5], [1.5, -1.5]], - ] - ), # left - ) - assert mesh.edges[1] == pytest.approx( + [-1.00000001, -1.00000001], # left + ), abs=1e-8) + assert mesh.edges_transformed[1] == pytest.approx( np.array( [ [[1.5, -0.5], [1.5, 0.5]], # top @@ -64,57 +59,6 @@ def test__edges(): ] ), # left ) - assert mesh.edges[5] == pytest.approx( - np.array( - [ - [[0.5, 0.5], [0.5, 1.5]], # top - [[0.5, 1.5], [-0.5, 1.5]], # right - [[-0.5, 1.5], [-0.5, 0.5]], # bottom - [[-0.5, 0.5], [0.5, 0.5]], - ] - ), # left - ) - -# def test__edges_transformed(): -# -# grid = np.array( -# [ -# [-1.5, -1.5], -# [-1.5, 0.0], -# [-1.5, 1.5], -# [0.0, -1.5], -# [0.0, 0.0], -# [0.0, 1.5], -# [1.5, -1.5], -# [1.5, 0.0], -# [1.5, 1.5], -# ] -# ) -# -# mesh = aa.Mesh2DRectangular.overlay_grid( -# shape_native=(3, 3), grid=grid, buffer=1e-8 -# ) -# -# assert mesh.edges_transformed[0] == pytest.approx( -# np.array( -# [ -# [[1.5, -1.5], [1.5, -0.5]], # top -# [[1.5, -0.5], [0.5, -0.5]], # right -# [[0.5, -0.5], [0.5, -1.5]], # bottom -# [[0.5, -1.5], [1.5, -1.5]], -# ] -# ), # left -# ) -# assert mesh.edges[1] == pytest.approx( -# np.array( -# [ -# [[1.5, -0.5], [1.5, 0.5]], # top -# [[1.5, 0.5], [0.5, 0.5]], # right -# [[0.5, 0.5], [0.5, -0.5]], # bottom -# [[0.5, -0.5], [1.5, -0.5]], -# ] -# ), # left -# ) def test__edge_pixel_list(): grid = np.array( From 216a8fc987fdbead338cb38f892d3a01459e9cdc Mon Sep 17 00:00:00 2001 From: Jammy2211 Date: Tue, 23 Sep 2025 14:53:32 +0100 Subject: [PATCH 14/17] test__edges_transformed --- autoarray/dataset/imaging/w_tilde.py | 10 ++--- .../imaging/inversion_imaging_util.py | 23 +++++++----- .../inversion/inversion/imaging/w_tilde.py | 25 ++++++++----- .../pixelization/mappers/mapper_util.py | 31 +++++++++++++--- .../pixelization/mappers/rectangular.py | 18 ++++++++- autoarray/plot/mat_plot/two_d.py | 10 ++--- autoarray/structures/arrays/kernel_2d.py | 2 - autoarray/structures/mesh/rectangular_2d.py | 1 + .../inversion/imaging/test_imaging.py | 8 +++- .../inversion/test_inversion_util.py | 1 + .../pixelization/mappers/test_mapper_util.py | 2 +- .../pixelization/mappers/test_rectangular.py | 37 +++++++++++++++++++ .../structures/mesh/test_rectangular.py | 36 ------------------ 13 files changed, 126 insertions(+), 78 deletions(-) diff --git a/autoarray/dataset/imaging/w_tilde.py b/autoarray/dataset/imaging/w_tilde.py index 7fc1fdf51..2cb5a5778 100644 --- a/autoarray/dataset/imaging/w_tilde.py +++ b/autoarray/dataset/imaging/w_tilde.py @@ -17,9 +17,9 @@ def __init__( curvature_preload: np.ndarray, indexes: np.ndim, lengths: np.ndarray, - noise_map : np.ndarray, - psf : np.ndarray, - mask : np.ndarray, + noise_map: np.ndarray, + psf: np.ndarray, + mask: np.ndarray, noise_map_value: float, ): """ @@ -102,5 +102,5 @@ def psf_operator_matrix_dense(self): self.mask.derive_indexes.native_for_slim ).astype("int"), native_shape=self.noise_map.shape_native, - correlate=False - ) \ No newline at end of file + correlate=False, + ) diff --git a/autoarray/inversion/inversion/imaging/inversion_imaging_util.py b/autoarray/inversion/inversion/imaging/inversion_imaging_util.py index 009db6526..285fb5eeb 100644 --- a/autoarray/inversion/inversion/imaging/inversion_imaging_util.py +++ b/autoarray/inversion/inversion/imaging/inversion_imaging_util.py @@ -8,6 +8,7 @@ import numpy as np + def psf_operator_matrix_dense_from( kernel_native: np.ndarray, native_index_for_slim_index: np.ndarray, # shape (N_pix, 2), native (y,x) coords of masked pixels @@ -41,7 +42,7 @@ def psf_operator_matrix_dense_from( ker = kernel_native if correlate else kernel_native[::-1, ::-1] # Padded index grid: -1 everywhere, slim index where masked - index_padded = -np.ones((Ny + 2*ph, Nx + 2*pw), dtype=np.int64) + index_padded = -np.ones((Ny + 2 * ph, Nx + 2 * pw), dtype=np.int64) for p, (y, x) in enumerate(native_index_for_slim_index): index_padded[y + ph, x + pw] = p @@ -103,31 +104,33 @@ def w_tilde_data_imaging_from( """ # 1) weight map = image / noise^2 (safe where noise==0) - weight_map = jnp.where(noise_map_native > 0.0, - image_native / (noise_map_native ** 2), - 0.0) + weight_map = jnp.where( + noise_map_native > 0.0, image_native / (noise_map_native**2), 0.0 + ) Ky, Kx = kernel_native.shape ph, pw = Ky // 2, Kx // 2 # 2) pad so neighbourhood gathers never go OOB - padded = jnp.pad(weight_map, ((ph, ph), (pw, pw)), mode="constant", constant_values=0.0) + padded = jnp.pad( + weight_map, ((ph, ph), (pw, pw)), mode="constant", constant_values=0.0 + ) # 3) build broadcasted neighbourhood indices for all requested pixels # shift pixel coords into the padded frame - ys = native_index_for_slim_index[:, 0] + ph # (N,) - xs = native_index_for_slim_index[:, 1] + pw # (N,) + ys = native_index_for_slim_index[:, 0] + ph # (N,) + xs = native_index_for_slim_index[:, 1] + pw # (N,) # kernel-relative offsets - dy = jnp.arange(Ky) - ph # (Ky,) - dx = jnp.arange(Kx) - pw # (Kx,) + dy = jnp.arange(Ky) - ph # (Ky,) + dx = jnp.arange(Kx) - pw # (Kx,) # broadcast to (N, Ky, Kx) Y = ys[:, None, None] + dy[None, :, None] X = xs[:, None, None] + dx[None, None, :] # 4) gather patches and correlate (no kernel flip) - patches = padded[Y, X] # (N, Ky, Kx) + patches = padded[Y, X] # (N, Ky, Kx) return jnp.sum(patches * kernel_native[None, :, :], axis=(1, 2)) # (N,) diff --git a/autoarray/inversion/inversion/imaging/w_tilde.py b/autoarray/inversion/inversion/imaging/w_tilde.py index a61b311d1..1cb7a2f6d 100644 --- a/autoarray/inversion/inversion/imaging/w_tilde.py +++ b/autoarray/inversion/inversion/imaging/w_tilde.py @@ -79,7 +79,7 @@ def w_tilde_data(self): image_native=self.data.native.array, noise_map_native=self.noise_map.native.array, kernel_native=self.psf.native.array, - native_index_for_slim_index=self.data.mask.derive_indexes.native_for_slim + native_index_for_slim_index=self.data.mask.derive_indexes.native_for_slim, ) @property @@ -362,11 +362,10 @@ def _curvature_matrix_x1_mapper(self) -> np.ndarray: """ return inversion_util.curvature_matrix_via_w_tilde_from( - w_tilde=self.w_tilde.w_matrix, - mapping_matrix=self.mapping_matrix + w_tilde=self.w_tilde.w_matrix, mapping_matrix=self.mapping_matrix ) -# return self._curvature_matrix_mapper_diag + # return self._curvature_matrix_mapper_diag @property def _curvature_matrix_multi_mapper(self) -> np.ndarray: @@ -518,7 +517,9 @@ def mapped_reconstructed_data_dict(self) -> Dict[LinearObj, Array2D]: for linear_obj in self.linear_obj_list: reconstruction = reconstruction_dict[linear_obj] - if isinstance(linear_obj, AbstractMapper) and self.has(cls=AbstractLinearObjFuncList): + if isinstance(linear_obj, AbstractMapper) and self.has( + cls=AbstractLinearObjFuncList + ): mapped_reconstructed_image = inversion_util.mapped_reconstructed_data_via_image_to_pix_unique_from( data_to_pix_unique=linear_obj.unique_mappings.data_to_pix_unique, @@ -535,12 +536,16 @@ def mapped_reconstructed_data_dict(self) -> Dict[LinearObj, Array2D]: values=mapped_reconstructed_image, mask=self.mask ) - elif isinstance(linear_obj, AbstractMapper) and not self.has(cls=AbstractLinearObjFuncList): + elif isinstance(linear_obj, AbstractMapper) and not self.has( + cls=AbstractLinearObjFuncList + ): - mapped_reconstructed_image = inversion_util.mapped_reconstructed_data_via_w_tilde_from( - w_tilde=self.w_tilde.psf_operator_matrix_dense, - mapping_matrix=self.mapping_matrix, - reconstruction=reconstruction, + mapped_reconstructed_image = ( + inversion_util.mapped_reconstructed_data_via_w_tilde_from( + w_tilde=self.w_tilde.psf_operator_matrix_dense, + mapping_matrix=self.mapping_matrix, + reconstruction=reconstruction, + ) ) mapped_reconstructed_image = Array2D( diff --git a/autoarray/inversion/pixelization/mappers/mapper_util.py b/autoarray/inversion/pixelization/mappers/mapper_util.py index f41aa07cf..47e183621 100644 --- a/autoarray/inversion/pixelization/mappers/mapper_util.py +++ b/autoarray/inversion/pixelization/mappers/mapper_util.py @@ -124,20 +124,17 @@ def create_transforms(traced_points): # stored in a (N, 2) array and return functions that # take in (N, 2) arrays and transform the values into # the range (0, 1) and the inverse transform - N = traced_points.shape[0] #// 2 + N = traced_points.shape[0] # // 2 t = jnp.arange(1, N + 1) / (N + 1) - sort_points = jnp.sort(traced_points, axis=0)#[::2] + sort_points = jnp.sort(traced_points, axis=0) # [::2] transform = partial(forward_interp, sort_points, t) inv_transform = partial(reverse_interp, t, sort_points) return transform, inv_transform -def adaptive_rectangular_transformed_grid_from( - source_plane_data_grid, - grid -): +def adaptive_rectangular_transformed_grid_from(source_plane_data_grid, grid): mu = source_plane_data_grid.mean(axis=0) scale = source_plane_data_grid.std(axis=0).min() source_grid_scaled = (source_plane_data_grid - mu) / scale @@ -150,6 +147,28 @@ def inv_full(U): return inv_full(grid) +def adaptive_rectangular_areas_from(source_grid_size, source_plane_data_grid): + + pixel_edges_1d = jnp.linspace(0, 1, source_grid_size + 1) + + mu = source_plane_data_grid.mean(axis=0) + scale = source_plane_data_grid.std(axis=0).min() + source_grid_scaled = (source_plane_data_grid - mu) / scale + + transform, inv_transform = create_transforms(source_grid_scaled) + + def inv_full(U): + return inv_transform(U) * scale + mu + + pixel_edges = inv_full(jnp.stack([pixel_edges_1d, pixel_edges_1d]).T) + + # lengths along each axis + pixel_lengths = jnp.diff(pixel_edges, axis=0).squeeze() # shape (N_source, 2) + + dy = pixel_lengths[:, 0] + dx = pixel_lengths[:, 1] + + return jnp.outer(dy, dx).flatten() def adaptive_rectangular_mappings_weights_via_interpolation_from( diff --git a/autoarray/inversion/pixelization/mappers/rectangular.py b/autoarray/inversion/pixelization/mappers/rectangular.py index 2ace32ba9..f9f398c69 100644 --- a/autoarray/inversion/pixelization/mappers/rectangular.py +++ b/autoarray/inversion/pixelization/mappers/rectangular.py @@ -122,6 +122,20 @@ def pix_sub_weights(self) -> PixSubWeights: weights=weights, ) + @cached_property + def areas_transformed(self): + """ + A class packing the ndarrays describing the neighbors of every pixel in the rectangular pixelization (see + `Neighbors` for a complete description of the neighboring scheme). + + The neighbors of a rectangular pixelization are computed by exploiting the uniform and symmetric nature of the + rectangular grid, as described in the method `mesh_util.rectangular_neighbors_from`. + """ + return mapper_util.adaptive_rectangular_areas_from( + source_grid_size=self.shape_native[0], + source_plane_data_grid=self.source_plane_data_grid.array, + ) + @cached_property def edges_transformed(self): """ @@ -138,5 +152,5 @@ def edges_transformed(self): return mapper_util.adaptive_rectangular_transformed_grid_from( source_plane_data_grid=self.source_plane_data_grid.array, - grid=edges_reshaped - ) \ No newline at end of file + grid=edges_reshaped, + ) diff --git a/autoarray/plot/mat_plot/two_d.py b/autoarray/plot/mat_plot/two_d.py index 1aa8e84d1..f9744803d 100644 --- a/autoarray/plot/mat_plot/two_d.py +++ b/autoarray/plot/mat_plot/two_d.py @@ -575,21 +575,21 @@ def _plot_rectangular_mapper( # bypass=True, # ) - norm = self.cmap.norm_from(array=pixel_values.array, use_log10=self.use_log10) + norm = self.cmap.norm_from( + array=pixel_values.array, use_log10=self.use_log10 + ) edges_transformed = mapper.edges_transformed edges_transformed_dense = np.moveaxis( - np.stack(np.meshgrid(*edges_transformed.T)), - 0, - 2 + np.stack(np.meshgrid(*edges_transformed.T)), 0, 2 ) plt.pcolormesh( edges_transformed_dense[..., 0], edges_transformed_dense[..., 1], pixel_values.array.reshape(shape_native), - shading='flat', + shading="flat", norm=norm, cmap=self.cmap.cmap, ) diff --git a/autoarray/structures/arrays/kernel_2d.py b/autoarray/structures/arrays/kernel_2d.py index 64f680c89..b2992278a 100644 --- a/autoarray/structures/arrays/kernel_2d.py +++ b/autoarray/structures/arrays/kernel_2d.py @@ -17,7 +17,6 @@ from autoarray.structures.arrays import array_2d_util - class Kernel2D(AbstractArray2D): def __init__( self, @@ -711,4 +710,3 @@ def convolve_mapping_matrix(self, mapping_matrix, mask, jax_method="direct"): return jax.vmap( self.convolve_image_no_blurring_for_mapping, in_axes=(1, None, None) )(mapping_matrix, mask, jax_method).T - diff --git a/autoarray/structures/mesh/rectangular_2d.py b/autoarray/structures/mesh/rectangular_2d.py index 59736c1e5..35e56003b 100644 --- a/autoarray/structures/mesh/rectangular_2d.py +++ b/autoarray/structures/mesh/rectangular_2d.py @@ -17,6 +17,7 @@ from autoarray.inversion.pixelization.mesh import mesh_util from autoarray.structures.grids import grid_2d_util + class Mesh2DRectangular(Abstract2DMesh): def __init__( diff --git a/test_autoarray/inversion/inversion/imaging/test_imaging.py b/test_autoarray/inversion/inversion/imaging/test_imaging.py index c77bd0862..46bfe55ac 100644 --- a/test_autoarray/inversion/inversion/imaging/test_imaging.py +++ b/test_autoarray/inversion/inversion/imaging/test_imaging.py @@ -144,7 +144,13 @@ def test__w_tilde_checks_noise_map_and_raises_exception_if_preloads_dont_match_n grid = aa.Grid2D.from_mask(mask=mask) w_tilde = WTildeImaging( - curvature_preload=None, indexes=None, lengths=None, noise_map_value=2.0, noise_map=None, psf=None, mask=mask + curvature_preload=None, + indexes=None, + lengths=None, + noise_map_value=2.0, + noise_map=None, + psf=None, + mask=mask, ) with pytest.raises(exc.InversionException): diff --git a/test_autoarray/inversion/inversion/test_inversion_util.py b/test_autoarray/inversion/inversion/test_inversion_util.py index e2eb96a3d..cd4132f1d 100644 --- a/test_autoarray/inversion/inversion/test_inversion_util.py +++ b/test_autoarray/inversion/inversion/test_inversion_util.py @@ -86,6 +86,7 @@ def test__reconstruction_positive_negative_from(): assert reconstruction == pytest.approx(np.array([1.0, -1.0, 3.0]), 1.0e-4) + def test__mapped_reconstructed_data_via_mapping_matrix_from(): mapping_matrix = np.array([[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]]) diff --git a/test_autoarray/inversion/pixelization/mappers/test_mapper_util.py b/test_autoarray/inversion/pixelization/mappers/test_mapper_util.py index a5b41a15f..2245dda72 100644 --- a/test_autoarray/inversion/pixelization/mappers/test_mapper_util.py +++ b/test_autoarray/inversion/pixelization/mappers/test_mapper_util.py @@ -13,7 +13,7 @@ def make_five_pixels(): return np.array([[0, 0], [0, 1], [1, 0], [1, 1], [1, 2]]) -def _test__sub_slim_indexes_for_pix_index(): +def test__sub_slim_indexes_for_pix_index(): pix_indexes_for_sub_slim_index = np.array( [[0, 4], [1, 4], [2, 4], [0, 4], [1, 4], [3, 4], [0, 4], [3, 4]] ).astype("int") diff --git a/test_autoarray/inversion/pixelization/mappers/test_rectangular.py b/test_autoarray/inversion/pixelization/mappers/test_rectangular.py index ef8123b99..68e38e4e6 100644 --- a/test_autoarray/inversion/pixelization/mappers/test_rectangular.py +++ b/test_autoarray/inversion/pixelization/mappers/test_rectangular.py @@ -1,4 +1,5 @@ import numpy as np +import pytest import autoarray as aa @@ -72,3 +73,39 @@ def test__pixel_signals_from__matches_util(grid_2d_sub_1_7x7, image_7x7): ) assert (pixel_signals == pixel_signals_util).all() + + +def test__edges_transformed(mask_2d_7x7): + + grid = aa.Grid2DIrregular( + [ + [-1.5, -1.5], + [-1.5, 0.0], + [-1.5, 1.5], + [0.0, -1.5], + [0.0, 0.0], + [0.0, 1.5], + [1.5, -1.5], + [1.5, 0.0], + [1.5, 1.5], + ], + ) + + mesh = aa.Mesh2DRectangular.overlay_grid( + shape_native=(3, 3), grid=grid, buffer=1e-8 + ) + + mapper_grids = aa.MapperGrids( + mask=mask_2d_7x7, + source_plane_data_grid=grid, + source_plane_mesh_grid=mesh, + ) + + mapper = aa.Mapper(mapper_grids=mapper_grids, regularization=None) + + assert mapper.edges_transformed[4] == pytest.approx( + np.array( + [1.5, 1.5], # left + ), + abs=1e-8, + ) diff --git a/test_autoarray/structures/mesh/test_rectangular.py b/test_autoarray/structures/mesh/test_rectangular.py index 42ab9d9ce..a63489733 100644 --- a/test_autoarray/structures/mesh/test_rectangular.py +++ b/test_autoarray/structures/mesh/test_rectangular.py @@ -23,42 +23,6 @@ def test__neighbors__compare_to_mesh_util(): assert (mesh.neighbors == neighbors_util).all() assert (mesh.neighbors.sizes == neighbors_sizes_util).all() -def test__edges_transformed(): - - grid = np.array( - [ - [-1.5, -1.5], - [-1.5, 0.0], - [-1.5, 1.5], - [0.0, -1.5], - [0.0, 0.0], - [0.0, 1.5], - [1.5, -1.5], - [1.5, 0.0], - [1.5, 1.5], - ] - ) - - mesh = aa.Mesh2DRectangular.overlay_grid( - shape_native=(3, 3), grid=grid, buffer=1e-8 - ) - - print(mesh.edges_transformed[0]) - - assert mesh.edges_transformed[0] == pytest.approx( - np.array( - [-1.00000001, -1.00000001], # left - ), abs=1e-8) - assert mesh.edges_transformed[1] == pytest.approx( - np.array( - [ - [[1.5, -0.5], [1.5, 0.5]], # top - [[1.5, 0.5], [0.5, 0.5]], # right - [[0.5, 0.5], [0.5, -0.5]], # bottom - [[0.5, -0.5], [1.5, -0.5]], - ] - ), # left - ) def test__edge_pixel_list(): grid = np.array( From 2ecd33b52f24596f58a1a73ebe03ab71d1c3a111 Mon Sep 17 00:00:00 2001 From: Jammy2211 Date: Tue, 23 Sep 2025 14:59:08 +0100 Subject: [PATCH 15/17] test__areas_transformed --- .../pixelization/mappers/test_rectangular.py | 37 +++++++++++++++++++ 1 file changed, 37 insertions(+) diff --git a/test_autoarray/inversion/pixelization/mappers/test_rectangular.py b/test_autoarray/inversion/pixelization/mappers/test_rectangular.py index 68e38e4e6..eecf4b529 100644 --- a/test_autoarray/inversion/pixelization/mappers/test_rectangular.py +++ b/test_autoarray/inversion/pixelization/mappers/test_rectangular.py @@ -75,6 +75,43 @@ def test__pixel_signals_from__matches_util(grid_2d_sub_1_7x7, image_7x7): assert (pixel_signals == pixel_signals_util).all() + + +def test__areas_transformed(mask_2d_7x7): + + grid = aa.Grid2DIrregular( + [ + [-1.5, -1.5], + [-1.5, 0.0], + [-1.5, 1.5], + [0.0, -1.5], + [0.0, 0.0], + [0.0, 1.5], + [1.5, -1.5], + [1.5, 0.0], + [1.5, 1.5], + ], + ) + + mesh = aa.Mesh2DRectangular.overlay_grid( + shape_native=(3, 3), grid=grid, buffer=1e-8 + ) + + mapper_grids = aa.MapperGrids( + mask=mask_2d_7x7, + source_plane_data_grid=grid, + source_plane_mesh_grid=mesh, + ) + + mapper = aa.Mapper(mapper_grids=mapper_grids, regularization=None) + + assert mapper.areas_transformed[4] == pytest.approx( + 4.0, + abs=1e-8, + ) + + + def test__edges_transformed(mask_2d_7x7): grid = aa.Grid2DIrregular( From 9c52d10998c61177bf2dd6d8497016b2691eb4ba Mon Sep 17 00:00:00 2001 From: Jammy2211 Date: Tue, 23 Sep 2025 19:25:24 +0100 Subject: [PATCH 16/17] all of inversion tests pass --- autoarray/__init__.py | 2 + autoarray/fixtures.py | 2 +- .../inversion/pixelization/mappers/factory.py | 12 +- .../pixelization/mappers/mapper_util.py | 2 - .../mappers/rectangular_uniform.py | 106 ++++++++++++++++++ .../inversion/pixelization/mesh/__init__.py | 1 + .../pixelization/mesh/rectangular.py | 2 +- .../pixelization/mesh/rectangular_uniform.py | 31 +++++ autoarray/plot/mat_plot/two_d.py | 99 ++++++++-------- autoarray/structures/mesh/rectangular_2d.py | 2 +- .../structures/mesh/rectangular_2d_uniform.py | 5 + .../imaging/test_inversion_imaging_util.py | 6 +- .../inversion/inversion/test_abstract.py | 4 +- .../inversion/inversion/test_factory.py | 10 +- .../pixelization/mappers/test_factory.py | 4 +- .../pixelization/mappers/test_mapper_util.py | 56 --------- .../pixelization/mappers/test_rectangular.py | 11 +- .../structures/mesh/test_rectangular.py | 16 +-- 18 files changed, 236 insertions(+), 135 deletions(-) create mode 100644 autoarray/inversion/pixelization/mappers/rectangular_uniform.py create mode 100644 autoarray/inversion/pixelization/mesh/rectangular_uniform.py create mode 100644 autoarray/structures/mesh/rectangular_2d_uniform.py diff --git a/autoarray/__init__.py b/autoarray/__init__.py index 789dac386..54dc7a84f 100644 --- a/autoarray/__init__.py +++ b/autoarray/__init__.py @@ -43,6 +43,7 @@ from .inversion.pixelization.mappers.rectangular import MapperRectangular from .inversion.pixelization.mappers.delaunay import MapperDelaunay from .inversion.pixelization.mappers.voronoi import MapperVoronoi +from .inversion.pixelization.mappers.rectangular_uniform import MapperRectangularUniform from .inversion.pixelization.image_mesh.abstract import AbstractImageMesh from .inversion.pixelization.mesh.abstract import AbstractMesh from .inversion.inversion.imaging.mapping import InversionImagingMapping @@ -75,6 +76,7 @@ from .operators.over_sampling.over_sampler import OverSampler from .structures.grids.irregular_2d import Grid2DIrregular from .structures.mesh.rectangular_2d import Mesh2DRectangular +from .structures.mesh.rectangular_2d_uniform import Mesh2DRectangularUniform from .structures.mesh.voronoi_2d import Mesh2DVoronoi from .structures.mesh.delaunay_2d import Mesh2DDelaunay from .structures.arrays.kernel_2d import Kernel2D diff --git a/autoarray/fixtures.py b/autoarray/fixtures.py index 75957966f..6c1f14264 100644 --- a/autoarray/fixtures.py +++ b/autoarray/fixtures.py @@ -421,7 +421,7 @@ def make_rectangular_mapper_7x7_3x3(): adapt_data=aa.Array2D.ones(shape_native=(3, 3), pixel_scales=0.1), ) - return aa.MapperRectangular( + return aa.MapperRectangularUniform( mapper_grids=mapper_grids, border_relocator=make_border_relocator_2d_7x7(), regularization=make_regularization_constant(), diff --git a/autoarray/inversion/pixelization/mappers/factory.py b/autoarray/inversion/pixelization/mappers/factory.py index 965213b03..689f35011 100644 --- a/autoarray/inversion/pixelization/mappers/factory.py +++ b/autoarray/inversion/pixelization/mappers/factory.py @@ -4,6 +4,7 @@ from autoarray.inversion.pixelization.border_relocator import BorderRelocator from autoarray.inversion.regularization.abstract import AbstractRegularization from autoarray.structures.mesh.rectangular_2d import Mesh2DRectangular +from autoarray.structures.mesh.rectangular_2d_uniform import Mesh2DRectangularUniform from autoarray.structures.mesh.delaunay_2d import Mesh2DDelaunay from autoarray.structures.mesh.voronoi_2d import Mesh2DVoronoi @@ -39,10 +40,19 @@ def mapper_from( from autoarray.inversion.pixelization.mappers.rectangular import ( MapperRectangular, ) + from autoarray.inversion.pixelization.mappers.rectangular_uniform import ( + MapperRectangularUniform, + ) from autoarray.inversion.pixelization.mappers.delaunay import MapperDelaunay from autoarray.inversion.pixelization.mappers.voronoi import MapperVoronoi - if isinstance(mapper_grids.source_plane_mesh_grid, Mesh2DRectangular): + if isinstance(mapper_grids.source_plane_mesh_grid, Mesh2DRectangularUniform): + return MapperRectangularUniform( + mapper_grids=mapper_grids, + border_relocator=border_relocator, + regularization=regularization, + ) + elif isinstance(mapper_grids.source_plane_mesh_grid, Mesh2DRectangular): return MapperRectangular( mapper_grids=mapper_grids, border_relocator=border_relocator, diff --git a/autoarray/inversion/pixelization/mappers/mapper_util.py b/autoarray/inversion/pixelization/mappers/mapper_util.py index 47e183621..81b562bd2 100644 --- a/autoarray/inversion/pixelization/mappers/mapper_util.py +++ b/autoarray/inversion/pixelization/mappers/mapper_util.py @@ -161,8 +161,6 @@ def inv_full(U): return inv_transform(U) * scale + mu pixel_edges = inv_full(jnp.stack([pixel_edges_1d, pixel_edges_1d]).T) - - # lengths along each axis pixel_lengths = jnp.diff(pixel_edges, axis=0).squeeze() # shape (N_source, 2) dy = pixel_lengths[:, 0] diff --git a/autoarray/inversion/pixelization/mappers/rectangular_uniform.py b/autoarray/inversion/pixelization/mappers/rectangular_uniform.py new file mode 100644 index 000000000..0f2461672 --- /dev/null +++ b/autoarray/inversion/pixelization/mappers/rectangular_uniform.py @@ -0,0 +1,106 @@ +import jax.numpy as jnp + +from autoconf import cached_property + +from autoarray.inversion.pixelization.mappers.rectangular import MapperRectangular +from autoarray.inversion.pixelization.mappers.abstract import PixSubWeights + +from autoarray.inversion.pixelization.mappers import mapper_util + + +class MapperRectangularUniform(MapperRectangular): + """ + To understand a `Mapper` one must be familiar `Mesh` objects and the `mesh` and `pixelization` packages, where + the four grids grouped in a `MapperGrids` object are explained (`image_plane_data_grid`, `source_plane_data_grid`, + `image_plane_mesh_grid`,`source_plane_mesh_grid`) + + If you are unfamliar withe above objects, read through the docstrings of the `pixelization`, `mesh` and + `mapper_grids` packages. + + A `Mapper` determines the mappings between the masked data grid's pixels (`image_plane_data_grid` and + `source_plane_data_grid`) and the mesh's pixels (`image_plane_mesh_grid` and `source_plane_mesh_grid`). + + The 1D Indexing of each grid is identical in the `data` and `source` frames (e.g. the transformation does not + change the indexing, such that `source_plane_data_grid[0]` corresponds to the transformed value + of `image_plane_data_grid[0]` and so on). + + A mapper therefore only needs to determine the index mappings between the `grid_slim` and `mesh_grid`, + noting that associations are made by pairing `source_plane_mesh_grid` with `source_plane_data_grid`. + + Mappings are represented in the 2D ndarray `pix_indexes_for_sub_slim_index`, whereby the index of + a pixel on the `mesh_grid` maps to the index of a pixel on the `grid_slim` as follows: + + - pix_indexes_for_sub_slim_index[0, 0] = 0: the data's 1st sub-pixel maps to the mesh's 1st pixel. + - pix_indexes_for_sub_slim_index[1, 0] = 3: the data's 2nd sub-pixel maps to the mesh's 4th pixel. + - pix_indexes_for_sub_slim_index[2, 0] = 1: the data's 3rd sub-pixel maps to the mesh's 2nd pixel. + + The second dimension of this array (where all three examples above are 0) is used for cases where a + single pixel on the `grid_slim` maps to multiple pixels on the `mesh_grid`. For example, a + `Delaunay` triangulation, where every `grid_slim` pixel maps to three Delaunay pixels (the corners of the + triangles) with varying interpolation weights . + + For a `Rectangular` mesh every pixel in the masked data maps to only one pixel, thus the second + dimension of `pix_indexes_for_sub_slim_index` is always of size 1. + + The mapper allows us to create a mapping matrix, which is a matrix representing the mapping between every + unmasked data pixel annd the pixels of a mesh. This matrix is the basis of performing an `Inversion`, + which reconstructs the data using the `source_plane_mesh_grid`. + + Parameters + ---------- + mapper_grids + An object containing the data grid and mesh grid in both the data-frame and source-frame used by the + mapper to map data-points to linear object parameters. + regularization + The regularization scheme which may be applied to this linear object in order to smooth its solution, + which for a mapper smooths neighboring pixels on the mesh. + """ + + @cached_property + def pix_sub_weights(self) -> PixSubWeights: + """ + Computes the following three quantities describing the mappings between of every sub-pixel in the masked data + and pixel in the `Rectangular` mesh. + + - `pix_indexes_for_sub_slim_index`: the mapping of every data pixel (given its `sub_slim_index`) + to mesh pixels (given their `pix_indexes`). + + - `pix_sizes_for_sub_slim_index`: the number of mappings of every data pixel to mesh pixels. + + - `pix_weights_for_sub_slim_index`: the interpolation weights of every data pixel's mesh + pixel mapping + + These are packaged into the class `PixSubWeights` with attributes `mappings`, `sizes` and `weights`. + + The `sub_slim_index` refers to the masked data sub-pixels and `pix_indexes` the mesh pixel indexes, + for example: + + - `pix_indexes_for_sub_slim_index[0, 0] = 2`: The data's first (index 0) sub-pixel maps to the Rectangular + mesh's third (index 2) pixel. + + - `pix_indexes_for_sub_slim_index[2, 0] = 4`: The data's third (index 2) sub-pixel maps to the Rectangular + mesh's fifth (index 4) pixel. + + The second dimension of the array `pix_indexes_for_sub_slim_index`, which is 0 in both examples above, is used + for cases where a data pixel maps to more than one mesh pixel (for example a `Delaunay` triangulation + where each data pixel maps to 3 Delaunay triangles with interpolation weights). The weights of multiple mappings + are stored in the array `pix_weights_for_sub_slim_index`. + + For a Rectangular pixelization each data sub-pixel maps to a single mesh pixel, thus the second + 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, 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=self.source_plane_data_grid.over_sampled, + ) + ) + + return PixSubWeights( + mappings=mappings, + sizes=4 * jnp.ones(len(mappings), dtype="int"), + weights=weights, + ) \ No newline at end of file diff --git a/autoarray/inversion/pixelization/mesh/__init__.py b/autoarray/inversion/pixelization/mesh/__init__.py index a14f53f69..28f35f116 100644 --- a/autoarray/inversion/pixelization/mesh/__init__.py +++ b/autoarray/inversion/pixelization/mesh/__init__.py @@ -1,4 +1,5 @@ from .abstract import AbstractMesh as Mesh from .rectangular import Rectangular +from .rectangular_uniform import RectangularUniform from .voronoi import Voronoi from .delaunay import Delaunay diff --git a/autoarray/inversion/pixelization/mesh/rectangular.py b/autoarray/inversion/pixelization/mesh/rectangular.py index e95a05ec6..318ef0436 100644 --- a/autoarray/inversion/pixelization/mesh/rectangular.py +++ b/autoarray/inversion/pixelization/mesh/rectangular.py @@ -1,5 +1,5 @@ import numpy as np -from typing import Dict, Optional, Tuple +from typing import Optional, Tuple from autoarray.structures.grids.uniform_2d import Grid2D diff --git a/autoarray/inversion/pixelization/mesh/rectangular_uniform.py b/autoarray/inversion/pixelization/mesh/rectangular_uniform.py new file mode 100644 index 000000000..18a092610 --- /dev/null +++ b/autoarray/inversion/pixelization/mesh/rectangular_uniform.py @@ -0,0 +1,31 @@ +from autoarray.inversion.pixelization.mesh.rectangular import Rectangular + +from typing import Optional + + +from autoarray.structures.grids.uniform_2d import Grid2D +from autoarray.structures.mesh.rectangular_2d_uniform import Mesh2DRectangularUniform + +class RectangularUniform(Rectangular): + + def mesh_grid_from( + self, + source_plane_data_grid: Optional[Grid2D] = None, + source_plane_mesh_grid: Optional[Grid2D] = None, + ) -> Mesh2DRectangularUniform: + """ + Return the rectangular `source_plane_mesh_grid` as a `Mesh2DRectangular` object, which provides additional + functionality for perform operatons that exploit the geometry of a rectangular pixelization. + + Parameters + ---------- + source_plane_data_grid + The (y,x) grid of coordinates over which the rectangular pixelization is overlaid, where this grid may have + had exterior pixels relocated to its edge via the border. + source_plane_mesh_grid + Not used for a rectangular pixelization, because the pixelization grid in the `source` frame is computed + by overlaying the `source_plane_data_grid` with the rectangular pixelization. + """ + return Mesh2DRectangularUniform.overlay_grid( + shape_native=self.shape, grid=source_plane_data_grid.over_sampled + ) \ No newline at end of file diff --git a/autoarray/plot/mat_plot/two_d.py b/autoarray/plot/mat_plot/two_d.py index f9744803d..07e3172eb 100644 --- a/autoarray/plot/mat_plot/two_d.py +++ b/autoarray/plot/mat_plot/two_d.py @@ -568,64 +568,71 @@ def _plot_rectangular_mapper( if pixel_values is not None: - # self.plot_array( - # array=pixel_values, - # visuals_2d=visuals_2d, - # auto_labels=auto_labels, - # bypass=True, - # ) - - norm = self.cmap.norm_from( - array=pixel_values.array, use_log10=self.use_log10 - ) + from autoarray.inversion.pixelization.mappers.rectangular_uniform import MapperRectangularUniform + from autoarray.inversion.pixelization.mappers.rectangular import MapperRectangular - edges_transformed = mapper.edges_transformed + if isinstance(mapper, MapperRectangularUniform): - edges_transformed_dense = np.moveaxis( - np.stack(np.meshgrid(*edges_transformed.T)), 0, 2 - ) + self.plot_array( + array=pixel_values, + visuals_2d=visuals_2d, + auto_labels=auto_labels, + bypass=True, + ) - plt.pcolormesh( - edges_transformed_dense[..., 0], - edges_transformed_dense[..., 1], - pixel_values.array.reshape(shape_native), - shading="flat", - norm=norm, - cmap=self.cmap.cmap, - ) + else: - if self.colorbar is not False: + norm = self.cmap.norm_from( + array=pixel_values.array, use_log10=self.use_log10 + ) - cb = self.colorbar.set( - units=self.units, - ax=ax, + edges_transformed = mapper.edges_transformed + + edges_transformed_dense = np.moveaxis( + np.stack(np.meshgrid(*edges_transformed.T)), 0, 2 + ) + + plt.pcolormesh( + edges_transformed_dense[..., 0], + edges_transformed_dense[..., 1], + pixel_values.array.reshape(shape_native), + shading="flat", norm=norm, - cb_unit=auto_labels.cb_unit, - use_log10=self.use_log10, + cmap=self.cmap.cmap, ) - self.colorbar_tickparams.set(cb=cb) - extent_axis = self.axis.config_dict.get("extent") + if self.colorbar is not False: - if extent_axis is None: - extent_axis = extent + cb = self.colorbar.set( + units=self.units, + ax=ax, + norm=norm, + cb_unit=auto_labels.cb_unit, + use_log10=self.use_log10, + ) + self.colorbar_tickparams.set(cb=cb) - self.axis.set(extent=extent_axis) + extent_axis = self.axis.config_dict.get("extent") - self.tickparams.set() - self.yticks.set( - min_value=extent_axis[2], - max_value=extent_axis[3], - units=self.units, - pixels=shape_native[0], - ) + if extent_axis is None: + extent_axis = extent - self.xticks.set( - min_value=extent_axis[0], - max_value=extent_axis[1], - units=self.units, - pixels=shape_native[1], - ) + self.axis.set(extent=extent_axis) + + self.tickparams.set() + self.yticks.set( + min_value=extent_axis[2], + max_value=extent_axis[3], + units=self.units, + pixels=shape_native[0], + ) + + self.xticks.set( + min_value=extent_axis[0], + max_value=extent_axis[1], + units=self.units, + pixels=shape_native[1], + ) if not isinstance(self.text, list): self.text.set() diff --git a/autoarray/structures/mesh/rectangular_2d.py b/autoarray/structures/mesh/rectangular_2d.py index 35e56003b..64f51379b 100644 --- a/autoarray/structures/mesh/rectangular_2d.py +++ b/autoarray/structures/mesh/rectangular_2d.py @@ -109,7 +109,7 @@ def overlay_grid( origin=origin, ) - return Mesh2DRectangular( + return cls( values=grid_slim, shape_native=shape_native, pixel_scales=pixel_scales, diff --git a/autoarray/structures/mesh/rectangular_2d_uniform.py b/autoarray/structures/mesh/rectangular_2d_uniform.py new file mode 100644 index 000000000..231f6a6ed --- /dev/null +++ b/autoarray/structures/mesh/rectangular_2d_uniform.py @@ -0,0 +1,5 @@ +from autoarray.structures.mesh.rectangular_2d import Mesh2DRectangular + +class Mesh2DRectangularUniform(Mesh2DRectangular): + + pass \ No newline at end of file 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 e5abd0b4b..120a1da9c 100644 --- a/test_autoarray/inversion/inversion/imaging/test_inversion_imaging_util.py +++ b/test_autoarray/inversion/inversion/imaging/test_inversion_imaging_util.py @@ -183,7 +183,7 @@ def test__data_vector_via_w_tilde_data_two_methods_agree(): psf = kernel - pixelization = aa.mesh.Rectangular(shape=(20, 20)) + pixelization = aa.mesh.RectangularUniform(shape=(20, 20)) # TODO : Use pytest.parameterize @@ -266,7 +266,7 @@ def test__curvature_matrix_via_w_tilde_two_methods_agree(): psf = kernel - pixelization = aa.mesh.Rectangular(shape=(20, 20)) + pixelization = aa.mesh.RectangularUniform(shape=(20, 20)) mapper_grids = pixelization.mapper_grids_from( mask=mask, @@ -313,7 +313,7 @@ def test__curvature_matrix_via_w_tilde_preload_two_methods_agree(): psf = kernel - pixelization = aa.mesh.Rectangular(shape=(20, 20)) + pixelization = aa.mesh.RectangularUniform(shape=(20, 20)) for sub_size in range(1, 2, 3): grid = aa.Grid2D.from_mask(mask=mask, over_sample_size=sub_size) diff --git a/test_autoarray/inversion/inversion/test_abstract.py b/test_autoarray/inversion/inversion/test_abstract.py index 84338c4ae..5f1a918cf 100644 --- a/test_autoarray/inversion/inversion/test_abstract.py +++ b/test_autoarray/inversion/inversion/test_abstract.py @@ -115,8 +115,8 @@ def test__curvature_matrix__via_w_tilde__identical_to_mapping(): grid = aa.Grid2D.from_mask(mask=mask, over_sample_size=1) - mesh_0 = aa.mesh.Rectangular(shape=(3, 3)) - mesh_1 = aa.mesh.Rectangular(shape=(4, 4)) + mesh_0 = aa.mesh.RectangularUniform(shape=(3, 3)) + mesh_1 = aa.mesh.RectangularUniform(shape=(4, 4)) mapper_grids_0 = mesh_0.mapper_grids_from( mask=mask, diff --git a/test_autoarray/inversion/inversion/test_factory.py b/test_autoarray/inversion/inversion/test_factory.py index 602001429..b8cdaf31f 100644 --- a/test_autoarray/inversion/inversion/test_factory.py +++ b/test_autoarray/inversion/inversion/test_factory.py @@ -69,7 +69,7 @@ def test__inversion_imaging__via_mapper( settings=aa.SettingsInversion(use_w_tilde=False), ) - assert isinstance(inversion.linear_obj_list[0], aa.MapperRectangular) + assert isinstance(inversion.linear_obj_list[0], aa.MapperRectangularUniform) assert isinstance(inversion, aa.InversionImagingMapping) assert inversion.log_det_curvature_reg_matrix_term == pytest.approx( 7.2571757082, 1.0e-4 @@ -85,7 +85,7 @@ def test__inversion_imaging__via_mapper( settings=aa.SettingsInversion(use_w_tilde=True), ) - assert isinstance(inversion.linear_obj_list[0], aa.MapperRectangular) + assert isinstance(inversion.linear_obj_list[0], aa.MapperRectangularUniform) assert isinstance(inversion, aa.InversionImagingWTilde) assert inversion.log_det_curvature_reg_matrix_term == pytest.approx( 7.257175708246, 1.0e-4 @@ -236,7 +236,7 @@ def test__inversion_imaging__via_linear_obj_func_and_mapper( ) assert isinstance(inversion.linear_obj_list[0], aa.m.MockLinearObj) - assert isinstance(inversion.linear_obj_list[1], aa.MapperRectangular) + assert isinstance(inversion.linear_obj_list[1], aa.MapperRectangularUniform) assert isinstance(inversion, aa.InversionImagingMapping) assert inversion.log_det_curvature_reg_matrix_term == pytest.approx( 7.2571757082469945, 1.0e-4 @@ -257,7 +257,7 @@ def test__inversion_imaging__via_linear_obj_func_and_mapper( ) assert isinstance(inversion.linear_obj_list[0], aa.m.MockLinearObj) - assert isinstance(inversion.linear_obj_list[1], aa.MapperRectangular) + assert isinstance(inversion.linear_obj_list[1], aa.MapperRectangularUniform) assert isinstance(inversion, aa.InversionImagingWTilde) assert inversion.log_det_curvature_reg_matrix_term == pytest.approx( 7.2571757082469945, 1.0e-4 @@ -495,7 +495,7 @@ def test__inversion_interferometer__via_mapper( settings=aa.SettingsInversion(use_w_tilde=False), ) - assert isinstance(inversion.linear_obj_list[0], aa.MapperRectangular) + assert isinstance(inversion.linear_obj_list[0], aa.MapperRectangularUniform) assert isinstance(inversion, aa.InversionInterferometerMapping) assert inversion.mapped_reconstructed_data == pytest.approx( 1.0 + 0.0j * np.ones(shape=(7,)), 1.0e-4 diff --git a/test_autoarray/inversion/pixelization/mappers/test_factory.py b/test_autoarray/inversion/pixelization/mappers/test_factory.py index 31a24ba2a..d8f68507d 100644 --- a/test_autoarray/inversion/pixelization/mappers/test_factory.py +++ b/test_autoarray/inversion/pixelization/mappers/test_factory.py @@ -24,7 +24,7 @@ def test__rectangular_mapper(): grid.over_sampled[0, 0] = -2.0 grid.over_sampled[0, 1] = 2.0 - mesh = aa.mesh.Rectangular(shape=(3, 3)) + mesh = aa.mesh.RectangularUniform(shape=(3, 3)) mapper_grids = mesh.mapper_grids_from( mask=mask, @@ -35,7 +35,7 @@ def test__rectangular_mapper(): mapper = aa.Mapper(mapper_grids=mapper_grids, regularization=None) - assert isinstance(mapper, aa.MapperRectangular) + assert isinstance(mapper, aa.MapperRectangularUniform) assert mapper.image_plane_mesh_grid == None assert mapper.source_plane_mesh_grid.geometry.shape_native_scaled == pytest.approx( diff --git a/test_autoarray/inversion/pixelization/mappers/test_mapper_util.py b/test_autoarray/inversion/pixelization/mappers/test_mapper_util.py index 2245dda72..fdbaabe88 100644 --- a/test_autoarray/inversion/pixelization/mappers/test_mapper_util.py +++ b/test_autoarray/inversion/pixelization/mappers/test_mapper_util.py @@ -13,62 +13,6 @@ def make_five_pixels(): return np.array([[0, 0], [0, 1], [1, 0], [1, 1], [1, 2]]) -def test__sub_slim_indexes_for_pix_index(): - pix_indexes_for_sub_slim_index = np.array( - [[0, 4], [1, 4], [2, 4], [0, 4], [1, 4], [3, 4], [0, 4], [3, 4]] - ).astype("int") - pix_pixels = 5 - pix_weights_for_sub_slim_index = np.array( - [ - [0.1, 0.9], - [0.2, 0.8], - [0.3, 0.7], - [0.4, 0.6], - [0.5, 0.5], - [0.6, 0.4], - [0.7, 0.3], - [0.8, 0.2], - ] - ) - - ( - sub_slim_indexes_for_pix_index, - sub_slim_sizes_for_pix_index, - sub_slim_weights_for_pix_index, - ) = aa.util.mapper.sub_slim_indexes_for_pix_index( - pix_indexes_for_sub_slim_index=pix_indexes_for_sub_slim_index, - pix_weights_for_sub_slim_index=pix_weights_for_sub_slim_index, - pix_pixels=pix_pixels, - ) - - assert ( - sub_slim_indexes_for_pix_index - == np.array( - [ - [0, 3, 6, -1, -1, -1, -1, -1], - [1, 4, -1, -1, -1, -1, -1, -1], - [2, -1, -1, -1, -1, -1, -1, -1], - [5, 7, -1, -1, -1, -1, -1, -1], - [0, 1, 2, 3, 4, 5, 6, 7], - ] - ) - ).all() - assert (sub_slim_sizes_for_pix_index == np.array([3, 2, 1, 2, 8])).all() - - assert ( - sub_slim_weights_for_pix_index - == np.array( - [ - [0.1, 0.4, 0.7, -1, -1, -1, -1, -1], - [0.2, 0.5, -1, -1, -1, -1, -1, -1], - [0.3, -1, -1, -1, -1, -1, -1, -1], - [0.6, 0.8, -1, -1, -1, -1, -1, -1], - [0.9, 0.8, 0.7, 0.6, 0.5, 0.4, 0.3, 0.2], - ] - ) - ).all() - - def test__mapping_matrix(three_pixels, five_pixels): pix_indexes_for_sub_slim_index = np.array([[0], [1], [2]]) slim_index_for_sub_slim_index = np.array([0, 1, 2]) diff --git a/test_autoarray/inversion/pixelization/mappers/test_rectangular.py b/test_autoarray/inversion/pixelization/mappers/test_rectangular.py index eecf4b529..f80c67cde 100644 --- a/test_autoarray/inversion/pixelization/mappers/test_rectangular.py +++ b/test_autoarray/inversion/pixelization/mappers/test_rectangular.py @@ -22,7 +22,7 @@ def test__pix_indexes_for_sub_slim_index__matches_util(): over_sample_size=1, ) - mesh_grid = aa.Mesh2DRectangular.overlay_grid( + mesh_grid = aa.Mesh2DRectangularUniform.overlay_grid( shape_native=(3, 3), grid=grid.over_sampled ) @@ -47,7 +47,7 @@ def test__pix_indexes_for_sub_slim_index__matches_util(): def test__pixel_signals_from__matches_util(grid_2d_sub_1_7x7, image_7x7): - mesh_grid = aa.Mesh2DRectangular.overlay_grid( + mesh_grid = aa.Mesh2DRectangularUniform.overlay_grid( shape_native=(3, 3), grid=grid_2d_sub_1_7x7.over_sampled ) @@ -75,8 +75,6 @@ def test__pixel_signals_from__matches_util(grid_2d_sub_1_7x7, image_7x7): assert (pixel_signals == pixel_signals_util).all() - - def test__areas_transformed(mask_2d_7x7): grid = aa.Grid2DIrregular( @@ -93,7 +91,7 @@ def test__areas_transformed(mask_2d_7x7): ], ) - mesh = aa.Mesh2DRectangular.overlay_grid( + mesh = aa.Mesh2DRectangularUniform.overlay_grid( shape_native=(3, 3), grid=grid, buffer=1e-8 ) @@ -111,7 +109,6 @@ def test__areas_transformed(mask_2d_7x7): ) - def test__edges_transformed(mask_2d_7x7): grid = aa.Grid2DIrregular( @@ -128,7 +125,7 @@ def test__edges_transformed(mask_2d_7x7): ], ) - mesh = aa.Mesh2DRectangular.overlay_grid( + mesh = aa.Mesh2DRectangularUniform.overlay_grid( shape_native=(3, 3), grid=grid, buffer=1e-8 ) diff --git a/test_autoarray/structures/mesh/test_rectangular.py b/test_autoarray/structures/mesh/test_rectangular.py index a63489733..f4838ae3d 100644 --- a/test_autoarray/structures/mesh/test_rectangular.py +++ b/test_autoarray/structures/mesh/test_rectangular.py @@ -12,7 +12,7 @@ def test__neighbors__compare_to_mesh_util(): # I8 I 9I10I11I # I12I13I14I15I - mesh = aa.Mesh2DRectangular.overlay_grid( + mesh = aa.Mesh2DRectangularUniform.overlay_grid( shape_native=(7, 5), grid=np.zeros((2, 2)), buffer=1e-8 ) @@ -39,7 +39,7 @@ def test__edge_pixel_list(): ] ) - mesh = aa.Mesh2DRectangular.overlay_grid( + mesh = aa.Mesh2DRectangularUniform.overlay_grid( shape_native=(3, 3), grid=grid, buffer=1e-8 ) @@ -61,7 +61,7 @@ def test__shape_native_and_pixel_scales(): ] ) - mesh = aa.Mesh2DRectangular.overlay_grid( + mesh = aa.Mesh2DRectangularUniform.overlay_grid( shape_native=(3, 3), grid=grid, buffer=1e-8 ) @@ -82,7 +82,7 @@ def test__shape_native_and_pixel_scales(): ] ) - mesh = aa.Mesh2DRectangular.overlay_grid( + mesh = aa.Mesh2DRectangularUniform.overlay_grid( shape_native=(5, 4), grid=grid, buffer=1e-8 ) @@ -91,7 +91,7 @@ def test__shape_native_and_pixel_scales(): grid = np.array([[2.0, 1.0], [4.0, 3.0], [6.0, 5.0], [8.0, 7.0]]) - mesh = aa.Mesh2DRectangular.overlay_grid( + mesh = aa.Mesh2DRectangularUniform.overlay_grid( shape_native=(3, 3), grid=grid, buffer=1e-8 ) @@ -114,7 +114,7 @@ def test__pixel_centres__3x3_grid__pixel_centres(): ] ) - mesh = aa.Mesh2DRectangular.overlay_grid( + mesh = aa.Mesh2DRectangularUniform.overlay_grid( shape_native=(3, 3), grid=grid, buffer=1e-8 ) @@ -148,7 +148,7 @@ def test__pixel_centres__3x3_grid__pixel_centres(): ] ) - mesh = aa.Mesh2DRectangular.overlay_grid( + mesh = aa.Mesh2DRectangularUniform.overlay_grid( shape_native=(4, 3), grid=grid, buffer=1e-8 ) @@ -179,7 +179,7 @@ def test__interpolated_array_from(): pixel_scales=1.0, ) - grid_rectangular = aa.Mesh2DRectangular( + grid_rectangular = aa.Mesh2DRectangularUniform( values=grid, shape_native=grid.shape_native, pixel_scales=grid.pixel_scales ) From 52ca2ac257eb7b56e200337bab305baaa8e71d09 Mon Sep 17 00:00:00 2001 From: Jammy2211 Date: Tue, 23 Sep 2025 19:31:01 +0100 Subject: [PATCH 17/17] all unit tests pass --- .../inversion/pixelization/mappers/rectangular_uniform.py | 2 +- .../inversion/pixelization/mesh/rectangular_uniform.py | 3 ++- autoarray/plot/mat_plot/two_d.py | 8 ++++++-- autoarray/structures/mesh/rectangular_2d_uniform.py | 3 ++- 4 files changed, 11 insertions(+), 5 deletions(-) diff --git a/autoarray/inversion/pixelization/mappers/rectangular_uniform.py b/autoarray/inversion/pixelization/mappers/rectangular_uniform.py index 0f2461672..3c58813cb 100644 --- a/autoarray/inversion/pixelization/mappers/rectangular_uniform.py +++ b/autoarray/inversion/pixelization/mappers/rectangular_uniform.py @@ -103,4 +103,4 @@ def pix_sub_weights(self) -> PixSubWeights: mappings=mappings, sizes=4 * jnp.ones(len(mappings), dtype="int"), weights=weights, - ) \ No newline at end of file + ) diff --git a/autoarray/inversion/pixelization/mesh/rectangular_uniform.py b/autoarray/inversion/pixelization/mesh/rectangular_uniform.py index 18a092610..acefb8288 100644 --- a/autoarray/inversion/pixelization/mesh/rectangular_uniform.py +++ b/autoarray/inversion/pixelization/mesh/rectangular_uniform.py @@ -6,6 +6,7 @@ from autoarray.structures.grids.uniform_2d import Grid2D from autoarray.structures.mesh.rectangular_2d_uniform import Mesh2DRectangularUniform + class RectangularUniform(Rectangular): def mesh_grid_from( @@ -28,4 +29,4 @@ def mesh_grid_from( """ return Mesh2DRectangularUniform.overlay_grid( shape_native=self.shape, grid=source_plane_data_grid.over_sampled - ) \ No newline at end of file + ) diff --git a/autoarray/plot/mat_plot/two_d.py b/autoarray/plot/mat_plot/two_d.py index 07e3172eb..340366151 100644 --- a/autoarray/plot/mat_plot/two_d.py +++ b/autoarray/plot/mat_plot/two_d.py @@ -568,8 +568,12 @@ def _plot_rectangular_mapper( if pixel_values is not None: - from autoarray.inversion.pixelization.mappers.rectangular_uniform import MapperRectangularUniform - from autoarray.inversion.pixelization.mappers.rectangular import MapperRectangular + from autoarray.inversion.pixelization.mappers.rectangular_uniform import ( + MapperRectangularUniform, + ) + from autoarray.inversion.pixelization.mappers.rectangular import ( + MapperRectangular, + ) if isinstance(mapper, MapperRectangularUniform): diff --git a/autoarray/structures/mesh/rectangular_2d_uniform.py b/autoarray/structures/mesh/rectangular_2d_uniform.py index 231f6a6ed..688b7c53b 100644 --- a/autoarray/structures/mesh/rectangular_2d_uniform.py +++ b/autoarray/structures/mesh/rectangular_2d_uniform.py @@ -1,5 +1,6 @@ from autoarray.structures.mesh.rectangular_2d import Mesh2DRectangular + class Mesh2DRectangularUniform(Mesh2DRectangular): - pass \ No newline at end of file + pass