From d91500db040c11d5c55d6aec740043566a41b87e Mon Sep 17 00:00:00 2001 From: Claude Date: Sat, 22 Aug 2026 13:33:41 +0000 Subject: [PATCH 1/2] fix: form the reconstruction covariance via Cholesky, not elementwise sqrt `reconstruction_noise_map_with_covariance` applied np.sqrt elementwise to the whole inverse of curvature_reg_matrix. The off-diagonals of a covariance matrix are covariances and are routinely negative, so every one was NaN by construction -- for any matrix, however well-conditioned -- and each call emitted "RuntimeWarning: invalid value encountered in sqrt". The property's docstring promised a matrix accounting for "the covariance of the noise between pixels"; the entries carrying that covariance were the broken ones. The 1D reconstruction_noise_map was NOT affected: np.sqrt is elementwise, so it commutes with taking the diagonal. Replaced with reconstruction_covariance_matrix, which returns C directly and forms it from a Cholesky factorization. An A/B against np.linalg.inv refuted two hypotheses and confirmed one: - inv gives negative diagonals on well-formed SPD matrices: REFUTED (0 across cond 1e3-1e15, n=400, 20 trials each) - inv is materially less accurate on the diagonal: REFUTED (matches cho_solve; at cond 1e15 inv was marginally better) - near-coincident mesh vertices degrade the inverse: REFUTED (regularization keeps the matrix PD, cond ~6.8e7 at duplicated columns) - inv returns asymmetric output: CONFIRMED (5.2e-7 at cond 1e12, vs 2.6e-16 for cho_solve) - inv silently succeeds on indefinite matrices: CONFIRMED The last is the substantive one. cho_factor raises LinAlgError on a matrix with a negative eigenvalue; inv succeeds and returns a plausible-looking covariance. At eigenvalue -1e-8 all 300 diagonals came back negative; at -1.0, zero did -- no NaN, no warning, no error, and wrong numbers. A covariance is only defined for a positive-definite matrix, so this adds a definiteness check the code lacked entirely. Both existing call sites already catch LinAlgError (inversion_plots.py:169 and :395), verified rather than assumed. reconstruction_noise_map now computes sqrt(diag(C)) directly. It was correct only incidentally before, via the elementwise sqrt; the invariant is now stated, and it cannot silently become a variance if the matrix changes. The old name stays as a deprecated alias. Its values do change -- diagonal from std-dev to variance, off-diagonals from NaN to covariances -- so the warning says so explicitly. Tests: the only prior assertion checked [0, 0], a diagonal element, which is why this shipped. Added off-diagonal finiteness under -W error::RuntimeWarning, exact symmetry, the sqrt(diag(C)) invariant, the indefinite-matrix raise, and the deprecation. Repointed two plotter monkeypatch sites to the new property. Full suite: 1150 passed, 1 skipped. The 3 test_transformer.py pynufft failures are pre-existing on a6b07cd and unrelated. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_0133X4XhMV91SFjzV2mK4Ejh --- autoarray/inversion/inversion/abstract.py | 84 +++++++++++++---- .../inversion/inversion/test_abstract.py | 91 ++++++++++++++++++- .../inversion/plot/test_inversion_plotters.py | 8 +- 3 files changed, 160 insertions(+), 23 deletions(-) diff --git a/autoarray/inversion/inversion/abstract.py b/autoarray/inversion/inversion/abstract.py index 8bc50d2bd..670e605bf 100644 --- a/autoarray/inversion/inversion/abstract.py +++ b/autoarray/inversion/inversion/abstract.py @@ -1,4 +1,5 @@ import copy +import warnings import numpy as np from typing import Dict, List, Optional, Type, Union @@ -836,27 +837,78 @@ def log_det_regularization_matrix_term(self) -> float: return self._log_det_symmetric_from(self.regularization_matrix_reduced) @property - def reconstruction_noise_map_with_covariance(self) -> np.ndarray: + def reconstruction_covariance_matrix(self) -> np.ndarray: """ - Returns the noise-map of the reconstruction as a two dimension matrix which accounts for the covariance - of the noise between pixels. + Returns the covariance matrix of the reconstruction, ``C = [F + reg_coeff*H]^-1``. + + This is the inverse of the curvature matrix with regularization -- the same matrix used to solve for the + reconstruction via the linear inversion. Its diagonal holds the variance of each reconstructed pixel and + its off-diagonal entries the covariances between pixels; the off-diagonals are routinely negative, which is + a property of the covariance and not an error. + + 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. + + 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 + on an exactly singular matrix, so an indefinite `curvature_reg_matrix` -- which the inversion does + encounter, hence `Settings.no_regularization_add_to_curvature_diag_value` -- previously returned a + plausible-looking covariance with no error and no warning. A covariance is only defined for a + positive-definite matrix, so failing here is correct and is handled by the callers in + `inversion/plot/inversion_plots.py`. + - `np.linalg.inv` is LU-based and exploits neither the symmetry nor the positive-definiteness this matrix + has. Its output drifts out of symmetry as conditioning worsens (measured at ~5e-7 absolute at + `cond ~ 1e12`, against ~3e-16 for the Cholesky solve), while a covariance matrix is symmetric by + definition. + + This property is NumPy-only: the input is coerced with `np.asarray`, so a JAX `curvature_reg_matrix` + forces a device-to-host transfer. It is a post-fit diagnostic, not part of the likelihood, so it is not on + the JIT path. + + Returns + ------- + The covariance matrix of the reconstruction, of shape [total_params, total_params]. + """ + from scipy.linalg import cho_factor, cho_solve + + matrix = np.asarray(self.curvature_reg_matrix) + + covariance = cho_solve( + cho_factor(matrix), np.eye(matrix.shape[0], dtype=matrix.dtype) + ) + + # cho_solve is accurate but not bitwise symmetric; a covariance matrix is symmetric by definition. + return 0.5 * (covariance + covariance.T) - The diagonal of this matrix is the noise-map of the reconstruction, which can be used for analysing the - reconstruction with noise properties that are representative of the fit and therefore should be used - for any scientific analysis (e.g. source reconstructions of strong lenses). + @property + def reconstruction_noise_map_with_covariance(self) -> np.ndarray: + """ + Deprecated alias of `reconstruction_covariance_matrix`. - This noise-map is defined as the RMS standard deviation of the noise in every pixel of the reconstruction. - This definition is identical to the `noise_map` attributes of dataset objects. + This property previously returned ``np.sqrt(np.linalg.inv(curvature_reg_matrix))`` -- an elementwise square + root of the whole covariance matrix. Because the off-diagonal entries of a covariance matrix are routinely + negative, every such entry was `NaN` by construction, for any input matrix however well-conditioned, and + each call emitted `RuntimeWarning: invalid value encountered in sqrt`. - It is computed as the square root of the inverse of the curvature matrix with regularization, which is the - same matrix used to solve for the reconstruction via the linear inversion. + 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`. Returns ------- - The noise-map of the reconstruction as a two dimension matrix which accounts for the covariance of the noise - between pixels. + The covariance matrix of the reconstruction (see `reconstruction_covariance_matrix`). """ - return np.sqrt(np.linalg.inv(self.curvature_reg_matrix)) + warnings.warn( + "`reconstruction_noise_map_with_covariance` is deprecated; use " + "`reconstruction_covariance_matrix` instead. Note the values have changed: it now returns the " + "covariance matrix, so its diagonal holds variances rather than standard deviations (the previous " + "elementwise square root made every off-diagonal NaN). For the RMS noise of each pixel use " + "`reconstruction_noise_map`.", + DeprecationWarning, + stacklevel=2, + ) + + return self.reconstruction_covariance_matrix @property def reconstruction_noise_map(self): @@ -870,15 +922,15 @@ def reconstruction_noise_map(self): The noise-map of the reconstruction is the RMS standard deviation of the noise in every pixel of the reconstruction. This definition is identical to the `noise_map` attributes of dataset objects. - It is computed as the square root of the diagonal of the `reconstruction_noise_map_with_covariance` matrix, - which is the same matrix used to solve for the reconstruction via the linear inversion. + 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. Returns ------- The noise-map of the reconstruction as a one dimensional ndarray, which does not account for the covariance of the noise between pixels. """ - return np.diagonal(self.reconstruction_noise_map_with_covariance) + return np.sqrt(np.diag(self.reconstruction_covariance_matrix)) @property def reconstruction_noise_map_dict(self) -> Dict[LinearObj, np.ndarray]: diff --git a/test_autoarray/inversion/inversion/test_abstract.py b/test_autoarray/inversion/inversion/test_abstract.py index 8cbee8161..a98632661 100644 --- a/test_autoarray/inversion/inversion/test_abstract.py +++ b/test_autoarray/inversion/inversion/test_abstract.py @@ -1,3 +1,5 @@ +import warnings + import numpy as np import pytest @@ -681,14 +683,97 @@ def test__reconstruction_noise_map__asymmetric_curvature_reg_matrix__correct_dia inversion = aa.m.MockInversion(curvature_reg_matrix=curvature_reg_matrix) - assert inversion.reconstruction_noise_map_with_covariance[0, 0] == pytest.approx( - np.sqrt(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 ) +def test__reconstruction_covariance_matrix__off_diagonals_are_finite_and_negative(): + """ + The off-diagonal entries of a covariance matrix are covariances and are routinely negative. + + `reconstruction_covariance_matrix` previously applied `np.sqrt` elementwise to the whole inverse, so every + negative off-diagonal became NaN by construction -- for any matrix, however well-conditioned -- while + emitting `RuntimeWarning: invalid value encountered in sqrt`. Only the [0, 0] diagonal element was asserted, + so nothing caught it. + """ + curvature_reg_matrix = np.array([[1.0, 1.0, 1.0], [1.0, 2.0, 1.0], [1.0, 1.0, 3.0]]) + + inversion = aa.m.MockInversion(curvature_reg_matrix=curvature_reg_matrix) + + with warnings.catch_warnings(): + warnings.simplefilter("error", RuntimeWarning) + covariance = inversion.reconstruction_covariance_matrix + + assert np.all(np.isfinite(covariance)) + + # this matrix has anti-correlated pixels, so the off-diagonals are genuinely negative + assert covariance[0, 1] < 0.0 + assert covariance == pytest.approx(np.linalg.inv(curvature_reg_matrix), 1.0e-8) + + +def test__reconstruction_covariance_matrix__is_exactly_symmetric(): + """A covariance matrix is symmetric by definition; `np.linalg.inv` drifts out of symmetry.""" + rng = np.random.default_rng(1234) + q, _ = np.linalg.qr(rng.standard_normal((25, 25))) + curvature_reg_matrix = (q * np.logspace(0, 9, 25)) @ q.T + curvature_reg_matrix = 0.5 * (curvature_reg_matrix + curvature_reg_matrix.T) + + inversion = aa.m.MockInversion(curvature_reg_matrix=curvature_reg_matrix) + + covariance = inversion.reconstruction_covariance_matrix + + assert covariance == pytest.approx(covariance.T, abs=1.0e-15) + + +def test__reconstruction_noise_map__is_sqrt_of_covariance_diagonal(): + """ + The invariant, asserted directly rather than via hand-computed values. + + `reconstruction_noise_map` used to be `np.diagonal(...)` of an already-square-rooted matrix, which was + correct only incidentally -- because `np.sqrt` is elementwise. It now takes the square root of the + covariance diagonal itself, so the relationship is stated rather than emergent. + """ + curvature_reg_matrix = np.array([[1.0, 1.0, 1.0], [1.0, 2.0, 1.0], [1.0, 1.0, 3.0]]) + + inversion = aa.m.MockInversion(curvature_reg_matrix=curvature_reg_matrix) + + assert inversion.reconstruction_noise_map == pytest.approx( + np.sqrt(np.diag(inversion.reconstruction_covariance_matrix)), 1.0e-12 + ) + + +def test__reconstruction_covariance_matrix__raises_on_a_non_positive_definite_matrix(): + """ + A covariance is only defined for a positive-definite matrix. + + `np.linalg.inv` raises only on an exactly singular matrix, so an indefinite `curvature_reg_matrix` returned + a plausible-looking covariance with no error and no warning. The Cholesky factorization rejects it, and the + plotting and CSV callers already catch `LinAlgError`. + """ + # symmetric, non-singular, but indefinite (eigenvalues +1 and -1) + curvature_reg_matrix = np.array([[0.0, 1.0], [1.0, 0.0]]) + + inversion = aa.m.MockInversion(curvature_reg_matrix=curvature_reg_matrix) + + assert np.isfinite(np.linalg.inv(curvature_reg_matrix)).all() # inv is silent here + + with pytest.raises(np.linalg.LinAlgError): + inversion.reconstruction_covariance_matrix + + +def test__reconstruction_noise_map_with_covariance__is_deprecated_alias(): + curvature_reg_matrix = np.array([[1.0, 1.0, 1.0], [1.0, 2.0, 1.0], [1.0, 1.0, 3.0]]) + + inversion = aa.m.MockInversion(curvature_reg_matrix=curvature_reg_matrix) + + 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) + + def test__max_pixel_list_from_and_centre__returns_top_pixels_and_brightest_centre(): source_plane_mesh_grid = aa.Grid2DIrregular( diff --git a/test_autoarray/inversion/plot/test_inversion_plotters.py b/test_autoarray/inversion/plot/test_inversion_plotters.py index 567c222b0..ccb0aa595 100644 --- a/test_autoarray/inversion/plot/test_inversion_plotters.py +++ b/test_autoarray/inversion/plot/test_inversion_plotters.py @@ -79,8 +79,8 @@ def test__inversion_subplot_of_mapper__singular_curvature_reg_matrix( monkeypatch.setattr( type(inversion), - "reconstruction_noise_map_with_covariance", - property(lambda self: np.sqrt(np.linalg.inv(np.zeros((params, params))))), + "reconstruction_covariance_matrix", + property(lambda self: np.linalg.inv(np.zeros((params, params)))), ) with pytest.raises(np.linalg.LinAlgError): @@ -107,8 +107,8 @@ def test__save_reconstruction_csv__singular_curvature_reg_matrix( monkeypatch.setattr( type(inversion), - "reconstruction_noise_map_with_covariance", - property(lambda self: np.sqrt(np.linalg.inv(np.zeros((params, params))))), + "reconstruction_covariance_matrix", + property(lambda self: np.linalg.inv(np.zeros((params, params)))), ) with pytest.raises(np.linalg.LinAlgError): From f0aefa8c1a4a90136fd0bdb640797cc41d2b8a5b Mon Sep 17 00:00:00 2001 From: Claude Date: Sat, 22 Aug 2026 13:50:36 +0000 Subject: [PATCH 2/2] fix: raise LinAlgError on non-finite input, symmetrize before Cholesky Addresses an adversarial review of the previous commit. 1. MAJOR regression, now fixed. scipy's cho_factor/cho_solve default to check_finite=True and raise ValueError -- not LinAlgError -- on a NaN/inf matrix. Both callers (inversion_plots.py:169 and :397) catch only LinAlgError, and the CSV writer's docstring explicitly promises that a failure there may not abort the enclosing model-fit. So a NaN-contaminated curvature matrix would have killed a fit that previously wrote a nan column and continued. Verified before and after: the old code returned [nan, nan]; the new code raised builtins.ValueError past the guards. Fixed by checking finiteness explicitly and raising LinAlgError with a message naming the cause, rather than by broadening the callers' except clauses -- the property's contract is now "raises LinAlgError for any input that has no covariance", which holds downstream too. check_finite=False is then passed to scipy, so the explicit check costs nothing net. 2. cho_factor reads only the upper triangle, so an asymmetric input was silently inverted as though its lower triangle matched. Verified: for [[2.0, 0.5], [0.1, 2.0]] it gave diag 0.5333 against the true inverse's 0.5063. The output symmetrization does not address this -- it symmetrizes the result, not the input -- so the input is now symmetrized too. curvature_reg_matrix is F + H and symmetric by construction, so this is defensive, but the silent wrong answer is not an acceptable failure mode. 3. Softened the "no value change" claim for reconstruction_noise_map. The two forms are algebraically identical but only numerically equivalent, since the covariance is now formed by Cholesky rather than LU: ~7e-15 relative at cond 1e3 rising to ~4e-5 at cond 1e13. Neither is the more correct result. 4. Tests. The symmetry test was tautological -- 0.5 * (C + C.T) is bitwise symmetric for any C, so it could only fail if the line were deleted. It now also asserts accuracy against an exactly-constructed ground truth (A = Q diag(w) Q.T), with the tolerance set from the measured error (~3e-9 at cond 1e9) and a comment recording that np.linalg.inv is marginally more accurate on that matrix, not less. Added tests for the non-finite raise and the asymmetric-input symmetrization. Fixed a stale test name that said "asymmetric_curvature_reg_matrix" over a symmetric matrix. Downstream check for the value-changing deprecated alias: no consumer of reconstruction_noise_map_with_covariance exists in PyAutoGalaxy (3ca31bf), PyAutoLens (87e5827) or autolens_workspace -- zero hits for "with_covariance" in any of them. reconstruction_noise_map IS used, in autolens_workspace source_science.py scripts computing signal_to_noise_map = reconstruction / reconstruction_noise_map, and its values are unchanged. Full suite: 1152 passed, 1 skipped. The 3 test_transformer.py pynufft failures are pre-existing on a6b07cd and unrelated. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_0133X4XhMV91SFjzV2mK4Ejh --- autoarray/inversion/inversion/abstract.py | 29 ++++++++++- .../inversion/inversion/test_abstract.py | 49 +++++++++++++++++-- 2 files changed, 73 insertions(+), 5 deletions(-) diff --git a/autoarray/inversion/inversion/abstract.py b/autoarray/inversion/inversion/abstract.py index 670e605bf..348fc8c60 100644 --- a/autoarray/inversion/inversion/abstract.py +++ b/autoarray/inversion/inversion/abstract.py @@ -862,6 +862,14 @@ def reconstruction_covariance_matrix(self) -> np.ndarray: `cond ~ 1e12`, against ~3e-16 for the Cholesky solve), while a covariance matrix is symmetric by definition. + Every failure mode raises `LinAlgError`, including a non-finite `curvature_reg_matrix`. That case is + checked explicitly because scipy would otherwise raise `ValueError`, which the callers -- here and + downstream -- do not catch; the previous implementation returned a silently all-`NaN` matrix instead. + + The matrix is symmetrized on input. `cho_factor` reads only the upper triangle, so an asymmetric input + would otherwise be inverted as though its lower triangle matched, silently and with no diagnostic. + `curvature_reg_matrix` is `F + H` and symmetric by construction, so this is defensive only. + This property is NumPy-only: the input is coerced with `np.asarray`, so a JAX `curvature_reg_matrix` forces a device-to-host transfer. It is a post-fit diagnostic, not part of the likelihood, so it is not on the JIT path. @@ -874,8 +882,21 @@ def reconstruction_covariance_matrix(self) -> np.ndarray: matrix = np.asarray(self.curvature_reg_matrix) + if not np.isfinite(matrix).all(): + raise np.linalg.LinAlgError( + "The curvature_reg_matrix contains non-finite entries (NaN or inf), so the reconstruction " + "covariance is undefined. Raised as LinAlgError so the plotting and CSV callers, which guard " + "on LinAlgError, degrade gracefully rather than aborting the model-fit." + ) + + # cho_factor reads only the upper triangle; symmetrize so an asymmetric input cannot be silently + # inverted as though its lower triangle matched its upper. + matrix = 0.5 * (matrix + matrix.T) + covariance = cho_solve( - cho_factor(matrix), np.eye(matrix.shape[0], dtype=matrix.dtype) + cho_factor(matrix, check_finite=False), + np.eye(matrix.shape[0], dtype=matrix.dtype), + check_finite=False, ) # cho_solve is accurate but not bitwise symmetric; a covariance matrix is symmetric by definition. @@ -925,6 +946,12 @@ 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. + 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 + conditioning-limited roundoff, measured at ~7e-15 relative at `cond ~ 1e3` rising to ~4e-5 at + `cond ~ 1e13`; neither result is the more correct one. + Returns ------- The noise-map of the reconstruction as a one dimensional ndarray, which does not account for the covariance diff --git a/test_autoarray/inversion/inversion/test_abstract.py b/test_autoarray/inversion/inversion/test_abstract.py index a98632661..e7ce2bd7b 100644 --- a/test_autoarray/inversion/inversion/test_abstract.py +++ b/test_autoarray/inversion/inversion/test_abstract.py @@ -678,7 +678,7 @@ def test__log_det_method__slogdet_is_finite_where_cholesky_fails_on_non_positive assert result == pytest.approx(np.linalg.slogdet(matrix)[1], 1.0e-8) -def test__reconstruction_noise_map__asymmetric_curvature_reg_matrix__correct_diagonal_noise_values(): +def test__reconstruction_noise_map__correct_diagonal_noise_values(): curvature_reg_matrix = np.array([[1.0, 1.0, 1.0], [1.0, 2.0, 1.0], [1.0, 1.0, 3.0]]) inversion = aa.m.MockInversion(curvature_reg_matrix=curvature_reg_matrix) @@ -713,20 +713,61 @@ def test__reconstruction_covariance_matrix__off_diagonals_are_finite_and_negativ assert covariance == pytest.approx(np.linalg.inv(curvature_reg_matrix), 1.0e-8) -def test__reconstruction_covariance_matrix__is_exactly_symmetric(): - """A covariance matrix is symmetric by definition; `np.linalg.inv` drifts out of symmetry.""" +def test__reconstruction_covariance_matrix__is_accurate_and_symmetric_when_ill_conditioned(): + """ + Ground truth is exact by construction: for `A = Q diag(w) Q.T` the inverse is `Q diag(1/w) Q.T`. + + The symmetry half only guards the symmetrization line -- `0.5 * (C + C.T)` is bitwise symmetric for any C -- + so the accuracy assertion against the constructed truth is what tests the factorization itself. + """ rng = np.random.default_rng(1234) q, _ = np.linalg.qr(rng.standard_normal((25, 25))) - curvature_reg_matrix = (q * np.logspace(0, 9, 25)) @ q.T + eigenvalues = np.logspace(0, 9, 25) + + curvature_reg_matrix = (q * eigenvalues) @ q.T curvature_reg_matrix = 0.5 * (curvature_reg_matrix + curvature_reg_matrix.T) + covariance_true = (q * (1.0 / eigenvalues)) @ q.T + inversion = aa.m.MockInversion(curvature_reg_matrix=curvature_reg_matrix) covariance = inversion.reconstruction_covariance_matrix + # cond ~ 1e9, so the achievable accuracy is eps * cond ~ 2e-7; the measured error is ~3e-9. This is not a + # claim that Cholesky beats LU here -- it does not, `np.linalg.inv` measures ~7e-10 on this matrix. + assert covariance == pytest.approx(covariance_true, abs=1.0e-7) assert covariance == pytest.approx(covariance.T, abs=1.0e-15) +def test__reconstruction_covariance_matrix__asymmetric_input_is_symmetrized_not_silently_upper_triangle(): + """ + `cho_factor` reads only the upper triangle, so an asymmetric input would be inverted as though its lower + triangle matched its upper -- silently, and differing from the true inverse. + """ + curvature_reg_matrix = np.array([[2.0, 0.5], [0.1, 2.0]]) + symmetrized = 0.5 * (curvature_reg_matrix + curvature_reg_matrix.T) + + inversion = aa.m.MockInversion(curvature_reg_matrix=curvature_reg_matrix) + + assert inversion.reconstruction_covariance_matrix == pytest.approx( + np.linalg.inv(symmetrized), 1.0e-8 + ) + + +def test__reconstruction_covariance_matrix__non_finite_matrix_raises_lin_alg_error(): + """ + scipy raises `ValueError` on a non-finite matrix, which the plotting and CSV callers do not catch -- they + guard on `LinAlgError`. The CSV writer explicitly promises not to abort the enclosing model-fit, so the + non-finite case is converted rather than allowed to escape. + """ + curvature_reg_matrix = np.array([[1.0, np.nan], [np.nan, 2.0]]) + + inversion = aa.m.MockInversion(curvature_reg_matrix=curvature_reg_matrix) + + with pytest.raises(np.linalg.LinAlgError, match="non-finite"): + inversion.reconstruction_covariance_matrix + + def test__reconstruction_noise_map__is_sqrt_of_covariance_diagonal(): """ The invariant, asserted directly rather than via hand-computed values.