#!/usr/bin/env python3 """Catalogue-level sensitivity fractions, candidate counts and 95% upper limits per (f, T). Hypothesis H(f, T): a fraction f of the galaxy's starlight is intercepted and re-radiated as a blackbody at temperature T. With g = f/(1-f) and the model ratios R_b(T, z) (s03_model.py) the waste-heat flux density in WISE band b of a galaxy observed with W1 flux F1 is S_b = g R_b F1 / (1 + g R_1) (the denominator removes the waste heat's own W1 contribution) Per galaxy and per hypothesis, using only catalogue quantities (AllWISE profile-fit fluxes): testable (1-eps) S_b >= K_S sigma_b in W3 or W4 -> the data could tell a host from an empty galaxy compatible F_b >= (1-eps) S_b - K_C sigma_b in W2, W3 AND W4 -> the galaxy is not excluded as a host candidate testable AND compatible eps = 0.2 is a tolerance for model/calibration error, K_S = 5, K_C = 2: round values, not optimised on the counts; data/limits_variants.csv shows how the limits move when each of them is changed (eps = 0.3 and 0.4 are the conservative choices if profile-fit band ratios were biased low; see s09_image_injection.py for the image-level test of that question). The Legacy Surveys DR9 forced photometry (unWISE coadds, optical-model shapes) is used only as a cross-check (variants 'legacy_w3_only' and 'allwise_and_legacy'): its W3/W4 errors are underestimated (see data/phot_checks.json), so they are inflated by the empirical factors LS_ERR_SCALE. Its usable subsample (fracflux_w3 < 0.5) is W3-bright selected, so 'allwise_on_legacy_subsample' gives the like-for-like AllWISE numbers for the same galaxies. Catalogue-level injection: the waste heat of H(f, T) (with 10% log-normal amplitude scatter) is added to the real catalogue fluxes of EVERY sample galaxy (each galaxy keeps its own depth, noise realisation and natural emission) and the same classification is run. NOTE (revision): because the injected flux and the threshold come from the same model, an injected host that is testable is recovered by construction (recov_eff is 0.996-1.000); the 'sensitivity fraction' sens_frac = recovered / all is therefore simply the fraction of galaxies in which the predicted signal exceeds K_S sigma. It is a model-based sensitivity correction, not an image-level completeness. Whether profile-fit photometry returns the modelled band ratios for an extended source is tested separately on images (s09_image_injection.py); the resulting response factors eta_b/eta_1 feed the variants 'psf_response_inject' (hosts injected with the measured response, thresholds unchanged) and 'psf_response_full' (thresholds also use the response; not used for any headline number). Limit: fraction of population galaxies hosting >= f at T <= CP95(N_cand ; N_testable) / (recovered/testable). Populations: all the full early-type sample -> "all outliers" limit, no assumption about hosts clean reliable WISE photometry AND no catalogue sign of a conventional power source (emission lines, radio, UV, post-starburst, AGN-like W1-W2) -> candidates here are by construction unexplained by the public catalogues; they are then inspected on images (s06_vetting.py) """ import json import os import numpy as np import pandas as pd import mirq_common as M from s03_model import Model K_S, K_C, EPS = 5.0, 2.0, 0.2 SDSS_MLIM = 17.77 RNG_SEED = 20261012 AGNCOL = 0.5 # W1-W2 (Vega) above which the hot-dust colour is called AGN-like NUVR_SF = 5.0 # NUV-r below which a galaxy is called UV-excess (residual star formation) # Empirical inflation of the Legacy Surveys W3/W4 flux errors. For line-less galaxies in the usable Legacy # subsample the spread of F/sigma below the median is 1.7 (from the 2.3rd percentile) to 1.9 (16th) in W3 and # 2.3 to 3.7 in W4 instead of ~1 (s04_baseline.py re-measures and records this in data/phot_checks.json as # empirical_error_scale_w3/4). Adopted: the larger of the two estimates per band. Revision: the first pass # hard-coded 2.6 and 7.2, which had been measured before the fracflux_w3 < 0.5 blending cut was introduced # and were never updated; the referee caught the mismatch. LS_ERR_SCALE = {2: 1.0, 3: 1.9, 4: 3.7} RESP_GRID = os.path.join(M.DATA, "psf_response_grid.csv") def load_master(): d = pd.read_csv(os.path.join(M.DATA, "master_slim.csv.gz"), dtype={"cc_flags": str, "ph_qual": str}, low_memory=False) d["cc_flags"] = d.cc_flags.fillna("9999").str.zfill(4) return d def phot(d, kind): """Dict of flux arrays (mJy) and validity mask for the chosen photometry.""" if kind == "psf": p = dict(F1=d.F1.values, F2=d.F2.values, F3=d.F3.values, F4=d.F4.values, e2=d.e2.values, e3=d.e3.values, e4=d.e4.values) ok = (d.sep_w.values <= 3.0) & (d.w1snr.values >= 5) elif kind == "ls": p = dict(F1=d.L1.values, F2=d.L2.values, F3=d.L3.values, F4=d.L4.values, e2=d.le2.values * LS_ERR_SCALE[2], e3=d.le3.values * LS_ERR_SCALE[3], e4=d.le4.values * LS_ERR_SCALE[4]) # DUP rows carry no flux; fracflux_w3 >= 0.5 means neighbours contribute at least half as much W3 flux # as the galaxy inside its own profile, where the forced fit is degenerate (checked on individual objects) ok = (d.sep_ls.values <= 1.5) & (d.ls_type.values != "DUP") & (d.fracflux_w3.values < 0.5) else: raise ValueError(kind) with np.errstate(invalid="ignore"): ok = ok & np.isfinite(p["F1"]) & (p["F1"] > 0) for b in (2, 3, 4): ok = ok & np.isfinite(p["F%d" % b]) & (p["e%d" % b] > 0) p["ok"] = ok return p class Engine: def __init__(self, d, use_ls=False, primary="psf", model="R", ks=K_S, kc=K_C, eps=EPS, test_bands=(3, 4), comp_bands=(2, 3, 4), sub_photosphere=False, sub_bands=(3, 4), resp=None): self.d = d self.m = Model(model) self.z = d.zhelio.values self.P = phot(d, primary) self.Q = phot(d, "ls" if primary == "psf" else "psf") if use_ls else None self.ks, self.kc, self.eps = ks, kc, eps self.test_bands, self.comp_bands, self.sub = test_bands, comp_bands, sub_photosphere # Bands in which the photosphere is subtracted when sub_photosphere is set. W2 is excluded by default: # the observed profile-fit W2/W1 of bright sample galaxies lies BELOW even the bluest Brown+14 template # for most of them (data/phot_checks.json, photosphere_template_check), so subtracting a template W2 # over-subtracts and would wrongly exclude real hosts of hot (400-600 K) waste heat at small f. self.sub_bands = sub_bands rng = np.random.default_rng(RNG_SEED) self._amp = np.exp(0.1 * rng.standard_normal(len(d))) # fixed injection amplitude scatter # profile-fit response to star-following emission, h_b = eta_b/eta_1 (s09_image_injection.py): # resp = None not used (primary) # resp = "inject" injected hosts are measured at h_b times the modelled band ratio; thresholds unchanged # resp = "full" the thresholds use h_b as well self.resp = resp self.h = {b: np.ones(len(d)) for b in (2, 3, 4)} if resp is not None: self.h = response_factors(d) def ratios(self, f, T): g = f / (1.0 - f) return g, {b: self.m.R_of(T, b, self.z) for b in (1, 2, 3, 4)} def _classify(self, p, g, R, inj=None): """inj = None: real fluxes. inj = amplitude array: inject the waste heat first.""" F1 = p["F1"] F = {b: p["F%d" % b] for b in (2, 3, 4)} if inj is not None: F = {b: F[b] + inj * g * R[b] * self.h[b] * p["F1"] for b in (2, 3, 4)} F1 = p["F1"] * (1 + inj * g * R[1]) den = 1.0 + g * R[1] F1star = F1 / den test = np.zeros(len(F1), bool) comp = np.ones(len(F1), bool) with np.errstate(invalid="ignore"): for b in (2, 3, 4): S = g * R[b] * F1star if self.resp == "full": S = S * self.h[b] e = p["e%d" % b] Fb = F[b] if self.sub and b in self.sub_bands: # remove the remaining stars (bluest Brown+14 elliptical) Fb = Fb - self.m.C_of(b, self.z, lower=True) * F1star if b in self.test_bands: test |= (1 - self.eps) * S >= self.ks * e if b in self.comp_bands: comp &= Fb >= (1 - self.eps) * S - self.kc * e return p["ok"] & test, p["ok"] & comp def classify(self, f, T, inject=False): g, R = self.ratios(f, T) inj = self._amp if inject else None test, comp = self._classify(self.P, g, R, inj) if self.Q is not None: _, comp_q = self._classify(self.Q, g, R, inj) comp = comp & (comp_q | ~self.Q["ok"]) return test, comp def run(self, f, T, sel=None): sel = np.ones(len(self.d), bool) if sel is None else sel test, comp = self.classify(f, T) test_i, comp_i = self.classify(f, T, inject=True) cand = test & comp & sel n, n_test = int(sel.sum()), int((test & sel).sum()) n_rec, n_test_i = int((test_i & comp_i & sel).sum()), int((test_i & sel).sum()) return dict(f=f, T=T, n=n, n_test=n_test, n_cand=int(cand.sum()), n_rec=n_rec, sens_frac=n_rec / n if n else np.nan, recov_eff=(n_rec / n_test_i) if n_test_i else np.nan), cand, test & sel def null_expectation(self, f, T, sel=None): """Expected number of candidates if every galaxy emitted ONLY the median old-star template in W3 and W4 (C_b F1) plus Gaussian noise of its catalogue sigma: the count produced by noise at the edge of testability. Real galaxies have more natural emission than the template, so this is a floor.""" from scipy.stats import norm sel = np.ones(len(self.d), bool) if sel is None else sel g, R = self.ratios(f, T) test, _ = self._classify(self.P, g, R) p = self.P F1star = p["F1"] / (1.0 + g * R[1]) prob = np.ones(len(F1star)) with np.errstate(invalid="ignore", divide="ignore"): for b in (3, 4): if b in self.comp_bands: thr = (1 - self.eps) * g * R[b] * F1star - self.kc * p["e%d" % b] prob = prob * norm.sf((thr - self.m.C_of(b, self.z) * p["F1"]) / p["e%d" % b]) return float(np.nansum(prob[test & sel])) def predicted_w12(self, f, T, z=0.08): """W1-W2 (Vega) a host of H(f,T) would show at redshift z (stellar template + waste heat).""" g = f / (1.0 - f) c2 = np.interp(z, self.m.z, self.m.C[0]) r1, r2 = (np.interp(z, self.m.z, self.m.R[int(np.where(np.isclose(self.m.T, T))[0][0]), b - 1]) for b in (1, 2)) return 2.5 * np.log10((c2 + g * r2) / (1 + g * r1) * M.WISE_F0[1] / M.WISE_F0[2]) def response_factors(d): """Per-galaxy h_b for b = 2, 3, 4: the factor by which AllWISE profile fitting is expected to return the waste heat of a host relative to the band ratio assumed by the signal model, for waste heat that follows the starlight. h_b = (eta_b / eta_1) * anchor(theta_e), where eta_b is the response of a background-limited profile fit to a de Vaucouleurs source with the galaxy's SDSS effective radius and axis ratio (data/psf_response_grid.csv) and the anchor corrects eta_1 to the W1 fraction that AllWISE actually returns for galaxies of that size (data/psf_w1_anchor.csv). This is the conservative ("anchored") version; see s09_image_injection.py.""" from scipy.interpolate import RegularGridInterpolator grid = pd.read_csv(RESP_GRID) res = np.array(sorted(grid.devrad.unique())) bas = np.array(sorted(grid.ba.unique())) devrad = d.theta_e.values / np.sqrt(d.ba.values) x = np.column_stack([np.log(np.clip(devrad, res[0], res[-1])), np.clip(d.ba.values, bas[0], bas[-1])]) eta = {} for b in (1, 2, 3, 4): tab = grid.pivot(index="devrad", columns="ba", values="eta%d" % b).loc[res, bas].values eta[b] = RegularGridInterpolator((np.log(res), bas), tab)(x) anc = pd.read_csv(os.path.join(M.DATA, "psf_w1_anchor.csv")) a = np.interp(np.log(d.theta_e.values), np.log(anc.theta_mid.values), anc.anchor.values) return {b: eta[b] / eta[1] * a for b in (2, 3, 4)} def limit(n_cand, n_test, recov_eff): if n_test <= 0 or not np.isfinite(recov_eff) or recov_eff <= 0: return np.nan return float(min(1.0, M.binom_ul(n_cand, n_test) / min(1.0, recov_eff))) def flags(d): """Catalogue evidence for conventional power sources / unreliable photometry (independent of f, T).""" cc = d.cc_flags out = pd.DataFrame(index=d.index) out["x_line"] = d.line_class.isin(["SF", "composite", "AGN", "nodata"]).values # ionised gas (or no data) with np.errstate(invalid="ignore"): out["x_agncol"] = ((d.w1mpro - d.w2mpro) >= AGNCOL).values out["x_uv"] = ((d.nuv_r < NUVR_SF) & (d.e_NUVmag < 0.3)).values # UV-excess: recent SF out["x_psb"] = (d.lick_hd_a >= 3.0).values # post-starburst out["x_radio"] = (d.first_n10.notna() | d.nvss_s14.notna()).values # FIRST (10") / NVSS (20") out["x_phot"] = ((cc.str[2] != "0") | (cc.str[3] != "0") | (d.n_allwise_6as != 1) | (d.nb != 1) | ~(d.sep_w <= 3.0)).values # artifact / blend / no match return out def clean_mask(fl, use_agncol=True): bad = fl.x_line | fl.x_uv | fl.x_psb | fl.x_radio | fl.x_phot if use_agncol: bad = bad | fl.x_agncol return ~bad.values def selection_factor(d, f): """Fraction of sample galaxies that would stay in the SDSS sample if dimmed by (1-f). Dividing a limit by this converts 'fraction of the observed sample' into 'fraction of galaxies of given pre-harvest luminosity'.""" dm = -2.5 * np.log10(1 - f) return float(np.mean(d.msel + dm <= SDSS_MLIM)), float(np.mean((d.msel + dm <= SDSS_MLIM) & (d.snr * (1 - f) >= 10))) def grid_table(eng, d, fl, fs, Ts, label="", extra_sel=None): rows, cands = [], {} for T in Ts: for f in fs: use_col = eng.predicted_w12(f, T) < 0.3 # the hypothesis itself must not predict AGN colours clean = clean_mask(fl, use_agncol=use_col) if extra_sel is not None: clean = clean & extra_sel sel_all = None if extra_sel is None else extra_sel r, cand, test = eng.run(f, T, sel_all) rc, cand_c, test_c = eng.run(f, T, clean) row = dict(label=label, T=T, f=f, n=r["n"], n_test=r["n_test"], sens_frac=r["sens_frac"], recov_eff=r["recov_eff"], n_cand=r["n_cand"], frac_compatible=r["n_cand"] / r["n_test"] if r["n_test"] else np.nan, ul_all=limit(r["n_cand"], r["n_test"], r["recov_eff"]), ul_zero=limit(0, r["n_test"], r["recov_eff"]), n_clean=rc["n"], n_test_clean=rc["n_test"], sens_frac_clean=rc["sens_frac"], n_cand_clean=rc["n_cand"], ul_clean=limit(rc["n_cand"], rc["n_test"], rc["recov_eff"]), ul_zero_clean=limit(0, rc["n_test"], rc["recov_eff"]), agncol_used=bool(use_col), n_null=eng.null_expectation(f, T, sel_all), n_null_clean=eng.null_expectation(f, T, clean)) for k in ("x_line", "x_agncol", "x_uv", "x_psb", "x_radio", "x_phot"): row["n_cand_" + k[2:]] = int((cand & fl[k].values).sum()) row["sel_mag"], row["sel_pess"] = selection_factor(d, f) rows.append(row) cands[(float(f), float(T))] = (np.where(cand)[0], np.where(cand_c)[0]) return pd.DataFrame(rows), cands SHOW = ["T", "f", "n_test", "sens_frac", "n_cand", "n_null", "ul_all", "ul_zero", "n_test_clean", "n_cand_clean", "n_null_clean", "ul_clean"] def main(): d = load_master() fl = flags(d) fs = np.unique(np.concatenate([M.F_GRID, M.F_FINE])) summary = {"K_S": K_S, "K_C": K_C, "EPS": EPS, "n_sample": len(d), "flag_counts": {k: int(fl[k].sum()) for k in fl.columns}, "n_clean_with_agncol": int(clean_mask(fl, True).sum()), "n_clean_without_agncol": int(clean_mask(fl, False).sum())} # ---------------- primary ---------------- eng = Engine(d) summary["n_valid_psf"] = int(eng.P["ok"].sum()) summary["n_valid_ls"] = int(phot(d, "ls")["ok"].sum()) tab, cands = grid_table(eng, d, fl, fs, M.T_GRID, label="primary") tab.to_csv(os.path.join(M.DATA, "limits_primary.csv"), index=False, float_format="%.5g") print(tab[tab.f.isin(M.F_GRID)][SHOW].to_string(index=False, float_format=lambda v: "%.3g" % v)) # membership table: every galaxy that is a candidate for any (f, T) of the brief grid rows = [] for (f, T), (idx, idx_c) in cands.items(): if np.isclose(M.F_GRID, f).any(): cset = set(idx_c.tolist()) rows += [(i, f, T, i in cset) for i in idx] cm = pd.DataFrame(rows, columns=["row", "f", "T", "clean"]) agg = cm.groupby("row").agg(f_max=("f", "max"), n_ft=("f", "size"), clean_any=("clean", "max")).reset_index() lst = cm.sort_values(["f", "T"]).groupby("row").apply( lambda x: ";".join("%.2f@%d%s" % (a, b, "c" if c else "") for a, b, c in zip(x.f, x["T"], x.clean)), include_groups=False) agg["ft_list"] = agg.row.map(lst) cd = d.loc[agg.row.values].copy() for c in ("f_max", "n_ft", "clean_any", "ft_list"): cd[c] = agg[c].values for c in fl.columns: cd[c] = fl.loc[agg.row.values, c].values cd = cd.sort_values(["f_max", "clean_any"], ascending=[False, False]) cd.to_csv(os.path.join(M.CACHE, "candidates_all.csv"), index=False, float_format="%.6g") summary["n_candidate_galaxies_any_ft"] = int(len(cd)) summary["n_candidate_galaxies_f>=0.1"] = int((cd.f_max >= 0.1).sum()) summary["n_clean_candidate_galaxies_f>=0.1"] = int(((cd.f_max >= 0.1) & cd.ft_list.str.contains(r"0\.[13]0@\d+c")).sum()) # ---------------- variants (brief grid plus f = 0.2, 0.25, 0.5) ---------------- variants = [] def add(label, eng_v, extra_sel=None, fl_v=fl, fgrid=M.F_VAR): t, _ = grid_table(eng_v, d, fl_v, fgrid, M.T_GRID, label=label, extra_sel=extra_sel) variants.append(t) add("primary", eng) add("allwise_and_legacy", Engine(d, use_ls=True)) add("legacy_w3_only", Engine(d, primary="ls", test_bands=(3,), comp_bands=(3,))) ls_ok = phot(d, "ls")["ok"] add("allwise_on_legacy_subsample", eng, extra_sel=ls_ok) # like-for-like for the two above add("allwise_w3_only_on_legacy_subsample", Engine(d, test_bands=(3,), comp_bands=(3,)), extra_sel=ls_ok) if os.path.exists(RESP_GRID): add("psf_response_inject", Engine(d, resp="inject")) add("psf_response_full", Engine(d, resp="full")) add("ks4", Engine(d, ks=4.0)) add("ks7", Engine(d, ks=7.0)) add("kc3", Engine(d, kc=3.0)) add("eps0", Engine(d, eps=0.0)) add("eps0.3", Engine(d, eps=0.3), fgrid=fs) # full fine grid: these feed the conservative curves add("eps0.4", Engine(d, eps=0.4), fgrid=fs) add("kbol_low", Engine(d, model="R_lo")) add("kbol_high", Engine(d, model="R_hi")) add("w3_only", Engine(d, test_bands=(3,), comp_bands=(3,))) add("w4_only", Engine(d, test_bands=(4,), comp_bands=(4,))) add("no_w2", Engine(d, comp_bands=(3, 4))) add("photosphere_subtracted", Engine(d, sub_photosphere=True), fgrid=fs) # W3 and W4 only # DIAGNOSTIC ONLY, not a valid limit: template photosphere also subtracted in W2 (see Engine.sub_bands) add("photosphere_subtracted_incl_w2_INVALID", Engine(d, sub_photosphere=True, sub_bands=(2, 3, 4))) add("point_like_ext0", eng, extra_sel=(d.ext_flg == 0).values) add("resolved_ext_gt0", eng, extra_sel=(d.ext_flg > 0).values) add("z_lt_0.06", eng, extra_sel=(d.z < 0.06).values) add("z_ge_0.06", eng, extra_sel=(d.z >= 0.06).values) add("strictly_lineless", eng, extra_sel=(d.line_class == "lineless").values) fl2 = fl.copy() fl2["x_radio"] = False add("clean_ignoring_radio", eng, fl_v=fl2) fl3 = fl.copy() fl3["x_uv"] = False add("clean_ignoring_uv", eng, fl_v=fl3) var = pd.concat(variants, ignore_index=True) var.to_csv(os.path.join(M.DATA, "limits_variants.csv"), index=False, float_format="%.5g") json.dump(summary, open(os.path.join(M.DATA, "limits_summary.json"), "w"), indent=1) print(json.dumps(summary, indent=1)) if __name__ == "__main__": main()