fix: harden scientific units, provenance, and BL19B2 resume

src/saxsabs/cli.py

@@ -46,14 +46,19 @@ def _die(message: str) -> None:
     raise SystemExit(1)
 
 
+_EXPECTED_INPUT_ERRORS = (OSError, UnicodeError, ValueError, TypeError, ImportError)
+
+
 def _clean_column_name(name: object) -> str:
     return "".join(ch for ch in str(name).strip().lower() if ch.isalnum())
 
 
 def _column_score(name: str, role: str) -> int:
+    raw_name = str(name)
+    name = _clean_column_name(name)
     if role == "q":
         exact = {"q", "chi", "radial", "2theta", "twotheta", "s", "x"}
-        prefixes = ("q", "chi", "radial", "twotheta")
+        prefixes = ("chi", "radial", "twotheta")
         suffixes = ("q",)
     else:
         exact = {"i", "intensity", "irel", "iabs", "signal", "count", "counts", "y"}
@@ -62,6 +67,8 @@ def _column_score(name: str, role: str) -> int:
 
     if name in exact:
         return 300
+    if role == "q" and q_axis_kind(raw_name) == "q":
+        return 200
     if any(name.startswith(prefix) and len(name) > len(prefix) for prefix in prefixes):
         return 200
     if any(name.endswith(suffix) and len(name) > len(suffix) for suffix in suffixes):
@@ -133,7 +140,7 @@ def _resolve_column(
     best_col = None
     best_score = 0
     for col in columns:
-        score = _column_score(_clean_column_name(col), role)
+        score = _column_score(col, role)
         if score > best_score:
             best_col = col
             best_score = score
@@ -775,13 +782,21 @@ def main() -> None:
         return
 
     if args.command == "parse-header":
-        header = json.loads(args.header_json.read_text(encoding="utf-8"))
-        exp, mon, trans = parse_header_values(header)
+        try:
+            header = json.loads(args.header_json.read_text(encoding="utf-8-sig"))
+            if not isinstance(header, dict):
+                raise ValueError("header JSON top level must be an object")
+            exp, mon, trans = parse_header_values(header)
+        except _EXPECTED_INPUT_ERRORS as exc:
+            _die(f"parse-header failed: {exc}")
         print(json.dumps({"exp_s": exp, "i0": mon, "trans": trans}, ensure_ascii=False))
         return
 
     if args.command == "parse-external1d":
-        result = read_external_1d_profile(args.input)
+        try:
+            result = read_external_1d_profile(args.input)
+        except _EXPECTED_INPUT_ERRORS as exc:
+            _die(f"parse-external1d failed: {exc}")
         print(
             json.dumps(
                 {
@@ -846,7 +861,7 @@ def main() -> None:
                     i_ref=profile_intensity(reference),
                     q_window=(args.qmin, args.qmax),
                 )
-        except ValueError as exc:
+        except _EXPECTED_INPUT_ERRORS as exc:
             _die(f"estimate-k failed: {exc}")
         print(
             json.dumps(
@@ -890,7 +905,7 @@ def main() -> None:
                 sample_profile=sample,
                 buffer_profile=buffer_profile,
             )
-        except ValueError as exc:
+        except _EXPECTED_INPUT_ERRORS as exc:
             _die(f"subtract-buffer failed: {exc}")
         print(
             json.dumps(
@@ -942,7 +957,7 @@ def main() -> None:
                 err_fluorescence=err_fluo,
                 fluorescence_profile=fluo_profile,
             )
-        except ValueError as exc:
+        except _EXPECTED_INPUT_ERRORS as exc:
             _die(f"subtract-fluorescence failed: {exc}")
         print(
             json.dumps(

src/saxsabs/core/buffer_subtraction.py

@@ -38,10 +38,12 @@ class BufferSubtractionResult:
         Propagated uncertainty.
     alpha : float
         Scaling factor applied to the buffer curve.
-    high_q_residual_mean : float
-        Mean intensity in the high-*q* diagnostic window (should be ≈0).
-    high_q_check_passed : bool
-        *True* if |mean| < 3 × σ in the diagnostic window.
+    high_q_residual_mean : float | None
+        Mean intensity in the high-*q* diagnostic window (should be ≈0), or
+        ``None`` when fewer than three points make the diagnostic undefined.
+    high_q_check_passed : bool | None
+        *True* if |mean| < 3 × σ in the diagnostic window, *False* when it
+        fails, or ``None`` when the diagnostic was not performed.
     err_statistical : np.ndarray
         Statistical component, excluding the uncertainty contribution from alpha.
     alpha_uncertainty : float | None
@@ -52,8 +54,8 @@ class BufferSubtractionResult:
     i_subtracted: np.ndarray
     err_subtracted: np.ndarray
     alpha: float
-    high_q_residual_mean: float = 0.0
-    high_q_check_passed: bool = True
+    high_q_residual_mean: float | None = None
+    high_q_check_passed: bool | None = None
     alpha_uncertainty: float | None = None
     err_statistical: np.ndarray | None = None
 
@@ -69,6 +71,16 @@ def _as_1d_float_array(name: str, values: np.ndarray | None, *, require_finite:
     return arr
 
 
+def _square_uncertainty(name: str, values: np.ndarray) -> np.ndarray:
+    """Square uncertainty values while rejecting finite overflow."""
+
+    with np.errstate(over="ignore", invalid="ignore"):
+        squared = np.square(values)
+    if np.any(np.isinf(squared)):
+        raise ValueError(f"{name} uncertainty propagation overflowed")
+    return squared
+
+
 def _prepare_source_grid(
     q_source: np.ndarray,
     y_source: np.ndarray,
@@ -123,7 +135,7 @@ def _prepare_variance_grid(
     """
     order = np.argsort(q_source)
     q_sorted = q_source[order]
-    variance_sorted = np.square(sigma_source[order])
+    variance_sorted = _square_uncertainty(label, sigma_source[order])
     uq, inv = np.unique(q_sorted, return_inverse=True)
     if uq.size < 2:
         raise ValueError(f"{label} q grid must contain at least 2 unique points")
@@ -134,7 +146,11 @@ def _prepare_variance_grid(
     for group in range(uq.size):
         group_variance = variance_sorted[inv == group]
         if np.all(np.isfinite(group_variance)):
-            variance_of_mean[group] = float(group_variance.sum() / group_variance.size**2)
+            with np.errstate(over="ignore", invalid="ignore", divide="ignore"):
+                value = group_variance.sum() / group_variance.size**2
+            if not np.isfinite(value):
+                raise ValueError(f"{label} uncertainty propagation overflowed")
+            variance_of_mean[group] = float(value)
     return uq, variance_of_mean
 
 
@@ -293,19 +309,36 @@ def subtract_buffer(
             q_s, q_b, e_b, label="buffer uncertainty"
         )
     else:
-        buffer_variance = np.square(e_b)
+        buffer_variance = _square_uncertainty("err_buffer", e_b)
 
     # Subtraction
-    i_sub = i_s - alpha * i_b
+    with np.errstate(over="ignore", invalid="ignore"):
+        i_sub = i_s - alpha * i_b
+    if not np.all(np.isfinite(i_sub)):
+        raise ValueError("buffer subtraction produced non-finite intensity")
 
     # Unknown input errors intentionally yield NaN, never an optimistic partial budget.
-    variance_statistical = np.square(e_s) + alpha**2 * buffer_variance
+    with np.errstate(over="ignore", invalid="ignore"):
+        alpha_squared = np.square(alpha)
+    if np.isinf(alpha_squared):
+        raise ValueError("buffer uncertainty propagation overflowed")
+    variance_sample = _square_uncertainty("err_sample", e_s)
+    with np.errstate(over="ignore", invalid="ignore"):
+        variance_statistical = variance_sample + alpha_squared * buffer_variance
+    if np.any(np.isinf(variance_statistical)):
+        raise ValueError("buffer uncertainty propagation overflowed")
     err_statistical = np.sqrt(variance_statistical)
     variance_sub = variance_statistical.copy()
     if alpha_uncertainty is None:
         variance_sub = variance_sub + np.full_like(i_b, np.nan)
     else:
-        variance_sub = variance_sub + np.square(i_b * alpha_uncertainty)
+        with np.errstate(over="ignore", invalid="ignore"):
+            alpha_term = i_b * alpha_uncertainty
+        alpha_variance = _square_uncertainty("alpha", alpha_term)
+        with np.errstate(over="ignore", invalid="ignore"):
+            variance_sub = variance_sub + alpha_variance
+        if np.any(np.isinf(variance_sub)):
+            raise ValueError("buffer uncertainty propagation overflowed")
     err_sub = np.sqrt(variance_sub)
 
     # High-q diagnostic
@@ -315,8 +348,8 @@ def subtract_buffer(
         residual_std = float(np.std(i_sub[mask]))
         check_ok = abs(residual_mean) < 3.0 * max(residual_std, 1e-30)
     else:
-        residual_mean = 0.0
-        check_ok = True  # not enough points for diagnostic
+        residual_mean = None
+        check_ok = None  # not enough points for diagnostic
 
     return BufferSubtractionResult(
         q=q_s,

src/saxsabs/core/calibration.py

@@ -304,7 +304,30 @@ def estimate_k_factor_robust(
         raise ValueError("q overlap with reference is insufficient")
 
     i_meas_interp = np.interp(q_ref_used, q_m, i_m)
-    valid = np.isfinite(i_meas_interp) & (i_meas_interp > positive_floor)
+    # ``np.interp`` linearly bridges a bad measured point.  That can turn a
+    # segment such as (1.0, -0.1, 1.0) into apparently positive values and
+    # invent a plausible K.  A reference point between measured samples is
+    # usable only when both bracketing measured endpoints pass the floor.  An
+    # exact source point needs only that source value itself.
+    raw_upper = np.searchsorted(q_m, q_ref_used, side="left")
+    exact_source = (raw_upper < q_m.size) & np.isclose(
+        q_m[np.clip(raw_upper, 0, q_m.size - 1)],
+        q_ref_used,
+        rtol=0.0,
+        atol=1e-14,
+    )
+    upper = np.clip(raw_upper, 1, q_m.size - 1)
+    lower = upper - 1
+    exact_indices = np.clip(raw_upper, 0, q_m.size - 1)
+    lower[exact_source] = exact_indices[exact_source]
+    upper[exact_source] = exact_indices[exact_source]
+    segment_valid = (
+        np.isfinite(i_m[lower])
+        & (i_m[lower] > positive_floor)
+        & np.isfinite(i_m[upper])
+        & (i_m[upper] > positive_floor)
+    )
+    valid = segment_valid & np.isfinite(i_meas_interp) & (i_meas_interp > positive_floor)
     if int(valid.sum()) < min_points:
         raise ValueError("measured signal too weak or non-positive in overlap region")