diff --git a/autoarray/dataset/interferometer/dataset.py b/autoarray/dataset/interferometer/dataset.py index 69ac03ec1..4835c529d 100644 --- a/autoarray/dataset/interferometer/dataset.py +++ b/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( diff --git a/autoarray/inversion/inversion/interferometer/inversion_interferometer_util.py b/autoarray/inversion/inversion/interferometer/inversion_interferometer_util.py index 8313d8f41..9c4345153 100644 --- a/autoarray/inversion/inversion/interferometer/inversion_interferometer_util.py +++ b/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 diff --git a/test_autoarray/dataset/interferometer/test_dataset.py b/test_autoarray/dataset/interferometer/test_dataset.py index c8823c7df..e6562fe94 100644 --- a/test_autoarray/dataset/interferometer/test_dataset.py +++ b/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, ):