import json import math from pathlib import Path import numpy as np from scipy.stats import linregress # Dissipation-Budgeted Nonreversible Sampling: minimal first-passage MVP. # Dynamics: dX = [-V'(X) + u] dt + sqrt(2D) dW, # V(x)=(x^2-1)^2/4, target is the right well x >= 0.8. SEED = 3131 DT = 0.002 TMAX = 12.0 D = 0.08 X0 = -1.0 TARGET = 0.8 N = 18000 SCALES = [0.0, 0.15, 0.30, 0.50, 0.80] BUDGETS = [0.02, 0.05, 0.15, 0.50] def force(x): # -V'(x) = x - x^3 return x - x**3 def simulate(u, budget=np.inf, seed=0): rng = np.random.default_rng(seed) nstep = int(round(TMAX / DT)) x = np.full(N, X0, dtype=np.float64) alive = np.ones(N, dtype=bool) hit = np.full(N, TMAX + DT, dtype=np.float64) q = np.zeros(N, dtype=np.float64) # For a constant applied drift, q increment is exactly u^2/(2D) dt. qrate = u * u / (2.0 * D) # Store survival at regular reporting times for an independent hazard check. report_every = 50 survival = [] times = [] for k in range(nstep): t = k * DT active = alive if np.any(active): # A cumulative budget is enforced before each step. Partial final # drift is represented by a multiplicative control gate. remaining = np.maximum(budget - q[active], 0.0) gate = np.minimum(1.0, remaining / (qrate * DT)) if qrate > 0 else np.ones(np.sum(active)) drift_u = u * gate q[active] += (drift_u * drift_u / (2.0 * D)) * DT noise = np.sqrt(2.0 * D * DT) * rng.standard_normal(np.sum(active)) xa = x[active] xnew = xa + (force(xa) + drift_u) * DT + noise x[active] = xnew entered = xnew >= TARGET inds = np.flatnonzero(active) if np.any(entered): ii = inds[entered] # Linear interpolation makes first-passage estimates less dt-dependent. prev = xa[entered] frac = np.clip((TARGET - prev) / (xnew[entered] - prev + 1e-30), 0, 1) hit[ii] = t + frac * DT alive[ii] = False if (k + 1) % report_every == 0: times.append((k + 1) * DT) survival.append(np.mean(alive)) return {"u": u, "budget": budget, "hit": hit, "q": q, "times": np.asarray(times), "survival": np.asarray(survival)} def summarize(r): hit = r["hit"] <= TMAX ht = r["hit"][hit] # Censored mean is not used as the primary speed metric. Estimate a late # approximately constant hazard from log survival, with a small floor. t = r["times"] s = np.maximum(r["survival"], 1.0 / N) use = (t >= 3.0) & (s > 0.03) & (s < 0.98) if np.sum(use) >= 3: fit = linregress(t[use], np.log(s[use])) rate = max(0.0, -fit.slope) r2 = fit.rvalue**2 else: rate, r2 = float("nan"), float("nan") # Empirical all-time cumulative hazard identity: H=-log S. s0 = max(r["survival"][0], 1.0 / N) send = max(r["survival"][-1], 1.0 / N) H = -math.log(send / s0) # q is pathwise integrated only until first hit, as required by stopping. mean_q_hit = float(np.mean(r["q"][hit])) if np.any(hit) else float("nan") return {"u": r["u"], "budget": r["budget"], "hit_fraction": float(np.mean(hit)), "mean_hit_time": float(np.mean(ht)) if np.any(hit) else None, "median_hit_time": float(np.median(ht)) if np.any(hit) else None, "hazard_tail": float(rate), "hazard_tail_r2": float(r2), "mean_q_hit": mean_q_hit, "cum_hazard": H, "log_survival_ratio": float(-math.log(send / s0))} def main(): out = Path("results.json") # Fixed seeds make every control reproducible. Independent seeds avoid # pretending that noisy first-passage paths are exactly paired. runs = [] for j, u in enumerate(SCALES): runs.append(simulate(u, np.inf, SEED + j)) for j, b in enumerate(BUDGETS): runs.append(simulate(0.5, b, SEED + 100 + j)) summaries = [summarize(r) for r in runs] # Verify S(t)=exp(-integral hazard) in the discrete empirical sense using # interval hazard h_k=-log(S_{k+1}/S_k)/dt. base = runs[0] ss = np.maximum(base["survival"], 1.0 / N) h = -np.diff(np.log(ss)) / np.diff(base["times"]) H_from_h = float(np.sum(h * np.diff(base["times"]))) H_direct = float(-np.log(ss[-1] / ss[0])) survival_identity_error = abs(H_from_h - H_direct) # Check Q formula on an uncapped constant-drift run: q should equal # u^2*T_hit/(2D) up to the final discretization/interpolation step. unc = runs[SCALES.index(0.5)] valid = unc["hit"] <= TMAX expected_q = 0.5**2 * unc["hit"][valid] / (2 * D) q_error = float(np.max(np.abs(unc["q"][valid] - expected_q))) if np.any(valid) else None # Falsifiable functional-form check: log hazard improvement against Q. sweep = [x for x in summaries[:len(SCALES)] if x["u"] > 0 and np.isfinite(x["hazard_tail"])] cal = [x for x in sweep if x["hazard_tail"] > 0 and x["mean_q_hit"] > 0] if len(cal) >= 3: xs = np.array([x["mean_q_hit"] for x in cal]) ys = np.log(np.array([x["hazard_tail"] for x in cal]) / max(summaries[0]["hazard_tail"], 1e-12)) lr = linregress(xs, ys) q_fit = {"slope_C": float(lr.slope), "intercept": float(lr.intercept), "r2": float(lr.rvalue**2), "n": len(cal)} else: q_fit = {"slope_C": None, "intercept": None, "r2": None, "n": len(cal)} result = {"config": {"N": N, "dt": DT, "Tmax": TMAX, "D": D, "target": TARGET, "seed": SEED}, "summaries": summaries, "math_checks": {"survival_identity_abs_error": survival_identity_error, "q_constant_drift_max_abs_error": q_error, "log_hazard_vs_q_fit": q_fit}, "interpretation": { "baseline": summaries[0], "uncapped_u_0.5": summaries[SCALES.index(0.5)], "budgeted": summaries[len(SCALES):]}} out.write_text(json.dumps(result, indent=2)) print(json.dumps(result, indent=2)) if __name__ == "__main__": main()