diff --git a/autoarray/inversion/inversion/abstract.py b/autoarray/inversion/inversion/abstract.py index 68b9be306..1a7e81e3e 100644 --- a/autoarray/inversion/inversion/abstract.py +++ b/autoarray/inversion/inversion/abstract.py @@ -418,19 +418,28 @@ def zeroed_ids_to_keep(self): - The concatenated pixel blocks occupy the **final** entries of the parameter vector, with total length: - total_pixels = sum(mapper.mesh.pixels for mapper in mappers) + total_pixels = sum(mapper.params for mapper in mappers) + + `mapper.params` is the length of the mapper's source-plane mesh grid — the + real size of its block. It is deliberately *not* `mapper.mesh.pixels`: a + `Delaunay` mesh's `pixels` is a description the caller supplied and can + disagree with the grid (every workspace script passed the appended grid + length until 2026-09), and offsetting by it mis-places every mapper's ring + but the last one's. --------------------------------------------------------------------------- Zeroed pixel convention --------------------------------------------------------------------------- For each mapper: - - `mapper.mesh.zeroed_pixels` must be a 1D array of **positive, mesh-local** - pixel indices in the range `[0, mapper.mesh.pixels - 1]`. + - `mapper.zeroed_pixels` is a 1D array of **positive, mesh-local** pixel + indices in the range `[0, mapper.params - 1]`, resolved by the mesh against + `mapper.params` (`AbstractMesh.zeroed_pixels_from`). - These indices identify pixels that should be **excluded** from the linear solve (e.g. edge pixels, masked regions, or padding pixels). - Indexing is defined purely within the mapper’s own pixelization (e.g. - row-major flattening for rectangular meshes). + row-major flattening for rectangular meshes, the last `zeroed_pixels` + vertices for a Delaunay mesh with an appended edge ring). This method converts all mesh-local zeroed pixel indices into **global parameter indices**, correctly offsetting for: @@ -461,21 +470,25 @@ def zeroed_ids_to_keep(self): n_total = int(self.total_params) - pixels_per_mapper = [int(m.mesh.pixels) for m in mapper_list] + # The real size of each mapper's block: the length of its mesh grid, never + # `mesh.pixels` (see the docstring). + pixels_per_mapper = [int(m.params) for m in mapper_list] total_pixels = int(sum(pixels_per_mapper)) # Global start index of concatenated pixel block pixel_start = n_total - total_pixels + zeros_local_list = [m.zeroed_pixels for m in mapper_list] + # Total number of zeroed pixels across all mappers (Python int => static) - total_zeroed = int(sum(len(m.mesh.zeroed_pixels) for m in mapper_list)) + total_zeroed = int(sum(len(zeros_local) for zeros_local in zeros_local_list)) n_keep = int(n_total - total_zeroed) # Build global indices-to-zero across all mappers zeros_global_list = [] offset = 0 - for m, n_pix in zip(mapper_list, pixels_per_mapper): - zeros_local = self._xp.asarray(m.mesh.zeroed_pixels, dtype=self._xp.int32) + for zeros_local, n_pix in zip(zeros_local_list, pixels_per_mapper): + zeros_local = self._xp.asarray(zeros_local, dtype=self._xp.int32) zeros_global_list.append(pixel_start + offset + zeros_local) offset += n_pix diff --git a/autoarray/inversion/mappers/abstract.py b/autoarray/inversion/mappers/abstract.py index 45bbcd602..e495c2f7d 100644 --- a/autoarray/inversion/mappers/abstract.py +++ b/autoarray/inversion/mappers/abstract.py @@ -104,6 +104,20 @@ def params(self) -> int: def pixels(self) -> int: return self.params + @property + def zeroed_pixels(self) -> np.ndarray: + """ + The **positive, mesh-local** indices of this mapper's parameters that the + inversion holds at zero (the edge pixels), resolved against the mapper's real + parameter count `params` rather than anything the mesh was constructed with. + + This is the only way the inversion should read a mesh's zeroed pixels: for a + `Delaunay` mesh the zeroed ring is the last `mesh.zeroed_pixels` points of the + source-plane mesh grid, which is well defined even when `mesh.pixels` does not + match that grid's length. + """ + return self.mesh.zeroed_pixels_from(pixels=self.params) + @property def mask(self): return self.source_plane_data_grid.mask diff --git a/autoarray/inversion/mesh/mesh/abstract.py b/autoarray/inversion/mesh/mesh/abstract.py index e21f14991..c8fb79e09 100644 --- a/autoarray/inversion/mesh/mesh/abstract.py +++ b/autoarray/inversion/mesh/mesh/abstract.py @@ -25,6 +25,34 @@ class AbstractMesh: def __eq__(self, other): return self.__dict__ == other.__dict__ and self.__class__ is other.__class__ + def zeroed_pixels_from(self, pixels: int) -> np.ndarray: + """ + Return the **positive** mesh-local indices of the pixels the inversion holds at + zero, for a mapper whose parameter block is `pixels` long. + + Meshes that know their edge pixels from their own geometry (the rectangular + family) expose them as a `zeroed_pixels` index array and ignore `pixels`; a + mesh that zeroes a *count* of appended edge vertices (`Delaunay`) overrides this + to resolve that count against `pixels`. A mesh with neither zeroes nothing. + + Parameters + ---------- + pixels + The number of linear parameters the mapper built from this mesh has, i.e. + the length of its source-plane mesh grid. + + Returns + ------- + np.ndarray + 1D array of positive mesh-local pixel indices to zero (possibly empty). + """ + zeroed_pixels = getattr(self, "zeroed_pixels", None) + + if zeroed_pixels is None: + return np.array([], dtype=int) + + return np.asarray(zeroed_pixels, dtype=int) + def relocated_grid_from( self, border_relocator: BorderRelocator, source_plane_data_grid: Grid2D, xp=np ) -> Grid2D: diff --git a/autoarray/inversion/mesh/mesh/delaunay.py b/autoarray/inversion/mesh/mesh/delaunay.py index cb2ea7374..50ef895ee 100644 --- a/autoarray/inversion/mesh/mesh/delaunay.py +++ b/autoarray/inversion/mesh/mesh/delaunay.py @@ -64,52 +64,72 @@ def __init__( - prevent poorly constrained boundary vertices from absorbing flux, - reduce edge artefacts in the reconstructed source. - Zeroed pixels are always placed at the **end of the mesh parameter vector** - and are not solved for; their values are fixed to zero. Internally, the - inversion accounts for these excluded parameters when constructing and - solving the linear system. + Zeroed pixels are always the **last `zeroed_pixels` vertices of the + source-plane mesh grid** the mapper is built from — the ring + `append_with_circle_edge_points` appends to the image-plane mesh grid — + and are not solved for; their values are fixed to zero. The ring is + resolved against the grid actually handed to the mapper + (`zeroed_pixels_from`), so it does not depend on `pixels` agreeing with + that grid's length: passing the appended grid length, or the interior + count, zeroes the same vertices. Internally, the inversion accounts for + these excluded parameters when constructing and solving the linear system. Parameters ---------- pixels : int - The number of active mesh vertices (linear parameters) used to represent - the source reconstruction. + The number of interior (active) mesh vertices, i.e. the number of + points drawn by the image mesh before any edge ring is appended. The + linear parameter count of a fit is set by the mesh grid the mapper + receives, so this is a description of the mesh rather than a control; + `total_pixels` adds the zeroed ring back on. areas_factor : float, optional The barycentric area of Delaunay triangles is used to weight the regularization matrix. This factor scales these areas, allowing for tuning of the regularization strength based on triangle size. zeroed_pixels : int, optional The number of edge mesh vertices to exclude from the inversion. These - are appended to the end of the mesh and fixed to zero. + are the last `zeroed_pixels` points of the mesh grid and are fixed to zero. """ - pixels = int(pixels) + zeroed_pixels - super().__init__() - self.pixels = pixels + self.pixels = int(pixels) self.areas_factor = areas_factor - self._zeroed_pixels = zeroed_pixels + self.zeroed_pixels = int(zeroed_pixels or 0) @property - def zeroed_pixels(self): + def total_pixels(self) -> int: + """ + The interior vertex count plus the zeroed edge ring — the length the mesh grid + is expected to have once `append_with_circle_edge_points` has run. + """ + return self.pixels + self.zeroed_pixels + + def zeroed_pixels_from(self, pixels: int) -> np.ndarray: """ - Return the **positive** mesh-local pixel indices to zero for a Delaunay mesh. + Return the **positive** mesh-local indices to zero for a mesh grid of `pixels` + vertices: the last `self.zeroed_pixels` of them. - For Delaunay meshes, `self.zeroed_pixels` is interpreted as a *count* of pixels - to be zeroed at the end of the pixel block. For example: - self.pixels = 780, self.zeroed_pixels = 30 - returns indices 750..779. + `pixels` is the mapper's real parameter count (`Mapper.params`, the length of + its source-plane mesh grid), not `self.pixels`, so the ring is the appended + edge points whatever count the mesh was constructed with. For example a grid of + 780 points and `zeroed_pixels = 30` gives indices 750..779. + + Parameters + ---------- + pixels + The number of vertices in the mesh grid the mapper was built from. Returns ------- np.ndarray 1D array of positive pixel indices to zero. """ - if self._zeroed_pixels <= 0: + if self.zeroed_pixels <= 0: return np.array([], dtype=int) - start = self.pixels - self._zeroed_pixels - return np.arange(start, self.pixels, dtype=int) + pixels = int(pixels) + + return np.arange(pixels - self.zeroed_pixels, pixels, dtype=int) @property def skip_areas(self): diff --git a/test_autoarray/inversion/inversion/test_abstract.py b/test_autoarray/inversion/inversion/test_abstract.py index 3b10a6588..c47dda43a 100644 --- a/test_autoarray/inversion/inversion/test_abstract.py +++ b/test_autoarray/inversion/inversion/test_abstract.py @@ -1194,3 +1194,63 @@ def test__mappings_from__pairs_each_clump_with_its_image_regions( assert len(mapping.source_contours) == mapping.pix_indexes.shape[0] assert len(mapping.image_regions) > 0 assert len(mapping.source_centre) == 2 + + +def _delaunay_mapper(pixels, grid_length, zeroed_pixels=5): + """ + A `MockMapper` over a Delaunay mesh whose grid is `grid_length` points long, the + last `zeroed_pixels` of them the appended edge ring. `pixels` is whatever the + caller told the mesh — the interior count, or (the pre-2026-09 workspace form) + the appended grid length. + """ + return aa.m.MockMapper( + mesh=aa.mesh.Delaunay(pixels=pixels, zeroed_pixels=zeroed_pixels), + parameters=grid_length, + regularization=aa.reg.Constant(), + ) + + +def _zeroed_ids(linear_obj_list): + inversion = aa.m.MockInversion( + linear_obj_list=linear_obj_list, + settings=aa.Settings(use_positive_only_solver=True, use_edge_zeroed_pixels=True), + ) + n_total = sum(linear_obj.params for linear_obj in linear_obj_list) + return sorted(set(range(n_total)) - set(inversion.zeroed_ids_to_keep.tolist())) + + +def test__zeroed_ids_to_keep__delaunay__ring_is_the_last_vertices_of_the_grid(): + """ + One mapper over a 25-point grid (20 interior + 5 ring): the ring is zeroed whether + the mesh was told `pixels=20` or `pixels=25`, with or without linear light profile + parameters ahead of the mesh block. + """ + for pixels in (20, 25): + assert _zeroed_ids([_delaunay_mapper(pixels, 25)]) == [20, 21, 22, 23, 24] + + light = aa.m.MockLinearObj(parameters=2, regularization=None) + + assert _zeroed_ids([light, _delaunay_mapper(pixels, 25)]) == [22, 23, 24, 25, 26] + + +def test__zeroed_ids_to_keep__two_delaunay_mappers__every_ring_zeroed(): + """ + Two mappers, each over a 25-point grid: both rings ([20..24] and [45..49]) are + zeroed under both call forms. Before the fix the first mapper's block was offset by + `mesh.pixels` rather than its real size, so `pixels=25` zeroed [15..19] — five + interior vertices — and left its ring live. + """ + expected = [20, 21, 22, 23, 24, 45, 46, 47, 48, 49] + + for pixels in (20, 25): + assert _zeroed_ids([_delaunay_mapper(pixels, 25), _delaunay_mapper(pixels, 25)]) == expected + + assert _zeroed_ids([_delaunay_mapper(20, 25), _delaunay_mapper(25, 25)]) == expected + + +def test__zeroed_ids_to_keep__delaunay_without_ring__nothing_zeroed(): + mapper = aa.m.MockMapper( + mesh=aa.mesh.Delaunay(pixels=9), parameters=9, regularization=aa.reg.Constant() + ) + + assert _zeroed_ids([mapper]) == [] diff --git a/test_autoarray/inversion/pixelization/interpolator/test_delaunay_nn.py b/test_autoarray/inversion/pixelization/interpolator/test_delaunay_nn.py index 5fe83ae53..daf673cad 100644 --- a/test_autoarray/inversion/pixelization/interpolator/test_delaunay_nn.py +++ b/test_autoarray/inversion/pixelization/interpolator/test_delaunay_nn.py @@ -18,8 +18,10 @@ def test__mesh_is_public_and_selects_natural_neighbor_interpolator(): areas_factor=0.4, ) - assert mesh.pixels == 108 - np.testing.assert_array_equal(mesh.zeroed_pixels, np.arange(100, 108)) + assert mesh.pixels == 100 + assert mesh.zeroed_pixels == 8 + assert mesh.total_pixels == 108 + np.testing.assert_array_equal(mesh.zeroed_pixels_from(pixels=108), np.arange(100, 108)) assert mesh.areas_factor == 0.4 assert mesh.max_cavity_triangles == 32 assert mesh.max_neighbors == 32 diff --git a/test_autoarray/inversion/pixelization/mesh/test_delaunay.py b/test_autoarray/inversion/pixelization/mesh/test_delaunay.py new file mode 100644 index 000000000..d0237bf43 --- /dev/null +++ b/test_autoarray/inversion/pixelization/mesh/test_delaunay.py @@ -0,0 +1,120 @@ +""" +The `Delaunay` mesh's zeroed edge ring. + +The ring is the last `zeroed_pixels` vertices of the mesh grid the mapper is built +from, resolved against that grid's length (`zeroed_pixels_from`) and never against +the `pixels` the mesh was constructed with. Until 2026-09 the constructor inflated +`pixels` by `zeroed_pixels` and derived the ring from the inflated count, so a +caller passing the appended grid length (every workspace script did) produced a +mesh whose `pixels` overstated the mapper's parameter count by the ring size — a +mismatch that cancelled for one mapper and mis-placed the ring for two, and that +broke the `to_dict` round-trip (the inflated `pixels` and an index-array +`zeroed_pixels` were both written to `tracer.json`). +""" + +import numpy as np +import pytest + +import autoarray as aa +from autonerves.dictable import from_dict, to_dict + + +def test__pixels_and_zeroed_pixels__stored_as_passed__not_inflated(): + mesh = aa.mesh.Delaunay(pixels=500, zeroed_pixels=30) + + assert mesh.pixels == 500 + assert mesh.zeroed_pixels == 30 + assert mesh.total_pixels == 530 + + mesh = aa.mesh.Delaunay(pixels=9) + + assert mesh.pixels == 9 + assert mesh.zeroed_pixels == 0 + assert mesh.total_pixels == 9 + + +def test__zeroed_pixels_from__last_n_vertices_of_the_grid_whatever_pixels_says(): + """ + Both call forms in the wild — `pixels=` (correct) and + `pixels=` (the workspace form) — must zero the same 30 + vertices: the appended ring at the end of a 530-point grid. + """ + expected = np.arange(500, 530) + + mesh = aa.mesh.Delaunay(pixels=500, zeroed_pixels=30) + assert (mesh.zeroed_pixels_from(pixels=530) == expected).all() + + mesh = aa.mesh.Delaunay(pixels=530, zeroed_pixels=30) + assert (mesh.zeroed_pixels_from(pixels=530) == expected).all() + + mesh = aa.mesh.Delaunay(pixels=25, zeroed_pixels=5) + assert (mesh.zeroed_pixels_from(pixels=40) == np.arange(35, 40)).all() + + +def test__zeroed_pixels_from__no_ring__empty(): + mesh = aa.mesh.Delaunay(pixels=9) + + zeroed = mesh.zeroed_pixels_from(pixels=9) + + assert zeroed.shape == (0,) + assert zeroed.dtype == int + + mesh = aa.mesh.Delaunay(pixels=9, zeroed_pixels=None) + + assert mesh.zeroed_pixels == 0 + assert mesh.zeroed_pixels_from(pixels=9).shape == (0,) + + +def test__zeroed_pixels_from__subclasses_inherit_the_count_semantics(): + for cls in (aa.mesh.DelaunayNN, aa.mesh.KNearestNeighbor, aa.mesh.KNNBarycentric): + mesh = cls(pixels=100, zeroed_pixels=10) + + assert mesh.pixels == 100 + assert mesh.zeroed_pixels == 10 + assert (mesh.zeroed_pixels_from(pixels=110) == np.arange(100, 110)).all() + + +def test__rectangular_mesh__zeroed_pixels_from__is_its_own_edge_ring(): + """ + The rectangular family knows its edge pixels from its shape; `pixels` is ignored. + """ + mesh = aa.mesh.RectangularUniform(shape=(4, 4)) + + assert (mesh.zeroed_pixels_from(pixels=16) == mesh.zeroed_pixels).all() + assert ( + mesh.zeroed_pixels_from(pixels=16) + == np.array([0, 1, 2, 3, 4, 7, 8, 11, 12, 13, 14, 15]) + ).all() + + +def test__to_dict__round_trips__pixels_and_zeroed_pixels_are_ints(): + mesh = aa.mesh.Delaunay(pixels=500, zeroed_pixels=30, areas_factor=0.25) + + mesh_dict = to_dict(mesh) + + assert mesh_dict["arguments"]["pixels"] == 500 + assert mesh_dict["arguments"]["zeroed_pixels"] == 30 + + mesh_reloaded = from_dict(mesh_dict) + + assert mesh_reloaded == mesh + assert mesh_reloaded.pixels == 500 + assert mesh_reloaded.zeroed_pixels == 30 + assert (mesh_reloaded.zeroed_pixels_from(pixels=530) == np.arange(500, 530)).all() + + +def test__mapper__zeroed_pixels__resolved_against_params(): + """ + A mapper over a 530-point grid zeroes the last 30 vertices whether the mesh was + told `pixels=500` or `pixels=530`. + """ + for pixels in (500, 530): + mapper = aa.m.MockMapper( + mesh=aa.mesh.Delaunay(pixels=pixels, zeroed_pixels=30), parameters=530 + ) + + assert (mapper.zeroed_pixels == np.arange(500, 530)).all() + + mapper = aa.m.MockMapper(mesh=aa.mesh.Delaunay(pixels=9), parameters=9) + + assert mapper.zeroed_pixels.shape == (0,)