""" Physical mixture model. Replaces the constant-σ Gaussian-bump ansatz of v3 with per-(p,q) widths σ_{p,q} = σ_0 · ν(p,q) where ν encodes how the resonance libration width scales with resonance order (= p − q for coprime p > q ≥ 1, with integer ratios p:1 treated as order = p − 1). We test three width-scaling laws: M0 ν(p,q) = 1 (constant; v3 baseline) M1 ν(p,q) = 1 / max(p − q, 1) (∝ 1 / order) M2 ν(p,q) = exp(−(p − q − 1)) (exponential suppression with order) For each model the two free parameters are α (Gauss-component weight) and σ_0 (absolute resonance scale). The likelihood-ratio test on M_k vs M_0 quantifies whether order-dependent broadening fits the data better. We also no longer collapse integer ratios (2:1, 3:1, 4:1, …) into a single center at {r}=0 — each gets its own bump with its own width, because they have different libration widths even though they share fractional part. Two further deliverables in this script: (a) The "favoured Farey" table: for each symbol, the posterior weight once the model is fit, identifying which low-order resonances carry the trap mass. (b) The Lithwick–Wu asymmetry test: fit *paired* Gaussians at (p/q − ε, p/q + ε) for first-order resonances (2:1, 3:2, 4:3, …) with one extra parameter ε. If ε > 0 significantly, this is the asymmetric pile-up + gap signature predicted by Lithwick & Wu (2012). """ import math import csv from collections import defaultdict, Counter import numpy as np import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt LN2 = math.log(2.0) # ----------------------- data load & helpers -------------------------- def farey_pairs(max_sum): out = [] for q in range(1, max_sum): for p in range(q + 1, max_sum - q + 1): if math.gcd(p, q) == 1: out.append((p, q)) return out def build_farey_alphabet(max_sum): """Per-(p,q) record: (p, q, fractional_part_in_[0,1), order).""" out = [] for p, q in farey_pairs(max_sum): order = p - q # = 1 for first-order MMRs (2:1, 3:2, 4:3, …) center = (p % q) / q if q > 1 else 0.0 out.append((p, q, center, order)) return out def load_ratios(path): systems = defaultdict(list) with open(path) as f: for row in csv.DictReader(f): try: p = float(row["pl_orbper"]) except (ValueError, TypeError): continue if math.isfinite(p) and p > 0: systems[row["hostname"]].append(p) out = [] for periods in systems.values(): if len(periods) < 2: continue periods = sorted(periods) for i in range(len(periods) - 1): r = periods[i + 1] / periods[i] if r > 1 and math.isfinite(r): out.append(r) return np.array(out) def gauss_density(x): return 1.0 / ((1.0 + x) * LN2) # ----------------------- width scaling laws -------------------------- WIDTH_MODELS = { "M0_constant": lambda p, q: 1.0, "M1_inv_order": lambda p, q: 1.0 / max(p - q, 1), "M2_exp_order": lambda p, q: math.exp(-(p - q - 1)), } def widths_for_model(alphabet, sigma_0, model_name): nu = WIDTH_MODELS[model_name] return np.array([sigma_0 * nu(p, q) for p, q, _, _ in alphabet]) # ----------------------- density (per-sigma version) ------------------ def farey_density_per_sigma(x, centers, weights, sigmas, n_wraps=3): """Vectorized: each center has its own sigma. x: (N,). Returns (N,).""" ks = np.arange(-n_wraps, n_wraps + 1) delta = x[:, None, None] - centers[None, :, None] + ks[None, None, :] # (N, K, 2W+1) inv_2s2 = (0.5 / (sigmas ** 2))[None, :, None] coef = (1.0 / (sigmas * math.sqrt(2 * math.pi)))[None, :, None] kernel = coef * np.exp(-(delta ** 2) * inv_2s2) per_center = kernel.sum(axis=2) return (per_center * weights[None, :]).sum(axis=1) def farey_normalizer_per_sigma(centers, weights, sigmas, n_grid=2001): xs = np.linspace(0, 1, n_grid) ys = farey_density_per_sigma(xs, centers, weights, sigmas) return np.trapezoid(ys, xs) def fit_physical_mixture(fracs, alphabet, alpha_grid, sigma_0_grid, model_name): centers = np.array([c for _, _, c, _ in alphabet]) weights = np.ones(len(alphabet)) / len(alphabet) base_vals = gauss_density(fracs) ll_grid = np.full((len(alpha_grid), len(sigma_0_grid)), -np.inf) for j, sigma_0 in enumerate(sigma_0_grid): sigmas = widths_for_model(alphabet, sigma_0, model_name) if np.any(sigmas <= 0): continue f = farey_density_per_sigma(fracs, centers, weights, sigmas) Z = farey_normalizer_per_sigma(centers, weights, sigmas) if Z <= 0: continue f_norm = f / Z for i, alpha in enumerate(alpha_grid): p = alpha * base_vals + (1 - alpha) * f_norm with np.errstate(divide="ignore", invalid="ignore"): ll_grid[i, j] = np.log(p).sum() idx = np.unravel_index(np.argmax(ll_grid), ll_grid.shape) return ll_grid, idx, ll_grid[idx], alpha_grid[idx[0]], sigma_0_grid[idx[1]] # ----------------------- Lithwick–Wu asymmetric test ----------------- def fit_lw_pair(fracs, centers_int, sigma_grid, eps_grid, weight=1.0): """For each first-order integer-adjacent resonance (2:1, 3:1, 4:1, …), centers_int is the fractional position (=0 in our convention for r near an integer). Fit *paired* Gaussians at (-ε, +ε) wrapped around 0: f(x) = 0.5 [N(x; −ε, σ) + γ N(x; +ε, σ)] where γ is the inside-vs-outside amplitude ratio. γ > 1 = excess inside, γ < 1 = excess outside; Lithwick–Wu predicts γ < 1 (pile-up just *below* the resonance ratio, i.e. at fractional part just under 1). """ pass # placeholder; LW analysis done in main() differently # ----------------------- main ---------------------------------------- def main(): rs = load_ratios("ps_dm.csv") fracs = rs - np.floor(rs) fracs = fracs[(fracs >= 0) & (fracs < 1)] print(f"n ratios: {len(rs)}; n fractional parts: {len(fracs)}") base_vals = gauss_density(fracs) ll_pure_gauss = np.log(base_vals).sum() print(f"pure Gauss LL: {ll_pure_gauss:.2f}") alphabet = build_farey_alphabet(max_sum=12) print(f"\nFarey-12 alphabet (22 symbols, per-(p,q) widths):") print(f" {'symbol':>8} {'frac':>7} {'order':>5} M0 M1 M2") for p, q, c, order in alphabet: m0 = WIDTH_MODELS["M0_constant"](p, q) m1 = WIDTH_MODELS["M1_inv_order"](p, q) m2 = WIDTH_MODELS["M2_exp_order"](p, q) sym = f"{p}:{q}" print(f" {sym:>8} {c:7.4f} {order:5d} {m0:.2f} {m1:.3f} {m2:.4f}") # ---------- fit each model ---------- alpha_grid = np.linspace(0.0, 1.0, 51) sigma_0_grid = np.exp(np.linspace(math.log(0.005), math.log(0.3), 41)) print(f"\n--- fitting {len(WIDTH_MODELS)} width-scaling models ---") results = {} for model_name in WIDTH_MODELS: ll_grid, idx, ll_best, a_best, s0_best = fit_physical_mixture( fracs, alphabet, alpha_grid, sigma_0_grid, model_name ) dll = ll_best - ll_pure_gauss results[model_name] = (ll_grid, a_best, s0_best, ll_best, dll) print(f" {model_name:>15s}: α={a_best:.3f} σ₀={s0_best:.4f} " f"LL={ll_best:.2f} ΔLL = {dll:+.2f}") # ---------- model selection (LRT) ---------- print("\n--- model selection (LRT pairwise) ---") lls = {m: r[3] for m, r in results.items()} print(f" 2·(LL[M1] − LL[M0]) = {2*(lls['M1_inv_order'] - lls['M0_constant']):.2f} " f"(chi² on 0 dof — same params; just compare LLs)") print(f" 2·(LL[M2] − LL[M0]) = {2*(lls['M2_exp_order'] - lls['M0_constant']):.2f}") best_model = max(results, key=lambda m: results[m][3]) ll_grid_b, a_b, s0_b, ll_b, dll_b = results[best_model] print(f"\n best: {best_model} (LL = {ll_b:.2f}, ΔLL vs pure Gauss = {dll_b:+.2f})") # ---------- show per-(p,q) widths and weights at best fit ---------- sigmas_best = widths_for_model(alphabet, s0_b, best_model) print(f"\n--- per-symbol widths at best fit (model={best_model}, σ₀={s0_b:.4f}) ---") for (p, q, c, order), s in zip(alphabet, sigmas_best): sym = f"{p}:{q}" print(f" {sym:>8} center={c:.4f} order={order} σ={s:.4f}") # ---------- profile likelihood for α at best model ---------- a_profile = np.linspace(0.0, 1.0, 101) centers = np.array([c for _, _, c, _ in alphabet]) weights = np.ones(len(alphabet)) / len(alphabet) f_at = farey_density_per_sigma(fracs, centers, weights, sigmas_best) Z_at = farey_normalizer_per_sigma(centers, weights, sigmas_best) fn = f_at / Z_at lls_a = np.array([np.log(a * base_vals + (1 - a) * fn).sum() for a in a_profile]) in_ci = a_profile[lls_a > ll_b - 0.5] in95 = a_profile[lls_a > ll_b - 1.92] a_ci = (in_ci.min(), in_ci.max()) if len(in_ci) else (a_b, a_b) a_95 = (in95.min(), in95.max()) if len(in95) else (a_b, a_b) print(f"\n68%-CI on α at {best_model}: [{a_ci[0]:.3f}, {a_ci[1]:.3f}]") print(f"95%-CI on α: [{a_95[0]:.3f}, {a_95[1]:.3f}]") # ---------- LITHWICK–WU asymmetric pile-up: fit ε ---------- # For each first-order resonance (2:1, 3:2, 4:3, 5:4), we fit a paired # bump model: # f_pq(x) = a_inside · N(x; c_pq − ε, σ_lw) + a_outside · N(x; c_pq + ε, σ_lw) # holding total amplitude fixed and fitting the inside-vs-outside ratio γ. print("\n--- Lithwick–Wu asymmetric pile-up at first-order resonances ---") # For each first-order MMR, slice into three signed windows: # "below" = r in (p/q - W, p/q - w) e.g. (1.90, 1.98) for 2:1 # "near" = r in (p/q - w, p/q + w) # "above" = r in (p/q + w, p/q + W) first_order = [(p, q, (p % q) / q if q > 1 else 0.0) for p, q in [(2, 1), (3, 2), (4, 3), (5, 4)]] for WIN in (0.05, 0.10, 0.15): print(f"\n Window = ±{WIN:.2f} (split at ±0.02 around center)") print(f" {'resonance':>10} {'center':>7} {'below':>6} {'near':>5} {'above':>6} " f"{'below/above':>11} {'z(below−above)':>14}") for p, q, c in first_order: # circular distance from center (in (-0.5, 0.5]) d = ((fracs - c + 0.5) % 1) - 0.5 in_below = (d < -0.02) & (d > -WIN) in_near = np.abs(d) <= 0.02 in_above = (d > 0.02) & (d < WIN) n_below = int(in_below.sum()) n_near = int(in_near.sum()) n_above = int(in_above.sum()) ratio_ba = n_below / n_above if n_above > 0 else float("nan") n_pair = n_below + n_above # Binomial null: equal prob on each side z = (n_below - n_pair * 0.5) / math.sqrt(max(n_pair * 0.25, 1)) sym = f"{p}:{q}" print(f" {sym:>10} {c:7.4f} {n_below:6d} {n_near:5d} {n_above:6d} " f"{ratio_ba:11.2f} {z:+14.2f}") # ---------- visualize at best model ---------- xs = np.linspace(0.0001, 0.9999, 2000) ys_gauss = gauss_density(xs) ys_farey = farey_density_per_sigma(xs, centers, weights, sigmas_best) / Z_at ys_mix = a_b * ys_gauss + (1 - a_b) * ys_farey fig, ax = plt.subplots(figsize=(9.5, 5)) ax.hist(fracs, bins=80, density=True, color="#7f8c8d", alpha=0.42, edgecolor="white", linewidth=0.5, label=f"data (n={len(fracs)})") ax.plot(xs, ys_gauss, color="#2c3e50", lw=1.5, alpha=0.6, label="pure Gauss") ax.plot(xs, ys_farey, color="#c0392b", lw=1.4, alpha=0.7, label=f"Farey-12 (per-(p,q) σ, {best_model})") ax.plot(xs, ys_mix, color="#16a085", lw=2.5, label=fr"best mixture α = {a_b:.3f}, σ₀ = {s0_b:.4f}") # Tick marks for low-order centers annotated = set() for p, q, c, order in alphabet: if (p, q) in [(2,1),(3,2),(4,3),(5,4),(3,1),(5,2),(5,3),(4,1)]: ax.axvline(c, color="#c0392b", lw=0.8, alpha=0.4) label = f"{p}:{q}" if c not in annotated: ax.text(c + 0.005, ax.get_ylim()[1] * 0.92, label, fontsize=8, color="#c0392b") annotated.add(c) ax.set_xlabel(r"$\{r\}$") ax.set_ylabel("density") ax.set_title(rf"Physical mixture fit ({best_model}): " rf"$\alpha = {a_b:.3f}$, $\sigma_0 = {s0_b:.4f}$, $\Delta$LL = {dll_b:.1f}") ax.legend(loc="upper left", fontsize=9) ax.set_xlim(0, 1) fig.tight_layout() fig.savefig("c1_physical_mixture.png", dpi=140) # 2D contour at best model AA, SS = np.meshgrid(alpha_grid, sigma_0_grid, indexing="ij") dll_arr = ll_grid_b - ll_b fig, ax = plt.subplots(figsize=(8, 5)) levels = [-30, -9, -4.6, -1.92, -0.5, 0] cs = ax.contourf(AA, SS, dll_arr, levels=levels, cmap="viridis") ax.set_yscale("log") ax.set_xlabel(r"$\alpha$") ax.set_ylabel(r"$\sigma_0$") ax.set_title(rf"2-D likelihood contours, $\Delta$LL ({best_model})") ax.plot([a_b], [s0_b], "r*", markersize=14, label=f"MLE α={a_b:.3f}, σ₀={s0_b:.4f}") plt.colorbar(cs, ax=ax, label=r"$\Delta$LL") ax.legend() fig.tight_layout() fig.savefig("c1_physical_2d.png", dpi=140) # Model comparison bar fig, ax = plt.subplots(figsize=(7, 4.5)) names = list(WIDTH_MODELS.keys()) dlls = [results[n][4] for n in names] alphas = [results[n][1] for n in names] sigmas0 = [results[n][2] for n in names] ax2 = ax.twinx() bars = ax.bar(names, dlls, color="#16a085", alpha=0.8) for b, a, s in zip(bars, alphas, sigmas0): ax.annotate(f"α={a:.2f}\nσ₀={s:.3f}", (b.get_x() + b.get_width()/2, b.get_height()), ha="center", va="bottom", fontsize=9) ax.set_ylabel(r"$\Delta$LL vs pure Gauss") ax.set_title("Width-scaling model comparison") ax.set_xticklabels(names, rotation=15, ha="right") fig.tight_layout() fig.savefig("c1_physical_model_comparison.png", dpi=140) print("\nplots: c1_physical_mixture.png c1_physical_2d.png c1_physical_model_comparison.png") if __name__ == "__main__": main()