Guard apply_sparse_operator against unequal real/imag noise sigma

autoarray/dataset/interferometer/dataset.py

@@ -266,13 +266,59 @@ def apply_sparse_operator(
         use_jax
             If `True`, JAX is used to accelerate the NUFFT precision matrix computation.
 
+        Precondition
+        ------------
+        Every visibility must have equal real and imaginary noise sigma
+        (`noise_map.real == noise_map.imag`). The sparse operator's precision operator
+        `W~ = Re(F^H W F)` is built from the real-part sigma alone (see
+        `psf_precision_operator_from`, which passes `noise_map_real` to
+        `nufft_precision_operator_from`), a reduction that is exact only under that
+        equality. With unequal sigmas the sparse curvature matrix silently disagrees with
+        the dense `InversionInterferometerMapping` path, so this method raises a
+        `DatasetException` rather than returning a wrong operator.
+
         Returns
         -------
         Interferometer
             A new `Interferometer` dataset with the precomputed `InterferometerSparseOperator` attached,
             enabling efficient pixelized source reconstruction via the sparse linear algebra formalism.
+
+        Raises
+        ------
+        exc.DatasetException
+            If any visibility has unequal real and imaginary noise sigma.
         """
 
+        noise_map_real = np.asarray(self.noise_map.real)
+        noise_map_imag = np.asarray(self.noise_map.imag)
+
+        if not np.allclose(noise_map_real, noise_map_imag):
+
+            unequal = ~np.isclose(noise_map_real, noise_map_imag)
+
+            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(
+                denominator == 0.0, 1.0, denominator
+            )
+
+            raise exc.DatasetException(
+                "The sparse operator cannot be applied to this interferometer dataset because its "
+                "noise-map has unequal real and imaginary sigma.\n\n"
+                "The sparse operator's precision operator `W~ = Re(F^H W F)` is built from the "
+                "real-part noise sigma only (see `psf_precision_operator_from`, which passes "
+                "`noise_map_real` to `nufft_precision_operator_from`). That reduction is exact only "
+                "when every visibility has equal real and imaginary sigma "
+                "(`sigma_real == sigma_imag`).\n\n"
+                f"This dataset has {int(np.count_nonzero(unequal))} of {noise_map_real.size} "
+                "visibilities where the real and imaginary sigma differ (maximum relative difference "
+                f"{np.max(relative_difference):.3e}), so the sparse curvature matrix would silently "
+                "disagree with the dense path.\n\n"
+                "Either equalise the real and imaginary noise sigma of every visibility, or fit "
+                "without calling `apply_sparse_operator()` — the dense "
+                "`InversionInterferometerMapping` path handles unequal real and imaginary sigmas "
+                "correctly."
+            )
+
         if nufft_precision_operator is None:
 
             logger.info(

autoarray/inversion/inversion/interferometer/inversion_interferometer_util.py

@@ -545,6 +545,16 @@ class InterferometerSparseOperator:
       - mask / rectangle definition is fixed
       - dtype is fixed
       - batch_size is fixed
+
+    Precondition
+    ------------
+    The `nufft_precision_operator` this state is built from encodes `W~ = Re(F^H W F)` for a
+    single real-valued noise weighting, computed from the real-part noise sigma alone. It is
+    therefore exact only for datasets where every visibility has equal real and imaginary
+    noise sigma (`sigma_real == sigma_imag`); with unequal sigmas the curvature matrix
+    assembled here silently disagrees with the dense `InversionInterferometerMapping` path.
+    `Interferometer.apply_sparse_operator` enforces this precondition and raises a
+    `DatasetException` when it is violated.
     """
 
     dirty_image: np.ndarray

test_autoarray/dataset/interferometer/test_dataset.py

@@ -67,7 +67,8 @@ def test__from_fits__raise_error_dft_visibilities_limit__threads_kwarg(
 ):
     """``from_fits`` must forward ``raise_error_dft_visibilities_limit`` to the
     ``Interferometer`` constructor so callers loading large DFT-based datasets can opt out
-    of the >10,000-visibility safety check (e.g. for profiling the JAX-traceable DFT path)."""
+    of the >10,000-visibility safety check (e.g. for profiling the JAX-traceable DFT path).
+    """
     from astropy.io import fits
 
     n_visibilities = 10_001
@@ -224,6 +225,55 @@ def test__apply_sparse_operator__dft_and_nufft_dirty_image_match(mask_2d_7x7):
     )
 
 
+def test__apply_sparse_operator__unequal_real_imag_noise__raises_exception(mask_2d_7x7):
+    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_array = np.ones((n_visibilities, 2), dtype=np.float64)
+    noise_map_array[2, 1] = 2.0
+
+    noise_map = aa.VisibilitiesNoiseMap(visibilities=noise_map_array)
+    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_2d_7x7,
+        transformer_class=transformer.TransformerDFT,
+    )
+
+    with pytest.raises(aa.exc.DatasetException):
+        dataset.apply_sparse_operator(use_jax=False)
+
+
+def test__apply_sparse_operator__non_uniform_but_equal_real_imag_noise__is_applied(
+    mask_2d_7x7,
+):
+    n_visibilities = 5
+    rng = np.random.default_rng(seed=0)
+    data = aa.Visibilities(
+        visibilities=rng.normal(size=(n_visibilities, 2)).astype(np.float64)
+    )
+
+    sigma = np.array([1.0, 2.0, 0.5, 3.0, 1.5], dtype=np.float64)
+    noise_map = aa.VisibilitiesNoiseMap(visibilities=np.stack((sigma, sigma), axis=-1))
+    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_2d_7x7,
+        transformer_class=transformer.TransformerDFT,
+    ).apply_sparse_operator(use_jax=False)
+
+    assert dataset.sparse_operator is not None
+
+
 def test__different_interferometer_without_mock_objects__customize_constructor_inputs(
     mask_2d_7x7,
 ):