From f139271bb97fb7f290032f04ec0df59b629e9dda Mon Sep 17 00:00:00 2001 From: Jammy2211 Date: Wed, 30 Sep 2026 10:34:57 +0100 Subject: [PATCH] fix(inversion): guard data_subtracted_dict on data=None, harden sparse noise check and data-term dtype (#590) Corrective follow-up to Discussion #13 phase 1 (#589), from the independent review on PyAutoLens#757: - AbstractInversion.data_subtracted_dict raises a clear InversionException (naming fit.inversion_with_data) when the inversion was built without data, instead of TypeError / {mapper: None}; subplot_of_mapper and subplot_mappings catch it and skip the panel. - check_noise_map_real_imag_equal uses atol=0.0 so tiny unequal sigmas are rejected rather than waved through by the default absolute tolerance. - apply_sparse_operator computes the cached data_term from complex128-promoted copies so one-shot equals sparse_terms_from_chunks for complex64 data. Tests for all four; red-checked against the reverted source. Co-Authored-By: Claude Fable 5.1 --- autoarray/dataset/interferometer/dataset.py | 8 +++- autoarray/inversion/inversion/abstract.py | 9 ++++ .../inversion_interferometer_util.py | 8 ++-- autoarray/inversion/plot/inversion_plots.py | 6 ++- .../dataset/interferometer/test_dataset.py | 47 +++++++++++++++++++ .../interferometer/test_interferometer.py | 12 +++++ .../test_inversion_interferometer_util.py | 10 ++++ .../inversion/plot/test_inversion_plotters.py | 34 ++++++++++++++ 8 files changed, 127 insertions(+), 7 deletions(-) diff --git a/autoarray/dataset/interferometer/dataset.py b/autoarray/dataset/interferometer/dataset.py index c7155a045..e81abd65a 100644 --- a/autoarray/dataset/interferometer/dataset.py +++ b/autoarray/dataset/interferometer/dataset.py @@ -388,9 +388,13 @@ def apply_sparse_operator( # `fast_chi_squared` (term 3) nor `FitInterferometer.noise_normalization` reduces over # the visibility arrays on every likelihood call. The expressions are the ones those # two reductions use, so the cached values are bit-for-bit what they would compute. + # Promote local copies to complex128 so the one-shot `data_term` equals + # `sparse_terms_from_chunks`, which promotes every chunk to complex128. + data_c128 = np.asarray(self.data.array, dtype=np.complex128) + noise_c128 = np.asarray(self.noise_map.array, dtype=np.complex128) data_term = float( - np.sum(self.data.array.real**2.0 / self.noise_map.array.real**2.0) - + np.sum(self.data.array.imag**2.0 / self.noise_map.array.imag**2.0) + np.sum(data_c128.real**2.0 / noise_c128.real**2.0) + + np.sum(data_c128.imag**2.0 / noise_c128.imag**2.0) ) noise_normalization = float( fit_util.noise_normalization_complex_from(noise_map=self.noise_map.array) diff --git a/autoarray/inversion/inversion/abstract.py b/autoarray/inversion/inversion/abstract.py index 6241a6a16..771f999d0 100644 --- a/autoarray/inversion/inversion/abstract.py +++ b/autoarray/inversion/inversion/abstract.py @@ -7,6 +7,7 @@ from autonerves import cached_property, is_test_mode +from autoarray import exc from autoarray.dataset.imaging.dataset import Imaging from autoarray.dataset.interferometer.dataset import Interferometer from autoarray.inversion.inversion.dataset_interface import DatasetInterface @@ -878,6 +879,14 @@ def data_subtracted_dict(self) -> Dict[LinearObj, Array2D]: A dictionary of the data subtracted by the reconstructed images of combinations of all but one of the linear objects the inversion. """ + if self.data is None: + raise exc.InversionException( + "`data_subtracted_dict` cannot be computed because this inversion was built " + "without data (`data=None` -- the sparse-operator path, which uses the " + "precomputed `data_term` instead of the visibilities). `data_subtracted_dict` " + "needs the fitted data, so use `fit.inversion_with_data` (or `fit.data`) " + "instead of `fit.inversion`." + ) data_subtracted_dict = {} diff --git a/autoarray/inversion/inversion/interferometer/inversion_interferometer_util.py b/autoarray/inversion/inversion/interferometer/inversion_interferometer_util.py index 44b1ef743..096698e36 100644 --- a/autoarray/inversion/inversion/interferometer/inversion_interferometer_util.py +++ b/autoarray/inversion/inversion/interferometer/inversion_interferometer_util.py @@ -32,7 +32,9 @@ def disable_jax(): def check_noise_map_real_imag_equal(noise_map) -> None: """ Raise a `DatasetException` unless every visibility of `noise_map` has equal real and - imaginary noise sigma (`noise_map.real == noise_map.imag`, to `np.allclose`). + imaginary noise sigma (`noise_map.real == noise_map.imag`, to `np.allclose` with + `atol=0.0` and the default relative tolerance, so the check is scale-free and tiny sigmas + such as 1e-9 vs 2e-9 are not waved through by the default absolute tolerance). This is the precondition of the sparse operator: its precision operator `W~ = Re(F^H W F)` is built from the real-part sigma alone (see @@ -55,10 +57,10 @@ def check_noise_map_real_imag_equal(noise_map) -> None: noise_map_real = np.asarray(np.real(noise_map)) noise_map_imag = np.asarray(np.imag(noise_map)) - if np.allclose(noise_map_real, noise_map_imag): + if np.allclose(noise_map_real, noise_map_imag, atol=0.0): return - unequal = ~np.isclose(noise_map_real, noise_map_imag) + unequal = ~np.isclose(noise_map_real, noise_map_imag, atol=0.0) denominator = np.maximum(np.abs(noise_map_real), np.abs(noise_map_imag)) relative_difference = np.abs(noise_map_real - noise_map_imag) / np.where( diff --git a/autoarray/inversion/plot/inversion_plots.py b/autoarray/inversion/plot/inversion_plots.py index 3c7445b5a..02a3f061e 100644 --- a/autoarray/inversion/plot/inversion_plots.py +++ b/autoarray/inversion/plot/inversion_plots.py @@ -6,6 +6,8 @@ from autonerves import conf +from autoarray import exc + from autoarray.inversion.mappers.abstract import Mapper from autoarray.plot.array import plot_array from autoarray.plot.utils import ( @@ -84,7 +86,7 @@ def subplot_of_mapper( positions=positions, lines=lines, ) - except (AttributeError, KeyError): + except (AttributeError, KeyError, exc.InversionException): pass # panels 1-3: reconstructed operated data (plain, log10, + mesh grid overlay) @@ -394,7 +396,7 @@ def subplot_mappings( region_alpha=region_alpha, region_labels=region_labels, ) - except (AttributeError, KeyError): + except (AttributeError, KeyError, exc.InversionException): pass # panel 1: reconstructed operated data diff --git a/test_autoarray/dataset/interferometer/test_dataset.py b/test_autoarray/dataset/interferometer/test_dataset.py index 7fcb99d77..a7d805ffe 100644 --- a/test_autoarray/dataset/interferometer/test_dataset.py +++ b/test_autoarray/dataset/interferometer/test_dataset.py @@ -511,6 +511,53 @@ def test__apply_sparse_operator__populates_data_term_and_noise_normalization( ) +def test__apply_sparse_operator__complex64_data__data_term_is_reduced_in_complex128( + mask_2d_7x7, +): + n_visibilities = 7 + rng = np.random.default_rng(seed=1) + uv_wavelengths = rng.normal(size=(n_visibilities, 2)) * 5.0e4 + + values = np.full(n_visibilities, 10001.0 + 0.0j) + sigma = np.full(n_visibilities, 1.0 + 1.0j) + + def dataset_from(dtype): + return aa.Interferometer( + data=aa.Visibilities(visibilities=np.asarray(values, dtype=dtype)), + noise_map=aa.VisibilitiesNoiseMap( + visibilities=np.asarray(sigma, dtype=dtype) + ), + uv_wavelengths=uv_wavelengths, + real_space_mask=mask_2d_7x7, + transformer_class=transformer.TransformerDFT, + ) + + dataset_c64 = dataset_from(np.complex64) + assert dataset_c64.data.array.dtype == np.complex64 + + data_term_c64 = dataset_c64.apply_sparse_operator( + use_jax=False + ).sparse_operator.data_term + data_term_c128 = dataset_from(np.complex128).apply_sparse_operator( + use_jax=False + ).sparse_operator.data_term + + assert data_term_c128 == 700140007.0 + assert data_term_c64 == data_term_c128 + + # The dataset itself is not promoted. + assert dataset_c64.data.array.dtype == np.complex64 + + pytest.importorskip("nufftax") + + terms = aa.util.inversion_interferometer.sparse_terms_from_chunks( + [(uv_wavelengths, dataset_c64.data.array, dataset_c64.noise_map.array)], + real_space_mask=mask_2d_7x7, + ) + + assert terms.data_term == pytest.approx(data_term_c64, rel=1.0e-12) + + def test__apply_sparse_operator_from_chunks__matches_apply_sparse_operator( interferometer_7_lop, mask_2d_7x7 ): diff --git a/test_autoarray/inversion/inversion/interferometer/test_interferometer.py b/test_autoarray/inversion/inversion/interferometer/test_interferometer.py index 423384122..6044299ae 100644 --- a/test_autoarray/inversion/inversion/interferometer/test_interferometer.py +++ b/test_autoarray/inversion/inversion/interferometer/test_interferometer.py @@ -1334,6 +1334,18 @@ def test__fast_chi_squared__data_none_without_a_data_term__raises(): inversion.fast_chi_squared +def test__data_subtracted_dict__data_none__raises_pointing_to_inversion_with_data(): + dataset_sparse, mapper = _sparse_interface_setup() + + inversion = aa.Inversion( + dataset=_interface_from(dataset_sparse, data=None), + linear_obj_list=[mapper], + ) + + with pytest.raises(aa.exc.InversionException, match="inversion_with_data"): + inversion.data_subtracted_dict + + def test__fast_chi_squared__data_none__jax_matches_numpy(): pytest.importorskip("jax") diff --git a/test_autoarray/inversion/inversion/interferometer/test_inversion_interferometer_util.py b/test_autoarray/inversion/inversion/interferometer/test_inversion_interferometer_util.py index bfc4dd78f..ae8b0758d 100644 --- a/test_autoarray/inversion/inversion/interferometer/test_inversion_interferometer_util.py +++ b/test_autoarray/inversion/inversion/interferometer/test_inversion_interferometer_util.py @@ -1252,6 +1252,16 @@ def test__sparse_terms_from_chunks__accepts_visibilities_and_real_two_column_arr assert terms.noise_normalization == terms_complex.noise_normalization +def test__check_noise_map_real_imag_equal__tiny_sigmas_use_a_relative_tolerance(): + check = aa.util.inversion_interferometer.check_noise_map_real_imag_equal + + # 1e-9 vs 2e-9 differ by a factor of two; the default `atol=1e-8` would wave it through. + with pytest.raises(aa.exc.DatasetException): + check(np.array([1.0e-9 + 2.0e-9j])) + + check(np.array([1.0e-9 + 1.0e-9j])) + + def test__sparse_terms_from_chunks__unequal_real_imag_noise_in_a_later_chunk__raises(): pytest.importorskip("nufftax") diff --git a/test_autoarray/inversion/plot/test_inversion_plotters.py b/test_autoarray/inversion/plot/test_inversion_plotters.py index c43b4e128..735ea1a2d 100644 --- a/test_autoarray/inversion/plot/test_inversion_plotters.py +++ b/test_autoarray/inversion/plot/test_inversion_plotters.py @@ -1,4 +1,5 @@ import autoarray.plot as aplt +from autoarray import exc from autoarray.inversion.mappers.abstract import Mapper from autoarray.inversion.plot.inversion_plots import save_reconstruction_csv @@ -96,6 +97,39 @@ def test__inversion_subplot_of_mapper__singular_curvature_reg_matrix( assert str(Path(plot_path) / "inversion_0.png") in plot_patch.paths +def test__inversion_subplots__data_subtracted_dict_raises_inversion_exception( + rectangular_inversion_7x7_3x3, + plot_path, + plot_patch, + monkeypatch, +): + inversion = rectangular_inversion_7x7_3x3 + + def _raise(self): + raise exc.InversionException("built without data") + + monkeypatch.setattr(type(inversion), "data_subtracted_dict", property(_raise)) + + with pytest.raises(exc.InversionException): + inversion.data_subtracted_dict + + aplt.subplot_of_mapper( + inversion=inversion, + mapper_index=0, + output_path=plot_path, + output_format="png", + ) + assert str(Path(plot_path) / "inversion_0.png") in plot_patch.paths + + aplt.subplot_mappings( + inversion=inversion, + pixelization_index=0, + output_path=plot_path, + output_format="png", + ) + assert str(Path(plot_path) / "mappings_0.png") in plot_patch.paths + + def test__save_reconstruction_csv__singular_curvature_reg_matrix( rectangular_inversion_7x7_3x3, tmp_path,