"""Convert candidate counts into 95% upper limits on the fraction of FGK stars dimmed by >= f. Expected number of flagged dimmed stars for a true population fraction p(f): N_exp = p * sum_c N_sens,c(f) * eps_c(f) * R_c(f) * C(f) * V N_sens,c clean stars of group c (survey x Teff class x main sequence / subgiant) whose local calibration scatter makes f a >=5 sigma effect eps_c MODELLED probability that a dimmed star is in the clean sample relative to an undimmed one (selection.py; cross-checked there with the model-free "twin" estimate) R_c probability that a dimmed star in the clean sample is flagged (injection into the real D distribution) C end-to-end correction to R (injection.py): ratio of the recovery when the dimming is put into the catalogue magnitudes and the calibration is redone, to the analytic recovery V probability that a truly dimmed, flagged star survives the vetting rules (measured on overlaps) Upper limit p95 = PoissonUL95(n_unresolved) / exposure. "Sample" limits omit eps (fraction of stars as they appear in the spectroscopic samples); "population" limits include it. Run once per astrometric cut set (DIM_STARS_CUTSET); products of the "orig" set go to data/origcuts/. Writes data/limits.csv, data/limits_by_class.csv, data/limits_scope.csv, data/vetting_efficiency.csv """ import os import warnings import numpy as np import pandas as pd import dimlib as L from analysis import CLASSES NSIG_THR = L.NSIG_THR from common import DATA, WORK, RAW, cache as cpath, dpkl, log, xmatch_upload from injection import correction warnings.filterwarnings("ignore") SURVEYS = ("apogee", "lamost", "gspspec") def vetting_efficiency(C): """Probability that a genuinely dimmed star flagged in survey A would be wrongly 'resolved'.""" rows = [] # (i) second-survey rule: P(D_B >= 2 sigma_B | D_A < 2 sigma_A) for stars clean in both surveys for a in SURVEYS: for b in SURVEYS: if a == b: continue m = C[a][["source_id", "D", "sigma_i"]].merge(C[b][["source_id", "D", "sigma_i"]], on="source_id", suffixes=("_a", "_b")) if len(m) < 100: continue ina = m.D_a < NSIG_THR * m.sigma_i_a pfd = np.mean(m.D_b[ina] >= NSIG_THR * m.sigma_i_b[ina]) rows.append(dict(rule="second survey", survey=a, other=b, n_overlap=len(m), frac_of_sample=len(m) / len(C[a]), p_false_dismiss=pfd)) return pd.DataFrame(rows) def neighbour_false_dismiss(C): """Fraction of clean stars that have a Gaia neighbour which would trigger the 'brighter neighbour' rule if the star were dimmed by f (control sample of 20,000 stars per survey).""" cache = cpath("control_neighbours.csv") ctrl = pd.concat([C[sv].sample(n=min(20000, len(C[sv])), random_state=7).assign(survey=sv) for sv in SURVEYS]) ids = " ".join(map(str, sorted(set(ctrl.source_id.values)))) if os.path.exists(cache) and os.path.exists(cache + ".ids") and open(cache + ".ids").read() == ids: nb = pd.read_csv(cache) else: # the control sample changes whenever a clean sample changes up = ctrl[["source_id", "ra", "dec"]].drop_duplicates("source_id").rename( columns={"source_id": "cand_id", "ra": "ra_c", "dec": "dec_c"}) log("XMatch: Gaia neighbours for %d control stars" % len(up)) nb = xmatch_upload(up, "ra_c", "dec_c", "vizier:I/355/gaiadr3", 5.0, cols2="Source,Gmag") nb.to_csv(cache, index=False) open(cache + ".ids", "w").write(ids) nb = nb[nb.Source != nb.cand_id] rows = [] for sv in SURVEYS: cs = ctrl[ctrl.survey == sv][["source_id", "G"]] x = nb.merge(cs.rename(columns={"source_id": "cand_id", "G": "G_c"}), on="cand_id") for f in L.F_GRID: dm = -2.5 * np.log10(1 - f) dG = x.Gmag - (x.G_c + dm) trig = ((x.angDist < 5.0) & (dG < 0.0)) | ((x.angDist < 2.5) & (dG < 1.0)) rows.append(dict(rule="brighter neighbour", survey=sv, f=f, p_false_dismiss=x.cand_id[trig].nunique() / len(cs))) return pd.DataFrame(rows) def run(): C = {} for sv in SURVEYS: s = pd.read_pickle(dpkl(sv)) C[sv] = s[s.clean].reset_index(drop=True) C[sv]["survey"] = sv del s eps = pd.read_csv(os.path.join(DATA, "selection_eff.csv")) # model-free "twin" estimate of eps (APOGEE and LAMOST only): first disc model = default, plus the range tw = pd.read_csv(os.path.join(DATA, "selection_eff_twin.csv")) tw0 = tw[tw.model == tw.model.iloc[0]].set_index(["survey", "cls", "f"]).eps_twin twlo = tw.groupby(["survey", "cls", "f"]).eps_twin.min() twhi = tw.groupby(["survey", "cls", "f"]).eps_twin.max() Ccorr = {(sv, f): correction(sv, f) for sv in SURVEYS for f in L.F_GRID} cand = pd.read_csv(os.path.join(DATA, "candidates_all.csv")) ve = vetting_efficiency(C) vn = neighbour_false_dismiss(C) pd.concat([ve, vn], ignore_index=True).to_csv(os.path.join(DATA, "vetting_efficiency.csv"), index=False) print(ve.round(4).to_string(index=False)) print(vn.round(4).to_string(index=False)) def V(sv, f): """Vetting survival probability for a truly dimmed star flagged in survey sv.""" p2 = 0.0 for _, r in ve[ve.survey == sv].iterrows(): p2 += r.frac_of_sample * r.p_false_dismiss if sv == "lamost": # repeat spectra: assume the same false-dismissal rate as LAMOST-vs-other surveys frac_rep = np.mean(C[sv].nspec > 1) pr = ve[ve.survey == sv].p_false_dismiss.mean() if (ve.survey == sv).any() else 0.1 p2 += frac_rep * pr p3 = float(vn[(vn.survey == sv) & (vn.f == f)].p_false_dismiss.iloc[0]) return max(0.0, 1.0 - p2 - p3) def exposure(tab, sv_of_row, f, require): """tab: clean stars (cls, ms, D, sigma_i, m_ms, m_cmd). Exposure pieces by survey, Teff class and main sequence (logg >= 4.1) / subgiant.""" sens, rec = L.flag_mask(tab, f, inject=True, require=require) out = [] for sv in np.unique(sv_of_row): for cname, _, _ in CLASSES: for ms in (True, False): q = (sv_of_row == sv) & (tab.cls.values == cname) & (tab.ms.values == ms) & sens n = int(q.sum()) if n == 0: continue e = eps[(eps.survey == sv) & (eps.cls == cname) & (eps.f == f)].iloc[0] k = (sv, cname, f) out.append(dict(survey=sv, cls=cname, ms=ms, f=f, n_sens=n, R=float(rec[q].mean()), C=Ccorr[(sv, f)], eps=e.eps, eps_lo=e.eps_lo, eps_hi=e.eps_hi, eps_smin=e.eps_sys_min, eps_smax=e.eps_sys_max, eps_twin=tw0.get(k, e.eps), eps_twin_min=twlo.get(k, e.eps_sys_min), eps_twin_max=twhi.get(k, e.eps_sys_max), V=V(sv, f))) return pd.DataFrame(out) rows, byc = [], [] TIERS = (("three-way", ("spec", "ms", "cmd"), "flag_"), ("spectroscopic-only", ("spec",), "flagS_")) # ------------------------------------------------------------ per survey for sv in SURVEYS: c = C[sv] for tier, req, col in TIERS: for f in L.F_GRID: ex = exposure(c, np.full(len(c), sv, dtype=object), f, req) ex["sample"], ex["tier"] = sv, tier byc.append(ex) cc = cand[(cand.survey == sv) & cand["%s%g" % (col, f)]] rows.append(summarise(sv, tier, f, ex, cc, len(c))) # ------------------------------------------------------------ combined (unique stars, best survey per star) allc = pd.concat([C[sv][["source_id", "survey", "cls", "D", "sigma_i", "m_ms", "m_cmd", "ms", "dist", "glat", "G", "feh"]] for sv in SURVEYS], ignore_index=True) allc = allc.sort_values("sigma_i").drop_duplicates("source_id").reset_index(drop=True) st = dict(n_unique=len(allc), frac_main_sequence=float(allc.ms.mean()), n_from_apogee=int((allc.survey == "apogee").sum()), n_from_lamost=int((allc.survey == "lamost").sum()), n_from_gspspec=int((allc.survey == "gspspec").sum())) for nm, v in (("dist", allc.dist.values), ("absb", np.abs(allc.glat.values)), ("G", allc.G.values)): for q in (5, 50, 95): st["%s_p%d" % (nm, q)] = float(np.percentile(v, q)) st["frac_absb_gt_17"] = float(np.mean(np.abs(allc.glat.values) > 17)) for cn, _, _ in CLASSES: st["frac_%s" % cn] = float(np.mean(allc.cls.values == cn)) pd.DataFrame([st]).to_csv(os.path.join(DATA, "combined_sample_stats.csv"), index=False) print(pd.DataFrame([st]).T.to_string()) key = set(zip(allc.source_id.values, allc.survey.values)) inkey = np.array([k in key for k in zip(cand.source_id.values, cand.survey.values)]) for tier, req, col in TIERS: for f in L.F_GRID: ex = exposure(allc, allc.survey.values, f, req) ex["sample"], ex["tier"] = "combined", tier byc.append(ex) cc = cand[cand["%s%g" % (col, f)].values & inkey] rows.append(summarise("combined", tier, f, ex, cc, len(allc))) # ------------------------------------------------------------ combined, main-sequence stars only (logg >= 4.1) ams = allc[allc.ms.values].reset_index(drop=True) cms = cand[cand.ms.values & inkey] for tier, req, col in TIERS: for f in L.F_GRID: ex = exposure(ams, ams.survey.values, f, req) ex["sample"], ex["tier"] = "combined_MS", tier byc.append(ex) rows.append(summarise("combined_MS", tier, f, ex, cms[cms["%s%g" % (col, f)].values], len(ams))) out = pd.DataFrame(rows) out.to_csv(os.path.join(DATA, "limits.csv"), index=False) byc = pd.concat(byc, ignore_index=True) byc.to_csv(os.path.join(DATA, "limits_by_class.csv"), index=False) # ------------------------------------------------------------ which stars carry the limit (combined, three-way) scope = [] ex3 = byc[(byc["sample"] == "combined") & (byc.tier == "three-way")] c3 = cand[inkey] for f in L.F_GRID: e = ex3[ex3.f == f].copy() e["E"] = e.n_sens * e.R * e.C * e.V * e.eps un = c3[c3["flag_%g" % f].values & (c3["status_%g" % f] == "unresolved").values] groups = [("%s %s" % (cn, "dwarfs" if ms else "subgiants"), (e.cls == cn) & (e.ms == ms), (un.cls == cn) & (un.ms == ms), (allc.cls == cn) & (allc.ms == ms)) for cn, _, _ in CLASSES for ms in (True, False)] groups += [("all dwarfs (logg >= %.1f)" % L.LOGG_MS, e.ms == True, un.ms == True, allc.ms == True), ("all subgiants (logg < %.1f)" % L.LOGG_MS, e.ms == False, un.ms == False, allc.ms == False)] groups += [("by survey: " + sv, e.survey == sv, un.survey == sv, allc.survey == sv) for sv in SURVEYS] for name, qe, qu, qa in groups: g = e[qe] E = float(g.E.sum()) scope.append(dict(f=f, group=name, n_stars=int(qa.sum()), n_sens=int(g.n_sens.sum()), mean_R=float((g.n_sens * g.R * g.C).sum() / max(g.n_sens.sum(), 1)), exposure=E, share_of_exposure=E / float(e.E.sum()), n_unresolved=int(qu.sum()), UL95_pop=L.poisson_ul(int(qu.sum())) / E if E > 0 else np.nan)) scope = pd.DataFrame(scope) scope.to_csv(os.path.join(DATA, "limits_scope.csv"), index=False) print(scope.to_string(index=False, float_format=lambda x: "%.3g" % x)) pd.set_option("display.width", 280) pd.set_option("display.max_columns", 40) print(out.to_string(index=False, float_format=lambda x: "%.3g" % x)) return out def summarise(name, tier, f, ex, cc, nclean): n_flag = len(cc) n_unres = int((cc["status_%g" % f] == "unresolved").sum()) n_conf2 = int(((cc["status_%g" % f] == "unresolved") & cc["confirmed2_%g" % f]).sum()) base = ex.n_sens * ex.R * ex.C # sensitive stars x end-to-end recovery E_samp = float((base * ex.V).sum()) E_pop = float((base * ex.V * ex.eps).sum()) E_pop_smin = float((base * ex.V * ex.eps_smin).sum()) E_pop_smax = float((base * ex.V * ex.eps_smax).sum()) E_novet_samp = float(base.sum()) E_novet_pop = float((base * ex.eps).sum()) ul = L.poisson_ul(n_unres) return dict(sample=name, tier=tier, f=f, dimming_mag=-2.5 * np.log10(1 - f), n_clean=nclean, n_sens=int(ex.n_sens.sum()), mean_R=float((ex.n_sens * ex.R).sum() / ex.n_sens.sum()), mean_C=float(base.sum() / (ex.n_sens * ex.R).sum()) if (ex.n_sens * ex.R).sum() > 0 else np.nan, mean_eps=float((ex.n_sens * ex.eps).sum() / ex.n_sens.sum()), mean_eps_twin=float((ex.n_sens * ex.eps_twin).sum() / ex.n_sens.sum()), mean_V=float((ex.n_sens * ex.V).sum() / ex.n_sens.sum()), n_flagged=n_flag, n_unresolved=n_unres, n_unresolved_measured_twice=n_conf2, frac_flagged=n_flag / max(ex.n_sens.sum(), 1), UL95_sample_vetted=ul / E_samp, UL95_pop_vetted=ul / E_pop, UL95_pop_vetted_eps_sysmax=ul / E_pop_smax, UL95_pop_vetted_eps_sysmin=ul / E_pop_smin, # same with the model-free twin estimate of eps for APOGEE and LAMOST (GSP-Spec keeps the model) UL95_pop_vetted_eps_twin=ul / float((base * ex.V * ex.eps_twin).sum()), UL95_pop_vetted_eps_twin_max=ul / float((base * ex.V * ex.eps_twin_max).sum()), UL95_pop_vetted_eps_twin_min=ul / float((base * ex.V * ex.eps_twin_min).sum()), UL95_pop_vetted_without_C=ul / float((ex.n_sens * ex.R * ex.V * ex.eps).sum()), UL95_sample_unvetted=L.poisson_ul(n_flag) / E_novet_samp, UL95_pop_unvetted=L.poisson_ul(n_flag) / E_novet_pop, UL95_pop_zero_candidates=L.poisson_ul(0) / E_pop) if __name__ == "__main__": run()