"""Revision step (referee minor point 11 and point (b)): injection-recovery test of the VETTING stage. What is tested. A single-element salt of S dex is injected into the real spectra of a random sample of stars drawn from the effective sample (the stars the limit actually rests on, including low-S/N ones and stars whose element is not detected), and the injected stars are passed through the same twin test and the same classification as the real candidates. This measures, for a salt of exactly S dex ADDED to a real star (with that star's real natural abundance, real noise, real twins), how often the spectrum test would wrongly refute it. How. The twin test only uses, per line, the integrated extra absorption of the target relative to its twins in the observed spectrum (d_obs) and in the pipeline's best-fit model (d_mod). The pipeline's own models provide the response of each line to the element's abundance: for every star already vetted whose catalogue abundance is high, d_mod / (10^x_cat - 10^x_twin) is the model equivalent width per unit abundance; a log-linear fit in Teff, logg and [Fe/H] gives the response R_line for any star. The injection adds A_line = R_line * 10^pred * (10^S - 1) to d_obs (the star now has S dex more of the element than its natural prediction) and, on the model side, assumes that the catalogue would report the total: x_cat' = log10(10^x_fit + 10^pred (10^S - 1)), capped at the edge of the model grid (the pipeline then raises flag 8), with the matching model absorption. What is NOT tested. The response of the GALAH DR4 neural-network fit itself to a single-element anomaly (whether the catalogue would in fact report x_cat', and whether flag_sp would stay 0) cannot be tested without rerunning the survey pipeline; that assumption remains a limitation. """ import csv, json import numpy as np import pandas as pd from common import * import s03_candidates as s3 import s04_vetting as s4 import s06_limits as s6 ELEMS = ['nd', 'la'] N_REC = 250 SALTS = [0.75, 1.0, 1.25, 1.5] SEED = 20261013 MIN_DN = 2.0 # rows with 10^x_cat - 10^x_twin above this calibrate the line response def line_response(vet_raw, d): """Per (element, line): robust log-linear fit of the model response d_mod / (10^x_cat - 10^x_twin) [mA per unit of 10^[X/Fe]] against Teff, logg, [Fe/H].""" out = {}; rows = [] for e in ELEMS: g = vet_raw[(vet_raw.element == e) & (vet_raw.status == 'ok')] dn = 10 ** g.x_cat.values - 10 ** g.x_twin.values idx = g.idx.values.astype(int) acc = {} for k in range(len(g)): if not np.isfinite(dn[k]) or dn[k] < MIN_DN: continue for w0, dob, sg, dmo in s4.parse_lines(g.lines.values[k]): if dmo > 0: acc.setdefault(w0, []).append((d['teff'][idx[k]], d['logg'][idx[k]], d['fe_h'][idx[k]], np.log10(dmo / dn[k]))) for w0, v in sorted(acc.items()): v = np.array(v, float) A = np.column_stack([np.ones(len(v)), (v[:, 0] - 5500.0) / 1000.0, v[:, 1] - 4.2, v[:, 2]]) y = v[:, 3]; keep = np.ones(len(v), bool) for _ in range(4): # sigma-clipped least squares b = np.linalg.lstsq(A[keep], y[keep], rcond=None)[0] res = y - A @ b; sc = 1.4826 * np.median(np.abs(res[keep] - np.median(res[keep]))) + 1e-6 keep = np.abs(res) < 3 * sc out[(e, w0)] = b rows.append(dict(element=e, line=w0, n=int(len(v)), const_log10_mA=b[0], per_1000K=b[1], per_dex_logg=b[2], per_dex_feh=b[3], scatter_dex=float(sc))) pd.DataFrame(rows).to_csv(DATA + 'injection_line_response.csv', index=False, float_format='%.4f') return out def main(verbose=True): d = s3.load(); c = s3.compute(d) vet, _ = s4.evaluate() # sets the per-line noise inflation infl = dict(s4.NOISE_INFLATION) raw = pd.read_csv(DATA + 'vetting_twin_test_raw.csv', keep_default_na=False, na_values=['nan']) resp = line_response(raw, d) rng = np.random.default_rng(SEED) todo = []; rec = {} for e in ELEMS: w = s3.eff_weight(d, c, e, s4.T_VET) pool = np.where(w > 0)[0] pick = np.sort(rng.choice(pool, N_REC, replace=False, p=w[pool] / w[pool].sum())) rec[e] = pick todo += [(int(i), e, 'receiver') for i in pick] rows = s4.measure_compact(d, c, todo, 'injection_receivers_raw.csv', verbose=verbose) # --- build the injected table inj = [] for r in rows: if r['status'] != 'ok': continue e = r['element']; i = r['idx']; edge = DR4_UPPER[e] a = s4.parse_lines(r['lines']) x_fit = r['x_cat'] if np.isfinite(r['x_cat']) else r['pred'] t = np.array([1.0, (d['teff'][i] - 5500.0) / 1000.0, d['logg'][i] - 4.2, d['fe_h'][i]]) R = np.array([10 ** float(resp[(e, w0)] @ t) if (e, w0) in resp else np.nan for w0 in a[:, 0]]) ok = np.isfinite(R) if ok.sum() == 0: continue a = a[ok]; R = R[ok] for S in [0.0] + SALTS: dn = 10 ** r['pred'] * (10 ** S - 1.0) # abundance added (units of 10^[X/Fe]) x_tot = np.log10(10 ** x_fit + dn) x_cat = min(x_tot, edge) # the catalogue cannot report more than the grid edge dn_mod = 10 ** x_cat - 10 ** x_fit lines = ';'.join('%.1f:%.2f:%.2f:%.2f' % (w0, dob + Rl * dn, sg, dmo + Rl * dn_mod) for (w0, dob, sg, dmo), Rl in zip(a, R)) inj.append(dict(kind='inj_%.2f' % S, element=e, sobject_id=r['sobject_id'], idx=i, x_cat=float(x_cat), ex_cat=float(x_cat - r['pred']), pred=r['pred'], x_twin=r['x_twin'], x_dr3=np.nan, n_twins=r['n_twins'], twin_mode=r['twin_mode'], n_lines=len(a), lines=lines, status='ok')) with open(DATA + 'injection_raw.csv', 'w', newline='') as fh: w = csv.DictWriter(fh, fieldnames=s4.RAW_COLS); w.writeheader() for r in inj: w.writerow({k: (f'{r[k]:.4f}' if isinstance(r[k], float) else r[k]) for k in s4.RAW_COLS}) v, _ = s4.evaluate('injection_raw.csv', infl=infl) idx = v.idx.values.astype(int) v['snr_px_ccd2'] = d['snr_px_ccd2'][idx].astype(float); v['teff'] = d['teff'][idx].astype(float) v['salt'] = [float(k.split('_')[1]) for k in v.kind] out = [] for e in ELEMS: sec_frac = float(c['sec'][e][rec[e]].mean()) for S in SALTS: g = v[(v.element == e) & (np.abs(v.salt - S) < 1e-6)] T = S - s3.MARGIN cand = g.ex_cat >= T # catalogued as a candidate under the stated assumption ver_true = np.array([s4.classify(r, T, s3.MARGIN) for _, r in g.iterrows()]) # refuted if excess < S excluded ver_add = np.array([s4.classify(r, T, s3.MARGIN - s6.M_REF) for _, r in g.iterrows()]) # ... < S - M_REF n = int(cand.sum()) k_true = int(((ver_true != 'refuted') & cand).sum()); k_add = int(((ver_add != 'refuted') & cand).sum()) k_sup = int(((ver_add == 'supported') & cand).sum()) hi = g.snr_px_ccd2 >= 50 out.append(dict(element=e, salt_added_dex=S, n_receivers=int(len(g)), frac_element_secure_in_catalogue=sec_frac, n_catalogued_as_candidate=n, not_refuted_at_S=k_true, eff_refute_at_S=k_true / max(n, 1), eff_refute_at_S_lower95=1 - binom_upper(n - k_true, n), not_refuted_at_S_minus_mref=k_add, eff_refute_at_S_minus_mref=k_add / max(n, 1), eff_refute_at_S_minus_mref_lower95=1 - binom_upper(n - k_add, n), frac_supported_3sigma=k_sup / max(n, 1), eff_refute_at_S_snr_ge50=float((ver_true[(cand & hi).values] != 'refuted').mean()) if (cand & hi).sum() else np.nan, eff_refute_at_S_minus_mref_snr_ge50=float((ver_add[(cand & hi).values] != 'refuted').mean()) if (cand & hi).sum() else np.nan, n_snr_ge50=int((cand & hi).sum()), analytic_efficiency_assumed=float(s6.norm.cdf(s6.M_REF / c['sig'][e])), median_rho=float(g.rho.median()), median_ex_spec=float(g.ex_spec.median()), median_ex_spec_minus_S=float((g.ex_spec - S).median()), rms_ex_spec_minus_S=float(1.4826 * (g.ex_spec - S - (g.ex_spec - S).median()).abs().median()))) # false-positive behaviour of the un-injected receivers (S = 0): how many look 'supported' at 3 sigma? g0 = v[(v.element == e) & (v.salt == 0.0)] z = (g0.dew_obs / g0.dew_sig) out.append(dict(element=e, salt_added_dex=0.0, n_receivers=int(len(g0)), n_catalogued_as_candidate=int((g0.ex_cat >= 0.75).sum()), frac_supported_3sigma=float((z > 3).mean()), median_ex_spec=float(g0.ex_spec.median()))) res = pd.DataFrame(out) res.to_csv(DATA + 'injection_recovery.csv', index=False, float_format='%.4f') # --- OPTIMISTIC variant, not used for any headline number: the injection says a real salt would # almost always be recovered at >= 3 sigma, so the 'ambiguous' candidates (not refuted, but below # 3 sigma) are unlikely to be salts. Counting only spectrum-supported candidates, with the # measured 3-sigma recovery fraction (95% lower bound) and the catalogue completeness as the # efficiency, shows how much of the limit is confusion from poor spectra. opt = [] try: key = json.load(open(DATA + 'key_numbers.json')); lim = pd.read_csv(DATA + 'limits_by_element.csv') for e in ELEMS: for S in (1.0, 1.25, 1.5): l = lim[(lim.element == e) & (np.abs(lim.salt_level_dex - S) < 1e-6)].iloc[0] r_ = res[(res.element == e) & (np.abs(res.salt_added_dex - S) < 1e-6)].iloc[0] n = int(r_.n_catalogued_as_candidate); ks = int(round(r_.frac_supported_3sigma * n)) eff3 = (1 - binom_upper(n - ks, n)) * (n / r_.n_receivers) effo = eff3 * key['catalogue_completeness_used'] opt.append(dict(element=e, salt_level_dex=S, n_supported=int(l.n_supported), N_eff=float(l.N_eff_no_extrapolation), efficiency=float(effo), f_ul95_supported_only=float(poisson_upper(int(l.n_supported)) / (l.N_eff_no_extrapolation * effo)))) pd.DataFrame(opt).to_csv(DATA + 'limits_supported_only_optimistic.csv', index=False, float_format='%.5g') except Exception as ex: print('optimistic variant skipped:', ex) v.drop(columns=['lines']).to_csv(DATA + 'injection_twin_test.csv', index=False, float_format='%.4f') json.dump(dict(n_receivers=N_REC, seed=SEED, m_ref=s6.M_REF, optimistic_supported_only=opt, table=res.replace({np.nan: None}).to_dict('records')), open(DATA + 'injection_summary.json', 'w'), indent=1) if verbose: pd.set_option('display.width', 250); pd.set_option('display.max_columns', 40) print(pd.read_csv(DATA + 'injection_line_response.csv').to_string(index=False)) print(res.to_string(index=False)) return res if __name__ == '__main__': main()