From f052e38191236def8ec714e729faf29e608d1e97 Mon Sep 17 00:00:00 2001 From: Jammy2211 Date: Thu, 15 Oct 2026 19:03:20 +0100 Subject: [PATCH 1/2] fix: exclude zero-weight Delaunay mappings from adaptive pixel signals adaptive_pixel_signals_from normalised each pixel's summed adapt signal by a count that included mapping entries with zero weight. For Delaunay tables those entries are not unique. A point lying exactly on a Delaunay vertex has barycentric weights (0, 0, 1) in any triangle around that vertex. A point on an edge has a zero weight in either neighbouring triangle. scipy's find_simplex and the JAX Delaunay walk legitimately pick different triangles, so the two backends counted different phantom neighbours. That gave different pixel signals, AdaptSplit regularization weights and regularization matrices H from identical inputs (PyAutoArray#627). On the 125x125 Aris setup this was the whole 118.76 log det H gap between NumPy and JAX. That setup is a stride-2 grid-subsample mesh with 1,375 vertex-coincident data points. The log-det factorizations agree with a 50-digit reference to <=2e-8 on each H. adaptive_pixel_signals_from gains a keyword exclude_zero_weight_mappings (default False, which keeps the previous result exactly). When True, the validity mask also requires a non-zero weight, so zero-weight entries add to neither the summed signal nor pixel_counts. The set of non-zero (index, weight) pairs is unique, so the result no longer depends on which simplex was chosen. The flag is a static Python bool, so the function stays jit-safe and differentiable under JAX. Scope: Delaunay tables only. The live human chose this on the issue ("Delaunay only (Recommended by review)"). The flag is the class attribute AbstractInterpolator.exclude_zero_weight_mappings: - False by default, which covers the rectangular interpolators. - True on InterpolatorDelaunay, inherited by InterpolatorDelaunayNN. - False again on InterpolatorKNearestNeighbor and so on KNNBarycentric. Mapper.pixel_signals_from passes the attribute through. KNN is excluded because KNNBarycentric clips negative barycentric weights to exactly zero for every data point outside its 3-nearest triangle. Those zeros cover whole regions, not rare coincidences. Excluding them would change KNN pixel signals broadly and would make the count, and so the likelihood, jump wherever a clipped weight reaches zero. On the 16x16 test geometry, 48% (vertex-coincident mesh) to 65% (jittered mesh) of KNNBarycentric's mapping slots are exact zeros. Delaunay and DelaunayNN have 0-2 such slots out of about 4,200. Behaviour change: pixel signals, adaptive regularization weights, H and the Bayesian evidence change on both backends, but only for Delaunay and DelaunayNN meshes that have data points exactly on vertices or edges. These are typically Overlay or grid-subsample image meshes. KNN, KNNBarycentric and rectangular meshes are unchanged bit for bit. So are Delaunay meshes without zero-weight entries. Tests: - scipy and JAX-walk Delaunay tables give identical pixel signals on a vertex-coincident grid. Without the fix they differed by up to 0.10. - A NumPy check that one vertex point described through two triangles gives identical signals. - AdaptSplit H parity at rtol 1e-12, and NumPy and JAX log det H agreeing with each other and with a high-precision reference to <1e-6. With the Delaunay flag forced off, the H check fails. - The default flag still counts zero-weight entries, checked against hand-computed values. - KNNBarycentric mapper signals equal the default (counting) util path, on a table that has clipped zero weights. - Delaunay and DelaunayNN mappers route the exclusion. Closes #627 Co-Authored-By: Claude Opus 5.5 Co-Authored-By: Claude Fable 5.1 --- autoarray/inversion/mappers/abstract.py | 1 + autoarray/inversion/mappers/mapper_util.py | 28 +- .../inversion/mesh/interpolator/abstract.py | 5 + .../inversion/mesh/interpolator/delaunay.py | 6 + autoarray/inversion/mesh/interpolator/knn.py | 5 + .../inversion/test_log_det_sparse.py | 117 ++++++++ .../pixelization/mappers/test_mapper_util.py | 254 ++++++++++++++++++ 7 files changed, 414 insertions(+), 2 deletions(-) diff --git a/autoarray/inversion/mappers/abstract.py b/autoarray/inversion/mappers/abstract.py index e495c2f7d..f224a9bcf 100644 --- a/autoarray/inversion/mappers/abstract.py +++ b/autoarray/inversion/mappers/abstract.py @@ -512,6 +512,7 @@ def pixel_signals_from(self, signal_scale: float, xp=np) -> np.ndarray: pix_size_for_sub_slim_index=self.pix_sizes_for_sub_slim_index, slim_index_for_sub_slim_index=self.over_sampler.slim_for_sub_slim, adapt_data=adapt_data, + exclude_zero_weight_mappings=self.interpolator.exclude_zero_weight_mappings, xp=xp, ) diff --git a/autoarray/inversion/mappers/mapper_util.py b/autoarray/inversion/mappers/mapper_util.py index ae0394f85..04cc10e2d 100644 --- a/autoarray/inversion/mappers/mapper_util.py +++ b/autoarray/inversion/mappers/mapper_util.py @@ -11,6 +11,7 @@ def adaptive_pixel_signals_from( pix_size_for_sub_slim_index: np.ndarray, slim_index_for_sub_slim_index: np.ndarray, adapt_data: np.ndarray, + exclude_zero_weight_mappings: bool = False, xp=np, ) -> np.ndarray: """ @@ -19,8 +20,11 @@ def adaptive_pixel_signals_from( The pixel signals are computed as follows: - 1) Divide by the number of mappe data points in the pixel, to ensure all pixels have the same + 1) Divide by the number of mapped data points in the pixel, to ensure all pixels have the same 'relative' signal (i.e. a pixel with 10 pixels doesn't have x2 the signal of one with 5). + By default every mapping entry within a data point's size counts, including entries whose + weight is zero. With `exclude_zero_weight_mappings=True` (set by Delaunay interpolators), + zero-weight entries are not counted -- see that parameter. 2) Divided by the maximum pixel-signal, so that all signals vary between 0 and 1. This ensures that the regularization weight_list are defined identically for any data quantity or signal-to-noise_map ratio. @@ -39,6 +43,20 @@ def adaptive_pixel_signals_from( A 1D array util every pixel on the grid to a pixel on the pixelization. adapt_data The image of the galaxy which is used to compute the weigghted pixel signals. + exclude_zero_weight_mappings + If `True`, mapping entries whose weight is exactly zero are not counted towards a pixel's number of + mapped data points. Delaunay interpolators set this (via their `exclude_zero_weight_mappings` class + attribute): a data point lying exactly on a Delaunay vertex has barycentric weights (0, 0, 1) in *any* + triangle around that vertex (a point on an edge has a zero weight in either neighbouring triangle), and + scipy's `find_simplex` and the JAX Delaunay walk legitimately choose different triangles. Counting those + zero-weight duplicates made the pixel signals, and through them the adaptive regularization matrix, + differ between the NumPy and JAX backends from identical inputs (PyAutoArray#627). The set of non-zero + (index, weight) pairs is unique, so excluding zero weights makes the result independent of that choice. + + It is `False` for every other interpolator. In particular `KNNBarycentric` clips negative barycentric + weights to exactly zero for every data point outside the triangle formed by its 3 nearest vertices -- + whole regions, not rare coincidences -- so excluding them there would change pixel signals broadly + and make the count (and hence the likelihood) jump wherever a clipped weight reaches zero. """ M_sub, B = pix_indexes_for_sub_slim_index.shape @@ -50,9 +68,15 @@ def adaptive_pixel_signals_from( # 2) Build a matching “parent‐slim” index for each flattened entry: I_sub = xp.repeat(xp.arange(M_sub), B) # (M_sub*B,) - # 3) Mask out any k >= pix_size_for_sub_slim_index[i] + # 3) Mask out any k >= pix_size_for_sub_slim_index[i] and, for Delaunay tables, any entry whose + # weight is zero (see `exclude_zero_weight_mappings` above, PyAutoArray#627). The flag is a + # static Python bool, so this is jit-safe; the mask is not differentiated, so the function + # stays differentiable under JAX. valid = (xp.arange(B)[None, :] < pix_size_for_sub_slim_index[:, None]).reshape(-1) + if exclude_zero_weight_mappings: + valid = valid & (flat_weights != 0) + flat_weights = xp.where(valid, flat_weights, 0.0) flat_pixidx = xp.where( valid, flat_pixidx, pixels diff --git a/autoarray/inversion/mesh/interpolator/abstract.py b/autoarray/inversion/mesh/interpolator/abstract.py index 5859bcb07..2373a2acc 100644 --- a/autoarray/inversion/mesh/interpolator/abstract.py +++ b/autoarray/inversion/mesh/interpolator/abstract.py @@ -3,6 +3,11 @@ class AbstractInterpolator: + # Whether zero-weight mapping entries are excluded when counting the data points mapped to each pixel for + # adaptive pixel signals (`mapper_util.adaptive_pixel_signals_from`). Only Delaunay tables set this: their + # zero-weight entries are arbitrary duplicates from the triangle choice (PyAutoArray#627). + exclude_zero_weight_mappings = False + def __init__( self, mesh, diff --git a/autoarray/inversion/mesh/interpolator/delaunay.py b/autoarray/inversion/mesh/interpolator/delaunay.py index bff723f69..97b4cdaed 100644 --- a/autoarray/inversion/mesh/interpolator/delaunay.py +++ b/autoarray/inversion/mesh/interpolator/delaunay.py @@ -737,6 +737,12 @@ def splitted_sizes(self): class InterpolatorDelaunay(AbstractInterpolator): + + # A data point on a Delaunay vertex (or edge) has zero barycentric weights whose vertices depend on which + # surrounding triangle is chosen, and scipy and the JAX walk choose differently, so zero-weight entries must + # not count towards adaptive pixel signals (PyAutoArray#627). + exclude_zero_weight_mappings = True + def __init__( self, mesh, diff --git a/autoarray/inversion/mesh/interpolator/knn.py b/autoarray/inversion/mesh/interpolator/knn.py index 313bec62a..9180c798a 100644 --- a/autoarray/inversion/mesh/interpolator/knn.py +++ b/autoarray/inversion/mesh/interpolator/knn.py @@ -145,6 +145,11 @@ def kernel_interpolate_points(points, query_chunk, values, k, radius_scale): class InterpolatorKNearestNeighbor(InterpolatorDelaunay): + # Not a Delaunay table: `KNNBarycentric` clips negative barycentric weights to exactly zero for every point + # outside its 3-nearest triangle, so its zero weights are region-wide, not choice-dependent duplicates. + # Excluding them would change pixel signals broadly and add a jump in the likelihood (PyAutoArray#627). + exclude_zero_weight_mappings = False + @property def dual_areas(self): """ diff --git a/test_autoarray/inversion/inversion/test_log_det_sparse.py b/test_autoarray/inversion/inversion/test_log_det_sparse.py index 65eaf693f..9e6e8d473 100644 --- a/test_autoarray/inversion/inversion/test_log_det_sparse.py +++ b/test_autoarray/inversion/inversion/test_log_det_sparse.py @@ -409,3 +409,120 @@ def raise_splu(*args, **kwargs): assert inversion.log_det_curvature_reg_matrix_term == pytest.approx( _dense_log_det(matrix), rel=1.0e-12 ) + + +def _vertex_coincident_adapt_split_mapper(regularization, xp): + """ + A Delaunay ``AdaptSplit`` mapper whose mesh is the stride-2 subsample of a 16x16 data grid + with a smooth ``0.01 * sin(7 * grid)`` distortion, so 64 data points sit exactly on mesh + vertices. scipy's ``find_simplex`` and the JAX Delaunay walk pick different (equally valid) + triangles for some of those points (PyAutoArray#627). + """ + n = 16 + + yy, xx = np.meshgrid(np.linspace(-1, 1, n), np.linspace(-1, 1, n), indexing="ij") + + data_grid = np.stack([yy.ravel(), xx.ravel()], axis=1) + data_grid = data_grid + 0.01 * np.sin(7 * data_grid[:, ::-1]) + + mesh_grid = data_grid.reshape(n, n, 2)[::2, ::2].reshape(-1, 2) + + adapt_data = np.exp(-np.sum(data_grid**2, axis=1) / 0.3) + + interpolator = aa.mesh.Delaunay(pixels=mesh_grid.shape[0]).interpolator_from( + source_plane_data_grid=aa.Grid2D.no_mask( + values=data_grid, shape_native=(data_grid.shape[0], 1), pixel_scales=1.0 + ), + source_plane_mesh_grid=aa.Grid2DIrregular(mesh_grid), + adapt_data=aa.Array2D.no_mask( + values=adapt_data, shape_native=(n, n), pixel_scales=0.1 + ), + xp=xp, + ) + + return aa.Mapper(interpolator=interpolator, regularization=regularization, xp=xp) + + +def _high_precision_log_det(matrix): + """ + A reference log determinant: a 50-digit mpmath determinant, or, if mpmath is unavailable, a + Cholesky carried out in ``np.longdouble`` (NumPy and SciPy linalg cast that back to fp64, + so the factorization is written out here). + """ + try: + import mpmath + except ImportError: + mpmath = None + + if mpmath is not None: + with mpmath.workdps(50): + return float(mpmath.log(mpmath.det(mpmath.matrix(matrix.tolist())))) + + lower = np.array(matrix, dtype=np.longdouble) + pixels = lower.shape[0] + + for j in range(pixels): + lower[j, j] = np.sqrt(lower[j, j] - np.sum(lower[j, :j] ** 2)) + lower[j + 1 :, j] = ( + lower[j + 1 :, j] - lower[j + 1 :, :j] @ lower[j, :j] + ) / lower[j, j] + + return float(2 * np.sum(np.log(np.diag(lower)))) + + +def test__log_det_regularization_matrix_term__vertex_coincident_adapt_split__numpy_matches_jax(): + """ + PyAutoArray#627: on a mesh whose vertices coincide with data points, the NumPy and JAX + backends used to build *different* ``AdaptSplit`` regularization matrices (max difference + 0.74 on max|H| 14.6 here, a log-determinant gap of 0.097), because the zero-weight mapping + entries of vertex-coincident points depend on which triangle each backend's point location + picks, and those entries were counted in the adaptive pixel signals. + + The log-determinant routes themselves were never the cause: given the same ``H`` both agree + with a 50-digit reference. That is pinned here only at moderate coefficients, where ``H`` is + well-posed (``cond ~ 2e9``). At high coefficients (``max(coefficient)^4 >~ 1e6``) the absolute + 1e-8 diagonal jitter puts ``H``'s constant-mode eigenvalue at round-off level, so the log + determinant is round-off-determined on either backend and no parity is asserted. + """ + jax = pytest.importorskip("jax") + jax.config.update("jax_enable_x64", True) + import jax.numpy as jnp + + regularization = aa.reg.AdaptSplit( + inner_coefficient=0.4, outer_coefficient=1.7, signal_scale=0.3 + ) + + mapper_np = _vertex_coincident_adapt_split_mapper(regularization, xp=np) + mapper_jax = _vertex_coincident_adapt_split_mapper(regularization, xp=jnp) + + matrix_np = np.asarray( + regularization.regularization_matrix_from(linear_obj=mapper_np) + ) + matrix_jax = regularization.regularization_matrix_from( + linear_obj=mapper_jax, xp=jnp + ) + + np.testing.assert_allclose( + np.asarray(matrix_jax), + matrix_np, + rtol=1.0e-12, + atol=1.0e-12 * np.abs(matrix_np).max(), + ) + + inversion_np = aa.m.MockInversion( + linear_obj_list=[mapper_np], regularization_matrix=matrix_np + ) + inversion_jax = aa.m.MockInversion( + linear_obj_list=[mapper_jax], regularization_matrix=matrix_jax + ) + inversion_jax.use_jax = True + + log_det_np = float(inversion_np.log_det_regularization_matrix_term) + log_det_jax = float(inversion_jax.log_det_regularization_matrix_term) + + assert abs(log_det_np - log_det_jax) < 1.0e-6 + + reference = _high_precision_log_det(matrix_np) + + assert abs(log_det_np - reference) < 1.0e-6 + assert abs(log_det_jax - reference) < 1.0e-6 diff --git a/test_autoarray/inversion/pixelization/mappers/test_mapper_util.py b/test_autoarray/inversion/pixelization/mappers/test_mapper_util.py index 9ce921a90..fc6839fb7 100644 --- a/test_autoarray/inversion/pixelization/mappers/test_mapper_util.py +++ b/test_autoarray/inversion/pixelization/mappers/test_mapper_util.py @@ -656,3 +656,257 @@ def test__mapping_matrix_from__no_mixed_precision__float64_on_both_backends(): assert jax_matrix.dtype == jnp.float64 assert np.asarray(jax_matrix) == pytest.approx(numpy_matrix, abs=1.0e-14) + + +# ---------------------------------------------------------------------------- +# Zero-weight mapping entries are not counted (PyAutoArray#627). +# +# A data point lying exactly on a Delaunay vertex has barycentric weights +# (0, 0, 1) in *every* triangle around that vertex, and a point on an edge has +# a zero weight in either neighbouring triangle. scipy's `find_simplex` and the +# JAX Delaunay walk legitimately choose different triangles for such points, so +# the two backends' mapping tables carry different zero-weight entries. +# `adaptive_pixel_signals_from` used to count those entries when normalising +# each pixel's signal, which made the pixel signals -- and through them the +# AdaptSplit regularization weights and matrix -- depend on the triangle chosen. +# Only non-zero-weight entries are counted now, and the set of non-zero +# (index, weight) pairs is unique. +# ---------------------------------------------------------------------------- + + +def _vertex_coincident_delaunay_setup(): + """ + A 16x16 data grid with a smooth `0.01 * sin(7 * grid)` distortion, and a mesh that is + the stride-2 subsample of the same points, so 64 data points sit exactly on mesh vertices + (the geometry of a grid-subsample / `Overlay` image mesh traced through a lens). + """ + n = 16 + + yy, xx = np.meshgrid(np.linspace(-1, 1, n), np.linspace(-1, 1, n), indexing="ij") + + data_grid = np.stack([yy.ravel(), xx.ravel()], axis=1) + data_grid = data_grid + 0.01 * np.sin(7 * data_grid[:, ::-1]) + + mesh_grid = data_grid.reshape(n, n, 2)[::2, ::2].reshape(-1, 2) + + adapt_data = np.exp(-np.sum(data_grid**2, axis=1) / 0.3) + + return data_grid, mesh_grid, adapt_data + + +def _delaunay_tables_from(delaunay_fn, data_grid, mesh_grid, xp): + from autoarray.inversion.mesh.interpolator.delaunay import ( + pixel_weights_delaunay_from, + ) + + mappings = np.asarray( + delaunay_fn(xp.asarray(mesh_grid), xp.asarray(data_grid), 0.5)[2] + ) + + sizes = (mappings >= 0).sum(axis=1) + + weights = np.asarray( + pixel_weights_delaunay_from( + data_grid=xp.asarray(data_grid), + mesh_grid=xp.asarray(mesh_grid), + pix_indexes_for_sub_slim_index=xp.asarray(mappings), + xp=xp, + ) + ) + + return mappings, sizes, weights + + +@requires_jax +def test__adaptive_pixel_signals_from__scipy_and_jax_delaunay_tables__agree(): + import jax + + jax.config.update("jax_enable_x64", True) + import jax.numpy as jnp + + from autoarray.inversion.mesh.interpolator.delaunay import ( + jax_delaunay, + scipy_delaunay, + ) + + data_grid, mesh_grid, adapt_data = _vertex_coincident_delaunay_setup() + + scipy_tables = _delaunay_tables_from(scipy_delaunay, data_grid, mesh_grid, xp=np) + jax_tables = _delaunay_tables_from(jax_delaunay, data_grid, mesh_grid, xp=jnp) + + # The precondition the test exists for: the two point-location routines pick different + # (equally valid) triangles for some of the vertex-coincident data points. + rows_differ = np.any( + np.sort(scipy_tables[0], axis=1) != np.sort(jax_tables[0], axis=1), axis=1 + ) + assert rows_differ.sum() > 0 + + def pixel_signals_from(tables): + mappings, sizes, weights = tables + + return aa.util.mapper.adaptive_pixel_signals_from( + pixels=mesh_grid.shape[0], + pixel_weights=weights, + signal_scale=1.0, + pix_indexes_for_sub_slim_index=mappings, + pix_size_for_sub_slim_index=sizes, + slim_index_for_sub_slim_index=np.arange(data_grid.shape[0]), + adapt_data=adapt_data, + exclude_zero_weight_mappings=True, + ) + + np.testing.assert_allclose( + pixel_signals_from(scipy_tables), + pixel_signals_from(jax_tables), + rtol=0.0, + atol=1.0e-12, + ) + + +def test__adaptive_pixel_signals_from__vertex_point_via_two_triangles__identical(): + """ + Pure NumPy: data point 0 lies exactly on vertex 0. Described through triangle (0, 1, 2) + or through triangle (0, 3, 4) its weights are (1, 0, 0); the other two data points are + ordinary interior points. Only the zero-weight neighbours differ, so the signals must not. + """ + kwargs = dict( + pixels=5, + signal_scale=1.0, + pix_size_for_sub_slim_index=np.array([3, 3, 3]), + slim_index_for_sub_slim_index=np.array([0, 1, 2]), + adapt_data=np.array([3.0, 1.0, 2.0]), + exclude_zero_weight_mappings=True, + ) + + weights = np.array([[1.0, 0.0, 0.0], [0.2, 0.3, 0.5], [0.6, 0.1, 0.3]]) + + signals_via_012 = aa.util.mapper.adaptive_pixel_signals_from( + pix_indexes_for_sub_slim_index=np.array([[0, 1, 2], [0, 1, 3], [0, 3, 4]]), + pixel_weights=weights, + **kwargs, + ) + signals_via_034 = aa.util.mapper.adaptive_pixel_signals_from( + pix_indexes_for_sub_slim_index=np.array([[0, 3, 4], [0, 1, 3], [0, 3, 4]]), + pixel_weights=weights, + **kwargs, + ) + + assert np.array_equal(signals_via_012, signals_via_034) + + # Pixel 0 is mapped by all three points with non-zero weight; the zero-weight entries of + # point 0 add nothing to pixels 1 and 2 (or 3 and 4) and are not counted there either. + pixel_sums = np.array([3.0 * 1.0 + 0.2 + 2.0 * 0.6, 0.3, 0.0, 0.5 + 0.2, 0.6]) + pixel_counts = np.array([3.0, 1.0, 1.0, 2.0, 1.0]) + expected = pixel_sums / pixel_counts + + assert signals_via_012 == pytest.approx(expected / expected.max(), rel=1.0e-14) + + +def test__adaptive_pixel_signals_from__default_flag__zero_weight_entries_still_counted(): + """ + Without `exclude_zero_weight_mappings` the util keeps its previous behaviour: every entry within + a data point's size counts towards its pixel, including zero-weight ones. Same table as the + two-triangle test above, described through triangle (0, 1, 2). + """ + signals = aa.util.mapper.adaptive_pixel_signals_from( + pixels=5, + signal_scale=1.0, + pixel_weights=np.array([[1.0, 0.0, 0.0], [0.2, 0.3, 0.5], [0.6, 0.1, 0.3]]), + pix_indexes_for_sub_slim_index=np.array([[0, 1, 2], [0, 1, 3], [0, 3, 4]]), + pix_size_for_sub_slim_index=np.array([3, 3, 3]), + slim_index_for_sub_slim_index=np.array([0, 1, 2]), + adapt_data=np.array([3.0, 1.0, 2.0]), + ) + + # Point 0's zero-weight entries on pixels 1 and 2 add nothing to the sums but are counted. + pixel_sums = np.array([3.0 * 1.0 + 0.2 + 2.0 * 0.6, 0.3, 0.0, 0.5 + 0.2, 0.6]) + pixel_counts = np.array([3.0, 2.0, 1.0, 2.0, 1.0]) + expected = pixel_sums / pixel_counts + + assert signals == pytest.approx(expected / expected.max(), rel=1.0e-14) + + +def _mapper_from(mesh, data_grid, mesh_grid, adapt_data): + interpolator = mesh.interpolator_from( + source_plane_data_grid=aa.Grid2D.no_mask( + values=data_grid, shape_native=(data_grid.shape[0], 1), pixel_scales=1.0 + ), + source_plane_mesh_grid=aa.Grid2DIrregular(mesh_grid), + adapt_data=aa.Array2D.no_mask( + values=adapt_data, shape_native=(adapt_data.shape[0], 1), pixel_scales=1.0 + ), + ) + + return aa.Mapper(interpolator=interpolator) + + +def _util_pixel_signals_from(mapper, adapt_data, exclude_zero_weight_mappings): + return aa.util.mapper.adaptive_pixel_signals_from( + pixels=mapper.pixels, + signal_scale=1.0, + pixel_weights=mapper.pix_weights_for_sub_slim_index, + pix_indexes_for_sub_slim_index=mapper.pix_indexes_for_sub_slim_index, + pix_size_for_sub_slim_index=mapper.pix_sizes_for_sub_slim_index, + slim_index_for_sub_slim_index=mapper.over_sampler.slim_for_sub_slim, + adapt_data=adapt_data, + exclude_zero_weight_mappings=exclude_zero_weight_mappings, + ) + + +def test__mapper_pixel_signals_from__knn_barycentric_keeps_counting_clipped_zero_weights(): + """ + `KNNBarycentric` clips negative barycentric weights to exactly zero for every data point outside + the triangle of its 3 nearest vertices, so its tables carry zero weights over whole regions. Those + entries keep being counted (PyAutoArray#627 scopes the exclusion to Delaunay tables), so its pixel + signals are exactly those of the default util path. + """ + data_grid, mesh_grid, adapt_data = _vertex_coincident_delaunay_setup() + mesh_grid = mesh_grid + np.random.default_rng(1).normal(0.0, 0.03, mesh_grid.shape) + + mapper = _mapper_from( + aa.mesh.KNNBarycentric(pixels=mesh_grid.shape[0]), + data_grid, + mesh_grid, + adapt_data, + ) + + weights = np.asarray(mapper.pix_weights_for_sub_slim_index) + in_size = ( + np.arange(weights.shape[1])[None, :] + < np.asarray(mapper.pix_sizes_for_sub_slim_index)[:, None] + ) + assert ((weights == 0.0) & in_size).sum() > 0 + + assert mapper.interpolator.exclude_zero_weight_mappings is False + + signals = mapper.pixel_signals_from(signal_scale=1.0) + + assert np.array_equal( + signals, + _util_pixel_signals_from( + mapper, adapt_data, exclude_zero_weight_mappings=False + ), + ) + assert not np.allclose( + signals, + _util_pixel_signals_from(mapper, adapt_data, exclude_zero_weight_mappings=True), + ) + + +def test__mapper_pixel_signals_from__delaunay_excludes_zero_weights(): + data_grid, mesh_grid, adapt_data = _vertex_coincident_delaunay_setup() + + for mesh in [ + aa.mesh.Delaunay(pixels=mesh_grid.shape[0]), + aa.mesh.DelaunayNN(pixels=mesh_grid.shape[0]), + ]: + mapper = _mapper_from(mesh, data_grid, mesh_grid, adapt_data) + + assert mapper.interpolator.exclude_zero_weight_mappings is True + + assert np.array_equal( + mapper.pixel_signals_from(signal_scale=1.0), + _util_pixel_signals_from( + mapper, adapt_data, exclude_zero_weight_mappings=True + ), + ) From da5f6cef8fdf898cdda22cb2e81db11a7fb09a44 Mon Sep 17 00:00:00 2001 From: Jammy2211 Date: Thu, 15 Oct 2026 20:48:03 +0100 Subject: [PATCH 2/2] test: skip the KNNBarycentric pixel-signal test on the no-JAX leg The KNN interpolator computes its weights with JAX, so the test added for PyAutoArray#627 failed with ModuleNotFoundError on the unittest-nojax leg (run 38795253065). importorskip("jax") matches the file's other JAX tests. Co-Authored-By: Claude Fable 5.1 --- .../inversion/pixelization/mappers/test_mapper_util.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/test_autoarray/inversion/pixelization/mappers/test_mapper_util.py b/test_autoarray/inversion/pixelization/mappers/test_mapper_util.py index fc6839fb7..9eae96393 100644 --- a/test_autoarray/inversion/pixelization/mappers/test_mapper_util.py +++ b/test_autoarray/inversion/pixelization/mappers/test_mapper_util.py @@ -859,7 +859,11 @@ def test__mapper_pixel_signals_from__knn_barycentric_keeps_counting_clipped_zero the triangle of its 3 nearest vertices, so its tables carry zero weights over whole regions. Those entries keep being counted (PyAutoArray#627 scopes the exclusion to Delaunay tables), so its pixel signals are exactly those of the default util path. + + The KNN interpolator computes its weights with JAX, so this test skips on the no-JAX leg. """ + pytest.importorskip("jax") + data_grid, mesh_grid, adapt_data = _vertex_coincident_delaunay_setup() mesh_grid = mesh_grid + np.random.default_rng(1).normal(0.0, 0.03, mesh_grid.shape)