Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
29 changes: 21 additions & 8 deletions autoarray/inversion/mesh/mesh_geometry/rectangular.py
Original file line number Diff line number Diff line change
Expand Up @@ -510,16 +510,29 @@ def areas_for_magnification(self):
@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 `rectangular_neighbors_from`.
The source-plane cell edges of every mesh pixel, transformed through
the adaptive CDF — one more edge than pixels per axis, shape
``(n + 1, 2)`` packing (y_edges, x_edges).

The interpolation mapper is node-based: reconstruction value
``(row r, col c)`` lives at node ``U_y = (n_y - r - 1)/(n_y - 3)``,
``U_x = (c - 1)/(n_x - 3)`` in unit CDF space (a guard ring occupies
the border; rows are flipped by the mapper's ``row = n - i``). The
cell edges are therefore the node midpoints — a uniform ``[0, 1]``
partition drew every cell up to ~1.5 mesh pixels away from where the
mapper scatters its flux (issue #372). Guard-node edges fall outside
the CDF domain and clamp to the data span in the inverse interp,
squashing the border guard cells to the region they actually cover.
"""

# edges defined in 0 -> 1 space, there is one more edge than pixel centers on each side
edges_y = self._xp.linspace(1, 0, self.shape_native[0] + 1)
edges_x = self._xp.linspace(0, 1, self.shape_native[1] + 1)
# Node-midpoint edges in unit CDF space (see docstring); the guard
# denominator keeps the degenerate n <= 3 meshes finite (their
# interpolator maps every point to the single interior node anyway).
n_y, n_x = self.shape_native
rows = self._xp.arange(n_y + 1)
cols = self._xp.arange(n_x + 1)
edges_y = (n_y - rows - 0.5) / max(n_y - 3, 1)
edges_x = (cols - 1.5) / max(n_x - 3, 1)

edges_reshaped = self._xp.stack([edges_y, edges_x]).T

Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -200,3 +200,53 @@ def test__edges_transformed(mask_2d_7x7):
),
abs=1e-8,
)


def test__edges_transformed__aligned_with_interpolation_node_convention():
"""
Regression test for issue #372: the pcolormesh plotting path draws value
(row r, col c) inside the cell bounded by edges_transformed — that cell
must be centred on the interpolation mapper's node for that value, not on
a uniform [0, 1] partition (which shifted plots by ~1.5 mesh pixels).

A delta function scattered through the mapper must land in plotted cells
whose weighted centroid matches the input point to sub-cell precision.
"""
from autoarray.inversion.mesh.interpolator.rectangular import (
adaptive_rectangular_mappings_weights_via_interpolation_from,
adaptive_rectangular_transformed_grid_from,
)

n = 10
rng = np.random.default_rng(0)
data_grid = rng.uniform(-1.0, 1.0, (5000, 2))
test_point = np.array([[0.3, -0.2]])

flat_indices, weights = (
adaptive_rectangular_mappings_weights_via_interpolation_from(
source_grid_size=n,
data_grid=data_grid,
data_grid_over_sampled=test_point,
)
)

# The node-midpoint unit-space edges (what edges_transformed now builds),
# pushed through the same CDF transform.
rows = np.arange(n + 1)
edges_y = (n - rows - 0.5) / (n - 3)
edges_x = (rows - 1.5) / (n - 3)
edges = np.stack([edges_y, edges_x]).T
edges_t = adaptive_rectangular_transformed_grid_from(data_grid, edges)
y_edges, x_edges = edges_t.T

centroid_y = 0.0
centroid_x = 0.0
for flat, weight in zip(flat_indices[0], weights[0]):
r, c = flat // n, flat % n
centroid_y += weight * 0.5 * (y_edges[r] + y_edges[r + 1])
centroid_x += weight * 0.5 * (x_edges[c] + x_edges[c + 1])

# Half a mesh cell in these units is ~0.15; the pre-fix uniform edges
# missed by ~0.4 in y.
assert centroid_y == pytest.approx(test_point[0, 0], abs=0.1)
assert centroid_x == pytest.approx(test_point[0, 1], abs=0.1)
Loading