Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
40 changes: 29 additions & 11 deletions autoarray/inversion/inversion/abstract.py
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,7 @@
import numpy as np
from typing import Dict, List, Optional, Type, Union

from autoconf import cached_property
from autoconf import cached_property, is_test_mode

from autoarray.dataset.imaging.dataset import Imaging
from autoarray.dataset.interferometer.dataset import Interferometer
Expand Down Expand Up @@ -712,13 +712,23 @@ def log_det_curvature_reg_matrix_term(self) -> float:
if not self.has(cls=AbstractRegularization):
return 0.0

return 2.0 * self._xp.sum(
self._xp.log(
self._xp.diag(
self._xp.linalg.cholesky(self.curvature_reg_matrix_reduced)
try:
return 2.0 * self._xp.sum(
self._xp.log(
self._xp.diag(
self._xp.linalg.cholesky(self.curvature_reg_matrix_reduced)
)
)
)
)
except np.linalg.LinAlgError:
if is_test_mode():
# Singular curvature-reg matrix from a fabricated test-mode model;
# this evidence term is discarded. Return a benign 0.0 so unguarded
# test-mode paths do not crash (matches the reconstruction guard in
# inversion_util). Normal runs re-raise unchanged. numpy-only: the
# JAX cholesky returns NaN rather than raising.
return 0.0
raise

@property
def log_det_regularization_matrix_term(self) -> float:
Expand All @@ -737,13 +747,21 @@ def log_det_regularization_matrix_term(self) -> float:
if not self.has(cls=AbstractRegularization):
return 0.0

return 2.0 * self._xp.sum(
self._xp.log(
self._xp.diag(
self._xp.linalg.cholesky(self.regularization_matrix_reduced)
try:
return 2.0 * self._xp.sum(
self._xp.log(
self._xp.diag(
self._xp.linalg.cholesky(self.regularization_matrix_reduced)
)
)
)
)
except np.linalg.LinAlgError:
if is_test_mode():
# Singular regularization matrix from a fabricated test-mode model;
# this evidence term is discarded. Benign 0.0 in test mode, unchanged
# (raise) in normal operation. numpy-only (JAX cholesky returns NaN).
return 0.0
raise

@property
def reconstruction_noise_map_with_covariance(self) -> np.ndarray:
Expand Down
20 changes: 19 additions & 1 deletion autoarray/inversion/inversion/inversion_util.py
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,8 @@

from typing import List, Optional, Type

from autoconf import is_test_mode

from autoarray.settings import Settings

from autoarray import exc
Expand Down Expand Up @@ -221,7 +223,19 @@ def reconstruction_positive_negative_from(
curvature_reg_matrix
The curvature_matrix plus regularization matrix, overwriting the curvature_matrix in memory.
"""
return xp.linalg.solve(curvature_reg_matrix, data_vector)
try:
return xp.linalg.solve(curvature_reg_matrix, data_vector)
except np.linalg.LinAlgError:
if is_test_mode():
# Test-mode fits run fabricated / under-converged models whose
# curvature-regularization matrix can be singular; the reconstruction
# value is discarded. Return a benign dummy so unguarded test-mode
# paths (search chaining, preloads, result construction) do not crash.
# Normal runs re-raise unchanged and resample via the likelihood's
# FitException guard — real-mode numerics are untouched. (The JAX path
# returns NaN rather than raising, so this branch is numpy-only.)
return xp.ones_like(data_vector)
raise


def reconstruction_positive_only_from(
Expand Down Expand Up @@ -350,6 +364,10 @@ def reconstruction_positive_only_from(
)

except (RuntimeError, np.linalg.LinAlgError, ValueError) as e:
if is_test_mode():
# See reconstruction_positive_negative_from: benign dummy in test mode,
# unchanged (raise InversionException) in normal operation.
return xp.ones_like(data_vector)
raise exc.InversionException() from e


Expand Down
53 changes: 53 additions & 0 deletions test_autoarray/inversion/inversion/test_inversion_util.py
Original file line number Diff line number Diff line change
Expand Up @@ -228,3 +228,56 @@ def test__preconditioner_matrix_via_mapping_matrix_from():
preconditioner_matrix
== np.array([[5.0, 2.0, 3.0], [4.0, 9.0, 6.0], [7.0, 8.0, 13.0]])
).all()


def test__reconstruction_positive_negative_from__singular_matrix_test_mode_returns_dummy(
monkeypatch,
):
# A rank-deficient curvature-regularization matrix makes np.linalg.solve
# raise a singular-matrix LinAlgError.
curvature_reg_matrix = np.array([[1.0, 1.0], [1.0, 1.0]])
data_vector = np.array([1.0, 2.0])

# Normal operation: the singular solve raises, so the search resamples.
monkeypatch.delenv("PYAUTO_TEST_MODE", raising=False)
with pytest.raises(np.linalg.LinAlgError):
aa.util.inversion.reconstruction_positive_negative_from(
data_vector=data_vector,
curvature_reg_matrix=curvature_reg_matrix,
)

# Test mode: the fabricated model's reconstruction is discarded, so a benign
# finite dummy is returned instead of crashing an unguarded test-mode path.
monkeypatch.setenv("PYAUTO_TEST_MODE", "1")
reconstruction = aa.util.inversion.reconstruction_positive_negative_from(
data_vector=data_vector,
curvature_reg_matrix=curvature_reg_matrix,
)
assert reconstruction.shape == data_vector.shape
assert np.all(np.isfinite(reconstruction))
assert (reconstruction == 1.0).all()


def test__reconstruction_positive_only_from__singular_matrix_test_mode_returns_dummy(
monkeypatch,
):
curvature_reg_matrix = np.array([[1.0, 1.0], [1.0, 1.0]])
data_vector = np.array([1.0, 2.0])

# Normal operation: the singular solve raises InversionException (resample).
monkeypatch.delenv("PYAUTO_TEST_MODE", raising=False)
with pytest.raises(aa.exc.InversionException):
aa.util.inversion.reconstruction_positive_only_from(
data_vector=data_vector,
curvature_reg_matrix=curvature_reg_matrix,
)

# Test mode: benign finite dummy instead of crashing.
monkeypatch.setenv("PYAUTO_TEST_MODE", "1")
reconstruction = aa.util.inversion.reconstruction_positive_only_from(
data_vector=data_vector,
curvature_reg_matrix=curvature_reg_matrix,
)
assert reconstruction.shape == data_vector.shape
assert np.all(np.isfinite(reconstruction))
assert (reconstruction == 1.0).all()
Loading