diff --git a/autoarray/inversion/inversion/interferometer/sparse.py b/autoarray/inversion/inversion/interferometer/sparse.py index b5faa8afb..4d5c66c93 100644 --- a/autoarray/inversion/inversion/interferometer/sparse.py +++ b/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, diff --git a/test_autoarray/inversion/inversion/interferometer/test_interferometer.py b/test_autoarray/inversion/inversion/interferometer/test_interferometer.py index 9da81b4bc..d91fb07e2 100644 --- a/test_autoarray/inversion/inversion/interferometer/test_interferometer.py +++ b/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