#!/usr/bin/env python3 """Collect every number quoted in report.md into data/report_numbers.txt, and write the machine-readable limit curves for the project's model: data/headline_limits.csv limit curve by covering fraction: 300 K, the envelope (weakest limit) over 150-400 K and over the full 100-600 K band, for the primary tolerance (0.2) and the conservative one (0.4), all outliers and clean subsample, with N_eff = 3/UL. EARLY-TYPE GALAXIES ONLY. data/neff_by_temperature.csv the same, resolved by temperature (for a model with a temperature prior) data/limit_crossings.csv covering fraction at which each limit falls to 10%, 1%, 0.1%, 1e-4 data/prior_work_comparison.csv matched-f comparison with published limits """ import json import os import numpy as np import pandas as pd import mirq_common as M # published limits used in the comparison (fraction of galaxies, 95%); see report section 7 for what each means PRIOR = [ ("Curtis+26 (129 galaxies, SED fits, 300 K)", 0.05, 0.146), ("Curtis+26 (129 galaxies, SED fits, 300 K)", 0.10, 0.101), ("Curtis+26 (129 galaxies, SED fits, 300 K) = project model value", 0.25, 0.026), ("Huang+26 (19,617 2MRS galaxies, fixed L*, 300 K, >=21%)", 0.20, 1.61e-4), ("Huang+26 (19,617 2MRS galaxies, fixed L*, 300 K, >=21%)", 0.25, 1.61e-4), ("Huang+26 (19,617 2MRS galaxies, fixed L*, 300 K, >=21%)", 0.30, 1.61e-4), ("Huang+26 (19,617 2MRS galaxies, fixed L*, 300 K, >=21%)", 0.50, 1.61e-4), ("Griffith+15 as used in the project model (N_eff 4.8e3)", 0.50, 3.0 / 4.8e3), ] def crossing(f, ul, level): """Smallest f at which the limit falls to `level` (log-log interpolation along the fine grid).""" f, ul = np.asarray(f, float), np.asarray(ul, float) ok = np.isfinite(ul) & (ul > 0) f, ul = f[ok], ul[ok] for i in range(1, len(f)): if ul[i - 1] > level >= ul[i]: t = (np.log(level) - np.log(ul[i - 1])) / (np.log(ul[i]) - np.log(ul[i - 1])) return float(np.exp(np.log(f[i - 1]) + t * (np.log(f[i]) - np.log(f[i - 1])))) return np.nan def envelope(t, Ts, cols): """Weakest (largest) limit over the temperatures Ts, per f.""" return t[t["T"].isin(Ts)].groupby("f")[cols].max() def main(): out = [] P = out.append J = lambda n: json.load(open(os.path.join(M.DATA, n))) # noqa: E731 fmt = lambda v: "%.3g" % v # noqa: E731 b = J("build_summary.json") P("== SAMPLE / DATA ==") P(json.dumps(b)) s = J("limits_summary.json") P(json.dumps(s)) P("\n== MODEL ==") tp = pd.read_csv(os.path.join(M.DATA, "model_templates.csv")) P("Brown+14 passive ellipticals: n=%d; L_bol/(nu L_nu)(W1) median %.2f (min %.2f, max %.2f, rms %.2f)" % (len(tp), tp.Lbol_over_nuLnu_W1.median(), tp.Lbol_over_nuLnu_W1.min(), tp.Lbol_over_nuLnu_W1.max(), tp.Lbol_over_nuLnu_W1.std())) for c in ("W1-W2", "W1-W3", "W1-W4"): P(" template %s: median %.2f (min %.2f, max %.2f)" % (c, tp[c].median(), tp[c].min(), tp[c].max())) cc = pd.read_csv(os.path.join(M.DATA, "model_colour_corrections.csv")) P("colour-correction validation: max |ratio-1| W3,W4 = %.4f; all bands T>=141K = %.4f" % (np.abs(cc[cc.band.isin(["W3", "W4"])].ratio - 1).max(), np.abs(cc[cc.spectrum != "B(100K)"].ratio - 1).max())) g = np.load(os.path.join(M.DATA, "model_grid.npz")) iz = int(np.argmin(np.abs(g["z"] - 0.08))) for i, T in enumerate(g["T"]): if T in M.T_GRID: P(" R_b(T=%d K, z=0.08): W1 %.3g W2 %.3g W3 %.3g W4 %.3g" % (T, *g["R"][i, :, iz])) P(" template colours C_b = F_b/F_1 at z=0.08: W2 %.3g W3 %.3g W4 %.3g ; bluest template: %.3g %.3g %.3g" % (*g["C"][:, iz], *g["C_min"][:, iz])) P(json.dumps(J("kbol_check.json"))) P("\n== PHOTOMETRY CHECKS ==") P(json.dumps(J("phot_checks.json"), indent=0)) if os.path.exists(os.path.join(M.DATA, "image_injection_summary.json")): P("\n== IMAGE-LEVEL PROFILE-FIT RESPONSE TEST (s09) ==") P(json.dumps(J("image_injection_summary.json"), indent=0)) P("\n== BASELINE ==") P(pd.read_csv(os.path.join(M.DATA, "baseline_by_class.csv")).to_string(index=False, float_format=fmt)) P("\nnatural emission as equivalent covering fraction (bright subsample):") P(pd.read_csv(os.path.join(M.DATA, "feq_percentiles.csv")).to_string(index=False, float_format=fmt)) P("\nper-galaxy caps f_max(T):") P(pd.read_csv(os.path.join(M.DATA, "fmax_percentiles.csv")).to_string(index=False, float_format=fmt)) st = pd.read_csv(os.path.join(M.DATA, "stack_limits.csv")) P("\nstacks (mean g = sum F_b / sum R_b F1; 95% bootstrap upper value), AllWISE:") P(st[st.phot == "psf"].to_string(index=False, float_format=lambda v: "%.4g" % v)) piv = st[st.phot == "psf"].copy() piv["f_ul95"] = piv.mean_g_ul95 / (1 + piv.mean_g_ul95) best = piv.groupby(["pop", "photosphere_subtracted", "T"]).f_ul95.min().unstack("T") P("\nstack: smaller of the W3 and W4 values, as a covering fraction (95% upper value):") P(best.to_string(float_format=lambda v: "%.4f" % v)) bestband = piv.loc[piv.groupby(["pop", "photosphere_subtracted", "T"]).f_ul95.idxmin(), ["pop", "photosphere_subtracted", "T", "band"]] P(bestband[bestband["pop"] == "all"].to_string(index=False)) P("\n== LIMITS (primary, tolerance eps = 0.2) ==") t = pd.read_csv(os.path.join(M.DATA, "limits_primary.csv")) cols = ["T", "f", "n_test", "sens_frac", "recov_eff", "n_cand", "n_null", "frac_compatible", "ul_all", "ul_zero", "n_clean", "n_test_clean", "sens_frac_clean", "n_cand_clean", "n_null_clean", "ul_clean", "ul_zero_clean", "sel_mag", "sel_pess", "agncol_used"] P(t[cols].to_string(index=False, float_format=fmt)) P("recovery efficiency of testable injected hosts (catalogue level, by construction ~1): min %.4f max %.4f" % (t.recov_eff.min(), t.recov_eff.max())) P("\ncandidate make-up on the brief grid (flags are not exclusive):") P(t[t.f.isin(M.F_GRID)][["T", "f", "n_cand", "n_cand_line", "n_cand_agncol", "n_cand_uv", "n_cand_psb", "n_cand_radio", "n_cand_phot"]].to_string(index=False)) t["floor_ratio_all"] = t.ul_all / t.ul_zero t["floor_ratio_clean"] = t.ul_clean / t.ul_zero_clean P("\nlimit / zero-candidate sensitivity floor:") P(t[t.f.isin([0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.5])].pivot(index="T", columns="f", values="floor_ratio_all").to_string(float_format=fmt)) P("clean:") P(t[t.f.isin([0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.5])].pivot(index="T", columns="f", values="floor_ratio_clean").to_string(float_format=fmt)) v = pd.read_csv(os.path.join(M.DATA, "limits_vetted.csv")) P("\n== LIMITS after image vetting of the clean candidates ==") P(v.to_string(index=False, float_format=fmt)) P(json.dumps(J("vetting_summary.json"), indent=0)) var = pd.read_csv(os.path.join(M.DATA, "limits_variants.csv")) e3 = var[var.label == "eps0.3"] e4 = var[var.label == "eps0.4"] ps = var[var.label == "photosphere_subtracted"] P("\n== f at which the 95% limit falls to a given level ==") rows = [] for T in M.T_GRID: s_ = t[t["T"] == T].sort_values("f") s4 = e4[e4["T"] == T].sort_values("f") row = dict(T=T) for lev in (1e-1, 1e-2, 1e-3, 1e-4): row["all_%g" % lev] = crossing(s_.f, s_.ul_all, lev) row["clean_%g" % lev] = crossing(s_.f, s_.ul_clean, lev) row["all_eps0.4_%g" % lev] = crossing(s4.f, s4.ul_all, lev) row["f_natural_50pct"] = crossing(s_.f, s_.frac_compatible, 0.5) row["f_natural_5pct"] = crossing(s_.f, s_.frac_compatible, 0.05) rows.append(row) cr = pd.DataFrame(rows) cr.to_csv(os.path.join(M.DATA, "limit_crossings.csv"), index=False, float_format="%.4g") P(cr.to_string(index=False, float_format=fmt)) # ---------------- temperature-resolved table ---------------- P("\n== TEMPERATURE-RESOLVED LIMITS AND N_eff = 3/UL (early-type galaxies only) ==") k = ["T", "f"] nt = t[k + ["n_test", "sens_frac", "n_cand", "ul_all", "ul_zero", "n_cand_clean", "ul_clean"]].merge( e3[k + ["ul_all", "ul_clean"]].rename(columns={"ul_all": "ul_all_eps0.3", "ul_clean": "ul_clean_eps0.3"}), on=k, how="left").merge( e4[k + ["n_cand", "ul_all", "ul_clean"]].rename(columns={"n_cand": "n_cand_eps0.4", "ul_all": "ul_all_eps0.4", "ul_clean": "ul_clean_eps0.4"}), on=k, how="left").merge( ps[k + ["ul_all", "ul_clean"]].rename(columns={"ul_all": "ul_all_photosphere_subtracted", "ul_clean": "ul_clean_photosphere_subtracted"}), on=k, how="left") for c in ("ul_all", "ul_clean", "ul_all_eps0.4"): nt["neff_" + c[3:]] = 3.0 / nt[c] nt.to_csv(os.path.join(M.DATA, "neff_by_temperature.csv"), index=False, float_format="%.4g") P(nt.to_string(index=False, float_format=fmt)) # ---------------- headline curves ---------------- P("\n== HEADLINE CURVES (early-type galaxies only). 'core' = weakest limit over 150-400 K; 'full' = weakest " "over the project model's warm band 100-600 K ==") env = pd.DataFrame(index=sorted(t.f.unique())) env.index.name = "f" s3 = t[t["T"] == 300.0].set_index("f") env["ul_all_300K"], env["ul_clean_300K"] = s3.ul_all, s3.ul_clean for tag, Ts in (("core150_400K", M.T_BAND_CORE), ("full100_600K", M.T_BAND_FULL)): e = envelope(t, Ts, ["ul_all", "ul_clean", "ul_zero"]) env["ul_all_" + tag], env["ul_clean_" + tag], env["ul_zero_" + tag] = e.ul_all, e.ul_clean, e.ul_zero ee = envelope(e4, Ts, ["ul_all", "ul_clean"]) env["ul_all_%s_eps0.4" % tag], env["ul_clean_%s_eps0.4" % tag] = ee.ul_all, ee.ul_clean ee3 = envelope(e3, Ts, ["ul_all"]) env["ul_all_%s_eps0.3" % tag] = ee3.ul_all env["T_weakest_all_" + tag] = [float(t[(t.f == f) & t["T"].isin(Ts)].sort_values("ul_all")["T"].iloc[-1]) for f in env.index] env["sens_frac_min_full100_600K"] = t.groupby("f").sens_frac.min() env["sel_mag"] = t.groupby("f").sel_mag.first() for c in [c for c in env.columns if c.startswith("ul_") and "zero" not in c]: env["neff_" + c[3:]] = 3.0 / env[c] env["population"] = "SDSS early-type galaxies only (0.02 1: this work is tighter) ==") rows = [] ev = env.set_index("f") for name, f, ul in PRIOR: r = ev.loc[f] rows.append(dict(reference=name, f=f, published_ul=ul, this_300K_all=r.ul_all_300K, this_core_all=r.ul_all_core150_400K, this_full_all=r.ul_all_full100_600K, this_300K_all_eps04=float(e4[(e4["T"] == 300) & (e4.f == f)].ul_all.iloc[0]), this_full_all_eps04=r["ul_all_full100_600K_eps0.4"], ratio_300K=ul / r.ul_all_300K, ratio_core=ul / r.ul_all_core150_400K, ratio_full=ul / r.ul_all_full100_600K, ratio_full_eps04=ul / r["ul_all_full100_600K_eps0.4"], this_300K_clean=r.ul_clean_300K)) pc = pd.DataFrame(rows) pc.to_csv(os.path.join(M.DATA, "prior_work_comparison.csv"), index=False, float_format="%.4g") P(pc.to_string(index=False, float_format=fmt)) # ---------------- low f: photosphere ---------------- P("\n== LOW f: primary vs photosphere-subtracted test (300 K and all T) ==") lo = t[t.f.isin([0.01, 0.02, 0.03, 0.05])][k + ["n_test", "sens_frac", "n_cand", "n_null", "frac_compatible", "ul_all", "ul_clean"]].merge( ps[k + ["n_cand", "frac_compatible", "ul_all", "ul_clean"]], on=k, suffixes=("", "_phsub")) P(lo.to_string(index=False, float_format=fmt)) # ---------------- variants ---------------- P("\n== VARIANTS: range of the all-outlier and clean limits across analysis variants ==") names_param = ["primary", "ks4", "ks7", "kc3", "eps0", "eps0.3", "eps0.4", "kbol_low", "kbol_high", "no_w2", "photosphere_subtracted", "allwise_and_legacy"] names_phot = names_param + ["legacy_w3_only"] for lab, names in (("model/threshold variants on AllWISE photometry", names_param), ("the same plus Legacy forced W3 photometry", names_phot)): sub = var[var.label.isin(names) & var.f.isin(M.F_VAR)] rng = sub.groupby(["T", "f"]).agg(ul_all_min=("ul_all", "min"), ul_all_max=("ul_all", "max"), ul_clean_min=("ul_clean", "min"), ul_clean_max=("ul_clean", "max")).reset_index() P(lab + ":") P(rng[rng.f.isin([0.1, 0.2, 0.3, 0.5])].to_string(index=False, float_format=fmt)) P("\nfull variant table for T = 100, 150, 200, 300 K:") P(var[var["T"].isin([100.0, 150.0, 200.0, 300.0]) & var.f.isin(M.F_VAR)][ ["label", "T", "f", "n", "n_test", "sens_frac", "n_cand", "ul_all", "n_test_clean", "n_cand_clean", "ul_clean"]] .to_string(index=False, float_format=fmt)) txt = "\n".join(out) open(os.path.join(M.DATA, "report_numbers.txt"), "w").write(txt + "\n") print(txt[:3000]) if __name__ == "__main__": main()