feat: BET sensitivity heatmap (BEaTmap-style) in Rouquerol tab

Author: Hj1308

app_bet.py

@@ -30,6 +30,7 @@
     diagnose_instrument_range,
     format_rouquerol_report,
     rouquerol_transform,
+    bet_sensitivity_heatmap,
 )
 
 # ════════════════════════════════════════════════════════════════════════════
@@ -273,6 +274,49 @@ def _plot_rouquerol_transform(p_rel, n, best_window) -> plt.Figure:
     return fig
 
 
+
+def _plot_bet_heatmap(heatmap_result, best_window) -> plt.Figure:
+    """Plot S_BET sensitivity heatmap (BEaTmap-style)."""
+    setup_plot_style()
+    s_bet = heatmap_result["s_bet"]
+    valid = heatmap_result["valid"]
+    p = heatmap_result["p_sorted"]
+    N = heatmap_result["n_points"]
+
+    s_masked = np.ma.masked_where(~valid | ~np.isfinite(s_bet), s_bet)
+
+    fig, ax = plt.subplots(figsize=(7, 5.5))
+    cmap = plt.cm.RdYlGn_r.copy()
+    cmap.set_bad(color="#e0e0e0")
+
+    im = ax.imshow(s_masked, aspect="auto", cmap=cmap,
+                   origin="lower", interpolation="nearest")
+
+    if best_window is not None:
+        p_lo = np.searchsorted(p, best_window.p_min)
+        p_hi = np.searchsorted(p, best_window.p_max)
+        if p_hi > p_lo:
+            rect = plt.Rectangle((p_lo, p_lo), p_hi - p_lo, p_hi - p_lo,
+                                 linewidth=2, edgecolor="blue",
+                                 facecolor="none", linestyle="--")
+            ax.add_patch(rect)
+
+    tick_step = max(1, N // 8)
+    tick_pos = np.arange(0, N, tick_step)
+    tick_labels = [f"{p[i]:.2f}" for i in tick_pos]
+    ax.set_xticks(tick_pos)
+    ax.set_xticklabels(tick_labels, fontsize=8, rotation=45)
+    ax.set_yticks(tick_pos)
+    ax.set_yticklabels(tick_labels, fontsize=8)
+
+    ax.set_xlabel("End point p/p₀")
+    ax.set_ylabel("Start point p/p₀")
+    plt.colorbar(im, ax=ax, label="S_BET (m² g⁻¹)")
+    ax.set_title("BET Sensitivity Heatmap", fontsize=10)
+    plt.tight_layout()
+    return fig
+
+
 def _match_instrument_window_by_pressure(p_ads, n_ads, bet_pts,
                                          start_pt, end_pt):
     """
@@ -415,6 +459,12 @@ def _match_instrument_window_by_pressure(p_ads, n_ads, bet_pts,
             )
         except Exception:
             instrument_window = None
+    heatmap_result = None
+    if rouquerol_result is not None:
+        try:
+            heatmap_result = bet_sensitivity_heatmap(p_ads, n_ads)
+        except Exception:
+            heatmap_result = None
 
 
 # ════════════════════════════════════════════════════════════════════════════
@@ -632,6 +682,19 @@ def _match_instrument_window_by_pressure(p_ads, n_ads, bet_pts,
             })
             st.dataframe(crit_df, use_container_width=True, hide_index=True)
 
+            if heatmap_result is not None:
+                st.divider()
+                st.markdown("**BET Sensitivity Heatmap**")
+                st.caption(
+                    "Each cell shows S_BET for a specific p/p₀ window "
+                    "(start × end). Colored = valid (Rouquerol PASS). "
+                    "Gray = invalid. Blue dashed = selected range."
+                )
+                fig_hm = _plot_bet_heatmap(heatmap_result, best)
+                st.pyplot(fig_hm, use_container_width=True)
+                plt.close(fig_hm)
+
+
             if instrument_window is not None:
                 st.divider()
                 st.markdown("**Instrument Range vs Rouquerol** (matched by p/p₀)")

rouquerol.py

@@ -407,3 +407,63 @@ def format_rouquerol_report(result: dict, sample_name: str = "Sample") -> str:
         ]
     )
     return "\n".join(lines)
+
+
+def bet_sensitivity_heatmap(
+    p_rel: np.ndarray,
+    n: np.ndarray,
+    *,
+    min_points: int = MIN_POINTS_DEFAULT,
+    criterion4_tol: float = CRITERION4_TOL_DEFAULT,
+) -> dict:
+    """Compute S_BET for ALL possible windows — BEaTmap-style sensitivity map.
+
+    For every contiguous window [i, j] with at least ``min_points`` points,
+    the BET fit is evaluated and S_BET, C, R² and Rouquerol validity are
+    stored in 2D matrices.  The result can be plotted as a heatmap to show
+    how sensitive S_BET is to the choice of pressure range.
+
+    Returns
+    -------
+    dict with keys:
+        s_bet    : (N, N) ndarray — S_BET for each window (NaN if too few
+                   points or non-physical)
+        c        : (N, N) ndarray — BET C constant
+        r2       : (N, N) ndarray — R² of the linear fit
+        valid    : (N, N) bool ndarray — Rouquerol-valid windows
+        p_sorted : (N,) ndarray — sorted p/p0 values (axis labels)
+        n_points : int — number of points after filtering
+    """
+    p_rel = np.asarray(p_rel, dtype=float)
+    n = np.asarray(n, dtype=float)
+    if p_rel.ndim != 1 or n.ndim != 1 or len(p_rel) != len(n):
+        raise ValueError("p_rel and n must be 1-D arrays of equal length.")
+
+    order = np.argsort(p_rel)
+    p_s, n_s = p_rel[order], n[order]
+    finite = np.isfinite(p_s) & np.isfinite(n_s) & (p_s > 0) & (p_s < 1) & (n_s > 0)
+    p_s, n_s = p_s[finite], n_s[finite]
+
+    N = len(p_s)
+    s_bet = np.full((N, N), np.nan)
+    c_val = np.full((N, N), np.nan)
+    r2_val = np.full((N, N), np.nan)
+    valid = np.full((N, N), False)
+
+    for i in range(N):
+        for j in range(i + min_points - 1, N):
+            win = evaluate_window(p_s, n_s, i, j, criterion4_tol=criterion4_tol)
+            if win is not None:
+                s_bet[i, j] = win.S_BET
+                c_val[i, j] = win.C
+                r2_val[i, j] = win.R2
+                valid[i, j] = win.valid
+
+    return {
+        "s_bet": s_bet,
+        "c": c_val,
+        "r2": r2_val,
+        "valid": valid,
+        "p_sorted": p_s,
+        "n_points": N,
+    }

tests/test_rouquerol.py

@@ -10,6 +10,7 @@
     monolayer_pressure_theory,
     rouquerol_transform,
     select_bet_range,
+    bet_sensitivity_heatmap,
 )
 
 
@@ -111,3 +112,20 @@ def test_bet_uncertainty_propagation():
     assert fit["sigma_C"] > 0
     S_true = Vm * N2_BET_FACTOR
     assert abs(fit["S_BET"] - S_true) < 3.0 * fit["sigma_S_BET"]
+
+
+def test_bet_sensitivity_heatmap_dimensions_and_stability():
+    """Heatmap returns NxN matrices; ideal isotherm gives stable S_BET."""
+    p, n, C, Vm = _ideal_bet_isotherm()
+    result = bet_sensitivity_heatmap(p, n)
+    N = result["n_points"]
+    assert result["s_bet"].shape == (N, N)
+    assert result["valid"].shape == (N, N)
+    assert result["r2"].shape == (N, N)
+    assert result["valid"].sum() > 0
+    S_true = Vm * N2_BET_FACTOR
+    valid_s = result["s_bet"][result["valid"]]
+    assert np.all(np.abs(valid_s - S_true) / S_true < 0.01)
+    for i in range(N):
+        for j in range(i):
+            assert np.isnan(result["s_bet"][i, j])