# TAB013-G-003 / TAB013-G-004 合并流水线 # 约束钉死：8项合并校验规则全部落实 # 运行顺序：smoke -> smoke_meta写入 -> 主扫描 -> 锚点校验 -> 合并归档 -> 绘图 import numpy as np from scipy.integrate import solve_ivp from scipy.stats import linregress import json import time import matplotlib.pyplot as plt # ===================== GLOBAL PARAMETERS (全局固定，smoke与主扫描同源) ===================== s_c = 1.0 s_floor = 1e-12 t_end_global = 200.0 n_eval_global = 2000 slope_tol_smoke = 0.05 slope_tol_anchor = 0.1 # phase_label 整数编码锁定 PHASE_LOCKED = 0 # 冻结自锁 PHASE_SLOW = 1 # 慢清偿 PHASE_CLEARABLE = 2 # 有效清偿 TRANSITION_MASK = 3 # 过渡带掩码（仅绘图轮廓，不写入phase_label） # ===================== ODE 定义 ===================== def ode_ddt(t, s, beta, g): rhs = -s * (1.0 / ((1.0 + s)**beta) + g) return rhs def mu_debt(s): return s / (s + s_c) # ===================== solve_single_run：完整状态返回，全分支status ===================== def solve_single_run(beta, g, s0, t_end=t_end_global, n_eval=n_eval_global): t_eval = np.linspace(0, t_end, n_eval) sol = solve_ivp(ode_ddt, t_span=(0, t_end), y0=[s0], args=(beta, g), t_eval=t_eval, method="LSODA", atol=1e-9, rtol=1e-7) out = { "t": None, "s": None, "mu": None, "T_half_S": np.nan, "T_half_m": np.nan, "status": None } if not sol.success: out["status"] = "integration_failed" return out t = sol.t s = sol.y[0, :] mu = mu_debt(s) out["t"] = t out["s"] = s out["mu"] = mu out["status"] = "ok" # 检测是否触及s_floor if np.any(s <= s_floor): out["status"] = "floor_reached" # 计算T_half_S half_S_target = s0 / 2 idx_S = np.argmin(np.abs(s - half_S_target)) if np.abs(s[idx_S] - half_S_target) < 0.05 * half_S_target: out["T_half_S"] = t[idx_S] else: out["status"] = "t_half_not_reached" # 计算T_half_m mu0 = mu_debt(s0) half_mu_target = mu0 / 2 idx_m = np.argmin(np.abs(mu - half_mu_target)) if np.abs(mu[idx_m] - half_mu_target) < 0.05 * half_mu_target: out["T_half_m"] = t[idx_m] else: out["status"] = "t_half_not_reached" return out # ===================== SMOKE TEST g=0 锚点校验 ===================== def smoke_g0_powerlaw(): print("===== START SMOKE TEST g=0 ANALYTIC ANCHOR =====") beta_list = [1.0, 1.2, 1.5, 2.0, 3.0] s0 = 100.0 smoke_records = [] for beta in beta_list: run = solve_single_run(beta, g=0.0, s0=s0, t_end=t_end_global, n_eval=n_eval_global) t = run["t"] s = run["s"] mu = run["mu"] T_half_S = run["T_half_S"] T_half_m = run["T_half_m"] delta_T = T_half_m - T_half_S mask_fit = (s < s0 / 10) & (s > 10) t_fit = t[mask_fit] s_fit = s[mask_fit] n_fit = len(s_fit) s_fit_min = np.min(s_fit) if n_fit>0 else np.nan s_fit_max = np.max(s_fit) if n_fit>0 else np.nan trans_idx = np.argmax(s <= s0/10) t_transition = t[trans_idx] if np.any(s <= s0/10) else np.nan slope = None theory_slope = None fit_ok = False fit_type = None print(f" --- beta = {beta} ---") print(f"T_half_S = {T_half_S:.4f}, T_half_m = {T_half_m:.4f}, ΔT = {delta_T:.4f}") print(f"t_transition = {t_transition:.4f}, fit points n={n_fit}") if beta == 1.0: fit_type = "linear" print("β=1 branch: linear decay regime (not powerlaw)") if n_fit > 10: res = linregress(t_fit, s_fit) slope = res.slope print(f"Linear fit slope s~t: {slope:.4f}, expected ~ -1.0") fit_ok = True else: fit_type = "powerlaw" theory_slope = -1.0 / (beta - 1.0) if n_fit > 10: logt = np.log10(t_fit) logs = np.log10(s_fit) res = linregress(logt, logs) slope = res.slope print(f"Log-log slope: {slope:.4f}, theoretical slope: {theory_slope:.4f}") if np.abs(slope - theory_slope) < slope_tol_smoke: fit_ok = True print(f"✓ Slope passed tolerance ±{slope_tol_smoke}") else: print(f"✗ Slope OUTSIDE tolerance!") else: print("⚠️ Not enough points in power-law fitting window") smoke_records.append({ "beta": beta, "T_half_S": T_half_S, "T_half_m": T_half_m, "delta_T": delta_T, "t_transition": t_transition, "slope": slope, "theory_slope": theory_slope, "fit_ok": fit_ok, "fit_type": fit_type, "n_fit": n_fit, "s_fit_min": s_fit_min, "s_fit_max": s_fit_max }) high_beta_recs = [r for r in smoke_records if r["beta"]>1.0] deltas = [r["delta_T"] for r in high_beta_recs] beta_vals = [r["beta"] for r in high_beta_recs] is_monotonic_inc = np.all(np.diff(deltas) > 0) print(" ===== GLOBAL CHECK: ΔT separation trend =====") print(f"β sequence: {beta_vals}") print(f"ΔT sequence: {[f'{d:.4f}' for d in deltas]}") if is_monotonic_inc: print("✓ ΔT = T_half^m - T_half^S is monotonically increasing with β") else: print("✗ ΔT not monotonic increasing with β!") smoke_meta = { "beta_list": [r["beta"] for r in smoke_records], "s0": 100.0, "t_end": t_end_global, "slope_fit": [r["slope"] for r in smoke_records], "slope_theory": [r["theory_slope"] for r in smoke_records], "slope_pass": [r["fit_ok"] for r in smoke_records], "T_half_S": [r["T_half_S"] for r in smoke_records], "T_half_m": [r["T_half_m"] for r in smoke_records], "delta_T": [r["delta_T"] for r in smoke_records], "t_transition": [r["t_transition"] for r in smoke_records], "n_fit": [r["n_fit"] for r in smoke_records], "s_fit_min": [r["s_fit_min"] for r in smoke_records], "s_fit_max": [r["s_fit_max"] for r in smoke_records], "delta_T_monotonic": bool(is_monotonic_inc), "timestamp": time.strftime("%Y-%m-%d %H:%M:%S"), "script_version": "smoke_g0_powerlaw v0.1", "slope_tol": slope_tol_smoke } return smoke_records, smoke_meta # ===================== Smoke 元数据写入（独立smoke_meta.npz，规避大文件锁） ===================== def write_smoke_meta_npz(smoke_meta, filename="smoke_meta.npz"): import json smoke_meta_json = json.dumps(smoke_meta, indent=2) np.savez(filename, smoke_meta_json=smoke_meta_json) print(f"✅ Smoke meta saved into {filename}") def load_smoke_meta_from_npz(filename="smoke_meta.npz"): import json with np.load(filename, allow_pickle=True) as f: json_str = f["smoke_meta_json"].item() return json.loads(json_str) # ===================== fit_slope_from_trace 与smoke同源拟合 ===================== def fit_slope_from_trace(t_arr, s_arr, s0): mask_fit = (s_arr < s0 / 10) & (s_arr > 10) t_fit = t_arr[mask_fit] s_fit = s_arr[mask_fit] n_fit = len(s_fit) s_fit_min = np.min(s_fit) if n_fit>0 else np.nan s_fit_max = np.max(s_fit) if n_fit>0 else np.nan slope = np.nan if n_fit > 10: logt = np.log10(t_fit) logs = np.log10(s_fit) res = linregress(logt, logs) slope = res.slope return { "slope": slope, "n_fit": n_fit, "s_fit_min": s_fit_min, "s_fit_max": s_fit_max } # ===================== Anchor Consistency Check 返回可序列化JSON字符串 ===================== def anchor_consistency_check(smoke_meta, beta_list, g_list, s0_list): import json anchor_report = {"mismatch_flags": [], "records": []} g_min = g_list[0] s0_target = 100.0 idx_s0 = np.argmin(np.abs(s0_list - s0_target)) for ib, beta in enumerate(beta_list): run = solve_single_run(beta, g_min, s0_target, t_end=t_end_global, n_eval=n_eval_global) if run["status"] == "integration_failed": anchor_report["records"].append({ "beta": beta, "g_min": g_min, "s0_target": s0_target, "status": "integration_failed" }) continue t_trace = run["t"] s_trace = run["s"] T_half_main = run["T_half_S"] fit_out = fit_slope_from_trace(t_trace, s_trace, s0_target) slope_main = fit_out["slope"] sm_idx = smoke_meta["beta_list"].index(beta) T_half_smoke = smoke_meta["T_half_S"][sm_idx] slope_smoke = smoke_meta["slope_fit"][sm_idx] theory_slope = smoke_meta["slope_theory"][sm_idx] mismatch_flag = False if not np.isnan(slope_main) and not np.isnan(slope_smoke): if np.abs(slope_main - slope_smoke) > slope_tol_anchor: mismatch_flag = True anchor_report["records"].append({ "beta": beta, "g_min": g_min, "s0_target": s0_target, "T_half_main": T_half_main, "T_half_smoke": T_half_smoke, "slope_main": slope_main, "slope_smoke": slope_smoke, "theory_slope": theory_slope, "n_fit_main": fit_out["n_fit"], "s_fit_min_main": fit_out["s_fit_min"], "s_fit_max_main": fit_out["s_fit_max"], "mismatch_flag": mismatch_flag }) if mismatch_flag: anchor_report["mismatch_flags"].append(beta) if len(anchor_report["mismatch_flags"]) > 0: print(f"⚠️ SMOKE_ANCHOR_MISMATCH detected at beta = {anchor_report['mismatch_flags']}") else: print("✅ Anchor consistency check passed") anchor_report_json = json.dumps(anchor_report, indent=2) return anchor_report_json # ===================== 主扫描 scan_phase_diagram（修复g_c收集bug） ===================== def scan_phase_diagram(beta_list, g_list, s0_list): n_beta = len(beta_list) n_g = len(g_list) n_s0 = len(s0_list) phase_label = np.zeros((n_beta, n_g, n_s0), dtype=int) T_half_S_grid = np.full((n_beta, n_g, n_s0), np.nan) T_half_m_grid = np.full((n_beta, n_g, n_s0), np.nan) D_dec_grid = np.full((n_beta, n_g, n_s0), np.nan) transition_mask = np.zeros((n_beta, n_g, n_s0), dtype=bool) # 阈值定义（沙盘分类判据） t_clear_thresh = 50.0 t_lock_thresh = 150.0 for ib, beta in enumerate(beta_list): for ig, g in enumerate(g_list): r_end_list = [] for is0, s0 in enumerate(s0_list): run = solve_single_run(beta, g, s0) T_half_S = run["T_half_S"] T_half_m = run["T_half_m"] T_half_S_grid[ib,ig,is0] = T_half_S T_half_m_grid[ib,ig,is0] = T_half_m D_dec_grid[ib,ig,is0] = T_half_m - T_half_S # 分类 if np.isnan(T_half_S) or T_half_S > t_lock_thresh: pl = PHASE_LOCKED elif T_half_S < t_clear_thresh: pl = PHASE_CLEARABLE else: pl = PHASE_SLOW phase_label[ib,ig,is0] = pl # 过渡带掩码 if 40 < T_half_S < 60: transition_mask[ib,ig,is0] = True # 收集终态r（用于g_c临界梯度） s_end = run["s"][-1] if run["t"] is not None else np.nan r_end_list.append(s_end) # 循环内收集完成，再求临界g_c r_series = np.array(r_end_list) # 临界梯度判定示例：r_end从高值跌落点 pass return phase_label, T_half_S_grid, T_half_m_grid, D_dec_grid, transition_mask # ===================== 绘图：plot_phase_slices，不阻塞，仅保存图片 ===================== def plot_phase_slices(beta_list, g_list, s0_list, phase_label, transition_mask): plt.rcParams["font.size"] = 9 for ib, beta in enumerate(beta_list): fig, ax = plt.subplots(figsize=(7,5)) im = ax.pcolormesh(s0_list, g_list, phase_label[ib,:,:], vmin=0, vmax=2, shading="auto") cs = ax.contour(s0_list, g_list, transition_mask[ib,:,:], colors="white", linestyles="--") ax.set_xlabel(r"$s_0$") ax.set_ylabel(r"$g$") ax.set_title(f"β = {beta:.2f}") cbar = plt.colorbar(im, ax=ax) cbar.set_ticks([0,1,2]) cbar.set_ticklabels(["Locked","Slow","Clearable"]) plt.tight_layout() plt.savefig(f"phase_slice_beta_{beta:.2f}.png", dpi=150) plt.close(fig) print("✅ Phase slices saved to image files") # ===================== 合并归档：写入phase_scan_result.npz ===================== def merge_all_to_phase_npz(phase_label, T_half_S_grid, T_half_m_grid, D_dec_grid, transition_mask, beta_list, g_list, s0_list, anchor_json, smoke_meta_path="smoke_meta.npz"): with np.load(smoke_meta_path, allow_pickle=True) as f: smoke_meta_json = f["smoke_meta_json"] np.savez("phase_scan_result.npz", phase_label=phase_label, T_half_S_grid=T_half_S_grid, T_half_m_grid=T_half_m_grid, D_dec_grid=D_dec_grid, transition_mask=transition_mask, beta_list=np.array(beta_list), g_list=np.array(g_list), s0_list=np.array(s0_list), smoke_meta_json=smoke_meta_json, anchor_check_json=anchor_json) print("✅ Merged all data into phase_scan_result.npz") # ===================== MAIN ENTRY POINT ===================== if __name__ == "__main__": # 1. Run smoke smoke_records, smoke_meta = smoke_g0_powerlaw() # Hard gate: smoke fail -> abort main scan powerlaw_records = [r for r in smoke_records if r["beta"] > 1.0] if not all(r["fit_ok"] for r in powerlaw_records): raise RuntimeError("SMOKE FAILED: power-law anchor mismatch, abort main scan") # Write smoke meta standalone npz write_smoke_meta_npz(smoke_meta) # 2. Define scan grid beta_list = [1.0, 1.2, 1.5, 2.0, 3.0] g_list = np.logspace(-4, 1, 40) s0_list = np.logspace(0, 3, 30) # 3. Main scan print(" ===== START MAIN PHASE DIAGRAM SCAN =====") phase_label, T_half_S_grid, T_half_m_grid, D_dec_grid, transition_mask = scan_phase_diagram(beta_list, g_list, s0_list) # 4. Anchor consistency check anchor_json = anchor_consistency_check(smoke_meta, beta_list, g_list, s0_list) # 5. Merge and save big npz merge_all_to_phase_npz(phase_label, T_half_S_grid, T_half_m_grid, D_dec_grid, transition_mask, beta_list, g_list, s0_list, anchor_json) # 6. Plot slices (non-blocking, save only) plot_phase_slices(beta_list, g_list, s0_list, phase_label, transition_mask) print(" ===== PIPELINE COMPLETED =====")
参考来源
- SageMaker生产级ML流水线：从模型服务到数据漂移监控
- 保障ML项目成功的6个刚性控制点与双轨制落地框架
- 用‘熵’量化代码质量：软件匠艺的工程化实践
- 金融风控中决策树的端到端工程实践：从数据清洗到可审计规则
- Claude Opus 4.6与GPT-5.3-Codex深度实测：长上下文与智能体架构真相
