From 051cfd3e53c65efa616ac31b03643a81c59902d4 Mon Sep 17 00:00:00 2001 From: Claude Date: Thu, 27 Aug 2026 14:33:24 +0000 Subject: [PATCH 1/2] fix: form the reconstruction covariance on the parameters the solve solved for `AbstractInversion.reconstruction` does not always solve the full linear system. Under `use_edge_zeroed_pixels` it subsets `curvature_reg_matrix` to `zeroed_ids_to_keep`, solves the reduced problem and scatters the answer back with exact zeros at the excluded pixels -- the mesh's poorly-constrained boundary vertices, zeroed precisely to keep the inversion stable. `reconstruction_covariance_matrix` inverted the FULL matrix regardless. So it re-admitted into an explicit inverse the very rows the solve dropped to stay stable, and `reconstruction_noise_map` reported a finite noise value for a pixel whose reconstruction reads exactly 0.0 because it was never solved for. The reconstruction and its noise map disagreed about which pixels had been estimated. Not opt-in, as it was previously scoped. `Delaunay.zeroed_pixels` is a count defaulting to 0, but `RectangularRTUAdaptDensity.zeroed_pixels` returns the whole edge ring unconditionally -- there is no setting to turn it off. Every `Rectangular*AdaptDensity` fit is affected: 108 of 784 parameters on a 28x28, 156 of 1600 on a 40x40. The structural fix, not just the local one ------------------------------------------ The predicate deciding whether the solve subsets the system lived inline in `reconstruction` and nowhere else, so the covariance had no way to know. It is now one property, `solve_ids_to_keep`, which both read. It returns None (not "every index") when the solve is full, so that path never enters the indexing and scatter-back code and stays byte-identical. It also reproduces the nesting deliberately: `use_edge_zeroed_pixels` is consulted only when `use_positive_only_solver` is on. That scoping was confirmed intentional in 84b9ed42 and the comment forbidding its "fix" is preserved and extended. Excluded entries are NaN, not zero ---------------------------------- Zero is a legitimate covariance value ("known exactly"), indistinguishable from a real result. These parameters were held at zero by construction and never estimated, so the honest report is that there is no number. Both consumers handle it, and neither needed changing: `norm_from` derives colour limits with `np.nanmax`, and `save_reconstruction_csv` already writes `nan` in this column when the covariance cannot be computed. Verified by rendering the mapper subplot in both linear and log10 scales and writing the CSV. A signal-to-noise map formed as `reconstruction / reconstruction_noise_map` gives `0.0 / NaN = NaN` silently, and `NaN >= threshold` is False, so these pixels fall outside a signal-to-noise cut rather than counting as significant. (Reporting 0 instead would have made that `0/0` -- the same NaN, plus a RuntimeWarning on every call.) This also changes the solved pixels, and needs a release note ------------------------------------------------------------- Inverting the submatrix is not the corresponding block of the full inverse. For symmetric positive-definite A, `[A^-1]_keep >= (A_keep)^-1` in the PSD ordering, so every kept pixel's variance is LOWER than before. The direction is guaranteed; a regression would show up as a kept pixel getting noisier, and is asserted as such. The MAGNITUDE was not measured on a real lens fit here -- this environment has no PyAutoLens. A structural proxy on 20x20/28x28/40x40 meshes puts the median ratio at 0.983-0.992 (max 0.995, never above 1), but that is exactly the kind of proxy that understated a related effect by an order of magnitude in the investigation behind this change, because a random mapping matrix spreads data support across the whole mesh while real ray tracing concentrates it in the arc. Measure on a real fit before quoting a number in the release note. This is unrelated to the NNLS active-set caveat on `reconstruction_noise_map`, where the shipped and free-set-restricted values bracket the truth. Here the excluded parameters are fixed at zero before the solver runs rather than chosen by the data, so conditioning on them is exact and there is no bracket. Also ---- - Corrects three comments that described `mapper_indices` as "ids of values which are on edge so zero-d and not solved for". It is not edge zeroing: it drops the linear objects carrying no regularization. Conflating the two made the reduced matrices and the zeroed solve look like one inconsistent treatment of a single index set when they are unrelated concerns. - `MockMapper` gains the `mesh` accessor for the `mesh` it already stored. Tests: 9 unit cases over a 4x4 mesh (4 interior parameters kept, 12 zeroed) and 1 end-to-end case through the real solver and mesh, asserting the invariant -- `reconstruction[i] == 0.0` exactly where `reconstruction_noise_map[i]` is NaN -- plus the full-system paths unchanged, the one-directional bound, and that a non-finite entry raises LinAlgError only when the solve actually used it. Full suite: 1177 passed, 55 skipped. Known, not fixed here (pre-existing, unreachable from this change): the `np.errstate(all="ignore")` guard in `plot/utils.py:norm_from` does not suppress `np.nanmax`'s All-NaN RuntimeWarning, which is issued through `warnings` rather than the floating-point error state. A per-mapper noise map cannot be all-NaN -- a rectangular mesh is at least 3x3 and keeps its interior, Delaunay zeroes a count strictly below its pixel total -- so this change cannot trigger it. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01Lkq5ww6eLEvJgPFGMgMU1C --- autoarray/inversion/inversion/abstract.py | 147 ++++++++++++- autoarray/inversion/mock/mock_mapper.py | 6 + .../inversion/inversion/test_abstract.py | 208 +++++++++++++++++- .../inversion/inversion/test_factory.py | 31 +++ 4 files changed, 378 insertions(+), 14 deletions(-) diff --git a/autoarray/inversion/inversion/abstract.py b/autoarray/inversion/inversion/abstract.py index 5163b6102..b651e9fdc 100644 --- a/autoarray/inversion/inversion/abstract.py +++ b/autoarray/inversion/inversion/abstract.py @@ -345,7 +345,10 @@ def regularization_matrix_reduced(self) -> Optional[np.ndarray]: if self.all_linear_obj_have_regularization: return self.regularization_matrix - # ids of values which are on edge so zero-d and not solved for. + # Restrict to the mapper parameters. This is NOT edge zeroing -- `mapper_indices` drops the linear + # objects that carry no regularization (light profiles and the like), which is a different index set + # from `zeroed_ids_to_keep` and applies for a different reason. The two were conflated by an earlier + # version of this comment. ids_to_keep = self.mapper_indices # Zero rows and columns in the matrix we want to ignore @@ -383,7 +386,10 @@ def curvature_reg_matrix_reduced(self) -> Optional[np.ndarray]: if self.all_linear_obj_have_regularization: return self.curvature_reg_matrix - # ids of values which are on edge so zero-d and not solved for. + # Restrict to the mapper parameters. This is NOT edge zeroing -- `mapper_indices` drops the linear + # objects that carry no regularization (light profiles and the like), which is a different index set + # from `zeroed_ids_to_keep` and applies for a different reason. The two were conflated by an earlier + # version of this comment. ids_to_keep = self.mapper_indices # Zero rows and columns in the matrix we want to ignore @@ -490,6 +496,44 @@ def zeroed_ids_to_keep(self): return keep_ids + @property + def solve_ids_to_keep(self) -> Optional[np.ndarray]: + """ + The global parameter indices the reconstruction actually solves for, or `None` when it solves the full + system. + + This is the single answer to "which parameters did the solve include?", and it exists so that every + quantity derived from the solve agrees with the solve about what it did. `reconstruction` subsets the + linear system by these indices and scatters its result back with exact zeros elsewhere; + `reconstruction_covariance_matrix` forms the covariance on the same submatrix and scatters back `NaN`. + + Before this property the predicate below lived inline in `reconstruction` and nowhere else, so the + covariance had no way to know the solve had been subset -- it inverted the full `curvature_reg_matrix` + and reported a noise value for pixels that were never solved for. Keep the two readers pointed here + rather than re-deriving the condition, or they can drift apart again. + + `None` rather than "every index" is deliberate: the full-system path must stay byte-identical to what it + was, and an `arange` would route it through indexing and scatter-back code it never used before. + + Note that `use_edge_zeroed_pixels` is consulted only when `use_positive_only_solver` is `True`, mirroring + the nesting in `reconstruction`. That scoping is deliberate -- see the comment there and + `Settings.use_edge_zeroed_pixels`. + + Returns + ------- + The global indices kept by the solve, or `None` if the full system was solved. + """ + if not self.settings.use_positive_only_solver: + return None + + if not self.settings.use_edge_zeroed_pixels: + return None + + if not self.has(cls=Mapper): + return None + + return self.zeroed_ids_to_keep + @cached_property def reconstruction(self) -> np.ndarray: """ @@ -513,13 +557,17 @@ def reconstruction(self) -> np.ndarray: # `use_positive_only_solver`: edge-zeroing is scoped to the positive-only solver, and the # positive-negative branch below solves the full system regardless of its value. This is # intended, not an oversight -- do not "fix" it by hoisting the check out of this branch. - if self.settings.use_edge_zeroed_pixels and self.has(cls=Mapper): + # `solve_ids_to_keep` encodes that nesting (it returns None unless BOTH settings are on), + # so it is safe to consult here and nowhere higher up. + ids_to_keep = self.solve_ids_to_keep + + if ids_to_keep is not None: # Use advanced indexing to select rows/columns - data_vector = self.data_vector[self.zeroed_ids_to_keep] - curvature_reg_matrix = self.curvature_reg_matrix[ - self.zeroed_ids_to_keep - ][:, self.zeroed_ids_to_keep] + data_vector = self.data_vector[ids_to_keep] + curvature_reg_matrix = self.curvature_reg_matrix[ids_to_keep][ + :, ids_to_keep + ] # Perform reconstruction via fnnls reconstruction_partial = ( @@ -536,11 +584,11 @@ def reconstruction(self) -> np.ndarray: # Scatter the partial solution back to the full shape if self._xp.__name__.startswith("jax"): - reconstruction = reconstruction.at[self.zeroed_ids_to_keep].set( + reconstruction = reconstruction.at[ids_to_keep].set( reconstruction_partial ) else: - reconstruction[self.zeroed_ids_to_keep] = reconstruction_partial + reconstruction[ids_to_keep] = reconstruction_partial return reconstruction @@ -570,7 +618,10 @@ def reconstruction_reduced(self) -> np.ndarray: if self.all_linear_obj_have_regularization: return self.reconstruction - # ids of values which are on edge so zero-d and not solved for. + # Restrict to the mapper parameters. This is NOT edge zeroing -- `mapper_indices` drops the linear + # objects that carry no regularization (light profiles and the like), which is a different index set + # from `zeroed_ids_to_keep` and applies for a different reason. The two were conflated by an earlier + # version of this comment. ids_to_keep = self.mapper_indices # Zero rows and columns in the matrix we want to ignore @@ -853,6 +904,30 @@ def reconstruction_covariance_matrix(self) -> np.ndarray: For the RMS standard deviation of each pixel (the quantity used for scientific analysis) use `reconstruction_noise_map`, which takes the square root of this matrix's diagonal. + Formed on the parameters the solve actually solved for + ----------------------------------------------------- + When `use_edge_zeroed_pixels` applies (see `solve_ids_to_keep`), `reconstruction` does not solve the + full system: it subsets `curvature_reg_matrix` to `zeroed_ids_to_keep`, solves the reduced problem and + scatters the answer back with **exact zeros** at the excluded pixels. Those pixels are the mesh's + poorly-constrained boundary vertices, zeroed precisely to keep the inversion stable. + + This matrix is formed on that same index set and scattered back the same way, so it describes the + estimator that was actually computed. The excluded rows and columns are `NaN`, not zero: zero is a + covariance value ("known exactly"), whereas these parameters were never estimated at all. The returned + shape is always `[total_params, total_params]`, so callers do not have to branch on the settings. + + Previously the full matrix was inverted regardless, which re-admitted into an explicit inverse the very + rows the solve dropped to stay stable, and reported a noise value for a pixel whose reconstruction reads + exactly `0` because it was never solved for. + + Note this also changes the values on the parameters that ARE solved. Inverting the submatrix is not the + corresponding block of the full inverse: for a symmetric positive-definite `A`, + `[A^-1]_keep >= (A_keep)^-1` in the positive-semidefinite ordering, so every kept pixel's variance is + lower here than it was. That is the correct quantity -- the excluded parameters are held at zero by + construction, so conditioning on them is exact, not an approximation. It is unrelated to the + NNLS active-set caveat documented on `reconstruction_noise_map`, where the pixels held at zero are + chosen by the data rather than fixed in advance. + The inverse is formed from a Cholesky factorization rather than `np.linalg.inv`, for two reasons: - `cho_factor` raises `LinAlgError` when the matrix is not positive-definite. `np.linalg.inv` raises only @@ -884,8 +959,22 @@ def reconstruction_covariance_matrix(self) -> np.ndarray: """ from scipy.linalg import cho_factor, cho_solve - matrix = np.asarray(self.curvature_reg_matrix) + full_matrix = np.asarray(self.curvature_reg_matrix) + + # Form the covariance on exactly the parameters the solve solved for. `solve_ids_to_keep` is None when + # the solve used the full system, in which case this path is unchanged. + ids_to_keep = self.solve_ids_to_keep + + if ids_to_keep is None: + matrix = full_matrix + else: + ids_to_keep = np.asarray(ids_to_keep) + matrix = full_matrix[ids_to_keep][:, ids_to_keep] + # The guard runs on the SUBMATRIX, not the full one. A non-finite entry in a row the solve excluded + # cannot reach the factorization, and the reconstruction does not fail on it either -- failing here + # would make the covariance stricter than the solve it describes. It must also run before the + # scatter-back below, which fills the excluded entries with NaN deliberately. if not np.isfinite(matrix).all(): raise np.linalg.LinAlgError( "The curvature_reg_matrix contains non-finite entries (NaN or inf), so the reconstruction " @@ -904,7 +993,18 @@ def reconstruction_covariance_matrix(self) -> np.ndarray: ) # cho_solve is accurate but not bitwise symmetric; a covariance matrix is symmetric by definition. - return 0.5 * (covariance + covariance.T) + covariance = 0.5 * (covariance + covariance.T) + + if ids_to_keep is None: + return covariance + + # Scatter back to the full parameter shape, so this property's shape does not depend on the settings. + # The excluded entries are NaN ("never estimated"), which is what the solve says about them -- their + # reconstruction is an exact structural zero, not a fitted value. + full_covariance = np.full(full_matrix.shape, np.nan, dtype=covariance.dtype) + full_covariance[np.ix_(ids_to_keep, ids_to_keep)] = covariance + + return full_covariance @property def reconstruction_noise_map_with_covariance(self) -> np.ndarray: @@ -919,6 +1019,10 @@ def reconstruction_noise_map_with_covariance(self) -> np.ndarray: It now returns the covariance matrix itself, so the values differ: the diagonal holds variances rather than standard deviations, and the off-diagonals hold covariances rather than `NaN`. + Note it also inherits `reconstruction_covariance_matrix`'s index set: under `use_edge_zeroed_pixels` the + rows and columns of parameters the solve excluded are `NaN` (never estimated), and the entries that + remain are the inverse of the submatrix rather than a block of the full inverse. See that property. + Returns ------- The covariance matrix of the reconstruction (see `reconstruction_covariance_matrix`). @@ -950,6 +1054,25 @@ def reconstruction_noise_map(self): It is computed as the square root of the diagonal of `reconstruction_covariance_matrix`, which is the inverse of the same matrix used to solve for the reconstruction via the linear inversion. + Pixels the solve never estimated are `NaN` + ------------------------------------------ + Under `use_edge_zeroed_pixels` (the shipped default, and unconditional for the + `Rectangular*AdaptDensity` mesh family, whose `zeroed_pixels` is the whole edge ring) the solve excludes + the mesh's boundary vertices and writes an exact `0.0` into `reconstruction` for them. This noise map + reports **`NaN`** at exactly those pixels, meaning "never estimated" -- they have no uncertainty because + they have no fitted value. + + So the two arrays agree on which pixels were solved: + `reconstruction[i] == 0.0` exactly at an excluded pixel, and `reconstruction_noise_map[i]` is `NaN` + there. Previously those pixels carried a finite noise value computed as though they had been solved. + + `NaN` propagates rather than raising, and the consumers handle it: the colour scales derive their + limits with `np.nanmax` (`plot/utils.py:norm_from`), and `save_reconstruction_csv` already writes `nan` + into this column when the covariance cannot be computed. A signal-to-noise map formed as + `reconstruction / reconstruction_noise_map` gives `0.0 / NaN = NaN` at these pixels, silently, and a + `NaN >= threshold` comparison is `False` -- so they fall outside any signal-to-noise cut rather than + being counted as significant. + This previously took the diagonal of an elementwise-square-rooted matrix. The two are algebraically identical -- `np.sqrt` is elementwise, so it commutes with taking the diagonal -- but only numerically equivalent, since the covariance is now formed by Cholesky rather than LU. The difference is diff --git a/autoarray/inversion/mock/mock_mapper.py b/autoarray/inversion/mock/mock_mapper.py index d9bc110e0..f73b0593b 100644 --- a/autoarray/inversion/mock/mock_mapper.py +++ b/autoarray/inversion/mock/mock_mapper.py @@ -57,6 +57,12 @@ def adapt_data(self): return super().adapt_data return self._adapt_data + @property + def mesh(self): + if self._mesh is None: + return super().mesh + return self._mesh + @property def mesh_geometry(self): if self._mesh_geometry is None: diff --git a/test_autoarray/inversion/inversion/test_abstract.py b/test_autoarray/inversion/inversion/test_abstract.py index e7ce2bd7b..e5c8f34de 100644 --- a/test_autoarray/inversion/inversion/test_abstract.py +++ b/test_autoarray/inversion/inversion/test_abstract.py @@ -683,7 +683,9 @@ def test__reconstruction_noise_map__correct_diagonal_noise_values(): inversion = aa.m.MockInversion(curvature_reg_matrix=curvature_reg_matrix) - assert inversion.reconstruction_covariance_matrix[0, 0] == pytest.approx(2.5, 1.0e-2) + assert inversion.reconstruction_covariance_matrix[0, 0] == pytest.approx( + 2.5, 1.0e-2 + ) assert inversion.reconstruction_noise_map == pytest.approx( np.sqrt(np.array([2.5, 1.0, 0.5])), 1.0e-3 ) @@ -812,7 +814,9 @@ def test__reconstruction_noise_map_with_covariance__is_deprecated_alias(): with pytest.warns(DeprecationWarning, match="reconstruction_covariance_matrix"): covariance = inversion.reconstruction_noise_map_with_covariance - assert covariance == pytest.approx(inversion.reconstruction_covariance_matrix, 1.0e-12) + assert covariance == pytest.approx( + inversion.reconstruction_covariance_matrix, 1.0e-12 + ) def test__max_pixel_list_from_and_centre__returns_top_pixels_and_brightest_centre(): @@ -883,3 +887,203 @@ def test__max_pixel_list_from__filter_neighbors__excludes_adjacent_pixels_from_t 0, 8, ] + + +def _zeroed_pixel_inversion( + curvature_reg_matrix, + use_positive_only_solver=True, + use_edge_zeroed_pixels=True, + with_mapper=True, +): + """ + A `MockInversion` over a 4x4 `RectangularUniform` mesh, whose `zeroed_pixels` is the edge ring. + + 16 parameters, of which the 12 edge pixels are zeroed and the 4 interior pixels [5, 6, 9, 10] are solved + for. A 4x4 kept block is the smallest one on which "the kept block equals the inverse of the submatrix" is + a real assertion rather than a scalar identity. + """ + if with_mapper: + linear_obj = aa.m.MockMapper( + mesh=aa.mesh.RectangularUniform(shape=(4, 4)), + parameters=16, + regularization=aa.reg.Constant(), + ) + else: + linear_obj = aa.m.MockLinearObj(parameters=16, regularization=aa.reg.Constant()) + + return aa.m.MockInversion( + linear_obj_list=[linear_obj], + curvature_reg_matrix=curvature_reg_matrix, + settings=aa.Settings( + use_positive_only_solver=use_positive_only_solver, + use_edge_zeroed_pixels=use_edge_zeroed_pixels, + ), + ) + + +def _spd_matrix(n=16, seed=0): + rng = np.random.default_rng(seed) + a = rng.normal(size=(n, n)) + return a @ a.T + n * np.eye(n) + + +def test__solve_ids_to_keep__none_unless_both_settings_and_a_mapper(): + """ + The single predicate for "did the solve subset the system?". + + It must reproduce the nesting in `reconstruction` exactly: `use_edge_zeroed_pixels` is consulted ONLY + when `use_positive_only_solver` is on, and only when a `Mapper` is present. That scoping is deliberate + (see the comment in `reconstruction`), so a change that made this property answer on + `use_edge_zeroed_pixels` alone would silently start subsetting the positive-negative solve. + """ + matrix = _spd_matrix() + + assert _zeroed_pixel_inversion(matrix).solve_ids_to_keep == pytest.approx( + np.array([5, 6, 9, 10]) + ) + + assert ( + _zeroed_pixel_inversion(matrix, use_edge_zeroed_pixels=False).solve_ids_to_keep + is None + ) + assert ( + _zeroed_pixel_inversion( + matrix, use_positive_only_solver=False + ).solve_ids_to_keep + is None + ) + assert _zeroed_pixel_inversion(matrix, with_mapper=False).solve_ids_to_keep is None + + +def test__reconstruction_covariance_matrix__formed_on_the_solved_indices(): + """ + The covariance must describe the estimator that was actually computed. + + `reconstruction` subsets `curvature_reg_matrix` to `zeroed_ids_to_keep` and scatters back exact zeros. + Previously this property inverted the FULL matrix regardless, re-admitting the poorly-constrained boundary + vertices the solve dropped to stay stable. Asserted structurally -- shape, which entries are NaN, and the + kept block against an independently computed inverse of the submatrix -- rather than against baked-in + numbers. + """ + matrix = _spd_matrix() + + covariance = _zeroed_pixel_inversion(matrix).reconstruction_covariance_matrix + + keep = np.array([5, 6, 9, 10]) + excluded = np.setdiff1d(np.arange(16), keep) + + # shape does not depend on the settings -- callers never branch on it + assert covariance.shape == (16, 16) + + assert np.isnan(covariance[excluded]).all() + assert np.isnan(covariance[:, excluded]).all() + + assert covariance[np.ix_(keep, keep)] == pytest.approx( + np.linalg.inv(matrix[np.ix_(keep, keep)]), 1.0e-8 + ) + + +def test__reconstruction_covariance_matrix__excluded_entries_are_nan_not_zero(): + """ + NaN means "never estimated"; zero would mean "known exactly", which is the opposite claim. + + Zero is a legitimate covariance value, so a consumer cannot tell it apart from a real result. These + parameters were held at zero by construction and never entered the solve, so the honest report is that + there is no number. + """ + covariance = _zeroed_pixel_inversion(_spd_matrix()).reconstruction_covariance_matrix + + assert not (covariance[0] == 0.0).any() + assert np.isnan(covariance[0]).all() + + +def test__reconstruction_noise_map__nan_at_the_pixels_the_solve_zeroed(): + """ + The invariant this task exists for: the reconstruction and its noise map agree on which pixels were solved. + + Also asserts no `RuntimeWarning` escapes. `np.sqrt` of NaN propagates silently (unlike `np.sqrt` of a + negative), so the NaN convention must not reintroduce the warning storm the elementwise-sqrt bug caused. + """ + matrix = _spd_matrix() + + with warnings.catch_warnings(): + warnings.simplefilter("error", RuntimeWarning) + noise_map = _zeroed_pixel_inversion(matrix).reconstruction_noise_map + + keep = np.array([5, 6, 9, 10]) + excluded = np.setdiff1d(np.arange(16), keep) + + assert np.isnan(noise_map[excluded]).all() + assert np.isfinite(noise_map[keep]).all() + + assert noise_map[keep] == pytest.approx( + np.sqrt(np.diag(np.linalg.inv(matrix[np.ix_(keep, keep)]))), 1.0e-8 + ) + + +def test__reconstruction_noise_map__kept_pixels_are_never_noisier_than_the_full_matrix(): + """ + Restricting the inverse is not a re-scaling of the excluded rows -- it changes the SOLVED pixels too. + + For a symmetric positive-definite `A`, `[A^-1]_keep >= (A_keep)^-1` in the positive-semidefinite ordering, + so every kept pixel's variance is lower here than it was. The direction is guaranteed by that inequality, + so it is asserted as a one-directional bound rather than as a magnitude: this is the value change that + needs a release note, and a regression would show up as a kept pixel getting NOISIER. + """ + matrix = _spd_matrix() + + restricted = _zeroed_pixel_inversion(matrix).reconstruction_noise_map + full = _zeroed_pixel_inversion( + matrix, use_edge_zeroed_pixels=False + ).reconstruction_noise_map + + keep = np.array([5, 6, 9, 10]) + + assert (restricted[keep] <= full[keep]).all() + + +@pytest.mark.parametrize( + "kwargs", + [ + {"use_edge_zeroed_pixels": False}, + {"use_positive_only_solver": False}, + {"with_mapper": False}, + ], +) +def test__reconstruction_covariance_matrix__full_system_path_is_unchanged(kwargs): + """ + Every route that does NOT subset the solve must return the full inverse exactly as before, with no NaN. + + `solve_ids_to_keep` returns None rather than "every index" precisely so this path never goes through the + indexing and scatter-back code at all. + """ + matrix = _spd_matrix() + + covariance = _zeroed_pixel_inversion( + matrix, **kwargs + ).reconstruction_covariance_matrix + + assert np.isfinite(covariance).all() + assert covariance == pytest.approx(np.linalg.inv(matrix), 1.0e-8) + + +def test__reconstruction_covariance_matrix__non_finite_entry_only_raises_when_the_solve_used_it(): + """ + The `LinAlgError` contract holds on the submatrix, which is what `inversion_plots.py` guards on. + + A NaN in a row the solve excluded cannot reach the factorization and does not stop the reconstruction + either, so failing on it would make the covariance stricter than the estimator it describes. A NaN in a + KEPT row still raises -- scipy would otherwise raise `ValueError`, which the plotting and CSV callers do + not catch, and the CSV writer promises not to abort the enclosing model-fit. + """ + excluded_nan = _spd_matrix() + excluded_nan[0, 0] = np.nan # index 0 is an edge pixel, so it is zeroed + + covariance = _zeroed_pixel_inversion(excluded_nan).reconstruction_covariance_matrix + assert np.isfinite(covariance[np.ix_([5, 6, 9, 10], [5, 6, 9, 10])]).all() + + kept_nan = _spd_matrix() + kept_nan[5, 5] = np.nan # index 5 is an interior pixel, so it is solved for + + with pytest.raises(np.linalg.LinAlgError, match="non-finite"): + _zeroed_pixel_inversion(kept_nan).reconstruction_covariance_matrix diff --git a/test_autoarray/inversion/inversion/test_factory.py b/test_autoarray/inversion/inversion/test_factory.py index e57001551..a04504166 100644 --- a/test_autoarray/inversion/inversion/test_factory.py +++ b/test_autoarray/inversion/inversion/test_factory.py @@ -184,6 +184,37 @@ def test__inversion_imaging__zeroed_pixels( assert inversion.reconstruction[4] > 0.0 +def test__inversion_imaging__zeroed_pixels__noise_map_agrees_with_the_reconstruction( + masked_imaging_7x7_no_blur, + rectangular_mapper_7x7_3x3, +): + """ + End-to-end: the reconstruction and its noise map agree on which pixels were solved. + + The unit tests in `test_abstract.py` drive this through `MockInversion` with a hand-built matrix; this + asserts the same invariant through the real solver, mesh and index bookkeeping, which is where the two + previously disagreed. The 3x3 mesh zeroes its 8 edge pixels, leaving only the centre solved for. + """ + inversion = aa.Inversion( + dataset=masked_imaging_7x7_no_blur, + linear_obj_list=[rectangular_mapper_7x7_3x3], + settings=aa.Settings( + use_positive_only_solver=True, use_edge_zeroed_pixels=True + ), + ) + + reconstruction = np.asarray(inversion.reconstruction) + noise_map = np.asarray(inversion.reconstruction_noise_map) + + assert inversion.solve_ids_to_keep == pytest.approx(np.array([4])) + + # a pixel is an exact structural zero in one array exactly when it is NaN in the other + assert np.array_equal(reconstruction == 0.0, np.isnan(noise_map)) + + assert np.isfinite(noise_map[4]) + assert inversion.reconstruction_covariance_matrix.shape == (9, 9) + + def test__inversion_imaging__via_linear_obj_func_and_mapper( masked_imaging_7x7_no_blur, rectangular_mapper_7x7_3x3, From a7703a2e1b4475aed35734e067727dda0704ae17 Mon Sep 17 00:00:00 2001 From: Claude Date: Thu, 27 Aug 2026 14:59:29 +0000 Subject: [PATCH 2/2] fix: correct the NaN/zero invariant, and record the measured value change Measuring this change on a real ray-traced fit (PyAutoLens + PyAutoGalaxy + PyAutoFit at main, against this branch) turned up an overclaim in the previous commit and supplied the magnitude it said was missing. The invariant was stated too strongly ------------------------------------- The docstring and the end-to-end test claimed a biconditional: a pixel reads exactly 0.0 in the reconstruction exactly where the noise map is NaN. That is false, and the difference is not marginal. The non-negative solver ALSO pins pixels it did solve for at exactly 0.0, wherever the fit wants no flux there, and those are fitted values carrying real error bars. On a 28x28 RectangularBilinearAdaptDensity fit, 784 parameters: reconstruction == 0.0 : 603 NaN in noise map : 108 <- only these were never estimated solved, but pinned at 0 by NNLS : 495 So `NaN => reconstruction == 0.0` holds, the converse does not, and `np.isnan(reconstruction_noise_map)` is the only way to identify the never-estimated pixels. The test passed only because the 3x3 fixture it runs on solves a single pixel, which makes the two sets coincide by accident. It now asserts the NaN set equals the excluded set exactly, and the one-way implication, with a comment saying why the converse is deliberately absent. This is also precisely the Defect 1 / Defect 2 boundary: 495 pixels pinned by the constraint are the NNLS active-set question (documented in #472, deferred), and 108 structurally zeroed ones are this change. Measured magnitude, at the evidence-optimal coefficient ------------------------------------------------------- The previous commit said the magnitude had not been measured on a real fit and that the structural proxy standing in for it should not be trusted. It should not have been: the proxy predicted a uniform 1-2% shift and the real answer is a much smaller median with a long tail. Rectangular 28x28 (108 of 784 zeroed), lambda* = 1 by log evidence on a 9-point grid, interior to the grid in every case: old/new p50 1.0002 p90 1.069 p99 1.44 max 1.81 21% of solved pixels move >1%, 9% move >10% Delaunay, set up as the workspace's own delaunay.py does (Overlay image mesh plus a ring of edge points, zeroed_pixels=30) -- which answers the "Delaunay untested" caveat this change inherited: old/new p50 1.0045 p99 1.90 max 2.37 With the default zeroed_pixels=0 the property is unchanged to the bit, and solve_ids_to_keep is None, confirming the no-op path. Downstream source flux and magnification through the workspace's S/N >= 5 cut moved by 0.00% in every case. The pixels whose noise changes most sit at the mesh edge, where the reconstruction has no flux to move: the counts of lit pixels surviving the cut are identical old and new (7/7, 24/24, 121/121). So the release note should say: a systematic, one-directional correction concentrated at the mesh edge, negligible in the median, up to ~2x on individual edge pixels, with no effect on published source flux or magnification. Full suite: 1177 passed, 55 skipped. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01Lkq5ww6eLEvJgPFGMgMU1C --- autoarray/inversion/inversion/abstract.py | 28 +++++++++++++++++-- .../inversion/inversion/test_factory.py | 18 +++++++++--- 2 files changed, 39 insertions(+), 7 deletions(-) diff --git a/autoarray/inversion/inversion/abstract.py b/autoarray/inversion/inversion/abstract.py index b651e9fdc..ba5b8beb1 100644 --- a/autoarray/inversion/inversion/abstract.py +++ b/autoarray/inversion/inversion/abstract.py @@ -928,6 +928,21 @@ def reconstruction_covariance_matrix(self) -> np.ndarray: NNLS active-set caveat documented on `reconstruction_noise_map`, where the pixels held at zero are chosen by the data rather than fixed in advance. + Measured on real ray-traced fits (Isothermal `einstein_radius=1.6` + shear, compact Sersic source, + `r=3.0"` mask, `over_sample_size_pixelization=4`, PSF and Poisson noise), at the regularization + coefficient the Bayesian evidence selects: + + - `RectangularBilinearAdaptDensity(28, 28)`, 108 of 784 parameters zeroed: the old value was + overstated by a median factor of 1.0002, but 21% of solved pixels move by more than 1%, 9% by more + than 10%, and the worst by 1.8x. The change is concentrated near the mesh edge rather than spread. + - `Delaunay` set up as the workspace's own example does (an `Overlay` image mesh plus a ring of edge + points, `zeroed_pixels=30`): median 1.0045, 99th percentile 1.90, worst 2.37. Delaunay is affected + too whenever `zeroed_pixels > 0`; with the default `zeroed_pixels=0` this property is unchanged to + the bit. + - Downstream source flux and magnification through the workspace's `S/N >= 5` cut moved by 0.00% in + every case: the pixels whose noise changes most are near the edge, where the reconstruction has no + flux to move. + The inverse is formed from a Cholesky factorization rather than `np.linalg.inv`, for two reasons: - `cho_factor` raises `LinAlgError` when the matrix is not positive-definite. `np.linalg.inv` raises only @@ -1062,9 +1077,16 @@ def reconstruction_noise_map(self): reports **`NaN`** at exactly those pixels, meaning "never estimated" -- they have no uncertainty because they have no fitted value. - So the two arrays agree on which pixels were solved: - `reconstruction[i] == 0.0` exactly at an excluded pixel, and `reconstruction_noise_map[i]` is `NaN` - there. Previously those pixels carried a finite noise value computed as though they had been solved. + So the two arrays agree on which pixels were solved. The implication runs one way, and the direction + matters: **`NaN` implies the reconstruction is exactly `0.0` there, but not the converse.** The + non-negative solver also pins pixels it DID solve for at exactly `0.0`, and those keep a finite, + meaningful noise value. On a representative fit -- `RectangularBilinearAdaptDensity(28, 28)`, 784 + parameters -- 603 pixels read `0.0` while only the 108 structurally excluded ones are `NaN`; the other + 495 were solved and pinned by the constraint. `np.isnan(reconstruction_noise_map)` is therefore the way + to identify the never-estimated pixels; `reconstruction == 0.0` is not, and conflates the two. + + Previously the excluded pixels carried a finite noise value computed as though they had been solved, + so there was no way to tell them apart from the pinned ones at all. `NaN` propagates rather than raising, and the consumers handle it: the colour scales derive their limits with `np.nanmax` (`plot/utils.py:norm_from`), and `save_reconstruction_csv` already writes `nan` diff --git a/test_autoarray/inversion/inversion/test_factory.py b/test_autoarray/inversion/inversion/test_factory.py index a04504166..aded1b5b1 100644 --- a/test_autoarray/inversion/inversion/test_factory.py +++ b/test_autoarray/inversion/inversion/test_factory.py @@ -206,12 +206,22 @@ def test__inversion_imaging__zeroed_pixels__noise_map_agrees_with_the_reconstruc reconstruction = np.asarray(inversion.reconstruction) noise_map = np.asarray(inversion.reconstruction_noise_map) - assert inversion.solve_ids_to_keep == pytest.approx(np.array([4])) + keep = np.asarray(inversion.solve_ids_to_keep) + excluded = np.setdiff1d(np.arange(reconstruction.shape[0]), keep) - # a pixel is an exact structural zero in one array exactly when it is NaN in the other - assert np.array_equal(reconstruction == 0.0, np.isnan(noise_map)) + assert keep == pytest.approx(np.array([4])) - assert np.isfinite(noise_map[4]) + # the NaN set is EXACTLY the set the solve excluded -- no more, no less + assert np.array_equal(np.flatnonzero(np.isnan(noise_map)), np.sort(excluded)) + + # NaN implies an exact structural zero. The converse is deliberately NOT asserted: the + # non-negative solver also pins pixels it DID solve for at exactly 0.0, and those keep a + # finite noise value. Asserting the biconditional would pass here only because this 3x3 + # fixture happens to solve a single pixel -- on a real fit 603 of 784 pixels read 0.0 while + # only the 108 excluded ones are NaN. + assert (reconstruction[excluded] == 0.0).all() + + assert np.isfinite(noise_map[keep]).all() assert inversion.reconstruction_covariance_matrix.shape == (9, 9)