"""Step 14: candidate bookkeeping and 95% upper limits on the fraction of disc galaxies with f >= 0.5/0.75/0.9. Candidate (classical Annis criterion): least-underluminous-band residual D < log10(1-f). Vetting stages: 1. catalogue-level flags, applied to EVERY sample galaxy (so they shrink the searchable sample): same-redshift SDSS DR12 spectroscopic neighbour within 2', second a.100 detection within 4', HI centroid > 40" from the optical galaxy, line-width error > 20 km/s, same-velocity SIMBAD galaxy within 5'; 2. object-level automated checks on candidates: (a) SDSS DR16 cone: spectroscopic neighbour in the beam or counterpart redshift mismatch; (b) ALFALFA-only galaxies: independent aperture photometry showing that the catalogue flux is too faint, accepted ONLY when the photometered segment is centred on the target (ap_valid, see s09); 3. visual classification (f >= 0.75 candidates only), with blind controls. PRIMARY RESULT (revision): stages 1-2 only. Stage 3 did not separate candidates from controls (it rejected about one third of both) and is not reproducible by code, so it is reported as a variant. N_eff = number of sample galaxies in which dimming by f WOULD have produced a surviving candidate = sum over catalogue-clean galaxies with D < 0 of w_sel (still inside the catalogue's magnitude AND surface-brightness range after dimming; s07) x eps2 (fraction of clean controls passing the DR16 cone checks) x p_ap (ALFALFA-only: probability that the aperture check would not wrongly excuse it, measured on control galaxies of the same apparent magnitude as the dimmed galaxy) [x eps3 (visual variant only: fraction of clean controls not visually rejected)]. Limits: classical one-sided Poisson, UL = mu_95(N_unresolved) / N_eff. No background is subtracted. Waste-heat temperature: the statistic uses the W1 (3.4 micron) band for part of the sample. Waste heat at temperature T adds nu L_nu = f L_bol B(3.4um, T) to that band and fills the deficit back in, so N_eff is recomputed as a function of T (limits_vs_temperature.csv). Outputs: data/limits.csv, data/limit_variants.csv, data/limits_vs_temperature.csv, data/candidates_final.csv, data/vetting_funnel.csv, data/efficiencies.csv, data/aperture_efficiency.csv """ import os import numpy as np import pandas as pd from scipy.stats import beta from common import DATA, poisson_ul, bb_frac, LAM_W1, BC_I, clopper_pearson FS = (0.5, 0.75, 0.9) BTF_CUT = -0.35 E_W_MAX = 20.0 ARTEFACT_CLASSES = frozenset({"C", "M", "S", "X"}) # visual classes treated as explained (see report) T_GRID = (50, 100, 150, 200, 250, 300, 350, 400, 500, 600) AP_MAG_BINS = (0.0, 15.5, 16.5, 99.0) # apparent i magnitude bins for the aperture-check efficiency AP_MIN_CONTROLS = 15 def load(): u = pd.read_csv(os.path.join(DATA, "union_sample.csv")) v = pd.read_csv(os.path.join(DATA, "vet_table.csv")) # hand-entered blind classifications (made from the contact sheets of s13); keyed by galaxy uid vis = pd.read_csv(os.path.join(DATA, "visual_classification.csv")) v = v.merge(vis[["uid", "code", "vis_class", "note"]], on="uid", how="left") v["oc_z_mismatch"] = v.oc_z_mismatch.fillna(False).astype(bool) v["n_spec_same"] = v.n_spec_same.fillna(0) for col in ("dmu", "D_rm", "in_a100", "in_alfalfa_gold"): v[col] = v.uid.map(u.set_index("uid")[col]) return u, v def catalogue_clean(t, simbad=True): ok = (t.n_spec_nb == 0) & (t.n_hi_nb == 0) & ~(t.hi_opt_off > 40.0) & (t.e_W <= E_W_MAX) return ok & (t.n_simbad_nb == 0) if simbad else ok def stage2_reason(r, thr, ap, off=0.0): """ap: 'centred' (baseline), 'any' (first-pass behaviour, no centring requirement) or 'none'. off: difference between the residual statistic in use and the baseline D (non-zero for the running-median variant).""" why = [] if r.n_spec_same > 0: why.append("DR16 spectroscopic neighbour in beam") if r.oc_z_mismatch: why.append("counterpart redshift mismatch") d_ap = r.D_ap if ap == "centred" else (r.D_ap_raw if ap == "any" else np.nan) if r.src == "alfalfa" and np.isfinite(d_ap) and d_ap + off > thr: why.append("catalogue flux too faint (aperture photometry)") return "; ".join(why) def aperture_controls(v): """ALFALFA-only non-candidates with an aperture measurement (300 main controls + 200 extra).""" k = v[(v.src == "alfalfa") & v.role.isin(["control", "control_ap"]) & v.phot_ok.fillna(False).astype(bool)] return k[["uid", "m_i", "dmag_med", "ap_valid", "D_ap_raw"]].copy() def aperture_survival(u, k, f, mode): """P(a truly dimmed ALFALFA-only galaxy is NOT excused by the aperture check), per union galaxy. A dimmed galaxy has residual r_i0 + log10(1-f) and is excused if r_i0 + 0.4 dmag > 0, where dmag is the catalogue-minus-aperture magnitude error. The distribution of dmag (and of 'segment centred') is taken from control galaxies whose catalogue magnitude matches that of the dimmed galaxy.""" p = np.ones(len(u)) if mode == "none": return p al = (u.src == "alfalfa").values mdim = (u.m_i - 2.5 * np.log10(1 - f)).values r0 = u.r_i.values kb = np.digitize(k.m_i.values, AP_MAG_BINS) gb = np.digitize(mdim, AP_MAG_BINS) usable = k.ap_valid.values if mode == "centred" else np.isfinite(k.dmag_med.values) excuse_thr = np.where(usable, 0.4 * k.dmag_med.values, -np.inf) # excused if excuse_thr > -r_i0 for b in np.unique(gb[al]): kk = kb == b if kk.sum() < AP_MIN_CONTROLS: kk = np.ones(len(k), bool) sel = al & (gb == b) & np.isfinite(r0) p[sel] = 1.0 - (excuse_thr[kk][None, :] > -r0[sel][:, None]).mean(axis=1) return p def w1_fill_shift(u, f, T): """Upward shift (dex) of the W1 residual caused by waste heat of temperature T from a fraction f of the stellar bolometric luminosity: log10[1 + f/(1-f) B(3.4um,T) L_bol / nuLnu(W1)], with L_bol = BC_I x nuLnu(i) and each galaxy's own W1/i ratio (catalogue median where a band is missing).""" ratio = u.w1_over_i.copy() for src in ("cf4", "alfalfa"): sel = u.src == src ratio[sel & ratio.isna()] = ratio[sel].median() return np.log10(1.0 + f / (1.0 - f) * bb_frac(LAM_W1, T) * BC_I / ratio.values) def detectable(u, f, T=None, dcol="D"): """Would dimming by f push the galaxy below the threshold? (D0 < 0, with W1 fill-in if T is given)""" if T is None: return (u[dcol] < 0).values rw = u.r_W1.values + w1_fill_shift(u, f, T) return np.fmax(u.r_i.values, rw) < 0 def evaluate(u, v, f, simbad=True, visual=False, artefacts=ARTEFACT_CLASSES, selection="w_sel", ap="centred", eps3_lower=False, dcol="D", T=None): """Run the vetting chain for one harvesting fraction; returns (candidate table, summary dict). selection: 'w_sel' (magnitude and surface brightness, 98th percentile; baseline), 'w_sel_p95', 'w_sel_p99.5', 'mag' (magnitude range only), 'none'.""" thr = np.log10(1 - f) use_visual = visual and f >= 0.75 # only the D < log10(0.25) candidates were inspected cu, cv = catalogue_clean(u, simbad), catalogue_clean(v, simbad) ctrl = v[(v.role == "control") & cv] s2_bad = (ctrl.n_spec_same > 0) | ctrl.oc_z_mismatch eps2 = 1 - s2_bad.mean() kv = ctrl[ctrl.vis_class.notna() & ~s2_bad] n_rej = int(kv.vis_class.isin(artefacts).sum()) eps3 = 1 - n_rej / len(kv) if eps3_lower: # one-sided 95% Clopper-Pearson lower bound on eps3 eps3 = beta.ppf(0.05, len(kv) - n_rej, n_rej + 1) # candidates: every union galaxy below the threshold; vetting columns from the vet table where available ids = u.uid[u[dcol] < thr] c = v[v.uid.isin(ids)].copy() miss = u[u.uid.isin(ids) & ~u.uid.isin(c.uid)].copy() # only possible for dcol != 'D' (never downloaded) if len(miss): miss["not_vetted"] = True c = pd.concat([c, miss], ignore_index=True) c["n_spec_same"] = c.n_spec_same.fillna(0) c["oc_z_mismatch"] = c.oc_z_mismatch.fillna(False).astype(bool) cvc = catalogue_clean(c, simbad) c["stage1"] = np.where(cvc, "", "catalogue flag (confusion / HI offset / width error)") off = (c[dcol] - c.D).fillna(0.0).values c["stage2"] = [stage2_reason(r, thr, ap, o) for r, o in zip(c.itertuples(), off)] c["stage3"] = np.where(c.vis_class.isin(artefacts), "visual: " + c.vis_class.astype(str), "") if use_visual else "" c["status"] = np.where(c.stage1 != "", "explained: " + c.stage1, np.where(c.stage2 != "", "explained: " + c.stage2, np.where(c.stage3 != "", "explained: " + c.stage3, "UNRESOLVED"))) c["f"] = f base = cu.values & detectable(u, f, T, dcol) if selection == "none": w = np.ones(len(u)) elif selection == "mag": w = u[f"inlim_{f}"].values.astype(float) else: w = u[f"{selection[:5]}_{f}{selection[5:]}"].values # w_sel_{f}, w_sel_{f}_p95, w_sel_{f}_p99.5 k = aperture_controls(v) p_ap = aperture_survival(u, k, f, ap) eps = eps2 * (eps3 if use_visual else 1.0) weight = base * w * p_ap * eps # per-galaxy contribution to N_eff return c, dict(f=f, thr=thr, base=base, weight=weight, eps=eps, eps2=eps2, eps3=eps3 if use_visual else 1.0, p_ap_mean=float(p_ap[base & (u.src == "alfalfa").values].mean()), n_ctrl_clean=len(ctrl), n_vis_ctrl=len(kv), n_vis_ctrl_rejected=n_rej, use_visual=use_visual) def row(f, case, n, neff, s, **extra): d = dict(f_min=f, case=case, N_unresolved=int(n), N_eff=round(float(neff), 1), UL95_fraction=poisson_ul(int(n)) / neff if neff > 0 else np.nan, eps_stage2=round(s["eps2"], 4), eps_aperture_alfalfa=round(s["p_ap_mean"], 4), eps_visual=round(s["eps3"], 4)) d.update(extra) return d def count(c, mask=None): un = c[c.status == "UNRESOLVED"] return len(un) if mask is None else int(mask(un).sum()) def main(): u, v = load() assert (catalogue_clean(u) == u.auto_clean).all() acc = pd.read_csv(os.path.join(DATA, "hi_parent_accounting.csv")).set_index("W50_min") gold = u.in_alfalfa_gold.values.astype(bool) isal = (u.src == "alfalfa").values rows, variants, funnel, final, effs, trows = [], [], [], [], [], [] # ---- aperture-check efficiency measured on controls (table for the report) k = aperture_controls(v) ap_rows = [] for lo, hi in zip(AP_MAG_BINS[:-1], AP_MAG_BINS[1:]): kk = k[(k.m_i >= lo) & (k.m_i < hi)] dv = kk.dmag_med[kk.ap_valid] ap_rows.append(dict(i_mag_lo=lo, i_mag_hi=hi, n_controls=len(kk), frac_segment_centred=kk.ap_valid.mean(), median_dmag=dv.median(), robust_sigma_dmag=1.4826 * (dv - dv.median()).abs().median(), frac_centred_and_dmag_gt_0p5=(kk.ap_valid & (kk.dmag_med > 0.5)).mean(), frac_centred_and_dmag_gt_1p0=(kk.ap_valid & (kk.dmag_med > 1.0)).mean())) pd.DataFrame(ap_rows).round(3).to_csv(os.path.join(DATA, "aperture_efficiency.csv"), index=False) print(" aperture-check controls (ALFALFA-only non-candidates) by catalogue i magnitude:") print(pd.DataFrame(ap_rows).round(3).to_string(index=False)) is_gold = lambda t: t.in_alfalfa_gold.astype(bool) # noqa: E731 for f in FS: c, s = evaluate(u, v, f) # PRIMARY: stages 1-2, no visual vetting w = s["weight"] neff = w.sum() unres = c[c.status == "UNRESOLVED"] unres_B = unres[unres.r_bar < BTF_CUT] cvis, svis = evaluate(u, v, f, visual=True) # variant: with visual vetting c["status_with_visual"] = cvis.status.values final.append(c) effs.append(dict(f=f, N_clean_and_D_lt_0=int((u.auto_clean & (u.D < 0)).sum()), mean_selection_weight=float(u[f"w_sel_{f}"][s["base"]].mean()), mean_selection_weight_mag_only=float(u[f"inlim_{f}"][s["base"]].mean()), n_controls_clean=s["n_ctrl_clean"], eps_stage2=s["eps2"], n_aperture_controls=len(k), eps_aperture_alfalfa_only=s["p_ap_mean"], N_eff=neff, N_eff_fraction_of_sample=neff / len(u), N_eff_cf4=w[~isal].sum(), N_eff_alfalfa_only=w[isal].sum(), N_eff_alfalfa_gold=w[gold].sum(), N_eff_V_ge_120=w[(u.V >= 120).values].sum(), n_visual_controls_clean=svis["n_vis_ctrl"], n_visual_controls_rejected=svis["n_vis_ctrl_rejected"], eps_visual=svis["eps3"], N_eff_with_visual=svis["weight"].sum())) s1 = c.stage1 == "" s2 = s1 & (c.stage2 == "") vis_rej = s2 & c.vis_class.isin(ARTEFACT_CLASSES) funnel.append(dict(f=f, threshold_dex=s["thr"], N_sample=len(u), N_catalogue_clean=int(u.auto_clean.sum()), N_raw=len(c), N_after_stage1=int(s1.sum()), N_after_stage2=int(s2.sum()), N_unresolved=len(unres), N_unresolved_baryon_deficient=len(unres_B), N_unresolved_without_a100_detection=int((unres.in_a100 == 0).sum()), N_after_visual=count(cvis) if f >= 0.75 else np.nan, N_visually_rejected=int(vis_rej.sum()) if f >= 0.75 else np.nan)) # ---------------- main table rows.append(row(f, "PRIMARY: waste heat cooler than ~250 K, any cause (stages 1-2)", len(unres), neff, s)) rows.append(row(f, "with visual vetting (stage 3)", count(cvis), svis["weight"].sum(), svis)) rows.append(row(f, "baryon-deficient candidates only", len(unres_B), (w * (u[f"rbar_dim_{f}"] < BTF_CUT).values).sum(), s)) for label, su, sc in (("CF4 members", ~isal, lambda t: t.src == "cf4"), ("ALFALFA-only members", isal, lambda t: t.src == "alfalfa"), ("ALFALFA gold subsample (uniform HI selection)", gold, is_gold), ("CF4 outside ALFALFA gold", ~gold, lambda t: ~is_gold(t)), ("V >= 120 km/s", (u.V >= 120).values, lambda t: t.V >= 120), ("V < 120 km/s", (u.V < 120).values, lambda t: t.V < 120), ("V >= 150 km/s", (u.V >= 150).values, lambda t: t.V >= 150)): rows.append(row(f, "stratum: " + label, count(c, sc), w[su].sum(), s)) # Closed-accounting bound for the HI-selected (ALFALFA gold) subsample: no optical selection term in # N_eff; instead add the fraction of the HI-selected parent that has no usable photometry at all # (where a disc dimmed out of the optical catalogue would have to be); see s11 and the report. # Each of the two terms is taken at 97.5% so that their sum is a bound at >= 95% (union bound). c0, s0 = evaluate(u, v, f, selection="none") ng = count(c0, is_gold) ne0 = s0["weight"][gold].sum() a_term = poisson_ul(ng, 0.975) / ne0 for wmin, lab, col in ((120, "all unusable sources counted, W50 >= 120 km/s parent", "N_possible_hidden"), (200, "all unusable sources counted, W50 >= 200 km/s parent", "N_possible_hidden"), (200, "by-eye ordinary sources removed, W50 >= 200 km/s parent", "N_possible_hidden_after_inspection")): b_term = poisson_ul(int(acc.loc[wmin, col]), 0.975) / acc.loc[wmin, "N_parent"] rows.append(dict(f_min=f, case=f"ALFALFA gold, closed-accounting bound ({lab})", N_unresolved=ng, N_eff=round(ne0, 1), UL95_fraction=a_term + b_term, eps_stage2=round(s0["eps2"], 4), eps_aperture_alfalfa=round(s0["p_ap_mean"], 4), eps_visual=1.0, term_TF=a_term, term_hidden=b_term)) # ---------------- limit as a function of waste-heat temperature has_w1 = u.r_W1.notna().values for T in T_GRID: cT, sT = evaluate(u, v, f, T=T) nT = sT["weight"].sum() col = f"f_warm_max_T{T}" # candidates whose W3/W4 fluxes allow a fraction >= f at this temperature (no WISE data counts as allowed) n_warm = int((~(unres[col] < f)).sum()) if col in unres.columns else np.nan trows.append(dict(f_min=f, T_K=T, N_eff=nT, N_eff_retained=nT / neff, frac_of_N_eff_using_W1=w[has_w1].sum() / neff, median_W1_fill_shift_dex=float(np.nanmedian(w1_fill_shift(u, f, T)[s["base"] & has_w1])), N_unresolved=len(unres), UL95_any_cause=poisson_ul(len(unres)) / nT, N_unresolved_allowed_by_W3W4=n_warm, UL95_waste_heat_at_T=poisson_ul(n_warm) / nT if n_warm == n_warm else np.nan)) # ---------------- variants. in_range = defensible analysis choices that define the quoted range def var(case, cs, in_range, sub=None, subc=None): cc, ss = cs ww = ss["weight"] if sub is None else ss["weight"][sub] variants.append(row(f, case, count(cc, subc), ww.sum(), ss, in_range=in_range)) var("baseline (primary)", (c, s), True) var("with visual vetting (classes C, M, S, X excused)", (cvis, svis), True) if f >= 0.75: var("with visual vetting, class C (companion) NOT excused", evaluate(u, v, f, visual=True, artefacts=frozenset({"M", "S", "X"})), True) var("with visual vetting, efficiency at its 95% lower bound", evaluate(u, v, f, visual=True, eps3_lower=True), True) var("selection range at the 95th percentile", evaluate(u, v, f, selection="w_sel_p95"), True) var("selection range at the 99.5th percentile", evaluate(u, v, f, selection="w_sel_p99.5"), True) var("no aperture-photometry check at all", evaluate(u, v, f, ap="none"), True) var("without the SIMBAD 5' neighbour flag", evaluate(u, v, f, simbad=False), True) var("running-median Tully-Fisher relation instead of clipped quadratic", evaluate(u, v, f, dcol="D_rm"), True) var("ALFALFA gold subsample only (uniform HI selection)", (c, s), True, gold, is_gold) var("CF4 outside ALFALFA gold only", (c, s), True, ~gold, lambda t: ~is_gold(t)) var("V >= 120 km/s only", (c, s), False, (u.V >= 120).values, lambda t: t.V >= 120) var("V < 120 km/s only", (c, s), False, (u.V < 120).values, lambda t: t.V < 120) var("waste heat at 300 K (W1 partly filled in)", evaluate(u, v, f, T=300), False) var("waste heat at 350 K (W1 partly filled in)", evaluate(u, v, f, T=350), False) var("OPTIMISTIC, not defensible: magnitude range only (no surface-brightness term; first-pass definition)", evaluate(u, v, f, selection="mag"), False) var("OPTIMISTIC, not defensible: no selection requirement on dimmed galaxies", evaluate(u, v, f, selection="none"), False) var("FIRST-PASS, for reference: aperture check without centring, visual vetting, magnitude range only", evaluate(u, v, f, visual=True, ap="any", selection="mag"), False) variants.append(dict(f_min=f, case="sensitivity floor: if every unresolved candidate is natural (N=0)", N_unresolved=0, N_eff=round(neff, 1), UL95_fraction=poisson_ul(0) / neff, eps_stage2=round(s["eps2"], 4), eps_aperture_alfalfa=round(s["p_ap_mean"], 4), eps_visual=1.0, in_range=False)) # raw, unvetted rate (the like-for-like comparison with Zackrisson et al. 2015: 11 / ~400) n_raw = int((u.D < s["thr"]).sum()) neff_raw = ((u.D < 0).values * u[f"w_sel_{f}"].values).sum() variants.append(dict(f_min=f, case="RAW, no vetting at all (all galaxies; compare Zackrisson+15: 11 / ~400)", N_unresolved=n_raw, N_eff=round(neff_raw, 1), UL95_fraction=poisson_ul(n_raw) / neff_raw, eps_stage2=1.0, eps_aperture_alfalfa=1.0, eps_visual=1.0, in_range=False, raw_fraction=n_raw / neff_raw)) # ---- symmetry check for f >= 0.5: faint tail vs mirrored bright tail (no limit is drawn from it). # D is the LEAST under-luminous band, so the mirror statistic is the LEAST over-luminous band. cl = u[u.auto_clean & u["inlim_0.5"]].copy() cl["Dmin"] = np.fmin(cl.r_i, cl.r_W1) for vmin in (80, 150): sx = cl[cl.V >= vmin] nf, nb = int((sx.D < -0.301).sum()), int((sx.Dmin > 0.301).sum()) print(f" symmetric-tail check V>={vmin}: N={len(sx)}, faint tail {nf}, mirrored bright tail {nb} (no limit drawn)") lim, var_t, tt = pd.DataFrame(rows), pd.DataFrame(variants), pd.DataFrame(trows) lim.to_csv(os.path.join(DATA, "limits.csv"), index=False) var_t.to_csv(os.path.join(DATA, "limit_variants.csv"), index=False) tt.round(5).to_csv(os.path.join(DATA, "limits_vs_temperature.csv"), index=False) pd.DataFrame(funnel).to_csv(os.path.join(DATA, "vetting_funnel.csv"), index=False) pd.DataFrame(effs).round(4).to_csv(os.path.join(DATA, "efficiencies.csv"), index=False) fin = pd.concat(final, ignore_index=True) keep = ["f", "uid", "src", "name", "PGC", "AGC", "ra", "dec", "cz", "W", "e_W", "in_a100", "inc", "V", "Dist", "D", "D_rm", "r_i", "r_W1", "D_other", "m_i", "m_W1", "ap_valid", "D_ap", "D_ap_raw", "dmag_med", "cen_off", "r_bar", "fgas", "logMs", "logMHI", "dmu", "in_alfalfa_gold", "n_spec_nb", "n_hi_nb", "hi_opt_off", "n_simbad_nb", "n_simbad_nb10", "simbad_nb_sep", "simbad_nb_name", "n_spec_same", "sep_spec_same", "self_dv", "nuLnu_W3_over_i", "nuLnu_W4_over_i", "f_warm_max", "f_warm_max_T100", "f_warm_max_T200", "f_warm_max_T300", "code", "vis_class", "note", "stage1", "stage2", "status", "status_with_visual"] fin[[x for x in keep if x in fin.columns]].round(4).to_csv(os.path.join(DATA, "candidates_final.csv"), index=False) # compact list of the unresolved f >= 0.75 candidates (the follow-up list of the report) un = fin[(fin.f == 0.75) & (fin.status == "UNRESOLVED")].copy() un["name"] = np.where(un.src == "alfalfa", "AGC " + un.AGC.astype(int).astype(str), un.name.astype(str)) un["also_f_ge_0p9"] = un.D < -1.0 un["baryon_deficient"] = un.r_bar < BTF_CUT un = un.sort_values("D")[["name", "uid", "ra", "dec", "cz", "W", "e_W", "inc", "V", "Dist", "D", "r_bar", "fgas", "in_a100", "in_alfalfa_gold", "also_f_ge_0p9", "baryon_deficient", "simbad_nb_sep", "f_warm_max", "vis_class", "status_with_visual", "note"]] un.round(3).to_csv(os.path.join(DATA, "unresolved_candidates.csv"), index=False) # ---- printed summary pd.set_option("display.width", 250) pd.set_option("display.max_colwidth", 110) vv = v[v.vis_class.notna()] print(" visual classes (rows: class; columns: role):") print(pd.crosstab(vv.vis_class, vv.role).to_string()) f75 = fin[fin.f == 0.75] s2ok = (f75.stage1 == "") & (f75.stage2 == "") e75 = effs[1] print(f" visual vetting: rejected {int((s2ok & f75.vis_class.isin(ARTEFACT_CLASSES)).sum())} of {int(s2ok.sum())} " f"stage-2 survivors (f >= 0.75) vs {e75['n_visual_controls_rejected']} of {e75['n_visual_controls_clean']} clean controls") lo, hi = clopper_pearson(e75["n_visual_controls_clean"] - e75["n_visual_controls_rejected"], e75["n_visual_controls_clean"]) print(f" visual efficiency {e75['eps_visual']:.2f} (68% interval {lo:.2f}-{hi:.2f})") print(pd.DataFrame(effs).round(3).T.to_string()) print(pd.DataFrame(funnel).T.to_string()) print(lim[["f_min", "case", "N_unresolved", "N_eff", "UL95_fraction"]].round(5).to_string(index=False)) print(var_t[["f_min", "case", "N_unresolved", "N_eff", "UL95_fraction", "in_range"]].round(5).to_string(index=False)) for f in FS: r = var_t[(var_t.f_min == f) & var_t.in_range.astype(bool)] print(f" f >= {f}: range of the 95% limit over defensible variants: " f"{100 * r.UL95_fraction.min():.2f}% - {100 * r.UL95_fraction.max():.2f}%") print(tt[["f_min", "T_K", "N_eff", "N_eff_retained", "median_W1_fill_shift_dex", "UL95_any_cause", "N_unresolved_allowed_by_W3W4", "UL95_waste_heat_at_T"]].round(4).to_string(index=False)) print(" status of the f >= 0.75 candidates (primary):") print(f75.status.value_counts().to_string()) un = f75[f75.status == "UNRESOLVED"] print(un[["uid", "name", "src", "in_a100", "V", "e_W", "D", "r_bar", "vis_class", "status_with_visual"]].round(2).to_string(index=False)) if __name__ == "__main__": main()