diff --git a/autoarray/inversion/inversion/factory.py b/autoarray/inversion/inversion/factory.py index f6fd4555e..addb1f643 100644 --- a/autoarray/inversion/inversion/factory.py +++ b/autoarray/inversion/inversion/factory.py @@ -266,6 +266,12 @@ def _use_interferometer_numba( FFT route is faster (see `Settings.interferometer_numba_nnz_per_source_max`). - `import numba` succeeds. + NaN mapping weights (an overflowing or degenerate `DelaunayNN` Sibson row) are + deliberately not a condition: the kernel propagates them into `F` exactly as the FFT + and mapping routes do, so every route returns the same NaN evidence and the sample is + rejected identically (pinned by the NaN-weight tests in + `test_interferometer_numba.py` and `test_factory.py`). + Parameters ---------- linear_obj_list diff --git a/autoarray/inversion/inversion/interferometer_numba/sparse.py b/autoarray/inversion/inversion/interferometer_numba/sparse.py index 58cc8b4be..cc5364a12 100644 --- a/autoarray/inversion/inversion/interferometer_numba/sparse.py +++ b/autoarray/inversion/inversion/interferometer_numba/sparse.py @@ -101,6 +101,9 @@ def _check_preconditions(self) -> None: These are raised, not worked around: each one would otherwise compute a different (and silently wrong) `F`, or none at all. The factory checks the same conditions up front and simply does not route here when one fails. + + NaN mapping weights are not a precondition: they propagate into `F` and give the + same NaN evidence as the FFT route, which is the intended rejection of the sample. """ if self.has(cls=AbstractLinearObjFuncList): raise exc.InversionException( diff --git a/test_autoarray/inversion/inversion/interferometer_numba/test_interferometer_numba.py b/test_autoarray/inversion/inversion/interferometer_numba/test_interferometer_numba.py index def6e4636..c341b5a5c 100644 --- a/test_autoarray/inversion/inversion/interferometer_numba/test_interferometer_numba.py +++ b/test_autoarray/inversion/inversion/interferometer_numba/test_interferometer_numba.py @@ -103,6 +103,38 @@ def _rectangular_mapper_from(mask, shape=(3, 3), over_sample_size=1): ) +def _delaunay_nn_mapper_from(mask, pixels=25, max_neighbors=None, seed=3): + """ + A `DelaunayNN` (Sibson) mapper on a seeded irregular mesh, so most masked pixels fall + strictly inside a triangle and carry wide natural-neighbour rows (5-8 vertices here, + against at most 3 for `Delaunay`). An Overlay mesh on this small mask would put the + mesh vertices on the pixel centres, giving width-1 rows that test nothing. + + `max_neighbors` lowers the mesh's class-level `SIBSON_MAX_NEIGHBORS` cap on the + instance, which makes the wider rows overflow and carry NaN weights. + """ + grid = aa.Grid2D.from_mask(mask=mask, over_sample_size=1) + + mesh = aa.mesh.DelaunayNN(pixels=pixels) + + if max_neighbors is not None: + mesh.max_neighbors = max_neighbors + + rng = np.random.default_rng(seed=seed) + + interpolator = mesh.interpolator_from( + source_plane_data_grid=grid, + source_plane_mesh_grid=aa.Grid2DIrregular( + rng.uniform(-4.0, 4.0, size=(pixels, 2)) + ), + ) + + return aa.Mapper( + interpolator=interpolator, + regularization=aa.reg.Constant(coefficient=1.0), + ) + + def _assert_numba_matches_sparse(dataset_sparse, mapper): """ Pin the numba `direct_conv` inversion against the NumPy FFT (`xp=np`) inversion, which @@ -172,6 +204,22 @@ def test__numba_inversion__rectangular__matches_sparse_numpy_inversion(): _assert_numba_matches_sparse(dataset_sparse=dataset_sparse, mapper=mapper) +def test__numba_inversion__delaunay_nn__matches_sparse_numpy_inversion(): + """ + Sibson rows are variable-width (up to `SIBSON_MAX_NEIGHBORS`), unlike the fixed + 3-wide Delaunay rows; `kernel_inputs_from` must drop each row's padding entries. + """ + mask, _, dataset_sparse = _dataset_from() + mapper = _delaunay_nn_mapper_from(mask=mask) + + pix_sizes = np.asarray(mapper.pix_sizes_for_sub_slim_index) + + assert pix_sizes.max() > 3 + assert len(np.unique(pix_sizes)) > 1 + + _assert_numba_matches_sparse(dataset_sparse=dataset_sparse, mapper=mapper) + + def test__numba_inversion__control__one_percent_scale_of_curvature_matrix_fails_the_pin(): """ The parity pin above is only meaningful if it can fail: a 1% rescaling of `F` — far @@ -334,6 +382,72 @@ def test__factory__routes_to_numba_below_the_gate_and_sparse_above_it(): assert not isinstance(inversion, InversionInterferometerSparseNumba) +def test__factory__delaunay_nn__routes_to_numba_with_default_settings(): + """ + Sibson rows carry more non-zeros than Delaunay ones, so the risk is that a + `DelaunayNN` mapper trips the nnz-per-source-column gate onto the FFT route. With the + packaged gate it must not, and the factory-built inversion matches the FFT one. + """ + mask, _, dataset_sparse = _dataset_from() + mapper = _delaunay_nn_mapper_from(mask=mask) + + assert numba_util.nnz_per_source_column_from( + mapper=mapper + ) <= aa.Settings().interferometer_numba_nnz_per_source_max + + inversion = aa.Inversion(dataset=dataset_sparse, linear_obj_list=[mapper]) + + assert isinstance(inversion, InversionInterferometerSparseNumba) + + inversion_sparse = InversionInterferometerSparse( + dataset=dataset_sparse, linear_obj_list=[mapper], xp=np + ) + + assert inversion.log_det_curvature_reg_matrix_term == pytest.approx( + inversion_sparse.log_det_curvature_reg_matrix_term, rel=1.0e-10 + ) + assert float(inversion.fast_chi_squared) == pytest.approx( + float(inversion_sparse.fast_chi_squared), rel=1.0e-10 + ) + + +def test__numba_inversion__delaunay_nn_nan_weights__give_nan_evidence_like_sparse(): + """ + A Sibson row that overflows `max_neighbors` (or is degenerate) carries NaN weights, + so the model sample is rejected. NaN weights are deliberately not a routing condition: + the FFT and mapping routes turn them into a NaN evidence, and the numba route must do + the same rather than return a finite number -- so the factory still routes here. + """ + mask, dataset, dataset_sparse = _dataset_from() + mapper = _delaunay_nn_mapper_from(mask=mask, max_neighbors=6) + + nan_rows = np.isnan(np.asarray(mapper.pix_weights_for_sub_slim_index)).any(axis=1) + + assert nan_rows.any() and not nan_rows.all() + + def log_evidence_terms_from(inversion): + return ( + float(inversion.fast_chi_squared) + + float(inversion.regularization_term) + + float(inversion.log_det_curvature_reg_matrix_term) + - float(inversion.log_det_regularization_matrix_term) + ) + + inversion_numba = aa.Inversion(dataset=dataset_sparse, linear_obj_list=[mapper]) + + assert isinstance(inversion_numba, InversionInterferometerSparseNumba) + + for inversion in ( + inversion_numba, + InversionInterferometerSparse( + dataset=dataset_sparse, linear_obj_list=[mapper], xp=np + ), + aa.Inversion(dataset=dataset, linear_obj_list=[mapper]), + ): + assert np.isnan(np.asarray(inversion.curvature_matrix)).any() + assert np.isnan(log_evidence_terms_from(inversion)) + + def test__factory__gate_of_zero_disables_the_numba_path(): mask, _, dataset_sparse = _dataset_from() mapper = _delaunay_mapper_from(mask=mask) diff --git a/test_autoarray/inversion/inversion/test_factory.py b/test_autoarray/inversion/inversion/test_factory.py index 870fac628..706ccbcba 100644 --- a/test_autoarray/inversion/inversion/test_factory.py +++ b/test_autoarray/inversion/inversion/test_factory.py @@ -5,6 +5,10 @@ import autoarray as aa +from autoarray.inversion.inversion.imaging_numba.sparse import ( + InversionImagingSparseNumba, +) + def test__inversion_imaging__via_linear_obj_func_list(masked_imaging_7x7_no_blur): mask = masked_imaging_7x7_no_blur.mask @@ -351,6 +355,107 @@ def test__inversion_imaging__compare_mapping_and_sparse_operator_values( ) +def _delaunay_nn_mapper_from(mask, pixels=25, max_neighbors=None, seed=3): + """ + A `DelaunayNN` (Sibson) mapper on a seeded irregular mesh, so the masked pixels carry + wide, variable-width natural-neighbour rows rather than the 3-wide Delaunay ones. + `max_neighbors` lowers the mesh's class-level cap on the instance so rows overflow + and carry NaN weights. + """ + mesh = aa.mesh.DelaunayNN(pixels=pixels) + + if max_neighbors is not None: + mesh.max_neighbors = max_neighbors + + rng = np.random.default_rng(seed=seed) + + interpolator = mesh.interpolator_from( + source_plane_data_grid=aa.Grid2D.from_mask(mask=mask, over_sample_size=1), + source_plane_mesh_grid=aa.Grid2DIrregular( + rng.uniform(-4.0, 4.0, size=(pixels, 2)) + ), + ) + + return aa.Mapper( + interpolator=interpolator, + regularization=aa.reg.Constant(coefficient=1.0), + ) + + +def test__inversion_imaging__delaunay_nn__numba_sparse_operator_matches_mapping( + masked_imaging_7x7, +): + pytest.importorskip("numba") + + mapper = _delaunay_nn_mapper_from(mask=masked_imaging_7x7.mask) + + pix_sizes = np.asarray(mapper.pix_sizes_for_sub_slim_index) + + assert pix_sizes.max() > 3 + assert len(np.unique(pix_sizes)) > 1 + + inversion_sparse_operator = aa.Inversion( + dataset=masked_imaging_7x7.apply_sparse_operator_cpu(), + linear_obj_list=[mapper], + ) + inversion_mapping = aa.Inversion( + dataset=masked_imaging_7x7, + linear_obj_list=[mapper], + ) + + assert isinstance(inversion_sparse_operator, InversionImagingSparseNumba) + assert isinstance(inversion_mapping, aa.InversionImagingMapping) + + np.testing.assert_allclose( + np.asarray(inversion_sparse_operator.curvature_matrix), + np.asarray(inversion_mapping.curvature_matrix), + rtol=1.0e-10, + atol=1.0e-10 * np.abs(np.asarray(inversion_mapping.curvature_matrix)).max(), + ) + np.testing.assert_allclose( + np.asarray(inversion_sparse_operator.reconstruction), + np.asarray(inversion_mapping.reconstruction), + rtol=1.0e-8, + ) + assert inversion_sparse_operator.log_det_curvature_reg_matrix_term == pytest.approx( + inversion_mapping.log_det_curvature_reg_matrix_term, rel=1.0e-10 + ) + assert inversion_sparse_operator.regularization_term == pytest.approx( + inversion_mapping.regularization_term, rel=1.0e-8 + ) + + +def test__inversion_imaging__delaunay_nn_nan_weights__numba_gives_nan_like_mapping( + masked_imaging_7x7, +): + """ + Overflowing Sibson rows carry NaN weights so the model sample is rejected. The mapping + route turns them into a NaN `log_det_curvature_reg_matrix_term`; the numba sparse + route must do the same, not silently return a finite one. + """ + pytest.importorskip("numba") + + mapper = _delaunay_nn_mapper_from(mask=masked_imaging_7x7.mask, max_neighbors=6) + + nan_rows = np.isnan(np.asarray(mapper.pix_weights_for_sub_slim_index)).any(axis=1) + + assert nan_rows.any() and not nan_rows.all() + + inversion_sparse_operator = aa.Inversion( + dataset=masked_imaging_7x7.apply_sparse_operator_cpu(), + linear_obj_list=[mapper], + ) + + assert isinstance(inversion_sparse_operator, InversionImagingSparseNumba) + + for inversion in ( + inversion_sparse_operator, + aa.Inversion(dataset=masked_imaging_7x7, linear_obj_list=[mapper]), + ): + assert np.isnan(np.asarray(inversion.curvature_matrix)).any() + assert np.isnan(inversion.log_det_curvature_reg_matrix_term) + + def test__inversion_imaging__linear_obj_func_and_non_func_give_same_terms( masked_imaging_7x7_no_blur, rectangular_mapper_7x7_3x3,