"""Core analysis functions: quality cuts, spectroscopic luminosity calibration, deficits.""" import itertools import os import numpy as np import pandas as pd from scipy.spatial import cKDTree # ------------------------------------------------------------------------------------------ # Pre-declared sample definition (identical for all surveys unless noted) # ------------------------------------------------------------------------------------------ TEFF_MIN, TEFF_MAX = 4000.0, 6500.0 # FGK LOGG_MIN = 3.5 # main sequence + subgiants LOGG_MS = 4.1 # logg >= LOGG_MS : "main-sequence" subsample; below: "subgiant" FEH_MIN, FEH_MAX = -1.0, 0.6 # Astrometric cuts. Two sets are carried through the whole pipeline: # "orig" parallax/error >= 20, RUWE < 1.4 : declared before looking at the data; # "final" parallax/error >= 50, RUWE < 1.2 : tightened AFTER a first look at the faint tail, whose # rate was found to be about 10 times higher at RUWE > 1.2 and at parallax/error < 50 than # elsewhere (data/astrometry_diagnostic.csv). This is a post-hoc choice, so every limit is # reported for both sets and the f = 0.25 and 0.5 limits are quoted as the range between them. CUTSETS = {"final": dict(poe_min=50.0, ruwe_max=1.2), "orig": dict(poe_min=20.0, ruwe_max=1.4)} CUTSET = os.environ.get("DIM_STARS_CUTSET", "final") POE_MIN = CUTSETS[CUTSET]["poe_min"] # parallax_over_error RUWE_MAX = CUTSETS[CUTSET]["ruwe_max"] SFD_MAX = 0.10 # E(B-V)_SFD of the full sightline DIST_MAX = 1000.0 / 1.0 # pc (parallax > 1 mas imposed at download) KNN_K = 50 KNN_SCALE = dict(teff=60.0, logg=0.08, feh=0.06) KNN_RMAX = 2.5 # max scaled distance to the K-th neighbour ("populated" parameter space) COLOUR_ZMAX = 4.0 # |colour residual| / local robust sigma, for BP-RP and G-Ks F_GRID = (0.25, 0.5, 0.75, 0.9) KNN_K_FLOOR = 200 # main-sequence neighbours for the logg-free / Teff-free references CMD_COLOUR_SCALE = 0.02 # mag of BP-RP per scaled unit in the (colour, [Fe/H]) neighbour search EDGE_MULT = 1.5 # floor = q95 + EDGE_MULT * (q95 - q84) of the main-sequence neighbours EDGE_MIN_MAG = 0.13 # mag; minimum distance of the floor below the main-sequence median NSIG_SENS = 5.0 # sensitive to f if |log10(1-f)| >= NSIG_SENS * sigma_i NSIG_THR = 2.0 # flagged if D < log10(1-f) + NSIG_THR * sigma_i (and below both floors) def cstar(bprp, excess): """Corrected BP/RP flux excess factor C* (Riello et al. 2021, eq. 6 / table 2).""" x = np.asarray(bprp, float) f = np.where(x < 0.5, 1.154545 + 0.033772 * x + 0.032277 * x ** 2, np.where(x < 4.0, 1.162004 + 0.011464 * x + 0.049255 * x ** 2 - 0.005879 * x ** 3, 1.057572 + 0.140537 * x)) return np.asarray(excess, float) - f def cstar_sigma(G): G = np.asarray(G, float) return 0.0059898 + 8.817481e-12 * G ** 7.618399 def common_cuts(d, poe_min=None, ruwe_max=None): """Ordered list of (name, boolean mask) for the cuts shared by every survey.""" POE_MIN_, RUWE_MAX_ = (POE_MIN if poe_min is None else poe_min), (RUWE_MAX if ruwe_max is None else ruwe_max) cs = cstar(d.BP - d.RP, d.bprp_excess) cuts = [ ("FGK dwarf/subgiant box: %g<=Teff<=%g, logg>=%g, %g<=[Fe/H]<=%g" % (TEFF_MIN, TEFF_MAX, LOGG_MIN, FEH_MIN, FEH_MAX), d.teff.between(TEFF_MIN, TEFF_MAX) & (d.logg >= LOGG_MIN) & d.feh.between(FEH_MIN, FEH_MAX)), ("Gaia 5/6-parameter solution, G/BP/RP present, not duplicated", d.solved.isin([31, 95]) & d.G.notna() & d.BP.notna() & d.RP.notna() & (d.dup == 0)), ("parallax > 1 mas and parallax/error >= %g" % POE_MIN_, (d.plx > 1.0) & (d.poe >= POE_MIN_)), ("low extinction sightline: E(B-V)_SFD < %g" % SFD_MAX, d.sfd < SFD_MAX), ("RUWE < %g" % RUWE_MAX_, d.ruwe < RUWE_MAX_), ("image/photometry clean: ipd_frac_multi_peak<=2, ipd_gof_harmonic_amplitude<0.1, |C*|<5 sigma", (d.ipd_fmp <= 2) & (d.ipd_gofha < 0.1) & (np.abs(cs) < 5 * cstar_sigma(d.G))), ("2MASS Ks: unique match <1 arcsec, Ks quality A/B, e_Ks<0.1", d.K.notna() & (d.e_K < 0.1) & (d.d2m < 1.0) & (d.n2m == 1) & d.qfl.astype(str).str[2:3].isin(["A", "B"])), ] return cuts def survey_cuts(d, survey): if survey == "apogee": return [("APOGEE: S/N>=50, no bad STARFLAG/ASPCAPFLAG bits, single CCF component", (d.snr >= 50) & d.flag_ok & (d.ncomp.fillna(1) <= 1))] if survey == "lamost": return [("LAMOST: snr_g>=30, no fibre offset, sigma(Teff)<150 K, sigma(logg)<0.25", (d.snr >= 30) & (d.offsets == 0) & (d.e_teff < 150) & (d.e_logg < 0.25))] if survey == "gspspec": fl = d.flags_gspspec.astype(str) good = np.ones(len(d), bool) for i in range(13): # first 13 quality-flag characters (Recio-Blanco et al. 2023) ch = fl.str[i:i + 1] lim = "1" if i < 7 else "0" # vbroad/vrad bias flags and flux-noise flag <= 1, the rest == 0 good &= (ch <= lim).values # (rvs_spec_sig_to_noise is published only for the ~1M sources with public RVS spectra, so the # flux-noise flag, which encodes the S/N-driven parameter uncertainty, is used instead of an S/N cut) return [("GSP-Spec: quality flags 1-7 <= 1 (incl. flux-noise flag), flags 8-13 == 0", good)] raise ValueError(survey) def apply_cuts(d, survey, verbose=True, poe_min=None, ruwe_max=None): rows = [] m = np.ones(len(d), bool) rows.append(("all stars in download box", int(m.sum()))) for name, c in survey_cuts(d, survey) + common_cuts(d, poe_min, ruwe_max): m &= np.asarray(c, bool) rows.append((name, int(m.sum()))) if verbose: for r in rows: print(" %9d %s" % (r[1], r[0])) return m, rows # ------------------------------------------------------------------------------------------ # Empirical "spectroscopic" prediction of a quantity y from (Teff, logg, [Fe/H]) # ------------------------------------------------------------------------------------------ def _design(T, g, f, deg=3): x, y, z = (T - 5500.0) / 1000.0, g - 4.4, f cols = [] for i, j, k in itertools.product(range(deg + 1), repeat=3): if i + j + k <= deg: cols.append(x ** i * y ** j * z ** k) return np.vstack(cols).T def poly_fit(T, g, f, y, deg=3, nit=4): A = _design(T, g, f, deg) ok = np.isfinite(y) for _ in range(nit): Ao = A[ok] coef = np.linalg.solve(Ao.T @ Ao, Ao.T @ y[ok]) # normal equations (20 well-scaled terms) r = y - A @ coef s = 1.4826 * np.nanmedian(np.abs(r[ok] - np.nanmedian(r[ok]))) ok = np.isfinite(y) & (np.abs(r) < 3 * s) return coef def poly_eval(coef, T, g, f, deg=3): return _design(T, g, f, deg) @ coef class SpecCalib: """Leave-one-out local-median calibration of observables on spectroscopic parameters. prediction(y) = cubic polynomial in (Teff, logg, [Fe/H]) + median residual of the K nearest neighbours in scaled parameter space (the star itself excluded). Because medians are used, a sub-per-cent population of dimmed stars cannot move the prediction. """ def __init__(self, d, k=KNN_K, scale=KNN_SCALE): self.k = k self.T, self.g, self.f = d.teff.values.astype(float), d.logg.values.astype(float), d.feh.values.astype(float) self.X = np.vstack([self.T / scale["teff"], self.g / scale["logg"], self.f / scale["feh"]]).T self.tree = cKDTree(self.X) dist, idx = self.tree.query(self.X, k=k + 1) # leave-one-out: drop the star itself. With tied parameter vectors (GSP-Spec values are quantised) the # star is not always returned first, so it is located explicitly; if it is not among the k+1 nearest # (more than k exact ties) the farthest neighbour is dropped instead. n = np.arange(len(self.X)) isself = idx == n[:, None] col = np.where(isself.any(axis=1), isself.argmax(axis=1), k) keep = np.ones(idx.shape, bool) keep[n, col] = False self.idx = idx[keep].reshape(len(n), k) self.rk = dist[:, -1] self.coef = {} self.res = {} def calibrate(self, name, y): y = np.asarray(y, float) coef = poly_fit(self.T, self.g, self.f, y) r = y - poly_eval(coef, self.T, self.g, self.f) nb = r[self.idx] med = np.nanmedian(nb, axis=1) sig = 1.4826 * np.nanmedian(np.abs(nb - med[:, None]), axis=1) self.coef[name] = coef self.res[name] = r return r - med, sig # leave-one-out residual (obs - pred) and local robust scatter def predict_external(self, name, T, g, f, scale=KNN_SCALE): """Prediction for stars that are not part of the calibration set (e.g. candidates' alternative parameters).""" X = np.vstack([np.asarray(T, float) / scale["teff"], np.asarray(g, float) / scale["logg"], np.asarray(f, float) / scale["feh"]]).T dist, idx = self.tree.query(X, k=self.k) med = np.nanmedian(self.res[name][idx], axis=1) return poly_eval(self.coef[name], np.asarray(T, float), np.asarray(g, float), np.asarray(f, float)) + med, dist[:, -1] def knn_floor(X, ref, ys, k): """Main-sequence reference from the k nearest *reference* stars (ref mask) in feature space X. For each observable y (an absolute magnitude) returns dev = y - median(y of neighbours) (positive = fainter than the reference median) e0 = q95 - q50, e1 = q95 - q84 (percentiles of y among the neighbours) The faint floor of the main-sequence band is then e0 + edge_mult * e1 below the median: the 95th percentile plus edge_mult times the 84-to-95 percentile blur (2.6 sigma for a Gaussian with edge_mult = 1.5). Leave-one-out. """ ref = np.asarray(ref, bool) tree = cKDTree(X[ref]) _, idx = tree.query(X, k=k + 1) refpos = np.full(len(X), -1) refpos[ref] = np.arange(ref.sum()) isself = idx[:, 0] == refpos # the star itself is its own nearest reference star idx = np.where(isself[:, None], idx[:, 1:], idx[:, :-1]) out = [] for y in ys: y = np.asarray(y, float) nb = y[ref][idx] # np.percentile is vectorised; np.nanpercentile loops over rows in Python (minutes for 5e5 stars) pct = np.nanpercentile if np.isnan(nb).any() else np.percentile q50, q84, q95 = pct(nb, [50, 84, 95], axis=1) out.append((y - q50, q95 - q50, q95 - q84)) return out FLOOR_COLS = ["fl_%s_%s%s" % (tag, q, b) for tag in ("ms", "cmd") for q in ("dev", "e0", "e1") for b in ("G", "K")] def set_floor(d, edge_mult=None): """(Re)compute D_ms, m_ms, D_cmd, m_cmd from the stored neighbour percentiles for a floor multiplier.""" edge_mult = EDGE_MULT if edge_mult is None else edge_mult for tag in ("ms", "cmd"): dg, dk = d["fl_%s_devG" % tag].values, d["fl_%s_devK" % tag].values eg = np.maximum(d["fl_%s_e0G" % tag].values + edge_mult * d["fl_%s_e1G" % tag].values, EDGE_MIN_MAG) ek = np.maximum(d["fl_%s_e0K" % tag].values + edge_mult * d["fl_%s_e1K" % tag].values, EDGE_MIN_MAG) d["D_" + tag] = np.maximum(-0.4 * dg, -0.4 * dk) d["m_" + tag] = np.maximum(-0.4 * (dg - eg), -0.4 * (dk - ek)) # < 0: below the floor in both bands d["floor_" + tag] = -0.4 * ek # floor relative to the median (Ks), dex return d def compute_deficits(d, k=KNN_K, scale=None, edge_mult=None, floors=None): """Add predicted magnitudes/colours and the deficits to a (cut) sample. Returns (d, calib). D deficit w.r.t. the full spectroscopic prediction M(Teff, logg, [Fe/H]) -> "spectroscopic luminosity" D_ms deficit w.r.t. the median main-sequence star (logg >= LOGG_MS) of the same (Teff, [Fe/H]); m_ms = D_ms - floor, negative when the star is below the faint edge of that main sequence D_cmd same with photometric colour instead of Teff: (BP-RP_0, [Fe/H]); m_cmd = D_cmd - floor All are log10(L_parallax / L_predicted); G and Ks are evaluated separately and the less deficient band is kept, so a grey deficit must show in both bands. floors: optionally a frame returned by an earlier call for the SAME stars in the same order; its main-sequence neighbour percentiles are reused (they do not depend on k, and a common rescaling of the neighbour scales leaves the 2-D neighbour sets unchanged). """ scale = scale or KNN_SCALE d = d.reset_index(drop=True).copy() cal = SpecCalib(d, k=k, scale=scale) for name in ("MG", "MK", "bprp0", "gk0"): r, s = cal.calibrate(name, d[name].values) d["r_" + name] = r d["s_" + name] = s d["rk"] = cal.rk d["D_G"] = -0.4 * d.r_MG # log10(L_parallax / L_spectroscopic), G band d["D_K"] = -0.4 * d.r_MK # same, Ks band d["D"] = np.maximum(d.D_G, d.D_K) # the less deficient band: a grey deficit must show in both d["z_bprp"] = d.r_bprp0 / np.maximum(d.s_bprp0, 0.01) d["z_gk"] = d.r_gk0 / np.maximum(d.s_gk0, 0.02) d["sig_D"] = 0.4 * np.maximum(d.s_MG, d.s_MK) if floors is not None: assert len(floors) == len(d) and np.array_equal(floors.source_id.values, d.source_id.values) for c in FLOOR_COLS: d[c] = floors[c].values else: # logg-free and Teff-free main-sequence references: nearest main-sequence stars in 2-D ref = (d.logg.values >= LOGG_MS) & np.isfinite(d.bprp0.values) X2 = np.vstack([d.teff.values / KNN_SCALE["teff"], d.feh.values / KNN_SCALE["feh"]]).T X3 = np.vstack([np.where(np.isfinite(d.bprp0.values), d.bprp0.values, 9.0) / CMD_COLOUR_SCALE, d.feh.values / KNN_SCALE["feh"]]).T for tag, X in (("ms", X2), ("cmd", X3)): (dg, e0g, e1g), (dk, e0k, e1k) = knn_floor(X, ref, (d.MG.values, d.MK.values), KNN_K_FLOOR) for nm, v in (("devG", dg), ("e0G", e0g), ("e1G", e1g), ("devK", dk), ("e0K", e0k), ("e1K", e1k)): d["fl_%s_%s" % (tag, nm)] = v set_floor(d, edge_mult) return d, cal def flag_mask(c, f, nsig_thr=None, nsig_sens=None, require=("spec", "ms", "cmd"), inject=False): """Stars flagged as dimmed by >= f. sensitive : sigma_i <= |log10(1-f)| / nsig_sens (f is at least an nsig_sens-sigma effect for this star) flagged : D < log10(1-f) + nsig_thr * sigma_i (deficit w.r.t. the spectroscopic luminosity) and m_ms < 0 (below the faint edge of the main sequence at its Teff and [Fe/H]; logg not used) and m_cmd < 0 (below the faint edge of the main sequence at its colour and [Fe/H]; Teff not used) With inject=True the grey dimming log10(1-f) is first added to all three statistics of every star (injection test); the returned mask then marks the recovered stars. Returns (sensitive, flagged). """ nsig_thr = NSIG_THR if nsig_thr is None else nsig_thr nsig_sens = NSIG_SENS if nsig_sens is None else nsig_sens Df = np.log10(1.0 - f) add = Df if inject else 0.0 sens = c.sigma_i.values <= abs(Df) / nsig_sens fl = sens.copy() if "spec" in require: fl &= (c.D.values + add) < Df + nsig_thr * c.sigma_i.values if "ms" in require: fl &= (c.m_ms.values + add) < 0.0 if "cmd" in require: fl &= (c.m_cmd.values + add) < 0.0 return sens, fl def robust_sigma(x): x = np.asarray(x, float) x = x[np.isfinite(x)] return 1.4826 * np.median(np.abs(x - np.median(x))) def poisson_ul(n, cl=0.95): """One-sided Poisson upper limit on the mean for n observed events.""" from scipy.stats import chi2 return 0.5 * chi2.ppf(cl, 2 * (n + 1))