diff --git a/CHANGELOG.md b/CHANGELOG.md index 577cecb..41b13e3 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,6 +2,16 @@ ## [2.0.0] - Unreleased +### Fluorescence subtraction (1D) + +- Add an opt-in absolute-scale fluorescence kernel + (`I_corr = I_abs − β F(q)`) with `constant`, `high_q_mean`, + `high_q_median`, and `measured_profile` methods, a `fluorescence` + correction-ledger token, CLI `subtract-fluorescence`, Workbench Tab 3/Tab 2 + 1-D hooks, and optional BL19B2 integrate1d post-processing. Default remains + off. Unknown `u(F0)`/`u(β)` keeps combined uncertainty NaN. Detector-space + NIST blank subtraction is unchanged. + ### Resource-handle pass (2026-08-17) - Added `saxsabs.io.detector_images` as the common copy-and-close detector diff --git a/README.md b/README.md index db2e374..04c2f56 100644 --- a/README.md +++ b/README.md @@ -83,7 +83,7 @@ a display server. | Route | Best for | Start here | | --- | --- | --- | -| **CLI utilities** | normalization, header and 1D parsing, gated K estimation, gated buffer subtraction | `saxsabs --help` | +| **CLI utilities** | normalization, header and 1D parsing, gated K estimation, gated buffer and fluorescence subtraction | `saxsabs --help` | | **SAXSAbs Workbench** | interactive K calibration, batch processing, external-1D scaling | `saxsabs-workbench --lang en` | | **Strict BL19B2 runner** | validated campaign inputs under current BL19B2 conventions | [batch runbook](docs/bl19b2_abs2d_batch_runbook.md) | | **Python API** | reusable scientific calculations and file I/O | [API reference](docs/api.md) | diff --git a/SASAbs.py b/SASAbs.py index 91c6e8c..c01c200 100644 --- a/SASAbs.py +++ b/SASAbs.py @@ -351,6 +351,14 @@ def _read_package_version() -> str: "lbl_t3_buffer_file": "Buffer 1D file:", "lbl_t3_alpha": "\u03b1 (scale):", "lbl_t3_buffer_status": "(not loaded)", + "lf_t3_fluo": "Fluorescence subtraction", + "cb_t3_fluo_enable": "Enable fluorescence subtraction", + "lbl_t3_fluo_method": "Method:", + "lbl_t3_fluo_f0": "F0 (cm\u207b\xb9):", + "lbl_t3_fluo_file": "Measured F(q) file:", + "lbl_t3_fluo_status": "(not loaded)", + "lf_t2_fluo": "Fluorescence subtraction (1D)", + "cb_t2_fluo_enable": "Enable fluorescence subtraction", "lbl_t2_alpha": "BG \u03b1-scale:", "cb_t2_buffer_enable": "Enable BG \u03b1-scaling", # --- Output format --- @@ -738,6 +746,14 @@ def _read_package_version() -> str: "lbl_t3_buffer_file": "缓冲液1D文件:", "lbl_t3_alpha": "\u03b1 (缩放):", "lbl_t3_buffer_status": "(未加载)", + "lf_t3_fluo": "荧光扣除", + "cb_t3_fluo_enable": "启用荧光扣除", + "lbl_t3_fluo_method": "方法:", + "lbl_t3_fluo_f0": "F0 (cm\u207b\xb9):", + "lbl_t3_fluo_file": "实测 F(q) 文件:", + "lbl_t3_fluo_status": "(未加载)", + "lf_t2_fluo": "荧光扣除（1D）", + "cb_t2_fluo_enable": "启用荧光扣除", "lbl_t2_alpha": "背景 \u03b1缩放:", "cb_t2_buffer_enable": "启用背景 \u03b1-缩放", # --- 输出格式 --- @@ -935,6 +951,13 @@ def _read_package_version() -> str: "Leave u(\u03b1) blank when unknown; combined uncertainty stays NaN and is " "never assumed to be zero." ), + "lbl_t3_fluo_f0_uncertainty": "u(F0), optional:", + "lbl_t3_fluo_beta": "\u03b2 (scale):", + "lbl_t3_fluo_window": "High-q window:", + "hint_t3_fluo": ( + "Opt-in additive correction on absolute cm^-1 after K/thickness and optional " + "buffer. High-q methods are valid only where elastic SAXS is negligible." + ), }) I18N["zh"].update({ "lbl_k_record_readonly": "\u53ea\u8bfb\uff1b\u7531 Tab1 \u751f\u6210\uff08\u987b\u901a\u8fc7\u9884\u68c0\uff09", @@ -1024,6 +1047,13 @@ def _read_package_version() -> str: "hint_t3_alpha_uncertainty": ( "u(\u03b1) \u672a\u77e5\u65f6\u7559\u7a7a\uff1b\u5408\u6210\u4e0d\u786e\u5b9a\u5ea6\u4fdd\u6301 NaN\uff0c\u7edd\u4e0d\u9ed8\u8ba4\u4e3a 0\u3002" ), + "lbl_t3_fluo_f0_uncertainty": "u(F0)\uff08\u53ef\u9009\uff09:", + "lbl_t3_fluo_beta": "\u03b2 (\u7f29\u653e):", + "lbl_t3_fluo_window": "\u9ad8 q \u7a97\u53e3:", + "hint_t3_fluo": ( + "\u7edd\u5bf9 cm^-1 \u4e0a\u7684\u53ef\u9009\u52a0\u6027\u6263\u9664\uff0c\u4f4d\u4e8e K/\u539a\u5ea6\u4e0e\u53ef\u9009 buffer \u4e4b\u540e\u3002" + "\u4ec5\u5f53\u9ad8 q \u7a97\u5185\u5f39\u6027 SAXS \u53ef\u5ffd\u7565\u65f6\uff0chigh-q \u65b9\u6cd5\u624d\u6210\u7acb\u3002" + ), }) try: @@ -1283,6 +1313,15 @@ def save_figure(fig, path, **_kwargs): except Exception: subtract_buffer = None +try: + from saxsabs.core.fluorescence_subtraction import ( + combine_sequential_standard_uncertainties, + subtract_fluorescence, + ) +except Exception: + subtract_fluorescence = None + combine_sequential_standard_uncertainties = None + try: from saxsabs.core.execution_policy import ( parse_run_policy, @@ -1325,11 +1364,13 @@ def should_skip_all_existing(existing_flags, policy): try: from saxsabs.core.intensity_state import ( require_absolute_input_for_buffer_subtraction, + require_absolute_input_for_fluorescence_subtraction, require_relative_input_for_absolute_scaling, serialize_correction_ledger, ) except ImportError: require_absolute_input_for_buffer_subtraction = None + require_absolute_input_for_fluorescence_subtraction = None require_relative_input_for_absolute_scaling = None serialize_correction_ledger = None @@ -3579,6 +3620,7 @@ def profile_operator_metadata( calibration_context=None, *, corrections_applied=None, + extra_corrections=None, k_factor=None, thickness_cm=None, thickness_source=None, @@ -3615,6 +3657,8 @@ def profile_operator_metadata( corrections.append("flat_field") else: corrections = list(corrections_applied) + if extra_corrections: + corrections = list(corrections) + list(extra_corrections) serialized_corrections = serialize_correction_ledger(corrections) if k_factor is None: k_var = getattr(self, "global_vars", {}).get("k_factor") @@ -3654,6 +3698,7 @@ def save_profile_table( run_policy=None, calibration_context=None, corrections_applied=None, + extra_corrections=None, combined_uncertainty=None, uncertainty_metadata=None, thickness_cm=None, @@ -3680,6 +3725,7 @@ def save_profile_table( profile_metadata = self.profile_operator_metadata( calibration_context, corrections_applied=corrections_applied, + extra_corrections=extra_corrections, thickness_cm=thickness_cm, thickness_source=thickness_source, ) @@ -3745,6 +3791,14 @@ def save_profile_table( "buffer_source_sha256", "buffer_alpha", "buffer_alpha_uncertainty", + "fluorescence_method", + "fluorescence_f0", + "fluorescence_f0_uncertainty", + "fluorescence_beta", + "fluorescence_beta_uncertainty", + "fluorescence_high_q_window", + "fluorescence_source_name", + "fluorescence_source_sha256", "uncertainty_model", "uncertainty_type", ): @@ -4225,6 +4279,16 @@ def init_tab2_batch(self): self.t2_instr_tol_pct = tk.DoubleVar(value=0.5) self.t2_alpha = tk.DoubleVar(value=1.0) self.t2_alpha_enabled = tk.BooleanVar(value=False) + self.t2_fluo_enabled = tk.BooleanVar(value=False) + self.t2_fluo_method = tk.StringVar(value="constant") + self.t2_fluo_f0 = tk.StringVar(value="") + self.t2_fluo_f0_uncertainty = tk.StringVar(value="") + self.t2_fluo_beta = tk.DoubleVar(value=1.0) + self.t2_fluo_beta_uncertainty = tk.StringVar(value="") + self.t2_fluo_qmin = tk.StringVar(value="") + self.t2_fluo_qmax = tk.StringVar(value="") + self.t2_fluo_path = tk.StringVar() + self.t2_fluo_status = tk.StringVar(value=self.tr("lbl_t3_fluo_status")) self.t2_output_format = tk.StringVar(value="tsv") self.t2_export_cal2d = tk.BooleanVar(value=False) self.t2_cal2d_dtype = tk.StringVar(value="float32") @@ -4575,6 +4639,38 @@ def _sync_pol_entry_state(): self._register_i18n_widget(lbl_a2, "lbl_t2_alpha") ttk.Entry(row_alpha_fmt, textvariable=self.t2_alpha, width=6).pack(side="left") + fluo2 = ttk.LabelFrame(c5, text=self.tr("lf_t2_fluo"), style="Group.TLabelframe") + fluo2.pack(fill="x", pady=(4, 0)) + self._register_i18n_widget(fluo2, "lf_t2_fluo") + cb_fluo2 = ttk.Checkbutton( + fluo2, text=self.tr("cb_t2_fluo_enable"), variable=self.t2_fluo_enabled + ) + cb_fluo2.pack(anchor="w") + self._register_i18n_widget(cb_fluo2, "cb_t2_fluo_enable") + row_fluo2 = ttk.Frame(fluo2) + row_fluo2.pack(fill="x") + ttk.Combobox( + row_fluo2, + textvariable=self.t2_fluo_method, + values=["constant", "high_q_mean", "high_q_median", "measured"], + width=16, + state="readonly", + ).pack(side="left", padx=3) + ttk.Entry(row_fluo2, textvariable=self.t2_fluo_f0, width=8).pack(side="left", padx=3) + ttk.Entry(row_fluo2, textvariable=self.t2_fluo_f0_uncertainty, width=8).pack( + side="left", padx=3 + ) + ttk.Entry(row_fluo2, textvariable=self.t2_fluo_beta, width=6).pack(side="left", padx=3) + ttk.Entry(row_fluo2, textvariable=self.t2_fluo_beta_uncertainty, width=6).pack( + side="left", padx=3 + ) + ttk.Entry(row_fluo2, textvariable=self.t2_fluo_qmin, width=6).pack(side="left", padx=2) + ttk.Entry(row_fluo2, textvariable=self.t2_fluo_qmax, width=6).pack(side="left", padx=2) + self.add_file_row( + fluo2, self.tr("lbl_t3_fluo_file"), self.t2_fluo_path, "*.dat *.txt *.csv *.xml" + ) + ttk.Label(fluo2, textvariable=self.t2_fluo_status, style="Hint.TLabel").pack(anchor="w") + row_fmt2 = ttk.Frame(c5) row_fmt2.pack(fill="x", pady=(2, 0)) lbl_ofmt2 = ttk.Label(row_fmt2, text=self.tr("lbl_output_format")) @@ -4711,6 +4807,15 @@ def _sync_pol_entry_state(): self.t2_instr_tol_pct, self.t2_alpha, self.t2_alpha_enabled, + self.t2_fluo_enabled, + self.t2_fluo_method, + self.t2_fluo_f0, + self.t2_fluo_f0_uncertainty, + self.t2_fluo_beta, + self.t2_fluo_beta_uncertainty, + self.t2_fluo_qmin, + self.t2_fluo_qmax, + self.t2_fluo_path, self.t2_output_format, self.t2_export_cal2d, self.t2_cal2d_dtype, @@ -4977,6 +5082,70 @@ def init_tab3_external_1d(self): self.add_hint(buf_frame, "hint_t3_alpha_uncertainty", wraplength=360) ttk.Label(buf_frame, textvariable=self.t3_buffer_status, style="Hint.TLabel").pack(anchor="w", padx=3) + self.t3_fluo_enabled = tk.BooleanVar(value=False) + self.t3_fluo_method = tk.StringVar(value="constant") + self.t3_fluo_f0 = tk.StringVar(value="") + self.t3_fluo_f0_uncertainty = tk.StringVar(value="") + self.t3_fluo_beta = tk.DoubleVar(value=1.0) + self.t3_fluo_beta_uncertainty = tk.StringVar(value="") + self.t3_fluo_qmin = tk.StringVar(value="") + self.t3_fluo_qmax = tk.StringVar(value="") + self.t3_fluo_path = tk.StringVar() + self.t3_fluo_status = tk.StringVar(value=self.tr("lbl_t3_fluo_status")) + + fluo_frame = ttk.LabelFrame(top, text=self.tr("lf_t3_fluo"), style="Group.TLabelframe") + fluo_frame.grid(row=2, column=0, columnspan=2, sticky="nsew", padx=5, pady=4) + self._register_i18n_widget(fluo_frame, "lf_t3_fluo") + cb_fluo = ttk.Checkbutton( + fluo_frame, text=self.tr("cb_t3_fluo_enable"), variable=self.t3_fluo_enabled + ) + cb_fluo.pack(anchor="w", padx=3, pady=2) + self._register_i18n_widget(cb_fluo, "cb_t3_fluo_enable") + row_fluo_m = ttk.Frame(fluo_frame) + row_fluo_m.pack(fill="x", pady=1) + lbl_fluo_m = ttk.Label(row_fluo_m, text=self.tr("lbl_t3_fluo_method")) + lbl_fluo_m.pack(side="left") + self._register_i18n_widget(lbl_fluo_m, "lbl_t3_fluo_method") + self.t3_fluo_method_combo = ttk.Combobox( + row_fluo_m, + textvariable=self.t3_fluo_method, + values=["constant", "high_q_mean", "high_q_median", "measured"], + width=16, + state="readonly", + ) + self.t3_fluo_method_combo.pack(side="left", padx=5) + lbl_fluo_f0 = ttk.Label(row_fluo_m, text=self.tr("lbl_t3_fluo_f0")) + lbl_fluo_f0.pack(side="left", padx=(8, 0)) + self._register_i18n_widget(lbl_fluo_f0, "lbl_t3_fluo_f0") + ttk.Entry(row_fluo_m, textvariable=self.t3_fluo_f0, width=8).pack(side="left", padx=5) + lbl_fluo_f0u = ttk.Label(row_fluo_m, text=self.tr("lbl_t3_fluo_f0_uncertainty")) + lbl_fluo_f0u.pack(side="left") + self._register_i18n_widget(lbl_fluo_f0u, "lbl_t3_fluo_f0_uncertainty") + ttk.Entry(row_fluo_m, textvariable=self.t3_fluo_f0_uncertainty, width=8).pack( + side="left", padx=5 + ) + row_fluo_b = ttk.Frame(fluo_frame) + row_fluo_b.pack(fill="x", pady=1) + lbl_fluo_b = ttk.Label(row_fluo_b, text=self.tr("lbl_t3_fluo_beta")) + lbl_fluo_b.pack(side="left") + self._register_i18n_widget(lbl_fluo_b, "lbl_t3_fluo_beta") + ttk.Entry(row_fluo_b, textvariable=self.t3_fluo_beta, width=8).pack(side="left", padx=5) + ttk.Entry(row_fluo_b, textvariable=self.t3_fluo_beta_uncertainty, width=8).pack( + side="left", padx=5 + ) + lbl_fluo_w = ttk.Label(row_fluo_b, text=self.tr("lbl_t3_fluo_window")) + lbl_fluo_w.pack(side="left", padx=(8, 0)) + self._register_i18n_widget(lbl_fluo_w, "lbl_t3_fluo_window") + ttk.Entry(row_fluo_b, textvariable=self.t3_fluo_qmin, width=8).pack(side="left", padx=2) + ttk.Entry(row_fluo_b, textvariable=self.t3_fluo_qmax, width=8).pack(side="left", padx=2) + self.add_file_row( + fluo_frame, self.tr("lbl_t3_fluo_file"), self.t3_fluo_path, "*.dat *.txt *.csv *.xml" + ) + self.add_hint(fluo_frame, "hint_t3_fluo", wraplength=720) + ttk.Label(fluo_frame, textvariable=self.t3_fluo_status, style="Hint.TLabel").pack( + anchor="w", padx=3 + ) + # ---- Output format selector ---- self.t3_output_format = tk.StringVar(value="tsv") fmt_row = ttk.Frame(buf_frame) @@ -5075,6 +5244,15 @@ def init_tab3_external_1d(self): self.t3_buffer_path, self.t3_alpha, self.t3_alpha_uncertainty, + self.t3_fluo_enabled, + self.t3_fluo_method, + self.t3_fluo_f0, + self.t3_fluo_f0_uncertainty, + self.t3_fluo_beta, + self.t3_fluo_beta_uncertainty, + self.t3_fluo_qmin, + self.t3_fluo_qmax, + self.t3_fluo_path, self.t3_output_format, ], ) @@ -5213,6 +5391,7 @@ def require_relative_external_profile_for_scaling( *, correction_mode=None, apply_buffer=False, + apply_fluorescence=False, ): """Validate the exact Tab3 operators before changing an external profile.""" if require_relative_input_for_absolute_scaling is None: @@ -5231,6 +5410,8 @@ def require_relative_external_profile_for_scaling( raise ValueError(f"{profile_name}: unknown correction mode: {mode}") if apply_buffer: corrections_to_apply.append("buffer") + if apply_fluorescence: + corrections_to_apply.append("fluorescence") assessment = require_relative_input_for_absolute_scaling( profile, profile_name=profile_name, @@ -5763,6 +5944,278 @@ def subtract_external_absolute_buffer( result.err_subtracted, ) + @staticmethod + def disabled_fluorescence_payload(): + return { + "enabled": False, + "method": None, + "path": "", + "sha256": None, + "f0": None, + "f0_uncertainty": None, + "beta": None, + "beta_uncertainty": None, + "high_q_window": None, + "profile": None, + } + + @staticmethod + def external_fluorescence_audit_payload(fluo_info): + return { + "enabled": bool(fluo_info.get("enabled", False)), + "method": fluo_info.get("method"), + "path": str(fluo_info.get("path") or ""), + "sha256": fluo_info.get("sha256"), + "f0": fluo_info.get("f0"), + "f0_uncertainty": fluo_info.get("f0_uncertainty"), + "beta": fluo_info.get("beta"), + "beta_uncertainty": fluo_info.get("beta_uncertainty"), + "high_q_window": fluo_info.get("high_q_window"), + } + + @staticmethod + def parse_optional_nonnegative_uncertainty(value, *, label): + text = str(value if value is not None else "").strip() + if not text: + return None + try: + uncertainty = float(text) + except (TypeError, ValueError) as exc: + raise ValueError(f"{label} must be numeric") from exc + if not np.isfinite(uncertainty) or uncertainty < 0: + raise ValueError(f"{label} must be finite and >= 0") + return uncertainty + + @staticmethod + def parse_optional_window(qmin_value, qmax_value): + qmin_text = str(qmin_value if qmin_value is not None else "").strip() + qmax_text = str(qmax_value if qmax_value is not None else "").strip() + if not qmin_text and not qmax_text: + return None + if not qmin_text or not qmax_text: + raise ValueError("high-q window requires both qmin and qmax") + qmin = float(qmin_text) + qmax = float(qmax_text) + if not np.isfinite(qmin) or not np.isfinite(qmax) or qmin >= qmax: + raise ValueError("high-q window must contain two finite increasing values") + return (qmin, qmax) + + def prepare_workbench_fluorescence( + self, + *, + source="t3", + pipeline_mode="scaled", + calibration_context=None, + k_factor=None, + require_scaled_pipeline=True, + ): + """Load optional fluorescence settings; default is disabled.""" + status_var = getattr(self, f"{source}_fluo_status", None) + + def set_status(value): + if status_var is not None: + try: + status_var.set(value) + except Exception: + pass + + enabled_var = getattr(self, f"{source}_fluo_enabled", None) + enabled = bool(enabled_var.get()) if enabled_var is not None else False + if not enabled: + set_status("Fluorescence: disabled") + return self.disabled_fluorescence_payload() + + try: + if require_scaled_pipeline and str(pipeline_mode or "").strip().lower() != "scaled": + raise ValueError( + "raw 流程下禁止荧光扣除；荧光必须作用在绝对强度标度上。" + ) + method = str(getattr(self, f"{source}_fluo_method").get()).strip().lower() + beta = float(getattr(self, f"{source}_fluo_beta").get()) + if not np.isfinite(beta) or beta <= 0: + raise ValueError("Fluorescence beta must be finite and > 0") + f0_text = str(getattr(self, f"{source}_fluo_f0").get()).strip() + f0 = float(f0_text) if f0_text else None + f0_uncertainty = self.parse_optional_nonnegative_uncertainty( + getattr(self, f"{source}_fluo_f0_uncertainty").get(), + label="Fluorescence F0 uncertainty", + ) + beta_uncertainty = self.parse_optional_nonnegative_uncertainty( + getattr(self, f"{source}_fluo_beta_uncertainty").get(), + label="Fluorescence beta uncertainty", + ) + high_q_window = self.parse_optional_window( + getattr(self, f"{source}_fluo_qmin").get(), + getattr(self, f"{source}_fluo_qmax").get(), + ) + path_text = str(getattr(self, f"{source}_fluo_path").get()).strip() + profile = None + digest = None + resolved_path = "" + if method == "measured": + if calibration_context is None: + raise ValueError("Measured fluorescence requires a valid CalibrationContext.") + active_k = float(k_factor) + if not np.isfinite(active_k) or active_k <= 0: + raise ValueError("Measured fluorescence requires finite positive active K") + if not path_text: + raise ValueError("已启用 measured 荧光扣除，但未选择荧光曲线。") + if require_absolute_input_for_fluorescence_subtraction is None: + raise RuntimeError("fluorescence-state validation is unavailable") + fluo_path = Path(path_text).expanduser().resolve() + profile = self.prepare_external_profile_axis( + fluo_path, self.read_external_1d_profile(fluo_path) + ) + require_absolute_input_for_fluorescence_subtraction( + profile, profile_name="Fluorescence" + ) + self.require_external_profile_operator_provenance( + profile, + calibration_context, + "Fluorescence", + require_full_context_fingerprint=True, + required_k_factor=active_k, + ) + digest = self._optional_file_sha256(fluo_path) + resolved_path = str(fluo_path) + set_status(f"Fluorescence loaded: {fluo_path.name} ({len(profile['x'])} points)") + else: + set_status(f"Fluorescence method: {method}") + return { + "enabled": True, + "method": method, + "path": resolved_path, + "sha256": digest, + "f0": f0, + "f0_uncertainty": f0_uncertainty, + "beta": beta, + "beta_uncertainty": beta_uncertainty, + "high_q_window": high_q_window, + "profile": profile, + } + except Exception as exc: + set_status(f"Fluorescence error: {exc}") + raise + + @staticmethod + def fluorescence_uncertainty_metadata(fluo_info, result): + model = "u_combined^2=u_sample^2+beta^2*u_F^2+F^2*u_beta^2" + unknown = ( + fluo_info.get("beta_uncertainty") is None + or ( + str(fluo_info.get("method") or "") not in {"measured", "measured_profile"} + and fluo_info.get("f0_uncertainty") is None + ) + ) + payload = { + "fluorescence_method": result.method, + "fluorescence_f0": repr(float(result.f0)), + "fluorescence_f0_uncertainty": ( + "unknown" + if result.f0_uncertainty is None + else repr(float(result.f0_uncertainty)) + ), + "fluorescence_beta": repr(float(result.beta)), + "fluorescence_beta_uncertainty": ( + "unknown" + if result.beta_uncertainty is None + else repr(float(result.beta_uncertainty)) + ), + "uncertainty_model": model, + "uncertainty_type": ( + "combined_standard_unknown_fluorescence" + if unknown + else "combined_standard" + ), + } + if result.high_q_window is not None: + payload["fluorescence_high_q_window"] = ( + f"{result.high_q_window[0]:.17g},{result.high_q_window[1]:.17g}" + ) + if fluo_info.get("path"): + payload["fluorescence_source_name"] = Path(fluo_info["path"]).name + payload["fluorescence_source_sha256"] = fluo_info.get("sha256") + return payload + + @staticmethod + def merge_uncertainty_metadata(existing, extra): + merged = dict(existing or {}) + extra = dict(extra or {}) + old_model = merged.get("uncertainty_model") + new_model = extra.get("uncertainty_model") + if old_model and new_model and old_model != new_model: + extra["uncertainty_model"] = f"{old_model};{new_model}" + previous_type = str(merged.get("uncertainty_type") or "") + extra_type = str(extra.get("uncertainty_type") or "") + merged.update(extra) + if "unknown" in previous_type or "unknown" in extra_type: + merged["uncertainty_type"] = "combined_standard_unknown" + return merged + + @staticmethod + def subtract_external_absolute_fluorescence( + sample_q, + sample_i, + sample_err, + fluo_info, + sample_profile=None, + ): + """Apply the single audited fluorescence kernel; never use a weaker fallback.""" + + if subtract_fluorescence is None: + raise RuntimeError("formal fluorescence subtraction kernel is unavailable") + if sample_profile is None: + sample_profile = { + "intensity_state": "absolute_cm^-1", + "intensity_unit": "1/cm", + "i_col": "I_abs_cm^-1", + "operator_provenance": { + "intensity_state": "absolute_cm^-1", + "corrections_applied": '["k","thickness"]', + }, + } + fluo_profile = fluo_info.get("profile") + q_f = i_f = e_f = None + if fluo_profile is not None: + q_f = np.asarray(fluo_profile["x"], dtype=np.float64) + i_f = _profile_intensity(fluo_profile) + e_f = _profile_uncertainty(fluo_profile) + result = subtract_fluorescence( + np.asarray(sample_q, dtype=np.float64), + np.asarray(sample_i, dtype=np.float64), + np.asarray(sample_err, dtype=np.float64), + sample_profile=sample_profile, + method=fluo_info["method"], + f0=fluo_info.get("f0"), + f0_uncertainty=fluo_info.get("f0_uncertainty"), + beta=fluo_info.get("beta") if fluo_info.get("beta") is not None else 1.0, + beta_uncertainty=fluo_info.get("beta_uncertainty"), + high_q_window=fluo_info.get("high_q_window"), + q_fluorescence=q_f, + i_fluorescence=i_f, + err_fluorescence=e_f, + fluorescence_profile=fluo_profile, + ) + if result.err_statistical is None: + raise RuntimeError("fluorescence kernel did not return statistical uncertainty") + return result + + def apply_workbench_fluorescence_to_profile(self, q, i_abs, i_err, fluo_info): + """Apply optional fluorescence on an absolute 1-D profile.""" + + if not fluo_info or not fluo_info.get("enabled"): + return i_abs, i_err, None, None, () + result = self.subtract_external_absolute_fluorescence( + q, i_abs, i_err, fluo_info + ) + return ( + result.i_subtracted, + result.err_statistical, + result.err_subtracted, + self.fluorescence_uncertainty_metadata(fluo_info, result), + ("fluorescence",), + ) + def _regularize_xy_triplet(self, x, y, e=None, min_points=3, name="profile"): x = np.asarray(x, dtype=np.float64) y = np.asarray(y, dtype=np.float64) @@ -6407,7 +6860,10 @@ def dry_run_external_1d(self): buffer_var = getattr(self, "t3_buffer_enabled", None) buffer_enabled = bool(buffer_var.get()) if buffer_var is not None else False + fluo_var = getattr(self, "t3_fluo_enabled", None) + fluorescence_enabled = bool(fluo_var.get()) if fluo_var is not None else False buffer_info = {"enabled": False, "profile": None} + fluorescence_gate_error = None try: buffer_info = self.prepare_external_buffer( pipeline_mode=pipeline_mode, @@ -6418,6 +6874,17 @@ def dry_run_external_1d(self): buffer_gate_error = str(exc) warnings.append(buffer_gate_error) + try: + self.prepare_workbench_fluorescence( + source="t3", + pipeline_mode=pipeline_mode, + calibration_context=active_calibration_context, + k_factor=k, + ) + except (OSError, RuntimeError, TypeError, ValueError) as exc: + fluorescence_gate_error = str(exc) + warnings.append(fluorescence_gate_error) + meta_map = {} bg_prof = None dark_prof = None @@ -6478,6 +6945,7 @@ def dry_run_external_1d(self): Path(fp).name, correction_mode=mode, apply_buffer=buffer_enabled, + apply_fluorescence=fluorescence_enabled, ) x_label = prof["x_label"] x_conversion = prof["x_conversion"] @@ -6582,6 +7050,7 @@ def dry_run_external_1d(self): k_trust_error, thickness_gate_error, buffer_gate_error, + fluorescence_gate_error, resume_gate_error, ) ): @@ -6672,6 +7141,15 @@ def run_external_1d_batch(self): k_factor=k, ) buffer_audit = self.external_buffer_audit_payload(buffer_info) + fluorescence_info = self.prepare_workbench_fluorescence( + source="t3", + pipeline_mode=pipeline_mode, + calibration_context=active_calibration_context, + k_factor=k, + ) + fluorescence_audit = self.external_fluorescence_audit_payload( + fluorescence_info + ) meta_map = {} bg_prof = None @@ -6748,6 +7226,8 @@ def run_external_1d_batch(self): x_conversion = "" operator_fingerprint = "" buffer_applied = False + fluorescence_applied = False + fluorescence_f0_applied = None scale_factor = scale_factor_global if pipeline_mode == "scaled" else np.nan thk_cm_used = fixed_thk_cm if pipeline_mode == "scaled" else np.nan norm_s = np.nan @@ -6766,6 +7246,7 @@ def run_external_1d_batch(self): Path(fp).name, correction_mode=corr_mode, apply_buffer=bool(buffer_info["enabled"]), + apply_fluorescence=bool(fluorescence_info["enabled"]), ) points = len(prof["x"]) x_label = prof["x_label"] @@ -6909,6 +7390,53 @@ def run_external_1d_batch(self): ), } + if fluorescence_info["enabled"]: + fluo_ledger = list(input_state.corrections_applied) + ["k"] + if corr_mode == "k_over_d": + fluo_ledger.append("thickness") + if buffer_info["enabled"]: + fluo_ledger.append("buffer") + fluo_sample_profile = { + "intensity_state": "absolute_cm^-1", + "intensity_unit": "1/cm", + "i_col": "I_abs_cm^-1", + "operator_provenance": { + "intensity_state": "absolute_cm^-1", + "corrections_applied": serialize_correction_ledger( + fluo_ledger + ), + }, + } + previous_statistical = err_abs + previous_combined = combined_uncertainty + result = self.subtract_external_absolute_fluorescence( + prof["x"], + i_abs, + previous_statistical, + fluorescence_info, + sample_profile=fluo_sample_profile, + ) + i_abs = result.i_subtracted + if combine_sequential_standard_uncertainties is None: + raise RuntimeError( + "fluorescence uncertainty composition helper is unavailable" + ) + err_abs, combined_uncertainty = ( + combine_sequential_standard_uncertainties( + previous_statistical, + previous_combined, + result.err_statistical, + result.err_subtracted, + ) + ) + fluorescence_f0_applied = result.f0 + uncertainty_metadata = self.merge_uncertainty_metadata( + uncertainty_metadata, + self.fluorescence_uncertainty_metadata( + fluorescence_info, result + ), + ) + output_corrections = list(input_state.corrections_applied) output_corrections.append("k") if corr_mode == "k_over_d": @@ -6919,6 +7447,8 @@ def run_external_1d_batch(self): ) if buffer_info["enabled"]: output_corrections.append("buffer") + if fluorescence_info["enabled"]: + output_corrections.append("fluorescence") corrections_applied_serialized = serialize_correction_ledger( output_corrections ) @@ -6939,6 +7469,8 @@ def run_external_1d_batch(self): ) if buffer_info["enabled"]: buffer_applied = True + if fluorescence_info["enabled"]: + fluorescence_applied = True status = "成功" outputs = written_path.name ok += 1 @@ -6985,6 +7517,23 @@ def run_external_1d_batch(self): "BufferCorrectionsApplied": serialize_correction_ledger( buffer_audit["corrections_applied"] ), + "FluorescenceEnabled": fluorescence_audit["enabled"], + "FluorescenceApplied": fluorescence_applied, + "FluorescenceMethod": fluorescence_audit["method"] or "", + "FluorescenceF0": ( + fluorescence_f0_applied + if fluorescence_f0_applied is not None + else ( + fluorescence_audit["f0"] + if fluorescence_audit["f0"] is not None + else np.nan + ) + ), + "FluorescenceBeta": ( + fluorescence_audit["beta"] + if fluorescence_audit["beta"] is not None + else np.nan + ), "PipelineMode": pipeline_mode, "CorrMode": corr_mode, "K": k, @@ -7038,6 +7587,7 @@ def run_external_1d_batch(self): "output_dir": str(out_dir), "report_csv": str(report_path), "buffer": buffer_audit, + "fluorescence": fluorescence_audit, "summary": {"success": ok, "skipped": skip, "failed": fail}, } meta_path = report_dir / f"external1d_meta_{stamp}.json" @@ -8624,6 +9174,15 @@ def load_data(path): i_err = np.full_like(i_abs, np.nan) else: i_err = np.asarray(res.sigma, dtype=np.float64) * scale_factor + ( + i_abs, + i_err, + combined_uncertainty, + fluo_meta, + extra_corrections, + ) = self.apply_workbench_fluorescence_to_profile( + res.radial, i_abs, i_err, context.get("fluorescence") + ) issue = self.profile_health_issue(i_abs) if issue: raise ValueError(issue) @@ -8635,6 +9194,9 @@ def load_data(path): "Q_A^-1", output_format=output_format, run_policy=run_policy, + extra_corrections=extra_corrections, + combined_uncertainty=combined_uncertainty, + uncertainty_metadata=fluo_meta, ) outputs.append(f"{mode}:{written_path.name}") mode_stats[mode]["ok"] += 1 @@ -8707,6 +9269,15 @@ def load_data(path): i_err = np.full_like(i_abs, np.nan) else: i_err = np.asarray(res.sigma, dtype=np.float64) * scale_factor + ( + i_abs, + i_err, + combined_uncertainty, + fluo_meta, + extra_corrections, + ) = self.apply_workbench_fluorescence_to_profile( + res.radial, i_abs, i_err, context.get("fluorescence") + ) issue = self.profile_health_issue(i_abs) if issue: raise ValueError(issue) @@ -8718,6 +9289,9 @@ def load_data(path): "Q_A^-1", output_format=output_format, run_policy=run_policy, + extra_corrections=extra_corrections, + combined_uncertainty=combined_uncertainty, + uncertainty_metadata=fluo_meta, ) written_disp = ( f"{written_path.parent.name}/{written_path.name}" @@ -8754,6 +9328,15 @@ def load_data(path): i_err = np.full_like(i_abs, np.nan) else: i_err = np.asarray(merge.sigma, dtype=np.float64) * scale_factor + ( + i_abs, + i_err, + combined_uncertainty, + fluo_meta, + extra_corrections, + ) = self.apply_workbench_fluorescence_to_profile( + merge.radial, i_abs, i_err, context.get("fluorescence") + ) issue = self.profile_health_issue(i_abs) if issue: raise ValueError(issue) @@ -8765,6 +9348,9 @@ def load_data(path): "Q_A^-1", output_format=output_format, run_policy=run_policy, + extra_corrections=extra_corrections, + combined_uncertainty=combined_uncertainty, + uncertainty_metadata=fluo_meta, ) outputs.append(f"1d_sector_sum:{written_path.name}") mode_stats[mode]["ok"] += 1 @@ -9233,6 +9819,13 @@ def run_batch(self): "resume": resume, "run_policy": run_policy, "bg_alpha": float(self.t2_alpha.get()) if self.t2_alpha_enabled.get() else 1.0, + "fluorescence": self.prepare_workbench_fluorescence( + source="t2", + pipeline_mode="scaled", + calibration_context=active_calibration_context, + k_factor=float(self.global_vars["k_factor"].get()), + require_scaled_pipeline=False, + ), "output_format": self.t2_output_format.get() if hasattr(self, "t2_output_format") else "tsv", "export_cal2d": export_cal2d, "cal2d_root": cal2d_root, @@ -9692,7 +10285,10 @@ def _t2_preflight_config(self): "t2_polarization", "t2_output_root", "t2_mask_path", "t2_flat_path", "t2_resume_enabled", "t2_overwrite", "t2_workers", "t2_strict_instrument", "t2_instr_tol_pct", "t2_alpha", - "t2_alpha_enabled", "t2_output_format", "t2_export_cal2d", + "t2_alpha_enabled", "t2_fluo_enabled", "t2_fluo_method", + "t2_fluo_f0", "t2_fluo_f0_uncertainty", "t2_fluo_beta", + "t2_fluo_beta_uncertainty", "t2_fluo_qmin", "t2_fluo_qmax", + "t2_fluo_path", "t2_output_format", "t2_export_cal2d", "t2_cal2d_dtype", "t2_cal2d_apply_flat", "t2_mode_full", "t2_mode_sector", "t2_mode_chi", "t2_sec_min", "t2_sec_max", "t2_sector_ranges_text", "t2_sector_save_each", @@ -9737,14 +10333,23 @@ def _t3_preflight_config(self): "t3_sample_exp", "t3_sample_i0", "t3_sample_t", "t3_bg_exp", "t3_bg_i0", "t3_bg_t", "t3_sync_bg_from_global", "t3_resume_enabled", "t3_overwrite", "t3_buffer_enabled", - "t3_buffer_path", "t3_alpha", "t3_alpha_uncertainty", "t3_output_format", + "t3_buffer_path", "t3_alpha", "t3_alpha_uncertainty", + "t3_fluo_enabled", "t3_fluo_method", "t3_fluo_f0", + "t3_fluo_f0_uncertainty", "t3_fluo_beta", "t3_fluo_beta_uncertainty", + "t3_fluo_qmin", "t3_fluo_qmax", "t3_fluo_path", "t3_output_format", ) config = { name: self._preflight_var_value(getattr(self, name, None)) for name in names } files = list(dict.fromkeys(str(item) for item in getattr(self, "t3_files", []))) - for name in ("t3_meta_csv_path", "t3_bg1d_path", "t3_dark1d_path", "t3_buffer_path"): + for name in ( + "t3_meta_csv_path", + "t3_bg1d_path", + "t3_dark1d_path", + "t3_buffer_path", + "t3_fluo_path", + ): config[f"{name}_identity"] = self._preflight_file_identity(config[name]) config.update({ "schema": "saxsabs-workbench-tab3-preflight-v1", diff --git a/docs/api.md b/docs/api.md index 4d2fbf0..c4110a8 100644 --- a/docs/api.md +++ b/docs/api.md @@ -93,6 +93,11 @@ subtract_buffer(q_sample, i_sample, err_sample, q_buffer, i_buffer, err_buffer, *, alpha_uncertainty: float | None = None, sample_profile: Mapping[str, object], buffer_profile: Mapping[str, object]) -> BufferSubtractionResult +subtract_fluorescence(q, i_abs, err_abs, *, sample_profile, method, + f0=None, f0_uncertainty=None, beta=1.0, beta_uncertainty=None, + high_q_window=None, q_fluorescence=None, i_fluorescence=None, + err_fluorescence=None, fluorescence_profile=None, + residual_window=None) -> FluorescenceSubtractionResult propagate_absolute_uncertainty(intensity: np.ndarray, *, statistical_standard_uncertainty=None, k_relative_standard_uncertainty=None, standard_relative_standard_uncertainty=None, @@ -108,6 +113,12 @@ conflicting evidence remains ambiguous. A unitless metadata label `absolute` is ambiguous and is not treated as cm$^{-1}$. `subtract_buffer` requires `sample_profile` and `buffer_profile` provenance, interpolates a buffer onto the sample q grid when necessary, and propagates supplied uncertainties. +`subtract_fluorescence` is an opt-in additive correction on absolute cm$^{-1}$ +data after K, thickness, and optional buffer. Methods are `constant`, +`high_q_mean`, `high_q_median`, and `measured_profile`. A high-q constant is +valid only where elastic SAXS is negligible in that window. Missing `u(F0)` or +`u(β)` keeps combined uncertainty NaN. Negative intensities are reported, not +clipped. This is not detector-dark, NIST-blank, or solvent subtraction. `propagate_absolute_uncertainty` combines statistical and supplied standard uncertainty components; relative inputs must be relative standard uncertainties. @@ -162,6 +173,10 @@ saxsabs estimate-k --meas [--ref ] [--q-col ] [--i-col [--intensity-state relative] [--thickness-cm ] saxsabs subtract-buffer --sample --buffer [--alpha ] [--alpha-uncertainty ] +saxsabs subtract-fluorescence --sample --method + [--f0 ] [--f0-uncertainty ] [--beta ] + [--beta-uncertainty ] [--qmin <Å^-1>] [--qmax <Å^-1>] + [--fluorescence ] saxsabs bl19b2-abs2d --input-root (--poni |--pydidas-cali-yaml ) (--mu |--sample-thickness-cm ) --monitor-mode [workflow options] @@ -178,6 +193,7 @@ The main commands are: | `parse-external1d` | profile path | parsed-profile summary JSON | | `estimate-k` | relative measured profile; optional reference (built-in SRM 3600 if omitted); optional `--thickness-cm` | K-factor result JSON | | `subtract-buffer` | absolute sample and buffer profiles with cm⁻¹ units | subtraction diagnostic JSON | +| `subtract-fluorescence` | absolute sample profile; constant, high-q window, or measured F(q) | fluorescence diagnostic JSON | | `bl19b2-abs2d` | explicit BL19B2 inputs and semantics | batch result JSON and requested files | | `bl19b2-abs2d-v1-legacy` | explicit migration choices | legacy-compatible batch result with documented assumptions | diff --git a/docs/architecture.md b/docs/architecture.md index d7813b5..157810e 100644 --- a/docs/architecture.md +++ b/docs/architecture.md @@ -17,6 +17,7 @@ │ saxsabs.core.intensity_state (1D ledger) │ │ saxsabs.core.workbench_preflight_gate │ │ saxsabs.core.buffer_subtraction (BioSAXS) │ +│ saxsabs.core.fluorescence_subtraction (1D) │ │ saxsabs.workflows.bl19b2_abs2d / integrate1d │ │ saxsabs.io.parsers (header + 1D I/O) │ │ saxsabs.io.writers (canSAS/NXcanSAS) │ @@ -31,17 +32,19 @@ - **`src/saxsabs/core`**: pure computation logic (normalization, robust K-factor estimation, fingerprinted NIST 30 keV material attenuation, xraydb/Elam diagnostic attenuation, 1D intensity/correction state, Workbench - preflight fingerprints, and buffer subtraction). No GUI side-effects — + preflight fingerprints, buffer subtraction, and optional 1D fluorescence + subtraction). No GUI side-effects — deterministic and testable. - **`src/saxsabs/io`**: robust input parsing plus standard-format writers (canSAS XML and NXcanSAS HDF5). -- **`src/saxsabs/cli.py`**: seven headless subcommands: five focused utilities +- **`src/saxsabs/cli.py`**: eight headless subcommands: six focused utilities (`norm-factor`, `parse-header`, `parse-external1d`, `estimate-k`, - `subtract-buffer`), the safety-first `bl19b2-abs2d` workflow, and the - explicit `bl19b2-abs2d-v1-legacy` migration entry. `estimate-k` and - `subtract-buffer` apply the intensity-state gates; the first three utilities - remain thin parsers/calculators. The legacy entry requires explicit monitor - and thickness semantics and never silently restores v1 defaults. + `subtract-buffer`, `subtract-fluorescence`), the safety-first `bl19b2-abs2d` + workflow, and the explicit `bl19b2-abs2d-v1-legacy` migration entry. + `estimate-k`, `subtract-buffer`, and `subtract-fluorescence` apply the + intensity-state gates; the first three utilities remain thin + parsers/calculators. The legacy entry requires explicit monitor and + thickness semantics and never silently restores v1 defaults. - **`src/saxsabs/constants.py`**: pluggable reference-standard registry (SRM 3600, water, custom curves). - **`src/saxsabs/workbench_launcher.py`**: packaged launcher used by @@ -99,6 +102,14 @@ provenance. The Workbench calls the shared core `subtract_buffer`; if it is unavailable, formal subtraction fails closed with no weaker local fallback. +- Fluorescence subtraction is opt-in and absolute-scale only. It is an additive + 1-D term `I_corr = I_abs − β F(q)` applied after K, thickness, and optional + buffer. Methods are `constant`, `high_q_mean`, `high_q_median`, and + `measured_profile`. High-q estimation is valid only where elastic SAXS is + negligible in the stated window. The ledger token is `fluorescence`. Unknown + `u(F0)` or `u(β)` keeps combined uncertainty NaN. The kernel never clips + negative intensities. Detector-space NIST blank subtraction is unchanged: + a q-independent 1-D constant must not be written back onto raw counts. - The NIST 30 keV GUI export is a material-attenuation provenance JSON. It is invalidated whenever source/energy/preset/composition/density/porosity input changes. Nominal identity is inferred from the edited composition, not copied @@ -139,9 +150,10 @@ estimation, NIST 30 keV material core, Elam diagnostic calculator, 1D intensity ledger, signed-in-memory Workbench preflight, fixed-thickness enforcement, disabled legacy/resume controls, exact K-only/Kd/buffer gates, - absolute-buffer validation, provenance-aware scrollable μ UI, disabled Tab 3 - raw mode, screen-aware startup, strict BL19B2 workflows, standard writers, - bilingual GUI, CLI, CI, and paper assets. K-only formal scaling requires both + absolute-buffer validation, optional absolute 1D fluorescence subtraction, + provenance-aware scrollable μ UI, disabled Tab 3 raw mode, screen-aware + startup, strict BL19B2 workflows, standard writers, bilingual GUI, CLI, CI, + and paper assets. K-only formal scaling requires both the inherited-thickness ledger entry and numeric/source provenance. - **Strict campaign ownership**: formal multi-folder and per-sample campaigns remain owned by the strict CLI/batch runner. The Workbench is an interactive diff --git a/docs/reviewer-faq.md b/docs/reviewer-faq.md index 4491d9b..0fc7fd2 100644 --- a/docs/reviewer-faq.md +++ b/docs/reviewer-faq.md @@ -7,6 +7,14 @@ entry point. Reusable scientific and I/O logic lives under `src/saxsabs/`; the GUI remains separate because the current Workbench and strict BL19B2 campaign runner have intentionally different ownership boundaries. +## Is fluorescence the same as empty-cell or buffer subtraction? + +No. Empty-cell / NIST-blank subtraction happens in detector space before +absolute scaling. Buffer subtraction removes solvent on the absolute +cm$^{-1}$ scale. Fluorescence subtraction is an optional additive 1-D term +applied after K, thickness, and optional buffer. It is off by default, does +not change $K$, and is not written onto 2D raw counts. + ## How can this be tested without GUI? Core logic is exposed as importable APIs and CLI commands. Tests run headlessly in CI. diff --git a/examples/manual-verification.md b/examples/manual-verification.md index 46f863c..d0a5709 100644 --- a/examples/manual-verification.md +++ b/examples/manual-verification.md @@ -165,7 +165,16 @@ front end to the strict BL19B2 campaign runner. Negative, non-finite, or malformed values must fail closed. Temporarily make the shared core kernel unavailable and confirm formal subtraction fails closed rather than using a weaker fallback. -9. Before packaging the repository, confirm `git status` contains no audit +9. Enable fluorescence subtraction on an absolute Tab 3 (or Tab 2 1-D) result. + Confirm it is refused on relative/raw profiles, refused if `fluorescence` is + already in `corrections_applied`, and applied only after K/thickness and + optional buffer. For `constant`, set a planted F0 and confirm + `I_corr = I_abs − F0`. Leave `u(F0)` blank and confirm combined uncertainty + is NaN; set a finite value and confirm it enters the combined column. For + `high_q_mean`, a window with fewer than 3 finite points must fail closed. + Confirm `saxsabs subtract-fluorescence --help` lists the same methods, and + that making the shared kernel unavailable fails closed with no GUI fallback. +10. Before packaging the repository, confirm `git status` contains no audit outputs, build caches, downloaded literature, or private drive roots. Keep manual evidence outside the repository and use only anonymized, portable fixtures for version-controlled examples and tests. diff --git a/src/saxsabs/__init__.py b/src/saxsabs/__init__.py index 87d42ad..b0cfd70 100644 --- a/src/saxsabs/__init__.py +++ b/src/saxsabs/__init__.py @@ -24,9 +24,16 @@ IntensityStateAssessment, assess_intensity_state, require_absolute_input_for_buffer_subtraction, + require_absolute_input_for_fluorescence_subtraction, require_relative_input_for_absolute_scaling, ) from .core.buffer_subtraction import BufferSubtractionResult, subtract_buffer +from .core.fluorescence_subtraction import ( + FluorescenceMethod, + FluorescenceSubtractionResult, + parse_fluorescence_method, + subtract_fluorescence, +) from .core.preflight import evaluate_preflight_gate, PreflightGateSummary from .core.execution_policy import ( RunPolicy, @@ -113,10 +120,15 @@ "IntensityStateAssessment", "assess_intensity_state", "require_absolute_input_for_buffer_subtraction", + "require_absolute_input_for_fluorescence_subtraction", "require_relative_input_for_absolute_scaling", # buffer subtraction "BufferSubtractionResult", "subtract_buffer", + "FluorescenceMethod", + "FluorescenceSubtractionResult", + "parse_fluorescence_method", + "subtract_fluorescence", # I/O "parse_header_values", "parse_header_values_with_meta", diff --git a/src/saxsabs/cli.py b/src/saxsabs/cli.py index 5cdb051..11b1e4a 100644 --- a/src/saxsabs/cli.py +++ b/src/saxsabs/cli.py @@ -1,6 +1,6 @@ """Command-line interface for headless SAXS calibration operations. -Provides seven subcommands: five small utilities plus the safe BL19B2 workflow +Provides eight subcommands: six small utilities plus the safe BL19B2 workflow and its explicit v1 migration entry. """ @@ -16,6 +16,7 @@ from . import __version__ from .core.buffer_subtraction import subtract_buffer +from .core.fluorescence_subtraction import subtract_fluorescence from .core.calibration import estimate_k_factor_robust from .core.intensity_state import require_relative_input_for_absolute_scaling from .core.normalization import compute_norm_factor @@ -384,6 +385,30 @@ def build_parser() -> argparse.ArgumentParser: p_sub.add_argument("--alpha", type=float, default=1.0) p_sub.add_argument("--alpha-uncertainty", type=float, default=None) + p_fluo = sub.add_parser( + "subtract-fluorescence", + help="Subtract additive fluorescence from an absolute sample profile", + ) + p_fluo.add_argument("--sample", required=True, type=Path) + p_fluo.add_argument( + "--method", + required=True, + choices=["constant", "high_q_mean", "high_q_median", "measured"], + help="How F(q) is obtained", + ) + p_fluo.add_argument("--f0", type=float, default=None, help="Constant F0 in cm^-1") + p_fluo.add_argument("--f0-uncertainty", type=float, default=None) + p_fluo.add_argument("--beta", type=float, default=1.0) + p_fluo.add_argument("--beta-uncertainty", type=float, default=None) + p_fluo.add_argument("--qmin", type=float, default=None, help="High-q window minimum") + p_fluo.add_argument("--qmax", type=float, default=None, help="High-q window maximum") + p_fluo.add_argument( + "--fluorescence", + type=Path, + default=None, + help="Measured additive fluorescence profile", + ) + p_bl = sub.add_parser( "bl19b2-abs2d", help="Process BL19B2 data with explicit monitor and thickness semantics", @@ -544,6 +569,58 @@ def main() -> None: ) return + if args.command == "subtract-fluorescence": + try: + sample = read_external_1d_profile(args.sample) + high_q_window = None + if args.qmin is not None or args.qmax is not None: + if args.qmin is None or args.qmax is None: + raise ValueError("--qmin and --qmax must be provided together") + high_q_window = (args.qmin, args.qmax) + q_fluo = i_fluo = err_fluo = None + fluo_profile = None + if args.fluorescence is not None: + fluo_profile = read_external_1d_profile(args.fluorescence) + q_fluo = fluo_profile["x"] + i_fluo = profile_intensity(fluo_profile) + err_fluo = profile_uncertainty(fluo_profile) + result = subtract_fluorescence( + sample["x"], + profile_intensity(sample), + profile_uncertainty(sample), + sample_profile=sample, + method=args.method, + f0=args.f0, + f0_uncertainty=args.f0_uncertainty, + beta=args.beta, + beta_uncertainty=args.beta_uncertainty, + high_q_window=high_q_window, + q_fluorescence=q_fluo, + i_fluorescence=i_fluo, + err_fluorescence=err_fluo, + fluorescence_profile=fluo_profile, + ) + except ValueError as exc: + _die(f"subtract-fluorescence failed: {exc}") + print( + json.dumps( + { + "points": int(result.q.size), + "method": result.method, + "beta": result.beta, + "beta_uncertainty": result.beta_uncertainty, + "f0": result.f0, + "f0_uncertainty": result.f0_uncertainty, + "high_q_residual_mean": result.high_q_residual_mean, + "high_q_check_passed": result.high_q_check_passed, + "high_q_points": result.high_q_points, + "negative_fraction": result.negative_fraction, + }, + ensure_ascii=False, + ) + ) + return + if args.command in {"bl19b2-abs2d", "bl19b2-abs2d-v1-legacy"}: from .workflows.bl19b2_abs2d import BL19B2Abs2DConfig, run_bl19b2_abs2d diff --git a/src/saxsabs/core/__init__.py b/src/saxsabs/core/__init__.py index c7ce875..fa80dec 100644 --- a/src/saxsabs/core/__init__.py +++ b/src/saxsabs/core/__init__.py @@ -38,10 +38,18 @@ assess_intensity_state, parse_correction_ledger, require_absolute_input_for_buffer_subtraction, + require_absolute_input_for_fluorescence_subtraction, require_relative_input_for_absolute_scaling, serialize_correction_ledger, ) from .buffer_subtraction import BufferSubtractionResult, subtract_buffer +from .fluorescence_subtraction import ( + FluorescenceMethod, + FluorescenceSubtractionResult, + parse_fluorescence_method, + subtract_fluorescence, + combine_sequential_standard_uncertainties, +) from .execution_policy import RunPolicy, parse_run_policy, should_skip_all_existing from .preflight import PreflightGateSummary, evaluate_preflight_gate from .reference_matching import ( @@ -108,10 +116,16 @@ "assess_intensity_state", "parse_correction_ledger", "require_absolute_input_for_buffer_subtraction", + "require_absolute_input_for_fluorescence_subtraction", "require_relative_input_for_absolute_scaling", "serialize_correction_ledger", "BufferSubtractionResult", "subtract_buffer", + "FluorescenceMethod", + "FluorescenceSubtractionResult", + "parse_fluorescence_method", + "subtract_fluorescence", + "combine_sequential_standard_uncertainties", "RunPolicy", "parse_run_policy", "should_skip_all_existing", diff --git a/src/saxsabs/core/fluorescence_subtraction.py b/src/saxsabs/core/fluorescence_subtraction.py new file mode 100644 index 0000000..a090105 --- /dev/null +++ b/src/saxsabs/core/fluorescence_subtraction.py @@ -0,0 +1,509 @@ +"""Additive fluorescence subtraction for absolute 1-D SAXS profiles. + +Implements an opt-in correction on the absolute cm^-1 scale: + + I_corr(q) = I_abs(q) − β × F(q) + σ_corr²(q) = σ_abs²(q) + β² × σ_F²(q) + F(q)² × σ_β² + +``F(q)`` is a non-negative additive term (sample X-ray fluorescence and any +q-independent inelastic background that remains after empty-cell subtraction). +It is applied after K, thickness, and optional buffer subtraction. Unknown +uncertainties stay NaN. + +This kernel is not detector-dark, NIST-blank, or solvent subtraction, and it +must not be applied in detector-count space without a solid-angle correction. +""" + +from __future__ import annotations + +from collections.abc import Mapping +from dataclasses import dataclass +from enum import Enum + +import numpy as np + +from saxsabs.core.intensity_state import require_absolute_input_for_fluorescence_subtraction + +DEFAULT_RESIDUAL_WINDOW = (0.15, 0.25) + + +class FluorescenceMethod(str, Enum): + """How the additive fluorescence term ``F(q)`` is obtained.""" + + CONSTANT = "constant" + HIGH_Q_MEAN = "high_q_mean" + HIGH_Q_MEDIAN = "high_q_median" + MEASURED_PROFILE = "measured_profile" + + +@dataclass(frozen=True) +class FluorescenceSubtractionResult: + """Container for fluorescence-subtracted SAXS data.""" + + q: np.ndarray + i_subtracted: np.ndarray + err_subtracted: np.ndarray + method: str + beta: float + f0: float + f_profile: np.ndarray + high_q_residual_mean: float = 0.0 + high_q_check_passed: bool = True + high_q_window: tuple[float, float] | None = None + high_q_points: int = 0 + negative_fraction: float = 0.0 + beta_uncertainty: float | None = None + f0_uncertainty: float | None = None + err_statistical: np.ndarray | None = None + + +def parse_fluorescence_method(method: object) -> FluorescenceMethod: + """Parse a method token; ``measured`` is accepted as measured_profile.""" + + token = str(method or "").strip().lower().replace("-", "_") + aliases = { + "constant": FluorescenceMethod.CONSTANT, + "high_q_mean": FluorescenceMethod.HIGH_Q_MEAN, + "high_qmean": FluorescenceMethod.HIGH_Q_MEAN, + "high_q_median": FluorescenceMethod.HIGH_Q_MEDIAN, + "high_qmedian": FluorescenceMethod.HIGH_Q_MEDIAN, + "measured": FluorescenceMethod.MEASURED_PROFILE, + "measured_profile": FluorescenceMethod.MEASURED_PROFILE, + } + parsed = aliases.get(token) + if parsed is None: + raise ValueError( + "fluorescence method must be constant, high_q_mean, high_q_median, " + f"or measured_profile; got {method!r}" + ) + return parsed + + +def _as_1d_float_array( + name: str, values: np.ndarray | None, *, require_finite: bool = True +) -> np.ndarray: + if values is None: + raise ValueError(f"{name} is required") + arr = np.asarray(values, dtype=np.float64) + if arr.ndim != 1: + raise ValueError(f"{name} must be a 1-D array") + if require_finite and not np.all(np.isfinite(arr)): + raise ValueError(f"{name} contains non-finite values") + return arr + + +def _optional_nonnegative_uncertainty(name: str, value: float | None) -> float | None: + if value is None: + return None + out = float(value) + if not np.isfinite(out) or out < 0: + raise ValueError(f"{name} must be finite and >= 0") + return out + + +def _validate_beta(beta: float) -> float: + value = float(beta) + if not np.isfinite(value) or value <= 0: + raise ValueError("Fluorescence scale factor beta must be finite and > 0") + return value + + +def _validate_f0(f0: float) -> float: + value = float(f0) + if not np.isfinite(value) or value < 0: + raise ValueError("f0 must be finite and >= 0") + return value + + +def _validate_window( + window: tuple[float, float] | None, *, name: str +) -> tuple[float, float] | None: + if window is None: + return None + try: + q_lo, q_hi = (float(value) for value in window) + except (TypeError, ValueError) as exc: + raise ValueError(f"{name} must contain two finite increasing values") from exc + if not np.isfinite(q_lo) or not np.isfinite(q_hi) or q_lo >= q_hi: + raise ValueError(f"{name} must contain two finite increasing values") + return (q_lo, q_hi) + + +def _prepare_source_grid( + q_source: np.ndarray, + y_source: np.ndarray, + *, + label: str, +) -> tuple[np.ndarray, np.ndarray]: + order = np.argsort(q_source) + q_sorted = q_source[order] + y_sorted = y_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") + if uq.size != q_sorted.size: + y_sum = np.zeros_like(uq, dtype=np.float64) + counts = np.zeros_like(uq, dtype=np.float64) + for idx, group in enumerate(inv): + y_sum[group] += y_sorted[idx] + counts[group] += 1.0 + y_sorted = y_sum / np.clip(counts, 1.0, None) + q_sorted = uq + return q_sorted, y_sorted + + +def _interpolate_on_grid( + q_target: np.ndarray, + q_source: np.ndarray, + y_source: np.ndarray, + *, + label: str, +) -> np.ndarray: + q_src, y_src = _prepare_source_grid(q_source, y_source, label=label) + tol = max( + 1e-12, + 1e-9 * max(abs(q_src[0]), abs(q_src[-1]), abs(q_target).max(initial=0.0)), + ) + if np.min(q_target) < q_src[0] - tol or np.max(q_target) > q_src[-1] + tol: + raise ValueError( + f"sample q grid extends outside {label} q range " + f"({q_src[0]:.6g} to {q_src[-1]:.6g})" + ) + return np.interp(q_target, q_src, y_src) + + +def _prepare_variance_grid( + q_source: np.ndarray, + sigma_source: np.ndarray, + *, + label: str, +) -> tuple[np.ndarray, np.ndarray]: + order = np.argsort(q_source) + q_sorted = q_source[order] + variance_sorted = np.square(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") + if uq.size == q_sorted.size: + return q_sorted, variance_sorted + + variance_of_mean = np.full(uq.shape, np.nan, dtype=np.float64) + 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 + ) + return uq, variance_of_mean + + +def _interpolate_variance_on_grid( + q_target: np.ndarray, + q_source: np.ndarray, + sigma_source: np.ndarray, + *, + label: str, +) -> np.ndarray: + q_src, variance_src = _prepare_variance_grid(q_source, sigma_source, label=label) + tol = max( + 1e-12, + 1e-9 * max(abs(q_src[0]), abs(q_src[-1]), abs(q_target).max(initial=0.0)), + ) + if np.min(q_target) < q_src[0] - tol or np.max(q_target) > q_src[-1] + tol: + raise ValueError( + f"sample q grid extends outside {label} q range " + f"({q_src[0]:.6g} to {q_src[-1]:.6g})" + ) + + upper = np.searchsorted(q_src, q_target, side="right") + upper = np.clip(upper, 1, q_src.size - 1) + lower = upper - 1 + span = q_src[upper] - q_src[lower] + weight_upper = (q_target - q_src[lower]) / span + weight_upper = np.clip(weight_upper, 0.0, 1.0) + weight_lower = 1.0 - weight_upper + + out = np.full(q_target.shape, np.nan, dtype=np.float64) + exact_lower = np.isclose(weight_upper, 0.0, rtol=0.0, atol=1e-14) + exact_upper = np.isclose(weight_upper, 1.0, rtol=0.0, atol=1e-14) + between = ~(exact_lower | exact_upper) + out[exact_lower] = variance_src[lower[exact_lower]] + out[exact_upper] = variance_src[upper[exact_upper]] + known = between & np.isfinite(variance_src[lower]) & np.isfinite(variance_src[upper]) + out[known] = ( + np.square(weight_lower[known]) * variance_src[lower[known]] + + np.square(weight_upper[known]) * variance_src[upper[known]] + ) + return out + + +def _window_mask(q: np.ndarray, i: np.ndarray, window: tuple[float, float]) -> np.ndarray: + q_lo, q_hi = window + return (q >= q_lo) & (q <= q_hi) & np.isfinite(i) + + +def _estimate_high_q_constant( + q: np.ndarray, + i_abs: np.ndarray, + window: tuple[float, float], + *, + use_median: bool, +) -> tuple[float, int]: + mask = _window_mask(q, i_abs, window) + n_points = int(mask.sum()) + if n_points < 3: + raise ValueError("high_q window must contain at least 3 finite intensity points") + values = i_abs[mask] + estimate = float(np.median(values) if use_median else np.mean(values)) + if not np.isfinite(estimate) or estimate < 0: + raise ValueError("estimated fluorescence constant must be finite and >= 0") + return estimate, n_points + + +def _residual_diagnostics( + q: np.ndarray, + i_corr: np.ndarray, + window: tuple[float, float], +) -> tuple[float, bool, int]: + mask = _window_mask(q, i_corr, window) + n_points = int(mask.sum()) + if n_points >= 3: + residual_mean = float(np.mean(i_corr[mask])) + residual_std = float(np.std(i_corr[mask])) + check_ok = abs(residual_mean) < 3.0 * max(residual_std, 1e-30) + return residual_mean, check_ok, n_points + return 0.0, True, n_points + + +def subtract_fluorescence( + q: np.ndarray, + i_abs: np.ndarray, + err_abs: np.ndarray | None, + *, + sample_profile: Mapping[str, object] | None = None, + method: str, + f0: float | None = None, + f0_uncertainty: float | None = None, + beta: float = 1.0, + beta_uncertainty: float | None = None, + high_q_window: tuple[float, float] | None = None, + q_fluorescence: np.ndarray | None = None, + i_fluorescence: np.ndarray | None = None, + err_fluorescence: np.ndarray | None = None, + fluorescence_profile: Mapping[str, object] | None = None, + residual_window: tuple[float, float] | None = None, +) -> FluorescenceSubtractionResult: + """Subtract an additive fluorescence term from an absolute SAXS curve. + + Parameters + ---------- + q, i_abs, err_abs + Sample profile on the absolute cm^-1 scale. Missing errors remain + unknown (NaN). + sample_profile + Intensity-state provenance. Required; unlabeled arrays are refused. + method + ``constant``, ``high_q_mean``, ``high_q_median``, or ``measured_profile``. + f0, f0_uncertainty + User constant and its standard uncertainty. Required for ``constant``. + Forbidden for ``high_q_*`` and ``measured_profile``. ``None`` + uncertainty keeps the combined result unknown. + beta, beta_uncertainty + Scale applied to ``F(q)``. ``None`` uncertainty keeps combined + uncertainty unknown; pass ``0.0`` only when β is treated as exact. + high_q_window + Required for ``high_q_*`` methods. + q_fluorescence, i_fluorescence, err_fluorescence, fluorescence_profile + Measured additive curve for ``measured_profile``. + residual_window + High-q diagnostic window. Defaults to ``high_q_window`` when present, + otherwise ``(0.15, 0.25)``. + """ + if sample_profile is None: + raise ValueError( + "subtract_fluorescence requires sample_profile with explicit " + "absolute_cm^-1 intensity_state and a cm^-1 intensity_unit" + ) + require_absolute_input_for_fluorescence_subtraction( + sample_profile, profile_name="sample" + ) + parsed_method = parse_fluorescence_method(method) + beta_value = _validate_beta(beta) + f0_uncertainty = _optional_nonnegative_uncertainty("f0_uncertainty", f0_uncertainty) + beta_uncertainty = _optional_nonnegative_uncertainty( + "beta_uncertainty", beta_uncertainty + ) + high_q_window = _validate_window(high_q_window, name="high_q_window") + residual_window = _validate_window(residual_window, name="residual_window") + + q_s = _as_1d_float_array("q", q) + i_s = _as_1d_float_array("i_abs", i_abs) + if q_s.shape != i_s.shape: + raise ValueError("q and i_abs shape mismatch") + e_s = ( + _as_1d_float_array("err_abs", err_abs, require_finite=False) + if err_abs is not None + else np.full_like(i_s, np.nan) + ) + if e_s.shape != i_s.shape: + raise ValueError("err_abs shape mismatch") + if np.any(np.isinf(e_s)): + raise ValueError("err_abs contains infinite values") + if np.any(np.isfinite(e_s) & (e_s < 0)): + raise ValueError("err_abs contains negative values") + + estimate_points = 0 + f_variance: np.ndarray + if parsed_method is FluorescenceMethod.CONSTANT: + if f0 is None: + raise ValueError("constant fluorescence method requires f0") + if high_q_window is not None: + raise ValueError("constant fluorescence method does not accept high_q_window") + if fluorescence_profile is not None or i_fluorescence is not None: + raise ValueError("constant fluorescence method does not accept a measured curve") + f0_value = _validate_f0(f0) + f_profile = np.full_like(i_s, f0_value) + if f0_uncertainty is None: + f_variance = np.full_like(i_s, np.nan) + else: + f_variance = np.full_like(i_s, f0_uncertainty**2) + elif parsed_method in { + FluorescenceMethod.HIGH_Q_MEAN, + FluorescenceMethod.HIGH_Q_MEDIAN, + }: + if f0 is not None: + raise ValueError("high_q fluorescence methods estimate f0 and refuse a user f0") + if high_q_window is None: + raise ValueError("high_q fluorescence methods require high_q_window") + if fluorescence_profile is not None or i_fluorescence is not None: + raise ValueError("high_q fluorescence methods do not accept a measured curve") + f0_value, estimate_points = _estimate_high_q_constant( + q_s, + i_s, + high_q_window, + use_median=parsed_method is FluorescenceMethod.HIGH_Q_MEDIAN, + ) + f_profile = np.full_like(i_s, f0_value) + if f0_uncertainty is None: + f_variance = np.full_like(i_s, np.nan) + else: + f_variance = np.full_like(i_s, f0_uncertainty**2) + else: + if f0 is not None: + raise ValueError("measured_profile refuses a scalar f0") + if high_q_window is not None: + raise ValueError("measured_profile does not accept high_q_window") + if f0_uncertainty is not None: + raise ValueError("measured_profile uses curve uncertainties, not f0_uncertainty") + if fluorescence_profile is None: + raise ValueError( + "measured_profile requires fluorescence_profile with explicit " + "absolute_cm^-1 intensity_state and a cm^-1 intensity_unit" + ) + require_absolute_input_for_fluorescence_subtraction( + fluorescence_profile, profile_name="fluorescence" + ) + q_f = _as_1d_float_array("q_fluorescence", q_fluorescence) + i_f = _as_1d_float_array("i_fluorescence", i_fluorescence) + if q_f.shape != i_f.shape: + raise ValueError("q_fluorescence and i_fluorescence shape mismatch") + if np.any(np.isfinite(i_f) & (i_f < 0)): + raise ValueError("i_fluorescence contains negative values") + e_f = ( + _as_1d_float_array("err_fluorescence", err_fluorescence, require_finite=False) + if err_fluorescence is not None + else np.full_like(i_f, np.nan) + ) + if e_f.shape != i_f.shape: + raise ValueError("err_fluorescence shape mismatch") + if np.any(np.isinf(e_f)): + raise ValueError("err_fluorescence contains infinite values") + if np.any(np.isfinite(e_f) & (e_f < 0)): + raise ValueError("err_fluorescence contains negative values") + if q_s.shape != q_f.shape or not np.allclose(q_s, q_f, rtol=0.0, atol=1e-8): + f_profile = _interpolate_on_grid( + q_s, q_f, i_f, label="fluorescence" + ) + f_variance = _interpolate_variance_on_grid( + q_s, q_f, e_f, label="fluorescence uncertainty" + ) + else: + f_profile = i_f + f_variance = np.square(e_f) + finite_f = f_profile[np.isfinite(f_profile)] + f0_value = float(np.mean(finite_f)) if finite_f.size else float("nan") + if not np.isfinite(f0_value) or f0_value < 0: + raise ValueError("measured fluorescence intensity must be finite and >= 0") + f0_uncertainty = None + + subtracted_term = beta_value * f_profile + i_corr = i_s - subtracted_term + + variance_statistical = np.square(e_s) + (beta_value**2) * f_variance + if parsed_method is not FluorescenceMethod.MEASURED_PROFILE and f0_uncertainty is None: + # Unknown u(F0) is a combined-budget gap, not a missing sample error. + variance_statistical = np.square(e_s) + err_statistical = np.sqrt(variance_statistical) + if beta_uncertainty is None or ( + parsed_method is not FluorescenceMethod.MEASURED_PROFILE and f0_uncertainty is None + ): + variance_combined = np.full_like(i_s, np.nan) + else: + variance_combined = variance_statistical + np.square(f_profile * beta_uncertainty) + err_combined = np.sqrt(variance_combined) + + diag_window = residual_window + if diag_window is None: + diag_window = high_q_window if high_q_window is not None else DEFAULT_RESIDUAL_WINDOW + residual_mean, check_ok, residual_points = _residual_diagnostics(q_s, i_corr, diag_window) + report_points = estimate_points if estimate_points else residual_points + finite = np.isfinite(i_corr) + if finite.any(): + negative_fraction = float(np.mean(i_corr[finite] < 0.0)) + else: + negative_fraction = float("nan") + + return FluorescenceSubtractionResult( + q=q_s, + i_subtracted=i_corr, + err_subtracted=err_combined, + method=parsed_method.value, + beta=beta_value, + f0=f0_value, + f_profile=f_profile, + high_q_residual_mean=residual_mean, + high_q_check_passed=check_ok, + high_q_window=high_q_window, + high_q_points=report_points, + negative_fraction=negative_fraction, + beta_uncertainty=beta_uncertainty, + f0_uncertainty=f0_uncertainty, + err_statistical=err_statistical, + ) + + +def combine_sequential_standard_uncertainties( + previous_statistical: np.ndarray, + previous_combined: np.ndarray | None, + next_statistical: np.ndarray, + next_combined: np.ndarray, +) -> tuple[np.ndarray, np.ndarray]: + """Compose two additive steps without mixing extras into the statistical term. + + ``next_statistical`` must have been computed using ``previous_statistical`` + as the incoming sample uncertainty. Any unknown combined term stays NaN. + """ + stat = np.asarray(next_statistical, dtype=np.float64) + nxt = np.asarray(next_combined, dtype=np.float64) + if previous_combined is None: + return stat, nxt + prev_stat = np.asarray(previous_statistical, dtype=np.float64) + prev_comb = np.asarray(previous_combined, dtype=np.float64) + extra_prev = np.square(prev_comb) - np.square(prev_stat) + extra_next = np.square(nxt) - np.square(stat) + extra_prev = np.clip(extra_prev, 0.0, None) + extra_next = np.clip(extra_next, 0.0, None) + combined = np.sqrt(np.square(stat) + extra_prev + extra_next) + unknown = ~np.isfinite(prev_comb) | ~np.isfinite(nxt) + combined = np.where(unknown, np.nan, combined) + return stat, combined diff --git a/src/saxsabs/core/intensity_state.py b/src/saxsabs/core/intensity_state.py index 402f434..b068fc4 100644 --- a/src/saxsabs/core/intensity_state.py +++ b/src/saxsabs/core/intensity_state.py @@ -46,6 +46,7 @@ def is_cm_inv_intensity_unit(value: object) -> bool: "polarization", "flat_field", "buffer", + "fluorescence", } ) @@ -95,6 +96,10 @@ def _canonical_correction(value: object) -> str: "flat": "flat_field", "flatfield": "flat_field", "buffer": "buffer", + "fluorescence": "fluorescence", + "fluo": "fluorescence", + "xrf": "fluorescence", + "fluorescencebackground": "fluorescence", } canonical = aliases.get(token) if canonical is None: @@ -329,3 +334,35 @@ def require_absolute_input_for_buffer_subtraction( if "buffer" in corrections: raise ValueError(f"{profile_name}: buffer profile is already buffer-subtracted") return assessment + + +def require_absolute_input_for_fluorescence_subtraction( + profile: Mapping[str, object], + *, + profile_name: str = "sample", +) -> IntensityStateAssessment: + """Require a traceable absolute profile before subtracting fluorescence.""" + + assessment = assess_intensity_state(profile) + if assessment.state is not IntensityState.ABSOLUTE_CM_INV: + raise ValueError( + f"{profile_name}: fluorescence subtraction requires explicit " + "absolute intensity in cm^-1" + ) + provenance = profile.get("operator_provenance") + provenance = provenance if isinstance(provenance, Mapping) else {} + raw_unit = profile.get("intensity_unit", provenance.get("intensity_unit", "")) + if not is_cm_inv_intensity_unit(raw_unit): + raise ValueError(f"{profile_name}: intensity_unit must be 1/cm") + corrections = set(assessment.corrections_applied) + missing = set(ABSOLUTE_CORRECTIONS) - corrections + if missing: + raise ValueError( + f"{profile_name}: absolute ledger is missing: " + + ", ".join(sorted(missing)) + ) + if "fluorescence" in assessment.protected_corrections: + raise ValueError( + f"{profile_name}: profile is already fluorescence-subtracted" + ) + return assessment diff --git a/src/saxsabs/io/parsers.py b/src/saxsabs/io/parsers.py index 49cabbf..a434f65 100644 --- a/src/saxsabs/io/parsers.py +++ b/src/saxsabs/io/parsers.py @@ -59,6 +59,14 @@ "buffersourcesha256": "buffer_source_sha256", "bufferalpha": "buffer_alpha", "bufferalphauncertainty": "buffer_alpha_uncertainty", + "fluorescencemethod": "fluorescence_method", + "fluorescencef0": "fluorescence_f0", + "fluorescencef0uncertainty": "fluorescence_f0_uncertainty", + "fluorescencebeta": "fluorescence_beta", + "fluorescencebetauncertainty": "fluorescence_beta_uncertainty", + "fluorescencehighqwindow": "fluorescence_high_q_window", + "fluorescencesourcename": "fluorescence_source_name", + "fluorescencesourcesha256": "fluorescence_source_sha256", "uncertaintymodel": "uncertainty_model", "uncertaintytype": "uncertainty_type", } diff --git a/src/saxsabs/io/writers.py b/src/saxsabs/io/writers.py index 3e104e2..2148ff2 100644 --- a/src/saxsabs/io/writers.py +++ b/src/saxsabs/io/writers.py @@ -52,6 +52,14 @@ "buffer_source_sha256", "buffer_alpha", "buffer_alpha_uncertainty", + "fluorescence_method", + "fluorescence_f0", + "fluorescence_f0_uncertainty", + "fluorescence_beta", + "fluorescence_beta_uncertainty", + "fluorescence_high_q_window", + "fluorescence_source_name", + "fluorescence_source_sha256", "uncertainty_model", "uncertainty_type", ) diff --git a/src/saxsabs/workflows/bl19b2_integrate1d.py b/src/saxsabs/workflows/bl19b2_integrate1d.py index 371b7bf..afed10d 100644 --- a/src/saxsabs/workflows/bl19b2_integrate1d.py +++ b/src/saxsabs/workflows/bl19b2_integrate1d.py @@ -43,6 +43,13 @@ class Integrate1DConfig: correct_solid_angle: bool = True polarization_factor: float | None = None resume: bool = True + fluorescence_method: str | None = None + fluorescence_f0: float | None = None + fluorescence_f0_uncertainty: float | None = None + fluorescence_beta: float = 1.0 + fluorescence_beta_uncertainty: float | None = None + fluorescence_qmin: float | None = None + fluorescence_qmax: float | None = None def output_root(self) -> Path: return Path(self.package_root) / "integration" @@ -246,6 +253,65 @@ def _write_new_or_verify(path: Path, data: bytes) -> None: temporary.replace(path) +def _absolute_1d_sample_profile() -> dict[str, object]: + return { + "intensity_state": "absolute_cm^-1", + "intensity_unit": "1/cm", + "operator_provenance": { + "intensity_state": "absolute_cm^-1", + "corrections_applied": '["k","thickness"]', + }, + } + + +def _fluorescence_config_payload(config: Integrate1DConfig) -> dict[str, Any]: + return { + "fluorescence_method": config.fluorescence_method, + "fluorescence_f0": config.fluorescence_f0, + "fluorescence_f0_uncertainty": config.fluorescence_f0_uncertainty, + "fluorescence_beta": config.fluorescence_beta, + "fluorescence_beta_uncertainty": config.fluorescence_beta_uncertainty, + "fluorescence_qmin": config.fluorescence_qmin, + "fluorescence_qmax": config.fluorescence_qmax, + } + + +def _apply_optional_fluorescence( + q: np.ndarray, + intensity: np.ndarray, + config: Integrate1DConfig, +) -> tuple[np.ndarray, dict[str, Any] | None]: + if not config.fluorescence_method: + return intensity, None + from saxsabs.core.fluorescence_subtraction import subtract_fluorescence + + window = None + if config.fluorescence_qmin is not None or config.fluorescence_qmax is not None: + if config.fluorescence_qmin is None or config.fluorescence_qmax is None: + raise ValueError("fluorescence_qmin and fluorescence_qmax must be provided together") + window = (float(config.fluorescence_qmin), float(config.fluorescence_qmax)) + result = subtract_fluorescence( + q, + intensity, + None, + sample_profile=_absolute_1d_sample_profile(), + method=config.fluorescence_method, + f0=config.fluorescence_f0, + f0_uncertainty=config.fluorescence_f0_uncertainty, + beta=config.fluorescence_beta, + beta_uncertainty=config.fluorescence_beta_uncertainty, + high_q_window=window, + ) + return result.i_subtracted, { + "method": result.method, + "f0": result.f0, + "beta": result.beta, + "negative_fraction": result.negative_fraction, + "high_q_residual_mean": result.high_q_residual_mean, + "high_q_check_passed": result.high_q_check_passed, + } + + def _profile_bytes(q: np.ndarray, intensity: np.ndarray) -> bytes: text = io.StringIO(newline="") writer = csv.writer(text, lineterminator="\n") @@ -462,6 +528,7 @@ def run_bl19b2_integrate1d(config: Integrate1DConfig) -> dict[str, Any]: "normalization_factor": 1.0, "do_not_repeat": sorted(DO_NOT_REPEAT), "pyFAI_version": _version("pyFAI"), + **_fluorescence_config_payload(config), } run_signature = _canonical_hash(signature_payload) out = config.output_root() @@ -504,7 +571,8 @@ def run_bl19b2_integrate1d(config: Integrate1DConfig) -> dict[str, Any]: "# BL19B2 absolute 1D integration\n\n" "The EDF inputs were already corrected for dark, background, monitor, transmission, " "thickness, and K. This step applies only the package mask and one solid-angle correction " - "during pyFAI CSR integration. No polarization correction is applied.\n" + "during pyFAI CSR integration. No polarization correction is applied. " + "Optional fluorescence subtraction, if configured, is applied to the absolute 1D curve.\n" ).encode("utf-8") _write_new_or_verify(out / "README.md", readme) @@ -541,6 +609,7 @@ def run_bl19b2_integrate1d(config: Integrate1DConfig) -> dict[str, Any]: else: image = _load_validate_edf(item, metadata, mask) q, intensity = _integrate(ai, image, mask, config) + intensity, fluorescence_diag = _apply_optional_fluorescence(q, intensity, config) profile_data = _profile_bytes(q, intensity) _write_new_or_verify(profile, profile_data) side_doc = { @@ -561,8 +630,12 @@ def run_bl19b2_integrate1d(config: Integrate1DConfig) -> dict[str, Any]: "dark": None, "flat": None, "normalization_factor": 1.0, - "do_not_repeat": sorted(DO_NOT_REPEAT), + "do_not_repeat": sorted( + set(DO_NOT_REPEAT) + | ({"fluorescence"} if fluorescence_diag else set()) + ), }, + "fluorescence": fluorescence_diag, } _write_new_or_verify(sidecar, _json_bytes(side_doc)) processed += 1 diff --git a/tests/test_bl19b2_integrate1d.py b/tests/test_bl19b2_integrate1d.py index d405961..9f91639 100644 --- a/tests/test_bl19b2_integrate1d.py +++ b/tests/test_bl19b2_integrate1d.py @@ -142,6 +142,58 @@ def integrate1d(self, image, npt, **kwargs): assert len(calls) == 1 +def test_integration_optional_constant_fluorescence_changes_1d_and_signature( + tmp_path: Path, monkeypatch: pytest.MonkeyPatch +): + package, _manifest = _build_package(tmp_path) + + class FakeIntegrator: + def integrate1d(self, _image, npt, **_kwargs): + return SimpleNamespace( + radial=np.linspace(0.001, 1.0, npt), + intensity=np.full(npt, 4.0), + ) + + monkeypatch.setattr(integration, "_load_integrator", lambda _path: FakeIntegrator()) + monkeypatch.setattr( + integration, + "_load_validate_edf", + lambda _item, _metadata, _mask: np.arange(4, dtype=np.float32).reshape(2, 2), + ) + monkeypatch.setattr(integration, "_version", lambda _name: "test-pyfai") + + result = integration.run_bl19b2_integrate1d( + integration.Integrate1DConfig( + package, + fluorescence_method="constant", + fluorescence_f0=1.5, + fluorescence_f0_uncertainty=0.0, + fluorescence_beta_uncertainty=0.0, + ) + ) + assert result["processed"] == 1 + profile = package / "integration" / "profiles" / "problem" / "sample_00001_abs1d_cm-1.csv" + rows = [line.split(",") for line in profile.read_text(encoding="utf-8").splitlines()[1:]] + intensities = [float(row[1]) for row in rows] + assert intensities[0] == pytest.approx(2.5) + sidecar = json.loads( + ( + package + / "integration" + / "metadata" + / "problem" + / "sample_00001_abs1d.json" + ).read_text(encoding="utf-8") + ) + assert sidecar["fluorescence"]["method"] == "constant" + assert sidecar["fluorescence"]["f0"] == pytest.approx(1.5) + signature = json.loads( + (package / "integration" / "config" / "run_signature.json").read_text(encoding="utf-8") + ) + assert signature["payload"]["fluorescence_method"] == "constant" + assert signature["payload"]["fluorescence_f0"] == pytest.approx(1.5) + + def test_resume_uses_canonical_selection_and_scientific_metadata_hash( tmp_path: Path, monkeypatch: pytest.MonkeyPatch ): diff --git a/tests/test_cli.py b/tests/test_cli.py index b675f1f..7ced143 100644 --- a/tests/test_cli.py +++ b/tests/test_cli.py @@ -396,6 +396,68 @@ def test_cli_subtract_buffer_refuses_unlabeled_profiles( assert "absolute" in capsys.readouterr().err.lower() +def test_cli_subtract_fluorescence( + tmp_path: Path, + capsys: pytest.CaptureFixture[str], + monkeypatch: pytest.MonkeyPatch, +): + sample = tmp_path / "sample.csv" + _write_absolute_profile(sample, "0.01,12.0\n0.02,11.0\n0.20,6.0\n") + monkeypatch.setattr( + sys, + "argv", + [ + "saxsabs", + "subtract-fluorescence", + "--sample", + str(sample), + "--method", + "constant", + "--f0", + "2.0", + "--f0-uncertainty", + "0.0", + "--beta-uncertainty", + "0.0", + ], + ) + main() + out = json.loads(capsys.readouterr().out) + assert out["points"] == 3 + assert out["method"] == "constant" + assert out["f0"] == pytest.approx(2.0) + assert out["negative_fraction"] == pytest.approx(0.0) + + +def test_cli_subtract_fluorescence_refuses_unlabeled_profiles( + tmp_path: Path, + capsys: pytest.CaptureFixture[str], + monkeypatch: pytest.MonkeyPatch, +): + sample = tmp_path / "sample.csv" + sample.write_text("q,i\n0.01,12.0\n0.02,11.0\n0.20,6.0\n", encoding="utf-8") + monkeypatch.setattr( + sys, + "argv", + [ + "saxsabs", + "subtract-fluorescence", + "--sample", + str(sample), + "--method", + "constant", + "--f0", + "2.0", + ], + ) + + with pytest.raises(SystemExit) as exc_info: + main() + + assert exc_info.value.code == 1 + assert "absolute" in capsys.readouterr().err.lower() + + def test_cli_version(capsys: pytest.CaptureFixture[str], monkeypatch: pytest.MonkeyPatch): from saxsabs import __version__ diff --git a/tests/test_fluorescence_subtraction.py b/tests/test_fluorescence_subtraction.py new file mode 100644 index 0000000..cd2044d --- /dev/null +++ b/tests/test_fluorescence_subtraction.py @@ -0,0 +1,417 @@ +"""Tests for absolute-scale fluorescence subtraction.""" + +import numpy as np +import pytest + +from saxsabs.core.fluorescence_subtraction import ( + FluorescenceSubtractionResult, + combine_sequential_standard_uncertainties, + parse_fluorescence_method, + subtract_fluorescence, +) + +ABS = { + "intensity_state": "absolute_cm^-1", + "intensity_unit": "1/cm", + "i_col": "I_abs_cm^-1", + "operator_provenance": { + "intensity_state": "absolute_cm^-1", + "corrections_applied": '["k","thickness"]', + }, +} + + +def _sub(*args, **kwargs): + kwargs.setdefault("sample_profile", ABS) + return subtract_fluorescence(*args, **kwargs) + + +def test_parse_fluorescence_method_aliases(): + assert parse_fluorescence_method("measured").value == "measured_profile" + assert parse_fluorescence_method("HIGH_Q_MEAN").value == "high_q_mean" + with pytest.raises(ValueError, match="fluorescence method"): + parse_fluorescence_method("compton") + + +def test_constant_subtraction_and_error(): + q = np.array([0.01, 0.02, 0.20]) + i = np.array([12.0, 11.0, 6.0]) + e = np.array([0.1, 0.1, 0.2]) + out = _sub( + q, + i, + e, + method="constant", + f0=2.0, + f0_uncertainty=0.0, + beta=1.0, + beta_uncertainty=0.0, + ) + assert isinstance(out, FluorescenceSubtractionResult) + np.testing.assert_allclose(out.i_subtracted, [10.0, 9.0, 4.0]) + np.testing.assert_allclose(out.err_subtracted, e) + np.testing.assert_allclose(out.err_statistical, e) + assert out.negative_fraction == 0.0 + assert out.method == "constant" + assert out.f0 == pytest.approx(2.0) + + +def test_beta_scales_constant_and_propagates_beta_uncertainty(): + q = np.array([0.01, 0.02, 0.03]) + out = _sub( + q, + np.full(3, 10.0), + np.full(3, 0.1), + method="constant", + f0=2.0, + f0_uncertainty=0.0, + beta=1.5, + beta_uncertainty=0.05, + ) + np.testing.assert_allclose(out.i_subtracted, 7.0) + expected = np.sqrt(0.1**2 + (2.0 * 0.05) ** 2) + np.testing.assert_allclose(out.err_subtracted, expected) + np.testing.assert_allclose(out.err_statistical, 0.1) + np.testing.assert_allclose(out.f_profile, 2.0) + + +def test_missing_f0_uncertainty_keeps_combined_unknown(): + q = np.array([0.01, 0.02, 0.20]) + out = _sub( + q, + np.ones(3) * 5.0, + np.ones(3) * 0.1, + method="constant", + f0=1.0, + ) + assert out.f0_uncertainty is None + np.testing.assert_allclose(out.err_statistical, 0.1) + assert np.all(np.isnan(out.err_subtracted)) + + +def test_f0_uncertainty_enters_statistical_and_combined(): + q = np.array([0.01, 0.02, 0.03]) + out = _sub( + q, + np.full(3, 8.0), + np.full(3, 0.2), + method="constant", + f0=1.0, + f0_uncertainty=0.3, + beta=2.0, + beta_uncertainty=0.0, + ) + expected = np.sqrt(0.2**2 + (2.0 * 0.3) ** 2) + np.testing.assert_allclose(out.err_statistical, expected) + np.testing.assert_allclose(out.err_subtracted, expected) + + +def test_high_q_mean_estimates_planted_constant(): + q = np.linspace(0.01, 0.30, 60) + i = np.exp(-q / 0.02) + 3.5 + out = _sub( + q, + i, + np.full_like(q, 0.01), + method="high_q_mean", + high_q_window=(0.25, 0.30), + f0_uncertainty=0.0, + beta_uncertainty=0.0, + ) + assert out.f0 == pytest.approx(3.5, rel=0.02) + np.testing.assert_allclose(out.i_subtracted, i - out.f0, rtol=1e-12) + assert out.high_q_points >= 3 + assert out.method == "high_q_mean" + + +def test_high_q_median_is_robust_to_one_spike(): + q = np.linspace(0.20, 0.30, 11) + i = np.full(q.shape, 4.0) + i[5] = 40.0 + mean_out = _sub( + q, i, np.full_like(q, 0.01), method="high_q_mean", + high_q_window=(0.20, 0.30), f0_uncertainty=0.0, beta_uncertainty=0.0, + ) + median_out = _sub( + q, i, np.full_like(q, 0.01), method="high_q_median", + high_q_window=(0.20, 0.30), f0_uncertainty=0.0, beta_uncertainty=0.0, + ) + assert mean_out.f0 > median_out.f0 + assert median_out.f0 == pytest.approx(4.0) + + +def test_measured_profile_interpolates_like_buffer(): + q_s = np.array([1.0]) + q_f = np.array([0.0, 2.0]) + out = _sub( + q_s, + np.array([10.0]), + np.array([0.4]), + method="measured_profile", + q_fluorescence=q_f, + i_fluorescence=np.array([2.0, 4.0]), + err_fluorescence=np.array([1.0, 3.0]), + fluorescence_profile=ABS, + beta=1.0, + beta_uncertainty=0.0, + ) + assert out.i_subtracted[0] == pytest.approx(7.0) + assert out.err_subtracted[0] == pytest.approx(np.sqrt(0.4**2 + 2.5)) + assert out.f0 == pytest.approx(3.0) + + +def test_measured_missing_curve_uncertainty_keeps_results_unknown(): + q = np.array([0.01, 0.02, 0.03]) + out = _sub( + q, + np.full(3, 10.0), + np.full(3, 0.1), + method="measured", + q_fluorescence=q, + i_fluorescence=np.full(3, 2.0), + err_fluorescence=None, + fluorescence_profile=ABS, + beta_uncertainty=0.0, + ) + assert np.all(np.isnan(out.err_statistical)) + assert np.all(np.isnan(out.err_subtracted)) + + +def test_refuses_unlabeled_and_negative_f0(): + q = np.array([0.01, 0.02, 0.03]) + with pytest.raises(ValueError, match="sample_profile"): + subtract_fluorescence(q, np.ones(3), np.ones(3), method="constant", f0=1.0) + with pytest.raises(ValueError, match="f0"): + _sub(q, np.ones(3) * 5, np.ones(3) * 0.1, method="constant", f0=-0.1) + + +def test_already_fluorescence_subtracted_is_refused(): + q = np.array([0.01, 0.02, 0.03]) + profile = { + **ABS, + "operator_provenance": { + "intensity_state": "absolute_cm^-1", + "corrections_applied": '["fluorescence","k","thickness"]', + }, + } + with pytest.raises(ValueError, match="already fluorescence-subtracted"): + subtract_fluorescence( + q, np.ones(3) * 5, np.ones(3) * 0.1, + sample_profile=profile, method="constant", f0=1.0, + ) + + +@pytest.mark.parametrize( + ("kwargs", "match"), + [ + ({"method": "nope", "f0": 1.0}, "fluorescence method"), + ({"method": "high_q_mean"}, "high_q_window"), + ( + { + "method": "high_q_mean", + "high_q_window": (0.29, 0.30), + }, + "at least 3", + ), + ({"method": "measured_profile"}, "fluorescence_profile"), + ({"method": "constant", "f0": 1.0, "beta": 0.0}, "beta"), + ({"method": "constant"}, "requires f0"), + ( + {"method": "high_q_mean", "high_q_window": (0.01, 0.03), "f0": 1.0}, + "refuse a user f0", + ), + ], +) +def test_invalid_method_arguments_raise(kwargs, match): + q = np.array([0.01, 0.02, 0.03]) + with pytest.raises(ValueError, match=match): + _sub(q, np.full(3, 5.0), np.full(3, 0.1), **kwargs) + + +def test_measured_profile_outside_q_range_raises(): + q_s = np.array([0.01, 0.20, 0.40]) + q_f = np.array([0.05, 0.10, 0.30]) + with pytest.raises(ValueError, match="outside fluorescence q range"): + _sub( + q_s, + np.ones(3) * 10.0, + np.ones(3) * 0.1, + method="measured", + q_fluorescence=q_f, + i_fluorescence=np.ones(3) * 2.0, + err_fluorescence=np.ones(3) * 0.05, + fluorescence_profile=ABS, + beta_uncertainty=0.0, + ) + + +def test_negative_err_abs_raises(): + q = np.array([0.01, 0.02, 0.03]) + with pytest.raises(ValueError, match="err_abs"): + _sub( + q, + np.ones(3) * 5.0, + np.array([0.1, -0.2, 0.1]), + method="constant", + f0=1.0, + ) + + +def test_infinite_err_abs_raises(): + q = np.array([0.01, 0.02, 0.03]) + err = np.array([0.1, np.inf, 0.1]) + with pytest.raises(ValueError, match="err_abs"): + _sub(q, np.ones(3) * 5.0, err, method="constant", f0=1.0) + + +def test_negative_fraction_is_reported_without_clipping(): + q = np.array([0.01, 0.02, 0.03]) + out = _sub( + q, + np.array([1.0, 0.4, 0.1]), + np.full(3, 0.01), + method="constant", + f0=0.5, + f0_uncertainty=0.0, + beta_uncertainty=0.0, + ) + np.testing.assert_allclose(out.i_subtracted, [0.5, -0.1, -0.4]) + assert out.negative_fraction == pytest.approx(2.0 / 3.0) + + +def test_buffer_then_fluorescence_ledger_is_allowed(): + profile = { + **ABS, + "operator_provenance": { + "intensity_state": "absolute_cm^-1", + "corrections_applied": '["buffer","k","thickness"]', + }, + } + q = np.array([0.01, 0.02, 0.03]) + out = subtract_fluorescence( + q, + np.full(3, 4.0), + np.full(3, 0.1), + sample_profile=profile, + method="constant", + f0=1.0, + f0_uncertainty=0.0, + beta_uncertainty=0.0, + ) + np.testing.assert_allclose(out.i_subtracted, 3.0) + + +def test_shipped_buffer_then_fluorescence_keeps_statistical_and_unknown_extras(): + """Drive the real Tab3 composition: buffer kernel then fluorescence kernel.""" + from saxsabs.core.buffer_subtraction import subtract_buffer + + q = np.array([0.01, 0.02, 0.03]) + i_sample = np.array([12.0, 11.0, 10.0]) + i_buffer = np.full(3, 2.0) + err_sample = np.full(3, 0.1) + err_buffer = np.full(3, 0.2) + alpha = 0.5 + u_alpha = 0.05 + f0 = 1.0 + u_f0 = 0.3 + beta = 1.0 + u_beta = 0.04 + + buffered = subtract_buffer( + q, + i_sample, + err_sample, + q, + i_buffer, + err_buffer, + alpha=alpha, + alpha_uncertainty=u_alpha, + sample_profile=ABS, + buffer_profile=ABS, + ) + after_buffer = { + **ABS, + "operator_provenance": { + "intensity_state": "absolute_cm^-1", + "corrections_applied": '["buffer","k","thickness"]', + }, + } + fluo = subtract_fluorescence( + buffered.q, + buffered.i_subtracted, + buffered.err_statistical, + sample_profile=after_buffer, + method="constant", + f0=f0, + f0_uncertainty=u_f0, + beta=beta, + beta_uncertainty=u_beta, + ) + stat, comb = combine_sequential_standard_uncertainties( + buffered.err_statistical, + buffered.err_subtracted, + fluo.err_statistical, + fluo.err_subtracted, + ) + + np.testing.assert_allclose(fluo.i_subtracted, i_sample - alpha * i_buffer - f0) + expected_stat = np.sqrt(err_sample**2 + (alpha * err_buffer) ** 2 + (beta * u_f0) ** 2) + np.testing.assert_allclose(stat, expected_stat) + expected_comb = np.sqrt( + expected_stat**2 + (i_buffer * u_alpha) ** 2 + (f0 * u_beta) ** 2 + ) + np.testing.assert_allclose(comb, expected_comb) + + unknown_alpha = subtract_buffer( + q, + i_sample, + err_sample, + q, + i_buffer, + err_buffer, + alpha=alpha, + sample_profile=ABS, + buffer_profile=ABS, + ) + fluo_after_unknown = subtract_fluorescence( + unknown_alpha.q, + unknown_alpha.i_subtracted, + unknown_alpha.err_statistical, + sample_profile=after_buffer, + method="constant", + f0=f0, + f0_uncertainty=u_f0, + beta=beta, + beta_uncertainty=u_beta, + ) + _stat2, comb2 = combine_sequential_standard_uncertainties( + unknown_alpha.err_statistical, + unknown_alpha.err_subtracted, + fluo_after_unknown.err_statistical, + fluo_after_unknown.err_subtracted, + ) + np.testing.assert_allclose(_stat2, expected_stat) + assert np.all(np.isnan(comb2)) + + +def test_combine_sequential_keeps_unknown_previous_combined(): + stat, comb = combine_sequential_standard_uncertainties( + np.array([0.2, 0.2]), + np.array([np.nan, np.nan]), + np.array([0.3, 0.3]), + np.array([0.4, 0.4]), + ) + np.testing.assert_allclose(stat, 0.3) + assert np.all(np.isnan(comb)) + + +def test_combine_sequential_adds_independent_extras(): + prev_stat = np.array([0.1]) + prev_comb = np.array([np.sqrt(0.1**2 + 0.3**2)]) + next_stat = np.array([np.sqrt(0.1**2 + 0.2**2)]) + next_comb = np.array([np.sqrt(0.1**2 + 0.2**2 + 0.4**2)]) + stat, comb = combine_sequential_standard_uncertainties( + prev_stat, prev_comb, next_stat, next_comb + ) + np.testing.assert_allclose(stat, next_stat) + np.testing.assert_allclose(comb, np.sqrt(0.1**2 + 0.2**2 + 0.3**2 + 0.4**2)) diff --git a/tests/test_intensity_state.py b/tests/test_intensity_state.py index d9df1ae..9ebcca9 100644 --- a/tests/test_intensity_state.py +++ b/tests/test_intensity_state.py @@ -7,6 +7,7 @@ assess_intensity_state, parse_correction_ledger, require_absolute_input_for_buffer_subtraction, + require_absolute_input_for_fluorescence_subtraction, require_relative_input_for_absolute_scaling, serialize_correction_ledger, ) @@ -212,3 +213,52 @@ def test_absolute_buffer_requires_unit_complete_ledger_and_no_prior_buffer(): } with pytest.raises(ValueError, match="absolute buffer ledger is missing"): require_absolute_input_for_buffer_subtraction(guard_only) + + +def test_fluorescence_aliases_are_canonical(): + assert parse_correction_ledger("fluo,xrf") == ("fluorescence",) + assert parse_correction_ledger(["fluorescence_background"]) == ("fluorescence",) + + +def test_absolute_fluorescence_requires_unit_complete_ledger_and_no_prior_fluorescence(): + valid = { + "intensity_state": "absolute_cm^-1", + "intensity_unit": "1/cm", + "operator_provenance": { + "intensity_state": "absolute_cm^-1", + "corrections_applied": '["buffer","k","thickness"]', + }, + } + assert require_absolute_input_for_fluorescence_subtraction(valid).is_absolute + + relative = { + "i_col": "I_rel", + "operator_provenance": {"intensity_state": "relative"}, + } + with pytest.raises(ValueError, match="absolute intensity"): + require_absolute_input_for_fluorescence_subtraction(relative) + + missing_unit = {**valid, "intensity_unit": ""} + with pytest.raises(ValueError, match="intensity_unit"): + require_absolute_input_for_fluorescence_subtraction(missing_unit) + + repeated = { + **valid, + "operator_provenance": { + "intensity_state": "absolute_cm^-1", + "corrections_applied": '["fluorescence","k","thickness"]', + }, + } + with pytest.raises(ValueError, match="already fluorescence-subtracted"): + require_absolute_input_for_fluorescence_subtraction(repeated) + + guard_only = { + "intensity_state": "absolute_cm^-1", + "intensity_unit": "1/cm", + "operator_provenance": { + "intensity_state": "absolute_cm^-1", + "do_not_repeat": '["k","thickness"]', + }, + } + with pytest.raises(ValueError, match="absolute ledger is missing"): + require_absolute_input_for_fluorescence_subtraction(guard_only) diff --git a/tests/test_public_exports.py b/tests/test_public_exports.py index ce174b2..9f2b62d 100644 --- a/tests/test_public_exports.py +++ b/tests/test_public_exports.py @@ -38,5 +38,6 @@ def test_top_level_reexports_scientific_safety_contracts(): assert callable(saxsabs.derive_fixed_thickness) assert callable(saxsabs.assess_intensity_state) assert callable(saxsabs.require_relative_input_for_absolute_scaling) + assert callable(saxsabs.require_absolute_input_for_fluorescence_subtraction) assert isinstance(saxsabs.XRAYDB_VERSION, str) assert saxsabs.XRAYDB_VERSION diff --git a/tests/test_workbench_scientific.py b/tests/test_workbench_scientific.py index 6eeb4ca..4ec86f5 100644 --- a/tests/test_workbench_scientific.py +++ b/tests/test_workbench_scientific.py @@ -1456,6 +1456,60 @@ def test_formal_buffer_subtraction_has_no_weaker_fallback(monkeypatch): ) +def test_formal_fluorescence_subtraction_has_no_weaker_fallback(monkeypatch): + module = _load_workbench_module() + app = module.SAXSAbsWorkbenchApp.__new__(module.SAXSAbsWorkbenchApp) + monkeypatch.setattr(module, "subtract_fluorescence", None) + with pytest.raises(RuntimeError, match="formal fluorescence subtraction kernel"): + app.subtract_external_absolute_fluorescence( + np.array([0.01, 0.02, 0.03]), + np.array([10.0, 9.0, 8.0]), + np.array([0.2, 0.2, 0.2]), + {"method": "constant", "f0": 1.0, "beta": 1.0}, + ) + + +def test_workbench_merge_keeps_unknown_when_either_step_is_unknown(): + module = _load_workbench_module() + merged = module.SAXSAbsWorkbenchApp.merge_uncertainty_metadata( + { + "uncertainty_model": "buffer-model", + "uncertainty_type": "combined_standard_unknown_alpha", + }, + { + "uncertainty_model": "fluo-model", + "uncertainty_type": "combined_standard", + "fluorescence_method": "constant", + }, + ) + assert merged["uncertainty_type"] == "combined_standard_unknown" + assert "buffer-model" in merged["uncertainty_model"] + assert "fluo-model" in merged["uncertainty_model"] + + +def test_workbench_constant_fluorescence_kernel_and_disabled_default(): + module = _load_workbench_module() + app = module.SAXSAbsWorkbenchApp.__new__(module.SAXSAbsWorkbenchApp) + app.t3_fluo_enabled = _Var(False) + payload = app.prepare_workbench_fluorescence(source="t3", pipeline_mode="scaled") + assert payload["enabled"] is False + + result = app.subtract_external_absolute_fluorescence( + np.array([0.01, 0.02, 0.03]), + np.array([10.0, 9.0, 8.0]), + np.array([0.2, 0.2, 0.2]), + { + "method": "constant", + "f0": 2.0, + "f0_uncertainty": 0.0, + "beta": 1.0, + "beta_uncertainty": 0.0, + "profile": None, + }, + ) + np.testing.assert_allclose(result.i_subtracted, [8.0, 7.0, 6.0]) + + @pytest.mark.parametrize( ("raw", "expected"), [ @@ -1661,6 +1715,117 @@ def counted_read(path): "corrections_applied": ["k", "thickness"], } + +def test_tab3_high_q_fluorescence_does_not_reuse_first_file_f0(tmp_path): + """Prepared high_q payload keeps f0=None for every queue file. + + Mutating the shared fluorescence_info['f0'] after file 1 makes file 2+ + hit the kernel's 'refuse a user f0' guard. + """ + module = _load_workbench_module() + app = module.SAXSAbsWorkbenchApp.__new__(module.SAXSAbsWorkbenchApp) + app.language = "en" + poni = tmp_path / "geometry.poni" + poni.write_text("poni", encoding="utf-8") + context = _calibration_context(module, poni) + fingerprint = context.fingerprint() + + def write_relative(path, values): + ledger = '["thickness"]' + path.write_text( + f"# calibration_context_fingerprint: {fingerprint}\n" + "# thickness_cm: 0.1\n" + "# thickness_source: upstream sample cell record\n" + "# intensity_state: relative\n" + "# intensity_unit: relative\n" + f"# corrections_applied: {ledger}\n" + f"# do_not_repeat: {ledger}\n" + "# q_A^-1 I_rel Error\n" + + "\n".join( + f"{q:.3f} {intensity:.3f} 0.1" + for q, intensity in zip((0.01, 0.02, 0.03), values) + ) + + "\n", + encoding="utf-8", + ) + + sample_a = tmp_path / "fluo_a.dat" + sample_b = tmp_path / "fluo_b.dat" + write_relative(sample_a, (10.0, 9.0, 8.0)) + write_relative(sample_b, (12.0, 11.0, 10.0)) + + app.t3_files = [str(sample_a), str(sample_b)] + app.global_vars = {"k_factor": _Var(2.5)} + app.t3_pipeline_mode = _Var("scaled") + app.t3_corr_mode = _Var("k_only") + app.t3_fixed_thk = _Var(1.0) + app.t3_buffer_enabled = _Var(False) + app.t3_buffer_path = _Var("") + app.t3_alpha = _Var(1.0) + app.t3_alpha_uncertainty = _Var("") + app.t3_buffer_status = _Var("") + app.t3_fluo_enabled = _Var(True) + app.t3_fluo_method = _Var("high_q_mean") + app.t3_fluo_f0 = _Var("") + app.t3_fluo_f0_uncertainty = _Var("0") + app.t3_fluo_beta = _Var(1.0) + app.t3_fluo_beta_uncertainty = _Var("0") + app.t3_fluo_qmin = _Var("0.01") + app.t3_fluo_qmax = _Var("0.03") + app.t3_fluo_path = _Var("") + app.t3_fluo_status = _Var("") + app.t3_output_root = _Var(str(tmp_path / "output")) + app.t3_resume_enabled = _Var(False) + app.t3_overwrite = _Var(False) + app.t3_output_format = _Var("tsv") + app.t3_x_mode = _Var("auto") + app.t3_wavelength_a = _Var("") + app.t3_meta_csv_path = _Var("") + app.t3_bg1d_path = _Var("") + app.t3_dark1d_path = _Var("") + app.t3_prog_bar = {} + app.root = SimpleNamespace(update_idletasks=lambda: None) + app.get_monitor_mode = lambda: "rate" + app.require_trusted_k_for_external = lambda *_args, **_kwargs: context + app.log = lambda _message: None + app.show_info = lambda *_args, **_kwargs: None + app.show_error = lambda _title, message: pytest.fail(message) + + app.t3_preflight_approval = module.approve_preflight( + app._t3_preflight_config(), "READY" + ) + prepared = app.prepare_workbench_fluorescence( + source="t3", + pipeline_mode="scaled", + calibration_context=context, + k_factor=2.5, + ) + assert prepared["enabled"] is True + assert prepared["f0"] is None + assert prepared["method"] == "high_q_mean" + + app.run_external_1d_batch() + + report_dir = tmp_path / "output" / "processed_external_1d_reports" + report_path = next(report_dir.glob("external1d_report_*.csv")) + report = __import__("pandas").read_csv(report_path) + assert report["Status"].tolist() == ["成功", "成功"] + assert report["FluorescenceApplied"].tolist() == [True, True] + assert report["FluorescenceF0"].tolist() == pytest.approx([22.5, 27.5]) + assert prepared["f0"] is None + + out_a = app.read_external_1d_profile( + tmp_path / "output" / "processed_external_1d_abs" / "fluo_a.dat" + ) + out_b = app.read_external_1d_profile( + tmp_path / "output" / "processed_external_1d_abs" / "fluo_b.dat" + ) + np.testing.assert_allclose(out_a["i_abs"], [2.5, 0.0, -2.5]) + np.testing.assert_allclose(out_b["i_abs"], [2.5, 0.0, -2.5]) + assert "fluorescence" in out_a["operator_provenance"]["corrections_applied"] + assert "fluorescence" in out_b["operator_provenance"]["corrections_applied"] + + def _make_complete_custom_record(module, tmp_path): files = {} for name in ("poni", "standard", "background", "dark", "reference"): @@ -1850,6 +2015,43 @@ def test_buffer_combined_uncertainty_and_provenance_roundtrip_all_formats( assert provenance["uncertainty_type"] == "combined_standard" +@pytest.mark.parametrize("output_format", ["tsv", "csv", "cansas_xml", "nxcansas_h5"]) +def test_fluorescence_provenance_roundtrip_all_formats(tmp_path, output_format): + if output_format == "nxcansas_h5": + pytest.importorskip("h5py") + module = _load_workbench_module() + app, _files = _make_complete_custom_record(module, tmp_path) + context = app.require_trusted_k_for_external(2.5) + + written = app.save_profile_table( + tmp_path / "fluo.dat", + np.array([0.01, 0.02, 0.03]), + np.array([10.0, 9.0, 8.0]), + np.array([0.1, 0.2, 0.3]), + "Q_A^-1", + output_format=output_format, + calibration_context=context, + corrections_applied=["fluorescence", "k", "thickness"], + extra_corrections=(), + combined_uncertainty=np.array([0.15, 0.25, 0.35]), + uncertainty_metadata={ + "fluorescence_method": "constant", + "fluorescence_f0": "2.0", + "fluorescence_f0_uncertainty": "0.0", + "fluorescence_beta": "1.0", + "fluorescence_beta_uncertainty": "0.0", + "uncertainty_model": "u_combined^2=u_sample^2+beta^2*u_F^2+F^2*u_beta^2", + "uncertainty_type": "combined_standard", + }, + ) + + profile = app.read_external_1d_profile(written) + provenance = profile["operator_provenance"] + assert provenance["fluorescence_method"] == "constant" + assert provenance["fluorescence_f0"] == "2.0" + assert "fluorescence" in provenance["corrections_applied"] + + @pytest.mark.parametrize("output_format", ["tsv", "csv", "cansas_xml", "nxcansas_h5"]) def test_unknown_buffer_alpha_uncertainty_stays_unknown_all_formats( tmp_path,