[#899] add_inversion() more robust vs partially-zero geometry matrices

tofu/data/_class10_compute.py

@@ -131,7 +131,8 @@ def compute_inversions(
     if hasattr(matrix, 'toarray'):
         mmm = matrix.toarray()
 
-    matnorm = np.mean(np.mean(mmm, axis=-1, where=mmm>0), axis=-1)
+    mmm = np.mean(mmm, axis=-1, where=mmm>0)
+    matnorm = np.mean(mmm, axis=-1, where=mmm>0)
     matrix_norm = matrix / matnorm[:, None, None] if m3d else matrix / matnorm
 
     # --------------------------------
@@ -181,6 +182,17 @@ def compute_inversions(
             ),
         )
 
+    # Safety check
+    if np.any(~np.isfinite(sol0)):
+        msg = (
+            "sol0 has non-finite values!\n"
+            f"\t- indok: {indok}\n"
+            f"\t- isnan(data[0, ...]): {np.any(~np.isfinite(data[0, ...]))}\n"
+            f"\t- isnan(mat0[0, ...]): {np.any(~np.isfinite(mat0[0, ...]))}\n"
+            f"\t- isnan(sol0): {np.any(~np.isfinite(sol0))}\n"
+        )
+        raise Exception(msg)
+
     if verb >= 1:
         # t1 = time.process_time()
         t1 = time.perf_counter()
@@ -798,6 +810,24 @@ def func_hess(x, mu=mu0, Tn=None, yn=None, TTn=TTn, Tyn=Tyn, R=R):
             else:
                 sol[ii, :] = sol[ii, :] + dcon['offset'][ic, :]
 
+        # safety check
+        if np.isnan(chi2n[ii]) or np.isnan(regularity[ii]):
+            lk1 = [
+                (sol0[indbsi], 'sol0[indbsi]'),
+                (Tni, 'Tni'),
+                (TTni, 'TTni'),
+                (Tyni, 'Tyni'),
+                (Ri, 'Ri'),
+                (yni, 'yni'),
+            ]
+            lstr = [f"\t- {v1}: {np.any(~np.isfinite(k1))}\n" for (k1, v1) in lk1]
+            msg = (
+                "non-finite inversion step (post):\n"
+                + "".join(lstr)
+                + f"\t- mu0: {mu0}\n"
+            )
+            raise Exception(msg)
+
         # post
         if chain:
             sol0[:] = sol[ii, :]

tofu/data/_class8_compute.py

@@ -1810,7 +1810,7 @@ def _concatenate_data_check(
     if len(dout) > 0:
         lstr = [f"\t- {k0}: {v0}" for k0, v0 in dout.items()]
         msg = (
-            f"The following data refers to no known camera in diag '{key}:\n"
+            f"The following data refers to no known camera in diag '{key}':\n"
             + "\n".join(lstr)
         )
         raise Exception(msg)