fig optimize

ice91 · ice91 · commit 3d75b3521ffc · 2025-09-19T13:04:03.000+08:00
diff --git a/scripts/fig_c1_alignment.py b/scripts/fig_c1_alignment.py
@@ -96,6 +96,18 @@ def run(config: Dict | None = None, which: str = "full") -> Dict[str, List[str]]
     ax.set_xlabel("alignment angle (rad)")
     ax.set_ylabel("PDF")
     ax.set_title("C1: alignment of $T_\\mu$ and $\\partial_\\mu\\epsilon$")
+    # ---- 容差線 & 註記統計 ----
+    ang_tol = float((config or {}).get("c1", {}).get("angle_tol", 1e-3))
+    ax.axvline(ang_tol, ls="--", lw=1.0, color="C1", label=f"tol={ang_tol:g}")
+    mean = float(np.mean(angles))
+    p95  = float(np.percentile(angles, 95))
+    ax.text(
+        0.98, 0.92,
+        f"mean={mean:.2e}\n95%={p95:.2e}",
+        transform=ax.transAxes, ha="right", va="top", fontsize=9
+    )
+    if ang_tol is not None:
+        ax.legend(frameon=False, fontsize=9)
     fig.tight_layout()
     pdf_path = paths["pdf"] / "fig2_c1_alignment.pdf"
     fig.savefig(pdf_path, bbox_inches="tight")
diff --git a/scripts/fig_c1_pure_trace.py b/scripts/fig_c1_pure_trace.py
@@ -112,6 +112,21 @@ def run(config: Dict | None = None, which: str = "full") -> Dict[str, List[str]]
     ax.set_yscale("log")
     ax.set_ylabel("magnitude (log)")
     ax.set_title("C1: torsion components")
+    # ---- 註記數值比 ----
+    eps = 1e-300
+    txt_lines = []
+    if axial > 0:
+        txt_lines.append(f"trace/axial ≈ {trace/max(axial, eps):.2e}")
+    else:
+        txt_lines.append("axial ≈ 0")
+    if tensor > 0:
+        txt_lines.append(f"trace/tensor ≈ {trace/max(tensor, eps):.2e}")
+    else:
+        txt_lines.append("tensor ≈ 0")
+    ax.text(
+        0.98, 0.95, "\n".join(txt_lines),
+        transform=ax.transAxes, ha="right", va="top", fontsize=9
+    )
     fig.tight_layout()
     pdf_path = paths["pdf"] / "fig1_c1_pure_trace.pdf"
     fig.savefig(pdf_path, bbox_inches="tight")
diff --git a/scripts/fig_c2_coeff_compare.py b/scripts/fig_c2_coeff_compare.py
@@ -92,6 +92,23 @@ def run(config: Dict | None = None, which: str = "full") -> Dict[str, List[str]]
     ax.set_ylabel(r"$\|\Delta\|$")
     ax.set_title("C2 residual (smoke)" if which == "smoke" else "C2 residual scan")
     ax.grid(True, ls=":", alpha=0.6)
+    
+    # ---- 目標殘差水平線 + 達成標註 ----
+    resid_goal = float((config or {}).get("equivalence", {}).get("resid_goal", 1e-10))
+    ax.axhline(resid_goal, ls="--", lw=1.0, color="C1", label=f"goal={resid_goal:.1e}")
+    try:
+        # 找到第一個低於目標的 IBP 容差
+        mask = resids <= resid_goal
+        if np.any(mask):
+            th_hit = thresholds[np.argmax(mask)]
+            ax.axvline(th_hit, ls=":", lw=1.0, color="C2")
+            ax.text(
+                th_hit, resid_goal, f"  hit @ {th_hit:.1e}",
+                va="bottom", ha="left", fontsize=8, color="C2"
+            )
+        ax.legend(frameon=False, fontsize=9)
+    except Exception:
+        pass
 
     pdf = paths["pdfdir"] / ("fig3_c2_coeff_compare.pdf" if which != "smoke" else "fig3_c2_coeff_compare_smoke.pdf")
     fig.tight_layout()
diff --git a/scripts/fig_c3_cT_heatmap.py b/scripts/fig_c3_cT_heatmap.py
@@ -148,6 +148,16 @@ def run(config: Dict | None = None, which: str = "full") -> Dict[str, List[str]]
     ax.set_title(r"C3: $|c_T-1|$ heatmap")
     cb = fig.colorbar(im, ax=ax)
     cb.set_label("|c_T-1|")
+    # ---- 等高線：|c_T - 1| = tol ----
+    ct_tol = float((config or {}).get("c3", {}).get("ct_tol", 1e-3))
+    try:
+        CS = ax.contour(
+            X, Y, np.abs(cT - 1.0),
+            levels=[ct_tol], colors="w", linewidths=1.2
+        )
+        ax.clabel(CS, fmt={ct_tol: f"tol={ct_tol:g}"}, inline=True, fontsize=8)
+    except Exception:
+        pass
 
     # optional locking overlay（若使用者提供 locking_curve）
     try:
diff --git a/scripts/fig_c3_degeneracy.py b/scripts/fig_c3_degeneracy.py
@@ -9,6 +9,7 @@
 
 import matplotlib.pyplot as plt
 import numpy as np
+from matplotlib.ticker import MaxNLocator
 
 REQUIRE_REAL = bool(int(os.environ.get("PALPT_REQUIRE_REAL_APIS", "0")))
 
@@ -78,6 +79,20 @@ def run(config: Dict | None = None, which: str = "full") -> Dict[str, List[str]]
     ax.set_title("DoF spectrum (smoke)" if which == "smoke" else "DoF spectrum")
     ax.grid(True, ls=":", alpha=0.5)
 
+    ax.xaxis.set_major_locator(MaxNLocator(integer=True))
+    # ---- 有效非零的閾值與計數 ----
+    deg_tol = float((config or {}).get("c3", {}).get("deg_tol", 1e-10))
+    ax.axhline(deg_tol, ls="--", lw=1.0, color="C1", label=f"tol={deg_tol:.1e}")
+    try:
+        n_phys = int(np.sum(eigs > deg_tol))
+        ax.text(
+            0.98, 0.92, f"#phys ≈ {n_phys}",
+            transform=ax.transAxes, ha="right", va="top", fontsize=9
+        )
+        ax.legend(frameon=False, fontsize=9)
+    except Exception:
+        pass
+
     pdf = paths["pdfdir"] / ("fig6_c3_degeneracy.pdf" if which != "smoke" else "fig6_c3_degeneracy_smoke.pdf")
     fig.tight_layout()
     fig.savefig(pdf, dpi=200)
diff --git a/scripts/fig_c3_dispersion.py b/scripts/fig_c3_dispersion.py
@@ -89,8 +89,22 @@ def run(config: Dict | None = None, which: str = "full") -> Dict[str, List[str]]
     ax.set_ylabel(r"$c_T(k)$")
     ax.set_title("Dispersion (smoke)" if which == "smoke" else "Dispersion")
     ax.grid(True, ls=":", alpha=0.6)
+    #ax.legend()
+    # ---- 容差帶 around 1 ----
+    ct_tol = float((config or {}).get("c3", {}).get("ct_tol", 1e-3))
+    ax.fill_between(k, 1.0-ct_tol, 1.0+ct_tol, alpha=0.12, label=f"±{ct_tol:g} band")
     ax.legend()
-
+    # ---- 註記最大偏差 ----
+    try:
+        d_locked = np.abs(cT_locked - 1.0)
+        d_unlk   = np.abs(cT_unlocked - 1.0)
+        note = (f"max|Δ_locked|={np.max(d_locked):.2e}\n"
+                f"max|Δ_unlocked|={np.max(d_unlk):.2e}")
+        ax.text(0.98, 0.05, note, transform=ax.transAxes,
+                ha="right", va="bottom", fontsize=9)
+    except Exception:
+        pass
+    
     pdf = paths["pdfdir"] / ("fig5_c3_dispersion.pdf" if which != "smoke" else "fig5_c3_dispersion_smoke.pdf")
     fig.tight_layout()
     fig.savefig(pdf, dpi=200)
diff --git a/scripts/fig_flux_ratio.py b/scripts/fig_flux_ratio.py
@@ -53,6 +53,9 @@ def build_from_config(config: dict | None):
     # 交給庫函式做理論值
     out = flux_ratio_FRW(R, {"flux": {"sigma": sigma, "c": c}})
     # out 期望含有：out["R"], out["R_DBI_CM"]
+    R = np.asarray(out["R"], dtype=float)
+    Q = np.asarray(out["R_DBI_CM"], dtype=float)
+    eps = np.finfo(float).eps
 
     # --- CSV ---
     csv_path = os.path.join(DATA_DIR, CSV_NAME)
@@ -65,11 +68,34 @@ def build_from_config(config: dict | None):
     # --- PDF ---
     pdf_path = os.path.join(PDF_DIR, PDF_NAME)
     plt.figure(figsize=(4.8, 3.2), dpi=180)
-    plt.plot(out["R"], out["R_DBI_CM"], lw=1.6, label=r"$\mathcal{R}_{X/Y}(R)$")
+    #plt.plot(out["R"], out["R_DBI_CM"], lw=1.6, label=r"$\mathcal{R}_{X/Y}(R)$")
+    plt.plot(R, Q, lw=1.6, label=r"$\mathcal{R}_{X/Y}(R)$")
     plt.axhline(1.0, ls="--", lw=1.0, label="unity")
     plt.xlabel(r"$R$")
     plt.ylabel(r"$\mathcal{R}_{X/Y}$")
     plt.title(r"Boundary flux ratio $\to 1$ with $R^{-\sigma}$ tail")
+    #plt.legend(frameon=False)
+    # ---- 擬合 R^{-sigma} 的收斂率 ----
+    try:
+        fit_Rmin = float(_cfg_get(config, ["flux", "fit", "Rmin"], R.min() + 0.5*(R.max()-R.min())))
+        mask = (R >= fit_Rmin) & (np.abs(Q-1.0) > 1e-14)
+        if np.any(mask):
+            X = np.log(R[mask])
+            Y = np.log(np.abs(Q[mask] - 1.0) + eps)
+            m, b = np.polyfit(X, Y, 1)     # Y ≈ m X + b
+            sigma_hat = -m
+            A_hat = np.exp(b)
+            # 以尾端平均號誌的符號還原 ±
+            sgn = np.sign(np.mean(Q[mask] - 1.0))
+            Qfit = 1.0 + sgn * A_hat * (R**(-sigma_hat))
+            plt.plot(R, Qfit, ls=":", lw=1.2, label=rf"fit: $1+A R^{{-\hat\sigma}}$, $\hat\sigma$≈{sigma_hat:.2f}")
+            plt.text(
+                0.02, 0.05, rf"fit window: $R\ge {fit_Rmin:g}$",
+                transform=plt.gca().transAxes, fontsize=8
+            )
+    except Exception:
+        pass
+
     plt.legend(frameon=False)
     plt.tight_layout()
     plt.savefig(pdf_path)
diff --git a/scripts/fig_gw_waveform_overlay.py b/scripts/fig_gw_waveform_overlay.py
@@ -118,6 +118,40 @@ def run(config: Dict | None = None, which: str = "full") -> Dict[str, List[str]]
     ax.set_ylabel("h")
     ax.set_title("GW waveform overlay")
     ax.legend()
+
+    # ---- 指標：max|Δh|、RMSE、相位差(以主頻估) ----
+    try:
+        h1 = d["gr"].astype(float)
+        h2 = d["model"].astype(float)
+        tt = d["t"].astype(float)
+        # 幅度差與 RMSE
+        max_abs = float(np.max(np.abs(h2 - h1)))
+        rmse = float(np.sqrt(np.mean((h2 - h1)**2)))
+        # 估時間位移：互相關最大點
+        s1 = h1 - np.mean(h1)
+        s2 = h2 - np.mean(h2)
+        corr = np.correlate(s1, s2, mode="full")
+        lag  = int(np.argmax(corr) - (len(s1) - 1))
+        dt   = float((tt[1] - tt[0]) * lag)
+        # 主頻：用 FFT 取最大頻率
+        n = len(h1)
+        dt0 = float(tt[1]-tt[0])
+        freq = np.fft.rfftfreq(n, d=dt0)
+        amp  = np.abs(np.fft.rfft(h1))
+        if len(freq) > 1:
+            ipeak = int(np.argmax(amp[1:]) + 1)
+            f0 = float(freq[ipeak])
+        else:
+            f0 = 0.0
+        phase = float(2*np.pi*f0*dt) if f0 > 0 else 0.0
+        note = (f"max|Δh|={max_abs:.2e}\n"
+                f"RMSE={rmse:.2e}\n"
+                f"Δt≈{dt:.2e},  Δφ≈{phase:.2e} rad")
+        ax.text(0.98, 0.05, note, transform=ax.transAxes,
+                ha="right", va="bottom", fontsize=9)
+    except Exception:
+        pass
+    
     fig.tight_layout()
     pdf_path = paths["pdf"] / "fig7_gw_waveform_overlay.pdf"
     fig.savefig(pdf_path, bbox_inches="tight")
diff --git a/scripts/fig_nlo_offsets.py b/scripts/fig_nlo_offsets.py
@@ -87,6 +87,23 @@ def _plot_pdf(out: dict) -> str:
         plt.loglog(kpos, ref, ls="--", lw=1.0, label=r"$\propto k^2$")
 
         plt.ylabel(r"$|\delta c_T^2|$")
+        # ---- 斜率擬合（log–log）與 Λ 估計 ----
+        try:
+            X = np.log(kpos)
+            Y = np.log(ypos)
+            m, b = np.polyfit(X, Y, 1)   # Y ≈ m X + b
+            # 從 δcT^2 = bcoef * k^2 / Λ^2 推回 Λ^2 的 robust 估計
+            # 需要配置裡的 bcoef（名為 Ricci_deps_deps_eff）
+            bcoef = float(_cfg_get(globals().get("CFG_CACHE", None), ["nlo", "coeffs", "Ricci_deps_deps_eff"], 0.0))
+            if bcoef > 0:
+                lam2_est = np.median(bcoef*(kpos**2)/ypos)
+                lam_est  = np.sqrt(lam2_est)
+                txt = rf"slope≈{m:.2f},  $\hat\Lambda$≈{lam_est:.2e}"
+            else:
+                txt = rf"slope≈{m:.2f}"
+            plt.text(0.03, 0.05, txt, transform=plt.gca().transAxes, fontsize=9)
+        except Exception:
+            pass
     else:
         # 全為 0 或只有 1 個正值：避免 log-scaling 警告，改線性 y 軸
         plt.plot(k, yabs, lw=1.6, label=r"$|\delta c_T^2(k)|\approx 0$")
