Skip to content

Commit 627576c

Browse files
authored
Merge pull request #406 from PyAutoLabs/claude/g-heart-green-fable-opus-uuiaqf
fix: NaN-fill the reconstruction CSV noise-map column on a singular curvature_reg matrix
2 parents ba8bde2 + f461005 commit 627576c

2 files changed

Lines changed: 65 additions & 2 deletions

File tree

autoarray/inversion/plot/inversion_plots.py

Lines changed: 20 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -369,6 +369,14 @@ def save_reconstruction_csv(
369369
One file is written per mapper: ``source_plane_reconstruction_{i}.csv``,
370370
with columns ``y``, ``x``, ``reconstruction``, ``noise_map``.
371371
372+
The reconstruction noise map inverts the curvature regularization matrix, which is
373+
singular for rank-deficient inversions (e.g. the reduced-iteration searches used by
374+
test profiles). The column schema is fixed, because consumers index the CSV by
375+
column name, so in that case the ``noise_map`` column is written as ``nan`` (with a
376+
logged warning) rather than omitted, and the file is still written: the
377+
reconstruction is a science product in its own right and must not be lost, nor may
378+
a failure here abort the enclosing model-fit.
379+
372380
Parameters
373381
----------
374382
inversion
@@ -383,10 +391,20 @@ def save_reconstruction_csv(
383391
y = mapper.source_plane_mesh_grid[:, 0]
384392
x = mapper.source_plane_mesh_grid[:, 1]
385393
reconstruction = inversion.reconstruction_dict[mapper]
386-
noise_map = inversion.reconstruction_noise_map_dict[mapper]
394+
395+
try:
396+
noise_map = inversion.reconstruction_noise_map_dict[mapper]
397+
except np.linalg.LinAlgError:
398+
logger.warning(
399+
f"save_reconstruction_csv: could not compute the reconstruction noise map for "
400+
f"mapper {i} (singular curvature_reg_matrix); writing the noise_map column of "
401+
f"source_plane_reconstruction_{i}.csv as nan."
402+
)
403+
noise_map = None
387404

388405
with open(output_path / f"source_plane_reconstruction_{i}.csv", mode="w", newline="") as f:
389406
writer = csv.writer(f)
390407
writer.writerow(["y", "x", "reconstruction", "noise_map"])
391408
for j in range(len(x)):
392-
writer.writerow([float(y[j]), float(x[j]), float(reconstruction[j]), float(noise_map[j])])
409+
noise_value = float("nan") if noise_map is None else float(noise_map[j])
410+
writer.writerow([float(y[j]), float(x[j]), float(reconstruction[j]), noise_value])

test_autoarray/inversion/plot/test_inversion_plotters.py

Lines changed: 45 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,8 @@
11
import autoarray.plot as aplt
22
from autoarray.inversion.mappers.abstract import Mapper
3+
from autoarray.inversion.plot.inversion_plots import save_reconstruction_csv
34

5+
import csv
46
import numpy as np
57
import pytest
68
from pathlib import Path
@@ -92,3 +94,46 @@ def test__inversion_subplot_of_mapper__singular_curvature_reg_matrix(
9294
)
9395

9496
assert str(Path(plot_path) / "inversion_0.png") in plot_patch.paths
97+
98+
99+
def test__save_reconstruction_csv__singular_curvature_reg_matrix(
100+
rectangular_inversion_7x7_3x3,
101+
tmp_path,
102+
monkeypatch,
103+
):
104+
inversion = rectangular_inversion_7x7_3x3
105+
106+
params = inversion.linear_obj_list[0].params
107+
108+
monkeypatch.setattr(
109+
type(inversion),
110+
"reconstruction_noise_map_with_covariance",
111+
property(lambda self: np.sqrt(np.linalg.inv(np.zeros((params, params))))),
112+
)
113+
114+
with pytest.raises(np.linalg.LinAlgError):
115+
inversion.reconstruction_noise_map_dict
116+
117+
save_reconstruction_csv(inversion=inversion, output_path=tmp_path)
118+
119+
csv_path = tmp_path / "source_plane_reconstruction_0.csv"
120+
121+
assert csv_path.exists()
122+
123+
with open(csv_path, mode="r") as f:
124+
reader = csv.reader(f)
125+
header_list = next(reader)
126+
row_list = [row for row in reader]
127+
128+
# The column schema is unchanged, because consumers index the CSV by column name.
129+
130+
assert header_list == ["y", "x", "reconstruction", "noise_map"]
131+
132+
mapper = inversion.cls_list_from(cls=Mapper)[0]
133+
134+
assert len(row_list) == len(mapper.source_plane_mesh_grid)
135+
136+
for row in row_list:
137+
assert np.isfinite(float(row[0]))
138+
assert np.isfinite(float(row[1]))
139+
assert np.isnan(float(row[3]))

0 commit comments

Comments
 (0)