diff --git a/autoarray/inversion/mesh/interpolator/rectangular.py b/autoarray/inversion/mesh/interpolator/rectangular.py index ab16ce294..1c0ebd6ed 100644 --- a/autoarray/inversion/mesh/interpolator/rectangular.py +++ b/autoarray/inversion/mesh/interpolator/rectangular.py @@ -544,11 +544,23 @@ def adaptive_rectangular_mappings_weights_via_interpolation_from( grid_over_sampled_transformed = transform_func(grid_over_sampled_scaled) grid_over_index = (source_grid_size - 3) * grid_over_sampled_transformed + 1 - # --- Step 4. Floor/ceil indices --- - ix_down = xp.floor(grid_over_index[:, 0]) - ix_up = xp.ceil(grid_over_index[:, 0]) - iy_down = xp.floor(grid_over_index[:, 1]) - iy_up = xp.ceil(grid_over_index[:, 1]) + # --- Step 4. Bracketing indices --- + # + # `ix_up` is ALWAYS `ix_down + 1`, never `ceil(grid_over_index)`. + # `transform()` ends in `xp.clip(F_q, 0.0, 1.0)`, so saturated points land + # on EXACTLY integer `grid_over_index` values. Under `ceil` the bracketing + # cell collapses there (`ix_up == ix_down`) and the interpolation weight is + # forced onto a single row; a 1-ULP change in the traced grid then moves + # such a point off the plateau and jumps its weight a whole mesh row. That + # is a discontinuity in the discretisation, and it is what made the eager + # and jitted likelihoods disagree by ~1.6e-3 on an otherwise smooth + # surface (autolens_workspace_test#279). Bracketing with `ix_down + 1` is + # continuous at integer coordinates, so the cell assignment no longer + # depends on round-off. Clamp so the `+ 1` cannot leave the mesh. + ix_down = xp.clip(xp.floor(grid_over_index[:, 0]), 0, source_grid_size - 2) + iy_down = xp.clip(xp.floor(grid_over_index[:, 1]), 0, source_grid_size - 2) + ix_up = ix_down + 1 + iy_up = iy_down + 1 # --- Step 5. Four corners --- idx_tl = xp.stack([ix_up, iy_down], axis=1) @@ -570,13 +582,20 @@ def flatten(idx, n): ) # --- 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) - - 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 + # + # `t_row` / `t_col` are fractional distances measured FROM the `down` node, + # so the `up` node carries `t` and the `down` node carries `1 - t`. The row + # weights were previously the other way round (`ix_up` carried `1 - t_row`, + # `ix_down` carried `t_row`), mirroring the interpolation in the row + # direction; the column weights were, and remain, correctly paired. No + # `+ 1e-12` guard is needed now the bracket is always exactly one cell wide. + t_row = grid_over_index[:, 0] - ix_down + t_col = grid_over_index[:, 1] - iy_down + + w_tl = t_row * (1 - t_col) + w_tr = t_row * t_col + w_bl = (1 - t_row) * (1 - t_col) + w_br = (1 - t_row) * t_col weights = xp.stack([w_tl, w_tr, w_bl, w_br], axis=1) return flat_indices, weights diff --git a/test_autoarray/inversion/pixelization/interpolator/test_rectangular.py b/test_autoarray/inversion/pixelization/interpolator/test_rectangular.py index fee895705..4cd5c1f28 100644 --- a/test_autoarray/inversion/pixelization/interpolator/test_rectangular.py +++ b/test_autoarray/inversion/pixelization/interpolator/test_rectangular.py @@ -293,6 +293,111 @@ def test__mappings_sizes_weights__shapes_and_weight_normalization(): assert np.allclose(w.sum(axis=1), 1.0, atol=1e-10) +def _index_space_nodes(idx, n): + """ + Recover each mapped pixel's (row, col) position in index space. + + Inverse of the module's ``flatten(ix, iy) = (n - ix) * n + iy``. + """ + return n - idx // n, idx % n + + +def test__mappings_sizes_weights__reproduces_the_query_position(): + """ + Bilinear interpolation must be exact for linear functions, which means the + four weights have to reconstruct the query itself: + + sum_i w_i * node_i == grid_over_index + + Partition of unity alone does NOT imply this — it is satisfied by any + consistent mis-pairing of corners to weights, which is exactly how the row + weights came to be mirrored (`ix_up` carrying `1 - t_row` instead of + `t_row`) and stayed that way from 2025-09-23 to 2026-08-26. The mirroring + was smooth, so gradient checks passed; only this property catches it. + See autolens_workspace_test#279. + """ + n = 16 + data_grid, over, weights = _seeded_inputs(seed=11) + + idx, w = adaptive_rectangular_mappings_weights_via_interpolation_from( + source_grid_size=n, + data_grid=data_grid, + data_grid_over_sampled=over, + mesh_weight_map=weights, + xp=np, + ) + + # Rebuild the index-space query the function discretises internally. + mu, scale = data_grid.mean(axis=0), data_grid.std(axis=0) + transform_func, _ = create_transforms( + (data_grid - mu) / scale, mesh_pixels=n, mesh_weight_map=weights, xp=np + ) + grid_over_index = (n - 3) * transform_func((over - mu) / scale) + 1 + + node_row, node_col = _index_space_nodes(idx, n) + + assert np.allclose(w.sum(axis=1), 1.0, atol=1e-10) + assert np.allclose((w * node_row).sum(axis=1), grid_over_index[:, 0], atol=1e-10) + assert np.allclose((w * node_col).sum(axis=1), grid_over_index[:, 1], atol=1e-10) + + +def test__mappings_sizes_weights__cell_assignment_is_continuous_at_integers(): + """ + ``transform()`` ends in ``clip(F_q, 0.0, 1.0)``, so saturated queries land on + EXACTLY integer ``grid_over_index`` values — systematically, not by chance. + Bracketing those with ``ceil`` collapsed the cell (``ix_up == ix_down``), so + a 1-ULP move off the plateau jumped a point's weight a whole mesh row. That + made the likelihood depend on floating-point association, which is how the + eager and jitted evaluations came to disagree by ~1.6e-3. + + Here the property is asserted directly on the interpolated value of a linear + ramp: approaching an integer row coordinate from either side must converge + to the value AT that coordinate. + """ + n = 16 + data_grid, _, weights = _seeded_inputs(seed=12) + mu, scale = data_grid.mean(axis=0), data_grid.std(axis=0) + + transform_func, inv = create_transforms( + (data_grid - mu) / scale, mesh_pixels=n, mesh_weight_map=weights, xp=np + ) + + def interpolated_ramp(over): + idx, w = adaptive_rectangular_mappings_weights_via_interpolation_from( + source_grid_size=n, + data_grid=data_grid, + data_grid_over_sampled=over, + mesh_weight_map=weights, + xp=np, + ) + node_row, node_col = _index_space_nodes(idx, n) + # A linear function of position; bilinear interpolation reproduces it + # exactly, so any discontinuity here is a cell-assignment jump. + return (w * (3.0 * node_row - 2.0 * node_col)).sum(axis=1) + + # Sweep the row coordinate densely across the whole data range. In index + # space that spans [1, n - 2], so the sweep crosses every interior integer + # boundary; a jump at any crossing is a cell-assignment discontinuity. + # This is deliberately placement-free — the inverse transform is a knot + # lookup and cannot land a query on an integer precisely enough to probe + # one boundary directly. + lo, hi = data_grid[:, 0].min(), data_grid[:, 0].max() + sweep = np.stack( + [np.linspace(lo, hi, 20001), np.full(20001, np.median(data_grid[:, 1]))], + axis=1, + ) + + values = interpolated_ramp(sweep) + steps = np.abs(np.diff(values)) + + # The ramp must actually vary, or continuity is vacuous. + assert values.max() - values.min() > 1.0 + + # A whole-row flip moves the ramp by ~3.0 (its row coefficient); a + # continuous scheme moves by ~the sweep resolution. + assert steps.max() < 0.05 + + # --------------------------------------------------------------------------- # Areas # ---------------------------------------------------------------------------