"""Step 10: robustness checks -> data/robustness.csv, data/recall_known_discs.csv. Variants (limits quoted for a 300 K swarm with f = 0.001, 0.01, 0.05 and 0.1): * detection threshold 4, 5 (baseline), 6 sigma; * dust cut E(B-V) < 0.10, 0.15 (baseline), 0.30; * a single photospheric anchor (Ks only / W1 only) instead of min(chi_Ks, chi_W1): numbers of positive and (mirror) negative outliers; * tolerance of the 'not excluded' test (TOL = 0, 0.2, 0.5) and the two-sided SED match; * vetting: without the manual image inspection; without any image information; with stars whose only identification is a published infrared-excess catalogue counted as unexplained (column UL95_strict_* of every vetted row, and its own row); * parent sample restricted to GCNS probability > 0.5; * chance-alignment estimate with the W4 trigger area capped at 6 instead of 8 arcsec. For variants that change the detection list (threshold, dust cut) the new detections are not vetted, so only the pure-sensitivity limit and the count-every-excess limit are given. The label-dependent limits (UL95_unexpl_*, UL95_strict_*) are divided by the vetting pass fraction, exactly as in data/limits.csv. External validation: recall of published WISE excess stars (Cotten & Song 2016 'Prime', Patel et al. 2014) that are in our searchable sample. """ import io import os import numpy as np import pandas as pd from common import RAW, DATA, cds_xmatch from swarm import swarm_response import importlib ex = importlib.import_module("03_excess") inj = importlib.import_module("04_injection") ctl = importlib.import_module("05_control") lm = importlib.import_module("09_limits") T_REF = 300 F_REF = (0.001, 0.01, 0.05, 0.1) CLS = {"A": ["A"], "FGK": ["FGK"], "M": ["M"], "all": ["A", "FGK", "M"]} def variant_rows(name, d, cal, thr, q3, q4, cand=None, tol=lm.TOL, vpass=None): """Detections, mirror negatives, recovered numbers and limits for one variant.""" d = d.copy() d["q3"], d["q4"] = q3, q4 d["det3"] = q3 & (d.chi3.values >= thr) d["det4"] = q4 & (d.chi4.values >= thr) d["det"] = d.det3 | d.det4 neg = (q3 & (d.chi3_neg.values <= -thr)) | (q4 & (d.chi4_neg.values <= -thr)) use, s = inj.star_arrays(d) resp = swarm_response(T_REF) c = d[d.det] if cand is not None: c = c.merge(cand[["sid", "explanation"]], on="sid", how="left") unexpl = ~c.explanation.fillna("unexplained").isin(lm.EXPLAINED).values strict = unexpl | c.explanation.isin(lm.CATALOGUE_ONLY).values vpass = vpass or lm.vetting_pass_fraction(cand) rows = [] recs, oks = {}, {} for f in F_REF: r = np.zeros(len(d), bool) r[use] = inj.inject(cal, s, f, resp, thr)[0] recs[f] = r oks[f] = lm.compatible(c, f, resp, thr, tol) for cls, members in CLS.items(): sel = d.cls.isin(members).values csel = c.cls.isin(members).values row = dict(variant=name, cls=cls, n_searchable=int(((q3 | q4) & sel).sum()), n_excess=int((d.det.values & sel).sum()), n_mirror_negative=int((neg & sel).sum())) for f in F_REF: n_rec = int(recs[f][sel].sum()) ok, okm = oks[f] tag = f"f{f:g}" row[f"N_rec_{tag}"] = n_rec row[f"UL95_sens_{tag}"] = lm.cp_upper(0, n_rec) row[f"k_all_{tag}"] = int((ok & csel).sum()) row[f"UL95_all_{tag}"] = lm.cp_upper(row[f"k_all_{tag}"], n_rec) row[f"k_all_sedmatch_{tag}"] = int((okm & csel).sum()) if cand is not None: row[f"k_unexpl_{tag}"] = int((ok & csel & unexpl).sum()) row[f"UL95_unexpl_{tag}"] = min(1.0, lm.cp_upper(row[f"k_unexpl_{tag}"], n_rec) / vpass[cls][0]) row[f"k_strict_{tag}"] = int((ok & csel & strict).sum()) row[f"UL95_strict_{tag}"] = min(1.0, lm.cp_upper(row[f"k_strict_{tag}"], n_rec) / vpass[cls][0]) rows.append(row) return rows def recall_known(d): """Published WISE excess stars inside our searchable A/FGK sample: how many do we flag?""" s = d[(d.q3 | d.q4) & d.cls.isin(["A", "FGK"])] cosd = np.cos(np.radians(s.DE_ICRS.values)) up = pd.DataFrame({"sid": s.sid.values, "ra": (s.RA_ICRS.values - s.pmRA.fillna(0).values * 16.0 / 3.6e6 / cosd) % 360.0, "dec": s.DE_ICRS.values - s.pmDE.fillna(0).values * 16.0 / 3.6e6}) rows = [] for name, cat in (("Cotten & Song 2016 Prime", "vizier:J/ApJS/225/15/table3"), ("Cotten & Song 2016 Reserved", "vizier:J/ApJS/225/15/table4"), ("Patel et al. 2014", "vizier:J/ApJS/212/10/catalog")): cache = os.path.join(RAW, "vet", "recall_" + cat.split("/")[-1] + "_" + cat.split("/")[-2] + ".csv") if os.path.exists(cache): txt = open(cache).read() else: txt = cds_xmatch(up.to_csv(index=False, float_format="%.6f"), cat, 5) with open(cache, "w") as f: f.write(txt) m = pd.read_csv(io.StringIO(txt), low_memory=False) ids = set(m.sid) x = s[s.sid.isin(ids)] rows.append(dict(catalogue=name, n_in_searchable_sample=len(x), n_flagged_excess=int(x.det.sum()), recall=float(x.det.mean()) if len(x) else np.nan, n_with_chi_gt3=int(((x.chi3 > 3) | (x.chi4 > 3)).sum()), median_max_chi_of_missed=float(np.fmax(x.chi3, x.chi4)[~x.det].median()) if len(x) else np.nan)) # Sgro & Song (2021): Gaia DR2 K/M dwarfs within 100 pc with W3/W4 excess. No VizieR table # exists, so the stars SIMBAD attaches to the paper are used (candidates, known excess stars # they re-analyse, and a few comparison objects such as Vega; matched by Gaia DR3 identifier). sgf = os.path.join(RAW, "vet", "simbad_sgro2021.csv") if os.path.exists(sgf): ids = set(pd.read_csv(sgf).sid) par = d[d.sid.isin(ids)] x = par[par.q3 | par.q4] rows.append(dict(catalogue="Sgro & Song 2021 (objects SIMBAD attaches to the paper; all classes)", n_in_parent_sample=len(par), n_in_searchable_sample=len(x), n_flagged_excess=int(x.det.sum()), recall=float(x.det.mean()) if len(x) else np.nan, n_with_chi_gt3=int(((x.chi3 > 3) | (x.chi4 > 3)).sum()), median_max_chi_of_missed=float(np.fmax(x.chi3, x.chi4)[~x.det].median()) if len(x) else np.nan)) out = pd.DataFrame(rows) out.to_csv(os.path.join(DATA, "recall_known_discs.csv"), index=False, float_format="%.4g") print(out.to_string()) return out def main(): d = inj.load_stars() cal = inj.load_cal() cand = pd.read_csv(os.path.join(DATA, "candidates.csv")) q3, q4 = d.q3.values.copy(), d.q4.values.copy() rows = [] rows += variant_rows("baseline (5 sigma, E(B-V) < 0.15, TOL 0.2)", d, cal, 5.0, q3, q4, cand) for thr in (4.0, 6.0): rows += variant_rows(f"threshold {thr:g} sigma", d, cal, thr, q3, q4) ebv0 = ex.EBV_MAX for ebv in (0.10, 0.30): ex.EBV_MAX = ebv _, a3, a4 = ex.quality(d) rows += variant_rows(f"E(B-V) < {ebv:.2f}", d, cal, 5.0, a3, a4) ex.EBV_MAX = ebv0 for tol in (0.0, 0.5): rows += variant_rows(f"tolerance TOL = {tol:g}", d, cal, 5.0, q3, q4, cand, tol=tol) # vetting variants c2 = cand.copy() c2.loc[c2.explanation == "background_visual", "explanation"] = "unexplained" rows += variant_rows("no manual image inspection", d, cal, 5.0, q3, q4, c2) c3 = c2.copy() c3.loc[c3.explanation == "background_image_offset", "explanation"] = "unexplained" rows += variant_rows("no image information at all (catalogues only)", d, cal, 5.0, q3, q4, c3) c4 = cand.copy() c4.loc[c4.explanation == "catalogued_excess_only", "explanation"] = "unexplained" rows += variant_rows("catalogue-only identifications counted as unexplained", d, cal, 5.0, q3, q4, c4, vpass=lm.vetting_pass_fraction(cand)) keep = (d.GCNSprob > 0.5).values rows += variant_rows("parent sample with GCNS probability > 0.5", d[keep].reset_index(drop=True), cal, 5.0, q3[keep], q4[keep], cand) rob = pd.DataFrame(rows) # single-anchor variants: counts only arows = [] for anc, lab in (("K", "Ks anchor only"), ("W", "W1 anchor only")): for cls, members in CLS.items(): sel = d.cls.isin(members).values pos = (q3 & (d[f"chi3{anc}"].values >= 5)) | (q4 & (d[f"chi4{anc}"].values >= 5)) neg = (q3 & (d[f"chi3{anc}"].values <= -5)) | (q4 & (d[f"chi4{anc}"].values <= -5)) arows.append(dict(variant=lab, cls=cls, n_searchable=int(((q3 | q4) & sel).sum()), n_excess=int((pos & sel).sum()), n_mirror_negative=int((neg & sel).sum()))) rob = pd.concat([rob, pd.DataFrame(arows)], ignore_index=True) rob.to_csv(os.path.join(DATA, "robustness.csv"), index=False, float_format="%.4g") pd.set_option("display.width", 250) print(rob[rob.cls.isin(["FGK", "M", "all"])][["variant", "cls", "n_searchable", "n_excess", "n_mirror_negative", "N_rec_f0.01", "k_all_f0.01", "UL95_all_f0.01", "k_strict_f0.01", "UL95_strict_f0.01", "k_unexpl_f0.01", "UL95_unexpl_f0.01", "k_all_f0.1", "UL95_all_f0.1", "k_strict_f0.1", "UL95_strict_f0.1", "k_unexpl_f0.1", "UL95_unexpl_f0.1"]].to_string(index=False)) # control estimate with a smaller W4 trigger area alt = ctl.main(r_max=ctl.R_MAX_ALT, verbose=False, d=d) alt.to_csv(os.path.join(DATA, "control_expectation_w4cap6.csv"), index=False, float_format="%.4g") print(alt[alt.band == "W3 or W4"].to_string(index=False)) recall_known(d) if __name__ == "__main__": main()