"""Step 9: upper limits on the fraction of stars hosting a partial Dyson swarm. For a swarm (f, T) and a stratum of stars: N_rec(f, T) = number of parent stars in which the injected swarm is recovered (step 4); k(f, T) = number of real excess stars that (i) are stars in which this swarm would be detectable, E_model >= THR * sigma in a searchable band, and (ii) are not excluded from hosting it: E_obs >= (1 - TOL) * E_model - 2 sigma in EVERY searchable band (one-sided: a star with MORE excess than the model could host the swarm plus something else, so it is counted - this is conservative); phi_95(f, T) = one-sided 95% Clopper-Pearson upper limit on k / N_rec. Four versions of k, from label-free to label-dependent: UL95_all every excess star counted, whatever it has been identified as. This is the headline: it uses no identification at all. UL95_strict unexplained stars PLUS stars whose only identification is an entry in a published infrared-excess catalogue (a swarm host would earn such an entry). UL95_unexplained unexplained stars only: every identification accepted (independent youth evidence, catalogue entries, background sources, companions). UL95_sensitivity k = 0: what the sample could have delivered. The two label-dependent limits are divided by the probability that a genuine on-star excess survives the background vetting (vetting_pass_fraction, about 0.94). Because the number of compatible stars changes steeply with f between 0.03 and 0.1 (the brightest unexplained M dwarfs are equivalent to f = 0.03-0.08), every limit is also computed with a looser tolerance (TOL_ALT = 0.5) and data/limit_ranges.csv gives, for every stratum and f, the range of each limit over the six temperatures and the two tolerances. Strata: A, FGK, M, FGK+M, all, and five M sub-ranges in absolute magnitude. A pooled stratum ('all', 'FGK+M', 'M') gives a limit on the fraction of its stars ONLY IF the swarm rate does not depend on spectral type inside it: its compatible stars and its recoverable stars are dominated by the brighter members. data/settled_fraction.csv therefore gives the maximum number of host systems per stratum (limit x number of parent stars) and their sum. The excess a host shows relative to its Ks-anchored photosphere: the host is dimmed by (1 - f) and the swarm adds a little flux at Ks, so with B = F_bol,obs / F_Ks,obs and c_b = P_b / F_Ks E_model,b = f F_bol,obs (r_b - c_b r_Ks) / ((1 - f) + f B r_Ks) (r_b = swarm response per unit absorbed bolometric flux, swarm.py). Outputs: data/limits.csv, data/limit_ranges.csv, data/settled_fraction.csv, data/completeness_by_subtype.csv, data/excess_incidence.csv, data/accounting.csv, data/unexplained.csv """ import json import os import numpy as np import pandas as pd from scipy.stats import beta from common import RAW, DATA, T_GRID, F_GRID from swarm import swarm_response import importlib ex = importlib.import_module("03_excess") inj = importlib.import_module("04_injection") im = importlib.import_module("07_images") cl = importlib.import_module("08_classify") DERIVED = os.path.join(RAW, "derived") TOL = 0.2 # tolerance on the model excess (bolometric-flux and colour calibration) TOL_ALT = 0.5 # looser tolerance used for the quoted ranges MATCH = 2.0 # factor for the stricter two-sided 'SED-matched' count EXPLAINED = ("young_star", "catalogued_excess_only", "other_known_type", "background_catalogued", "background_image_offset", "background_visual", "comoving_companion") CATALOGUE_ONLY = ("catalogued_excess_only",) # strata: name -> (classes, M_G lower bound, M_G upper bound); the M sub-ranges use the Mamajek # M_G of M0V (8.16), M2V (9.29), M4V (11.21), M5V (12.45), M6.5V (14.40) and L0V (16.6) STRATA = {"A": (["A"], -99, 99), "FGK": (["FGK"], -99, 99), "M": (["M"], -99, 99), "FGK+M": (["FGK", "M"], -99, 99), "all": (["A", "FGK", "M"], -99, 99), "M0-M1.5": (["M"], 8.16, 9.29), "M2-M3.5": (["M"], 9.29, 11.21), "M4-M4.5": (["M"], 11.21, 12.45), "M5-M6": (["M"], 12.45, 14.40), "M6.5-M9": (["M"], 14.40, 99)} CLASSES = {k: v[0] for k, v in STRATA.items() if v[1] == -99} # the unsplit strata (used by step 10) M_SUB = [k for k, v in STRATA.items() if v[1] != -99] def stratum_mask(df, name): members, lo, hi = STRATA[name] return (df.cls.isin(members) & (df.M_G >= lo) & (df.M_G < hi)).values def cp_upper(k, n, cl=0.95): """One-sided Clopper-Pearson upper limit on a binomial fraction.""" if n <= 0: return 1.0 if k >= n: return 1.0 return float(beta.ppf(cl, k + 1, n - k)) def vetting_pass_fraction(cand=None): """Probability that a genuine on-star excess (e.g. a swarm) survives the background vetting, and its uncertainty. p_img : probability that the image-offset rule rejects an on-star source. A flag in EITHER band rejects the star, so from the null stars p = 1 - (1 - p_W3)(1 - p_W4), where p_b = (1 - NULL_PCT/100) + the fraction of null stars with a lost centroid. The rule also flags a few per cent of the stars with a youth or published-disc identification (whose excess is presumably on the star); the larger of the two rates is used. p_cat : rate at which the catalogue background criteria fire on matched non-excess stars. """ out = {} p_img, p_err = 0.0, 0.0 try: cal = json.load(open(os.path.join(DATA, "image_centroid_calib.json"))) pb = [1 - im.NULL_PCT / 100.0 + 1 - cal[b]["n_null_usable"] / cal[b]["n_null"] for b in ("W3", "W4")] p_img = 1 - (1 - pb[0]) * (1 - pb[1]) p_err = float(np.sqrt(sum(p * (1 - p) / cal[b]["n_null"] for p, b in zip(pb, ("W3", "W4"))))) except (OSError, KeyError): pass if cand is not None and "flag_offset" in cand: nat = cand[cand.explanation.isin(["young_star", "catalogued_excess_only"])] if len(nat) >= 20: p_emp = float(nat.flag_offset.fillna(False).astype(bool).mean()) if p_emp > p_img: p_img, p_err = p_emp, float(np.sqrt(p_emp * (1 - p_emp) / len(nat))) try: br = pd.read_csv(os.path.join(DATA, "background_base_rates.csv")) br = br[br.criterion == "bg_any"].set_index("cls").rate_comparison except OSError: br = pd.Series(dtype=float) for name, (members, _, _) in STRATA.items(): p_cat = max([float(br.get(m, 0.0)) for m in members]) out[name] = ((1 - p_img) * (1 - p_cat), p_err * (1 - p_cat)) return out def model_excess(c, f, resp): """Model excess (Jy) in W3, W4 for hypothetical hosts with the observed properties of c.""" B = c.fbol.values / c.FKs.values den = (1.0 - f) + f * B * resp["Ks"] out = {} for b in ("3", "4"): cb = c[f"P{b}K"].values / c.FKs.values out[b] = f * c.fbol.values * (resp["W" + b] - cb * resp["Ks"]) / den return out def compatible(c, f, resp, thr=ex.THR, tol=TOL): m = model_excess(c, f, resp) detectable = np.zeros(len(c), bool) not_excl = np.ones(len(c), bool) match = np.ones(len(c), bool) for b in ("3", "4"): q = c["q" + b].values E, s = c[f"E{b}K"].values, c[f"sig{b}K"].values detectable |= q & (m[b] >= thr * s) not_excl &= ~q | (E >= (1 - tol) * m[b] - 2 * s) match &= ~q | ((E >= m[b] / MATCH - 2 * s) & (E <= m[b] * MATCH + 2 * s)) return detectable & not_excl, detectable & match def over_luminous_clean(df): """Youth indicator: more than 1 mag above the main-sequence ridge AND Gaia photometry that is not blended (a blended pair has a corrupted G-RP colour and only mimics over-luminosity).""" return (df.dM.values < -1.0) & cl.phot_clean(df) def accounting(d, cand): """One row per class: what the excess stars are, and how many chance alignments were expected.""" ctrl = pd.read_csv(os.path.join(DATA, "control_expectation.csv")) ctrl = ctrl[ctrl.band == "W3 or W4"].set_index("cls") try: direct = pd.read_csv(os.path.join(DATA, "control_direct_count.csv")) direct = direct[direct.band == "W3 or W4"] except OSError: direct = None try: br = pd.read_csv(os.path.join(DATA, "background_base_rates.csv")) br = br[br.criterion == "ev_gaia_neighbour_4as"].set_index("cls") except OSError: br = None cand = cand.drop(columns=[x for x in ("BPmag", "RPmag") if x in cand]).merge(d[["sid", "BPmag", "RPmag"]], on="sid") rows = [] for name in ("A", "FGK", "M", "all"): c = cand if name == "all" else cand[cand.cls == name] dd = d if name == "all" else d[d.cls == name] r = dict(cls=name, n_parent=len(dd), n_searchable=int((dd.q3 | dd.q4).sum()), n_excess=len(c), n_excess_W3=int(c.det3.sum()), n_excess_W4=int(c.det4.sum())) for e in EXPLAINED + ("unexplained",): r["n_" + e] = int((c.explanation == e).sum()) natural = c.explanation.isin(cl.NATURAL) r["n_not_young_or_known"] = int((~natural).sum()) r["expected_chance_alignments"] = float(ctrl.loc[name, "expected_false_positives"]) r["expected_chance_alignments_err"] = float(ctrl.loc[name, "expected_fp_err"]) if direct is not None: for cap in sorted(direct.cap_arcsec.unique()): r[f"direct_count_cap{cap:g}"] = int(direct[(direct.cls == name) & (direct.cap_arcsec == cap)].n_direct.iloc[0]) # efficiency of the image-offset test, measured on catalogued background sources bgc = c[c.explanation == "background_catalogued"] n_flag = int(c[~natural].flag_offset.sum()) r["n_offset_flagged_not_young"] = n_flag r["n_offset_flagged_young_or_catalogued"] = int(c[natural].flag_offset.sum()) if len(bgc) >= 10: eps = bgc.flag_offset.mean() r["offset_test_efficiency"] = eps r["offset_test_efficiency_err"] = np.sqrt(eps * (1 - eps) / len(bgc)) r["n_background_implied_by_offsets"] = n_flag / eps r["n_background_implied_err"] = (n_flag / eps) * np.sqrt(1 / max(n_flag, 1) + (1 - eps) / (eps * len(bgc))) # youth indicator: over-luminous with clean Gaia photometry (see over_luminous_clean) un = c[c.explanation == "unexplained"] yk = c[c.explanation == "young_star"] srch = dd[(dd.q3 | dd.q4)] r["unexplained_over_1mag_above_MS_any_photometry"] = int((un.dM < -1.0).sum()) r["unexplained_overluminous_clean"] = int(over_luminous_clean(un).sum()) r["unexplained_clean_photometry"] = int(cl.phot_clean(un).sum()) r["young_star_overluminous_clean"] = int(over_luminous_clean(yk).sum()) r["young_star_clean_photometry"] = int(cl.phot_clean(yk).sum()) cs = cl.phot_clean(srch) r["frac_searchable_clean_overluminous"] = float(over_luminous_clean(srch).sum() / max(cs.sum(), 1)) r["frac_searchable_over_1mag_above_MS_any_photometry"] = float((srch.dM < -1.0).mean()) # secondary indicators on the unexplained stars r["unexplained_in_sgro_song_2021"] = int(un.in_sgro_song_2021.fillna(False).astype(bool).sum()) r["excess_in_sgro_song_2021"] = int(c.in_sgro_song_2021.fillna(False).astype(bool).sum()) r["unexplained_gaia_neighbour_4as"] = int(un.ev_gaia_neighbour_4as.fillna(False).astype(bool).sum()) if br is not None and name in br.index: r["comparison_rate_gaia_neighbour_4as"] = float(br.loc[name, "rate_comparison"]) r["unexplained_gaia_neighbour_in_excess_of_base_rate"] = float( br.loc[name, "n_unexplained_in_excess_of_base_rate"]) rows.append(r) acc = pd.DataFrame(rows) acc.to_csv(os.path.join(DATA, "accounting.csv"), index=False, float_format="%.4g") print(acc.T.to_string()) return acc def blackbody_equivalent(d, cand): """Single-temperature blackbody (T, f) that reproduces the W3 and W4 excess of each unexplained star (only where both excesses exceed 2 sigma), and the secondary indicators that were NOT used to classify it -> data/unexplained.csv.""" cols = ["sid", "explanation", "simbad_id", "simbad_otype", "sp_type", "nW3", "nW4", "W1mW2", "JmKs", "GCNSprob", "IPDfmp", "bprp_excess", "gaia_phot_clean", "gaia_neighbour_minsep", "in_sgro_song_2021"] drop = [x for x in cols if x in d.columns and x != "sid"] c = d[d.det].drop(columns=drop).merge(cand[cols], on="sid") c = c[c.explanation == "unexplained"].copy() # 99th percentile of W1-W2 among searchable non-excess M dwarfs (hot-dust / AGN-colour flag) sm = d[(d.q3 | d.q4) & ~d.det & (d.cls == "M")] w12_99 = float((sm.w1mpro - sm.w2mpro).quantile(0.99)) Ts = np.arange(60.0, 1500.1, 5.0) r3 = np.array([swarm_response(T)["W3"] for T in Ts]) r4 = np.array([swarm_response(T)["W4"] for T in Ts]) ratio = r4 / r3 # decreases monotonically with T Tbb, fbb, flags = [], [], [] for _, s in c.iterrows(): if s.E3K > 2 * s.sig3K and s.E4K > 2 * s.sig4K: T = float(np.interp(s.E4K / s.E3K, ratio[::-1], Ts[::-1])) a = float(np.interp(T, Ts, r3)) x = s.E3K / (s.fbol * a) # = f / (1 - f) for a cool swarm Tbb.append(T) fbb.append(x / (1 + x)) else: Tbb.append(np.nan) fbb.append(np.nan) fl = [] if s.GCNSprob < 0.5: fl.append(f"GCNS probability {s.GCNSprob:.2f} (astrometry may be spurious)") if s.cls == "M" and (s.JmKs > 1.5 or s.W1mW2 > 1.0): fl.append(f"J-Ks = {s.JmKs:.1f} and W1-W2 = {s.W1mW2:.1f} are not the colours of a normal M dwarf") elif s.cls == "M" and s.W1mW2 > w12_99: fl.append(f"W1-W2 = {s.W1mW2:.2f} is above the 99th percentile of non-excess M dwarfs ({w12_99:.2f}): " "dust hotter than about 300 K or an AGN-coloured contaminant") if not bool(s.gaia_phot_clean): fl.append(f"Gaia photometry flagged as possibly blended (IPDfmp {s.IPDfmp:.0f}, BP/RP flux excess " f"{s.bprp_excess:.2f}): the G-RP colour and the height above the main sequence are unreliable") elif s.dM < -1.0: fl.append(f"{-s.dM:.1f} mag above the main sequence with clean photometry (youth candidate)") if np.isfinite(s.gaia_neighbour_minsep): fl.append(f"Gaia source at {s.gaia_neighbour_minsep:.1f} arcsec not shown to be co-moving") if bool(s.in_sgro_song_2021): fl.append("attached by SIMBAD to Sgro & Song 2021 (published excess-candidate list)") if isinstance(s.simbad_otype, str) and s.simbad_otype in ("Em*", "pA?", "pA*"): fl.append(f"SIMBAD type {s.simbad_otype}") flags.append("; ".join(fl)) out = pd.DataFrame({"gaia_dr3_source_id": c.sid.values, "simbad_id": c.simbad_id.values, "cls": c.cls.values, "sp_type": c.sp_type.values, "dist_pc": c.dist_pc.round(1).values, "Gmag": c.Gmag.round(2).values, "dM_above_MS": c.dM.round(2).values, "chi3": c.chi3.round(1).values, "chi4": c.chi4.round(1).values, "W3_excess_over_photosphere": (c.E3K / c.P3K).round(2).values, "W4_excess_over_photosphere": (c.E4K / c.P4K).round(2).values, "bb_T_K": np.round(Tbb, 0), "bb_f": [float(f"{v:.2g}") if np.isfinite(v) else np.nan for v in fbb], "offset_norm_W3": c.nW3.round(2).values, "offset_norm_W4": c.nW4.round(2).values, "W1mW2": c.W1mW2.round(2).values, "JmKs": c.JmKs.round(2).values, "GCNSprob": c.GCNSprob.round(2).values, "IPDfmp": c.IPDfmp.values, "bprp_excess": c.bprp_excess.round(2).values, "gaia_phot_clean": c.gaia_phot_clean.values, "gaia_neighbour_arcsec": c.gaia_neighbour_minsep.round(1).values, "in_sgro_song_2021": c.in_sgro_song_2021.values, "secondary_indicators": flags}) out = out.sort_values(["cls", "chi3"], ascending=[True, False]) out.to_csv(os.path.join(DATA, "unexplained.csv"), index=False) return out def main(thr=ex.THR, write=True, cand=None, stars=None, rec=None): d = inj.load_stars() if stars is None else stars if rec is None: z = np.load(os.path.join(DERIVED, "injection.npz")) rec, f_fine = z["rec"], z["f_fine"] else: rec, f_fine = rec if cand is None: cand = pd.read_csv(os.path.join(DATA, "candidates.csv")) c = d[d.det].merge(cand[["sid", "explanation"]], on="sid", how="left") c["explanation"] = c.explanation.fillna("unexplained") unexpl = ~c.explanation.isin(EXPLAINED).values strict = unexpl | c.explanation.isin(CATALOGUE_ONLY).values vpass = vetting_pass_fraction(cand) dmask = {name: stratum_mask(d, name) for name in STRATA} cmask = {name: stratum_mask(c, name) for name in STRATA} rows = [] for it, T in enumerate(T_GRID): resp = swarm_response(T) for f in F_GRID: jf = int(np.argmin(np.abs(f_fine - f))) ok, okm = compatible(c, f, resp, thr) ok2, _ = compatible(c, f, resp, thr, TOL_ALT) for name in STRATA: sel, csel = dmask[name], cmask[name] vp = vpass[name][0] n_par, n_rec = int(sel.sum()), int(rec[sel, it, jf].sum()) r = dict(cls=name, T_K=T, f=f, N_parent=n_par, N_rec=n_rec, completeness=n_rec / n_par, k_all=int((ok & csel).sum()), k_strict=int((ok & csel & strict).sum()), k_unexplained=int((ok & csel & unexpl).sum()), k_all_sedmatch=int((okm & csel).sum()), k_unexplained_sedmatch=int((okm & csel & unexpl).sum()), k_all_tol50=int((ok2 & csel).sum()), k_strict_tol50=int((ok2 & csel & strict).sum()), k_unexplained_tol50=int((ok2 & csel & unexpl).sum())) r["UL95_sensitivity"] = cp_upper(0, n_rec) r["UL95_all"] = cp_upper(r["k_all"], n_rec) r["UL95_strict"] = min(1.0, cp_upper(r["k_strict"], n_rec) / vp) r["UL95_unexplained"] = min(1.0, cp_upper(r["k_unexplained"], n_rec) / vp) r["UL95_all_tol50"] = cp_upper(r["k_all_tol50"], n_rec) r["UL95_strict_tol50"] = min(1.0, cp_upper(r["k_strict_tol50"], n_rec) / vp) r["UL95_unexplained_tol50"] = min(1.0, cp_upper(r["k_unexplained_tol50"], n_rec) / vp) r["UL95_unexplained_sedmatch"] = min(1.0, cp_upper(r["k_unexplained_sedmatch"], n_rec) / vp) r["UL95_all_sedmatch"] = cp_upper(r["k_all_sedmatch"], n_rec) r["vetting_pass_fraction"] = vp r["vetting_pass_fraction_err"] = vpass[name][1] r["limited_by"] = ("no sensitivity" if n_rec < 3 or n_rec / n_par < 0.01 else "excess sources" if r["UL95_all"] > 2 * r["UL95_sensitivity"] else "sensitivity") rows.append(r) lim = pd.DataFrame(rows) if not write: return lim lim.to_csv(os.path.join(DATA, "limits.csv"), index=False, float_format="%.4g") accounting(d, cand) blackbody_equivalent(d, cand) # --- ranges over temperature (100-600 K) and tolerance (TOL, TOL_ALT) for every stratum and f rrows = [] for (name, f), g in lim.groupby(["cls", "f"], sort=False): r = dict(cls=name, f=f, completeness_min=g.completeness.min(), completeness_max=g.completeness.max()) g3 = g[g.T_K == 300] for v in ("all", "strict", "unexplained"): both = np.r_[g[f"UL95_{v}"].values, g[f"UL95_{v}_tol50"].values] r[f"UL95_{v}_300K"] = float(g3[f"UL95_{v}"].iloc[0]) r[f"UL95_{v}_300K_tol50"] = float(g3[f"UL95_{v}_tol50"].iloc[0]) r[f"UL95_{v}_min"], r[f"UL95_{v}_max"] = float(both.min()), float(both.max()) r[f"k_{v}_min"] = int(min(g[f"k_{v}"].min(), g[f"k_{v}_tol50"].min())) r[f"k_{v}_max"] = int(max(g[f"k_{v}"].max(), g[f"k_{v}_tol50"].max())) rrows.append(r) rng = pd.DataFrame(rrows) rng.to_csv(os.path.join(DATA, "limit_ranges.csv"), index=False, float_format="%.4g") # --- completeness and limits by M sub-range (and the unsplit classes for reference) cs = lim[lim.cls.isin(["FGK", "M"] + M_SUB) & lim.f.isin([1e-3, 1e-2, 0.1])][ ["cls", "T_K", "f", "N_parent", "N_rec", "completeness", "k_all", "k_strict", "k_unexplained", "UL95_sensitivity", "UL95_all", "UL95_strict", "UL95_unexplained", "limited_by"]] cs.to_csv(os.path.join(DATA, "completeness_by_subtype.csv"), index=False, float_format="%.4g") # --- fraction of stars that could be 'settled' if every settled system built a swarm of # covering fraction 0.1%, 1% or 10% at temperature T, and the corresponding maximum # NUMBER of host systems per stratum (limit x number of parent stars in the stratum; # this extrapolates to the stars that are not searchable). srows = [] for name in STRATA: for T in T_GRID: for f in (0.001, 0.01, 0.1): r = lim[(lim.cls == name) & (lim.T_K == T) & np.isclose(lim.f, f)].iloc[0] row = dict(stratum=name, T_K=T, f=f, N_parent=r.N_parent, N_rec=r.N_rec, completeness=r.completeness, k_all=r.k_all, k_strict=r.k_strict, k_unexplained=r.k_unexplained, sensitivity_floor=r.UL95_sensitivity) for v in ("all", "strict", "unexplained"): row[f"max_settled_fraction_{v}"] = r[f"UL95_{v}"] row[f"max_systems_{v}"] = r[f"UL95_{v}"] * r.N_parent row[f"max_systems_{v}_tol50"] = r[f"UL95_{v}_tol50"] * r.N_parent srows.append(row) st = pd.DataFrame(srows) # sums of per-stratum maxima: conservative totals that do not assume a type-independent rate sums = [] for label, parts in (("sum A+FGK+M", ["A", "FGK", "M"]), ("sum A+FGK+M sub-ranges", ["A", "FGK"] + M_SUB)): for (T, f), g in st[st.stratum.isin(parts)].groupby(["T_K", "f"]): row = dict(stratum=label, T_K=T, f=f, N_parent=g.N_parent.sum(), N_rec=g.N_rec.sum(), completeness=g.N_rec.sum() / g.N_parent.sum(), k_all=g.k_all.sum(), k_strict=g.k_strict.sum(), k_unexplained=g.k_unexplained.sum()) for v in ("all", "strict", "unexplained"): row[f"max_systems_{v}"] = g[f"max_systems_{v}"].sum() row[f"max_systems_{v}_tol50"] = g[f"max_systems_{v}_tol50"].sum() row[f"max_settled_fraction_{v}"] = row[f"max_systems_{v}"] / row["N_parent"] sums.append(row) st = pd.concat([st, pd.DataFrame(sums)], ignore_index=True) st.to_csv(os.path.join(DATA, "settled_fraction.csv"), index=False, float_format="%.4g") # --- observed incidence of excesses (for comparison with debris-disc surveys) irows = [] for name, members in CLASSES.items(): sel = d.cls.isin(members) for b in ("3", "4"): q = sel & d["q" + b] for rmin in (0.1, 0.15, 0.25, 0.5, 1.0): # stars in which an excess of rmin x photosphere would be a >= THR detection sens = q & (rmin * d[f"P{b}K"] >= thr * d[f"sig{b}K"]) & (rmin * d[f"P{b}W"] >= thr * d[f"sig{b}W"]) hit = sens & d["det" + b] & (d[f"E{b}K"] >= rmin * d[f"P{b}K"]) n, k = int(sens.sum()), int(hit.sum()) hit_c = c[c.sid.isin(d.sid[hit])] k_un = int((~hit_c.explanation.isin(EXPLAINED)).sum()) k_yng = int(hit_c.explanation.isin(["young_star", "catalogued_excess_only"]).sum()) irows.append(dict(cls=name, band="W" + b, min_excess_over_photosphere=rmin, n_sensitive=n, n_excess=k, n_young_or_catalogued_excess=k_yng, n_unexplained=k_un, incidence=k / n if n else np.nan, lo68=float(beta.ppf(0.16, k, n - k + 1)) if k > 0 else 0.0, hi68=float(beta.ppf(0.84, k + 1, n - k)) if n > k else 1.0)) inc = pd.DataFrame(irows) inc.to_csv(os.path.join(DATA, "excess_incidence.csv"), index=False, float_format="%.4g") pd.set_option("display.width", 250) pd.set_option("display.max_columns", 40) print("vetting pass fraction:", {k: (round(v[0], 4), round(v[1], 4)) for k, v in vpass.items()}) for name in ("FGK", "M", "all"): print(f"--- {name}: 95% upper limits at T = 300 K") print(lim[(lim.cls == name) & (lim.T_K == 300)][["f", "N_rec", "completeness", "k_all", "k_strict", "k_unexplained", "k_all_tol50", "k_unexplained_tol50", "UL95_sensitivity", "UL95_all", "UL95_strict", "UL95_unexplained", "limited_by"]].to_string(index=False)) print(rng[rng.cls.isin(["FGK", "M", "all"]) & rng.f.isin([1e-3, 1e-2, 0.02, 0.03, 0.05, 0.07, 0.1, 0.5])][ ["cls", "f", "UL95_all_300K", "UL95_all_min", "UL95_all_max", "UL95_strict_min", "UL95_strict_max", "UL95_unexplained_300K", "UL95_unexplained_min", "UL95_unexplained_max"]].to_string(index=False)) print(cs[(cs.T_K == 300) & np.isclose(cs.f, 0.01)].to_string(index=False)) print(st[st.T_K.isin([300]) & np.isclose(st.f, 0.01)][["stratum", "N_parent", "completeness", "k_all", "k_strict", "k_unexplained", "max_systems_all", "max_systems_strict", "max_systems_unexplained"]].to_string(index=False)) return lim if __name__ == "__main__": main()