fix(inversion): guard data_subtracted_dict on data=None, harden sparse noise check and data-term dtype

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)

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 = {}
 

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(

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

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
 ):

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")
 

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")
 

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,