Merge pull request #315 from PyAutoLabs/feature/datacube-sparse-operator

fix(interferometer-sparse): guard against Delaunay mappers (issue #314)

Author: Jammy2211

autoarray/inversion/inversion/interferometer/sparse.py

@@ -7,6 +7,7 @@
     AbstractInversionInterferometer,
 )
 from autoarray.inversion.linear_obj.linear_obj import LinearObj
+from autoarray.inversion.mesh.mesh.delaunay import Delaunay
 from autoarray.settings import Settings
 from autoarray.inversion.mappers.abstract import Mapper
 from autoarray.structures.visibilities import Visibilities
@@ -101,6 +102,27 @@ def curvature_matrix_diag(self) -> np.ndarray:
 
         mapper = self.cls_list_from(cls=Mapper)[0]
 
+        # The interferometer sparse-operator curvature path
+        # (``InterferometerSparseOperator.curvature_matrix_via_sparse_operator_from``)
+        # has only been validated against ``Rectangular*`` meshes (single source
+        # pixel per image pixel, weight=1). When given a ``Delaunay`` mapper
+        # (three source pixels per image pixel via barycentric interpolation,
+        # weights summing to 1) the returned curvature matrix disagrees with the
+        # mapping path by ~34% Frobenius and the regularized matrix loses
+        # positive-definiteness, raising a numpy ``LinAlgError`` at the Cholesky
+        # call site in ``Inversion.log_det_curvature_reg_matrix_term``. Guard
+        # rather than silently mis-computing.
+        if isinstance(mapper.mesh, Delaunay):
+            raise NotImplementedError(
+                "Interferometer.apply_sparse_operator() is not implemented for "
+                "Delaunay-mesh pixelizations: the sparse curvature math has only "
+                "been validated against Rectangular meshes (Pmax=1, weight=1) "
+                "and is structurally wrong for barycentric-interpolated mappers "
+                "(Pmax=3). For Delaunay interferometer fits, use the plain DFT "
+                "or NUFFT path (i.e. omit the apply_sparse_operator step). "
+                "Tracking issue: https://github.com/PyAutoLabs/PyAutoArray/issues/314"
+            )
+
         return self.dataset.sparse_operator.curvature_matrix_via_sparse_operator_from(
             pix_indexes_for_sub_slim_index=mapper.pix_indexes_for_sub_slim_index,
             pix_weights_for_sub_slim_index=mapper.pix_weights_for_sub_slim_index,

test_autoarray/inversion/inversion/interferometer/test_interferometer.py

@@ -62,3 +62,73 @@ def test__fast_chi_squared(
     chi_squared = aa.util.fit.chi_squared_complex_from(chi_squared_map=chi_squared_map)
 
     assert inversion.fast_chi_squared == pytest.approx(chi_squared, 1.0e-4)
+
+
+def test__apply_sparse_operator__delaunay_mapper__raises_not_implemented():
+    """``InversionInterferometerSparse.curvature_matrix_diag`` must raise a
+    ``NotImplementedError`` rather than silently mis-computing when paired
+    with a Delaunay mapper.
+
+    The interferometer sparse-operator curvature path
+    (``InterferometerSparseOperator.curvature_matrix_via_sparse_operator_from``)
+    was only validated against ``Rectangular*`` meshes (single source pixel per
+    image pixel, weight 1). On a Delaunay mapper (three source pixels per image
+    pixel via barycentric interpolation, weights summing to 1) the returned
+    curvature matrix disagrees with the mapping path by ~34% Frobenius norm and
+    the regularized matrix loses positive-definiteness, raising a numpy
+    ``LinAlgError`` deep inside ``Inversion.log_det_curvature_reg_matrix_term``.
+    The guard at the entry point catches the mis-use early with a clear message.
+    """
+    mask = aa.Mask2D(
+        mask=[
+            [True, True, True, True, True, True, True],
+            [True, True, True, True, True, True, True],
+            [True, True, True, False, True, True, True],
+            [True, True, False, False, False, True, True],
+            [True, True, True, False, True, True, True],
+            [True, True, True, True, True, True, True],
+            [True, True, True, True, True, True, True],
+        ],
+        pixel_scales=2.0,
+    )
+
+    grid = aa.Grid2D.from_mask(mask=mask, over_sample_size=1)
+
+    mesh = aa.mesh.Delaunay(pixels=9)
+    image_mesh = aa.image_mesh.Overlay(shape=(3, 3))
+    image_mesh_grid = image_mesh.image_plane_mesh_grid_from(
+        mask=mask, adapt_data=None
+    )
+
+    interpolator = mesh.interpolator_from(
+        source_plane_data_grid=grid,
+        source_plane_mesh_grid=image_mesh_grid,
+    )
+    mapper = aa.Mapper(interpolator=interpolator)
+
+    n_visibilities = 5
+    rng = np.random.default_rng(seed=0)
+    data = aa.Visibilities(
+        visibilities=rng.normal(size=(n_visibilities, 2)).astype(np.float64)
+    )
+    noise_map = aa.VisibilitiesNoiseMap(
+        visibilities=np.ones((n_visibilities, 2), dtype=np.float64)
+    )
+    uv_wavelengths = rng.normal(size=(n_visibilities, 2)).astype(np.float64)
+
+    dataset = aa.Interferometer(
+        data=data,
+        noise_map=noise_map,
+        uv_wavelengths=uv_wavelengths,
+        real_space_mask=mask,
+        transformer_class=aa.TransformerDFT,
+    )
+    dataset_sparse = dataset.apply_sparse_operator(use_jax=False)
+
+    inversion = aa.Inversion(
+        dataset=dataset_sparse,
+        linear_obj_list=[mapper],
+    )
+
+    with pytest.raises(NotImplementedError, match=r"Delaunay"):
+        _ = inversion.curvature_matrix