"""Natural-pattern prediction, excess statistic, effective samples and raw candidate lists for GALAH DR4 FGK dwarfs/subgiants. Natural model: every neutron-capture abundance is a non-negative mixture of the solar s- and r-process patterns, 10^[X/Fe] = s_X * f_s + r_X * f_r (f_s = f_r = 1 is the scaled-solar mix). f_r follows the disc trend [Eu/Fe] = clip(-0.4 [Fe/H], 0, 0.4) (for probes other than Eu a secure Eu detection replaces it). f_s comes from the anchor element: heavy-peak probes (La, Ce, Nd, Sm, Eu): anchor Ba (secure detection required); light-peak probes (Rb, Sr, Zr, Mo, Ru): anchor Y (secure detection required); Ba probe: anchor Nd; Y probe: anchor Zr. For these two the anchor value is the secure detection if there is one, otherwise the star's own 90% detection threshold for the anchor (an upper bound), so the prediction is an upper bound and the excess a lower bound. Excess Ex_X = [X/Fe] - pred_X - zp_X(Teff); zp is the median residual, in 250 K bins, of ALL fitted values (detections and non-detections alike) in stars whose detection threshold is below +0.1 dex, which avoids the detection-threshold bias; it is applied to the anchored (not the bound) probes. """ import numpy as np, csv from scipy.stats import norm from common import * ANCHOR = {'la': 'ba', 'ce': 'ba', 'nd': 'ba', 'sm': 'ba', 'eu': 'ba', 'rb': 'y', 'sr': 'y', 'zr': 'y', 'mo': 'y', 'ru': 'y', 'ba': 'nd', 'y': 'zr'} BOUND_PROBES = ('ba', 'y') MARGIN = 0.25 # dex: a star salted by T+MARGIN must show an observed excess >= T LEVELS = [0.25, 0.5, 0.75, 1.0, 1.25] # thresholds T on the observed excess def load(): """DR4 cache plus per-star detection thresholds. The catalogue stores NaN for abundances that ran into the upper edge of the neural-network label range (flag bit 8, 'above limit'); these are exactly the stars a search for strong enhancements must not lose, so they are given the edge value (a lower bound on the true abundance) here.""" d = dict(np.load(CACHE + 'dr4.npz')) d.update(np.load(CACHE + 'dr4_x90.npz')) for e in NCAP: f = d['flag_' + e + '_fe']; x = d[e + '_fe'].astype('f4') fill = ((f & 8) > 0) & ~np.isfinite(x) x[fill] = DR4_UPPER[e] d[e + '_fe'] = x d['edgefill_' + e] = fill return d def secure(d, e): """Secure detection: reported as a detection (flag 0) or above the model range (flag bit 8), and the value is above the star's own 90% detection threshold.""" f = d['flag_' + e + '_fe']; x = d[e + '_fe'].astype(float) return ((f == 0) | ((f & 8) > 0)) & np.isfinite(x) & (x >= d['x90_' + e]) def at_edge(d, e): f = d['flag_' + e + '_fe']; x = d[e + '_fe'].astype(float) return ((f & 8) > 0) | (np.isfinite(x) & (x >= DR4_UPPER[e] - 0.03) & ((f & 3) == 0)) def base_mask(d): teff = d['teff'].astype(float); logg = d['logg'].astype(float); feh = d['fe_h'].astype(float) return (logg >= LOGG_MIN) & (teff < TEFF_MAX) & (teff >= TEFF_MIN) & np.isfinite(feh) def compute(d, sfrac=None, tol_fr=0.0): S = dict(S_FRAC if sfrac is None else sfrac); R = {e: 1 - S[e] for e in S} teff = d['teff'].astype(float); feh = d['fe_h'].astype(float) n = len(teff); base = base_mask(d) X = {e: d[e + '_fe'].astype(float) for e in NCAP} sec = {e: secure(d, e) for e in NCAP} eu_disc = np.clip(-0.4 * feh, 0.0, 0.4) + tol_fr fr_disc = 10 ** eu_disc fr_eu = np.where(sec['eu'], 10 ** np.clip(X['eu'], -0.5, 2.0), fr_disc) pred = {}; anchor_ok = {}; is_bound = {} for e in NCAP: a = ANCHOR[e] fr = fr_disc if e == 'eu' else fr_eu if e in BOUND_PROBES: val = np.where(sec[a], X[a], d['x90_' + a].astype(float)) ok = np.isfinite(val) is_bound[e] = ok & ~sec[a] else: val = X[a]; ok = sec[a]; is_bound[e] = np.zeros(n, bool) fs = np.maximum((10 ** np.where(ok, val, 0.0) - R[a] * fr) / S[a], 0.05) pred[e] = np.where(ok, np.log10(S[e] * fs + R[e] * fr), np.nan) anchor_ok[e] = ok # zero points zp = {}; tb = np.arange(TEFF_MIN, TEFF_MAX + 1, 250.0); ex = {}; sig = {} good = base & (d['flag_sp'] == 0) & (feh > -1) for e in NCAP: raw = X[e] - pred[e] f = d['flag_' + e + '_fe'] z = np.zeros(n); zrow = [] if e in BOUND_PROBES: zrow = [0.0] * (len(tb) - 1); sig[e] = np.nan else: cal = good & anchor_ok[e] & (d['x90_' + e] <= 0.1) & np.isin(f, [0, 1]) & np.isfinite(raw) glob = float(np.nanmedian(raw[cal])) if cal.sum() >= 30 else 0.0 for t0 in tb[:-1]: m = cal & (teff >= t0) & (teff < t0 + 250) med = float(np.nanmedian(raw[m])) if m.sum() >= 30 else glob z[(teff >= t0) & (teff < t0 + 250)] = med; zrow.append(med) r = (raw - z)[cal] sig[e] = float(1.4826 * np.nanmedian(np.abs(r - np.nanmedian(r)))) if cal.sum() >= 30 else np.nan zp[e] = zrow; ex[e] = raw - z; pred[e] = pred[e] + z return dict(base=base, X=X, sec=sec, pred=pred, anchor_ok=anchor_ok, ex=ex, zp=zp, tb=tb, sig=sig, is_bound=is_bound) def eff_weight(d, c, e, T, margin=None, x90=None): """Weight of each star in the effective sample for element e at observed-excess threshold T: the probability that the star, if its element e were raised by T + margin above the natural prediction, would be reported as a secure detection with observed excess >= T. Requirements: clean star (flag_sp == 0), anchor element securely detected, and the star's own 90% detection threshold x90 below the abundance it would then have (or below the upper edge of the DR4 label range, beyond which the pipeline raises flag 8 / clips). Stars with a secure measurement of e count 1 if their present excess is >= -margin; the others count Phi(margin / sigma_core), the probability that catalogue scatter does not pull them below T.""" margin = MARGIN if margin is None else margin x90 = d['x90_' + e].astype(float) if x90 is None else x90 pred = c['pred'][e]; edge = DR4_UPPER[e] detectable = c['base'] & (d['flag_sp'] == 0) & c['anchor_ok'][e] & np.isfinite(pred) & \ (x90 <= np.minimum(pred + T, edge - 0.03)) ex0 = c['ex'][e] known = c['sec'][e] & np.isfinite(ex0) sg = c['sig'][e] if np.isfinite(c['sig'][e]) else 0.25 w = np.zeros(len(x90)) w[detectable & known & (ex0 >= -margin)] = 1.0 if e not in BOUND_PROBES: w[detectable & ~known] = norm.cdf(margin / sg) return w def candidates(d, c, e, T): ok = c['base'] & (d['flag_sp'] == 0) & c['sec'][e] & c['anchor_ok'][e] ex = c['ex'][e] return ok & ((ex >= T) | (at_edge(d, e) & (ex >= min(T, 0.5)))) def main(): d = load(); c = compute(d) base = c['base'] & (d['flag_sp'] == 0) print('FGK dwarfs/subgiants:', c['base'].sum(), ' with flag_sp==0:', base.sum()) print('zero points (dex) per 250 K bin from %d K; core scatter sigma:' % TEFF_MIN) for e in NCAP: print(' ', e, ' '.join(f'{v:+.2f}' for v in c['zp'][e]), f' sigma={c["sig"][e]:.2f}') print('element: [N_eff, n_cand] at T =', LEVELS) for e in NCAP: row = [] for T in LEVELS: row.append(f'{eff_weight(d, c, e, T).sum():9.0f} {candidates(d, c, e, T).sum():6d}') print(f' {e:2s}', ' | '.join(row)) np.savez(CACHE + 'dr4_ex.npz', **{('ex_' + e): c['ex'][e].astype('f4') for e in NCAP}, **{('pred_' + e): c['pred'][e].astype('f4') for e in NCAP}, **{('sec_' + e): c['sec'][e] for e in NCAP}, **{('anc_' + e): c['anchor_ok'][e] for e in NCAP}) if __name__ == '__main__': main()