#!/usr/bin/env python3 """Individual-galaxy search: FP / Faber-Jackson fits, residual tails, selection-corrected upper limits. For a galaxy whose starlight is dimmed by a factor (1-f) at fixed sigma and R_e: dI = log I_e - (c0 + c1 log sigma + c2 log R_e + c3 z) shifts by log10(1-f) dFJ = log L - (d0 + d1 log sigma + d2 z) shifts by log10(1-f) (the classical FP residual in log R_e shifts by -b*log10(1-f); it carries the same information). Upper limit for a population in which a fraction q of early-type galaxies is dimmed by exactly f: n_obs(f) = number of sample galaxies with residual <= t_f, t_f = log10(1-f) (a-priori rule) M(f) = sum over sample galaxies i of survive_i(f) * P(dimmed twin of i has residual <= t_f) q_95(f) = PoissonUL95(n_obs) / M(f) Every observed outlier is treated as a potential signal (no background subtraction) -> conservative. "survive" = the dimmed twin would still satisfy every flux-dependent cut that defines the sample. """ import json import os import numpy as np import pandas as pd from scipy import optimize, special, stats import fpk_common as F F_GRID = F.F_GRID SDSS_SAFE_MLIM = 17.60 # conservative MGS limit applied to dimmed twins (17.77 nominal; 17.60 in early chunks) SDSS_MU50_LIM = 23.0 # MGS half-light surface-brightness limit below which targeting was complete SNR_MIN = {"sdss": 10.0, "6df": 5.0} R_MIN_6DF = 8.0 KAPPAS = (-1.0, 0.0, 1.0, 2.0) # pre-specified thresholds t = log10(1-f) + kappa * MAD SIGMA_AGREE = 0.06 # dex, SDSS quality cut: pipeline vs Portsmouth dispersion ELOGS_MAX_Q = 0.0434 # dex (10 per cent), SDSS quality cut # ---------------------------------------------------------------------------- # empirical spectral S/N versus fibre magnitude (used to predict the S/N of a dimmed twin) # ---------------------------------------------------------------------------- def snr_model(m, K, msky, smax): fl = 10 ** (-0.4 * (m - msky)) s = K * fl / np.sqrt(fl + 1.0) return 1.0 / np.sqrt(1.0 / s ** 2 + 1.0 / smax ** 2) def fit_snr_model(mfib, snr, lo, hi, step, p0): bins = np.arange(lo, hi + 1e-9, step) xs, ys = [], [] for a, b in zip(bins[:-1], bins[1:]): m = (mfib >= a) & (mfib < b) if m.sum() >= 30: xs.append(0.5 * (a + b)) ys.append(np.median(snr[m])) p, _ = optimize.curve_fit(snr_model, np.array(xs), np.array(ys), p0=p0, maxfev=40000) return [float(v) for v in p] def fibre_mag_6df(df): """J magnitude inside the 6.7-arcsec 6dF fibre from the catalogued Sersic index and effective radius.""" n = df.sersic_n.values bn = 2 * n - 1.0 / 3 + 4.0 / (405 * n) frac = special.gammainc(2 * n, bn * (3.35 / df.theta_e.values) ** (1.0 / n)) return df.msel.values - 2.5 * np.log10(frac) # ---------------------------------------------------------------------------- # would galaxy i still be in the sample if dimmed by (1-f)? returns (mask, snr_ratio g) # ---------------------------------------------------------------------------- def snr_ratio(df, f, cfg, variant): dm = -2.5 * np.log10(1 - f) if variant == "snr_linear": # sky-limited worst case: S/N proportional to flux return (1 - f) * np.ones(len(df)) p = cfg["snr_par"] return snr_model(df.mfib.values + dm, *p) / snr_model(df.mfib.values, *p) def survive(df, f, cfg, variant="default", quality=False): name = cfg["name"] dm = -2.5 * np.log10(1 - f) g = snr_ratio(df, f, cfg, variant) ok = df.mu_obs.values + dm <= cfg["mu_lim"] # stay inside the surface-brightness range actually sampled if name == "6df": ok &= df.msel.values + dm <= df.mlim.values ok &= df.snr.values * g >= SNR_MIN["6df"] ok &= df.Rtd.values * g ** 0.54 >= R_MIN_6DF # empirical R ~ (S/N)^0.54 in the 6dFGSv table elif name == "sdss": mlim = 17.77 if variant == "mlim1777" else SDSS_SAFE_MLIM ok &= df.msel.values + dm <= mlim ok &= df.mu50.values + dm <= SDSS_MU50_LIM ok &= df.snr.values * g >= SNR_MIN["sdss"] # the colour window of the sample definition slides with absolute magnitude mr = df.Mr.values + dm ok &= df.gr.values >= 0.63 - 0.02 * (mr + 20) if variant != "no_red_edge": # "no_red_edge": sample defined without the red colour limit, ok &= df.gr.values <= 1.03 - 0.02 * (mr + 20) # so that dimming which also reddens is not excluded if variant == "ew_boost": # emission lines not dimmed: EW grows as 1/(1-f); demand EW<3A ew = np.nan_to_num(df.ew_ha_port.values, nan=0.0) ok &= (ew / (1 - f)) < 3.0 if quality: ok &= np.abs(df.logs.values - df.logs_port.values) / g < SIGMA_AGREE ok &= df.e_logs.values / g < ELOGS_MAX_Q elif name == "howlett": ok &= df.msel.values + dm <= df.mlim.values return ok, g # ---------------------------------------------------------------------------- def fit_relations(df, verbose=False): """Truncated-regression fits of the I-direction FP and of Faber-Jackson; classical FP for reference.""" s0, r0, z0 = 2.25, 0.5, float(np.median(df.zdist)) y = df.logi.values X = np.column_stack([np.ones(len(df)), df.logs.values - s0, df.logr.values - r0, df.zdist.values - z0]) ylim = y - 0.4 * (df.mlim.values - df.msel.values) pre = F.fit_truncated_plane(y, X, ylim=ylim) cs = pre["coef"][1] yerr2 = (cs * df.e_logs.values) ** 2 fit = F.fit_truncated_plane(y, X, ylim=ylim, yerr2=yerr2, verbose=verbose) logL = df.logi.values + 2 * df.logr.values + np.log10(2 * np.pi) + 6.0 XJ = np.column_stack([np.ones(len(df)), df.logs.values - s0, df.zdist.values - z0]) ylimJ = logL - 0.4 * (df.mlim.values - df.msel.values) preJ = F.fit_truncated_plane(logL, XJ, ylim=ylimJ) fj = F.fit_truncated_plane(logL, XJ, ylim=ylimJ, yerr2=(preJ["coef"][1] * df.e_logs.values) ** 2) fit_nt = F.fit_truncated_plane(y, X, ylim=None, yerr2=yerr2) classic = F.weighted_orthogonal_fp(df.logr.values, df.logs.values, df.logi.values) info = dict( n=int(len(df)), s0=s0, r0=r0, z0=z0, fpI_coef=[float(v) for v in fit["coef"]], fpI_coef_err=[float(v) for v in fit["coef_err"]], fpI_s_int=fit["s_int"], fpI_rms_clipped=float(np.std(fit["resid"][fit["fitmask"]])), fpI_mad=float(F.mad_sigma(fit["resid"])), fpI_nfit=int(fit["fitmask"].sum()), fpI_coef_untruncated=[float(v) for v in fit_nt["coef"]], fj_coef=[float(v) for v in fj["coef"]], fj_s_int=fj["s_int"], fj_mad=float(F.mad_sigma(fj["resid"])), classic_direct={a: float(b) for a, b in classic["direct"].items()}, classic_orth={a: float(b) for a, b in classic["orthogonal"].items()}, ) res = pd.DataFrame({ "dI": fit["resid"], "dI_sd": fit["sd"], "dFJ": fj["resid"], "dFJ_sd": fj["sd"], "dR_orth": classic["resid_r_orth"], "dR_direct": classic["resid_r_direct"], }, index=df.index) return info, res, float(cs), float(fj["coef"][1]) def limits(df, res, cs, cj, cfg, label, variant="default", kappa=0.0, quality=False, f_grid=F_GRID, n_obs_override=None): """Counting upper limits for each dimming fraction f. kappa moves the threshold: t = shift + kappa*MAD.""" rows = [] for rel, col, c_sig in (("FP", "dI", cs), ("FJ", "dFJ", cj)): d = res[col].values mad = F.mad_sigma(d) for f in f_grid: shift = np.log10(1 - f) t = shift + kappa * mad ok, g = survive(df, f, cfg, variant, quality) # extra velocity-dispersion noise of the dimmed twin (spectrum S/N lower by the factor g) extra = np.abs(c_sig) * df.e_logs.values * np.sqrt(np.clip(1.0 / np.maximum(g, 1e-3) ** 2 - 1.0, 0, None)) with np.errstate(divide="ignore", invalid="ignore"): p = np.where(extra > 1e-6, stats.norm.cdf((t - shift - d) / np.maximum(extra, 1e-6)), (d + shift <= t).astype(float)) M = float(np.sum(p[ok])) n_obs = int(np.sum(d <= t)) if n_obs_override is not None and rel == "FP": n_obs = int(n_obs_override(t)) ul = float(F.poisson_ul(n_obs)) ulb = float(F.poisson_ul(n_obs, 1 - 0.05 / len(KAPPAS))) rows.append(dict(sample=label, variant=variant, relation=rel, kappa=kappa, f=float(f), shift_dex=shift, threshold_dex=t, N=int(len(df)), N_eff=int(ok.sum()), M_eff=M, n_obs=n_obs, n_mirror=int(np.sum(d >= -t)), ul95_n=ul, q95=ul / M if M > 0 else np.nan, q_bonf=ulb / M if M > 0 else np.nan)) return pd.DataFrame(rows) def best_of_thresholds(lim): """Valid 95% limit when the threshold is chosen among the len(KAPPAS) pre-specified ones: each threshold is evaluated at confidence 1-0.05/len(KAPPAS) (Bonferroni) and the smallest limit is kept.""" keys = ["sample", "variant", "relation", "f"] lim = lim[lim.kappa.isin(KAPPAS)] lim = lim.assign(_q=lim.q_bonf.fillna(np.inf)).reset_index(drop=True) idx = lim.groupby(keys)["_q"].idxmin() best = lim.loc[idx, keys + ["kappa", "threshold_dex", "N", "N_eff", "M_eff", "n_obs", "q_bonf"]] return best.rename(columns={"q_bonf": "q95_best"}).reset_index(drop=True) def sweep(df, res, cs, cj, cfg, label, **kw): return pd.concat([limits(df, res, cs, cj, cfg, label, kappa=k, **kw) for k in KAPPAS], ignore_index=True) def sdss_quality_mask(ds): """Data-quality cuts for the SDSS sample (all applied to the dimmed twins as well): * two independent velocity-dispersion codes (SDSS pipeline, Portsmouth pPXF) agree within 0.06 dex * 'clean' photometry flag; axis ratio >= 0.5 (removes inclined discs with rotation-boosted dispersions) * de Vaucouleurs and Petrosian magnitudes consistent (|deV - Petro + 0.08| < 0.3 mag) * velocity-dispersion error < 10 per cent """ with np.errstate(invalid="ignore"): agree = np.abs(ds.logs.values - ds.logs_port.values) < SIGMA_AGREE q = agree & (ds.clean.values == 1) & (ds.ba.values >= 0.5) q &= np.abs(ds.dmag_dev_petro.values + 0.08) < 0.3 q &= ds.e_logs.values < ELOGS_MAX_Q return q def load_samples(): """Read the three standardised samples and attach survey configuration.""" out = {} d6 = pd.read_csv(os.path.join(F.DATA, "sample_6df.csv")) d6 = d6[d6.in_sample == 1].reset_index(drop=True) d6["mfib"] = fibre_mag_6df(d6) cfg6 = dict(name="6df", mu_lim=float(np.percentile(d6.mu_obs, 99.0)), snr_par=fit_snr_model(d6.mfib.values, d6.snr.values, 11.5, 14.7, 0.4, [3.0, 16.5, 200.0])) out["6df"] = (d6, cfg6) ds = pd.read_csv(F.sample_path("sdss")) ds["mfib"] = ds.fibermag_r cfgs = dict(name="sdss", mu_lim=float(np.percentile(ds.mu_obs, 99.0)), snr_par=fit_snr_model(ds.mfib.values, ds.snr.values, 15.25, 19.5, 0.25, [20.0, 18.5, 60.0])) out["sdss"] = (ds, cfgs) dh = pd.read_csv(os.path.join(F.DATA, "sample_howlett.csv")) dh["e_logs"] = dh.e_logs.fillna(0.0) dh["mfib"] = dh.msel cfgh = dict(name="howlett", mu_lim=float(np.percentile(dh.mu_obs, 99.0)), snr_par=[1.0, 30.0, 1e9]) out["howlett"] = (dh, cfgh) return out def main(): S = load_samples() all_info, all_lim = {}, [] # ------------------------------------------------------------------ 6dFGSv (J band) d6, cfg6 = S["6df"] info, res6, cs, cj = fit_relations(d6) info["snr_par"] = cfg6["snr_par"] info["mu_lim"] = cfg6["mu_lim"] all_info["6df"] = info for v in ("default", "snr_linear"): all_lim.append(sweep(d6, res6, cs, cj, cfg6, "6dFGSv-J (unclipped)", variant=v)) r6 = pd.concat([d6, res6], axis=1) r6.to_csv(os.path.join(F.DATA, "resid_6df.csv"), index=False, float_format="%.8g") q = (d6.Mt.values <= 2.0) all_lim.append(sweep(d6[q].reset_index(drop=True), res6[q].reset_index(drop=True), cs, cj, cfg6, "6dFGSv-J E/S0 only (Mt<=2)")) c = (d6.in_6dFGSv_J.values == 1) all_lim.append(sweep(d6[c].reset_index(drop=True), res6[c].reset_index(drop=True), cs, cj, cfg6, "6dFGSv-J as published (chi2-clipped)")) # ------------------------------------------------------------------ Howlett+2022 (clipped; cross-check) dh, cfgh = S["howlett"] info, resh, csh, cjh = fit_relations(dh) all_info["howlett"] = info all_lim.append(sweep(dh, resh, csh, cjh, cfgh, "SDSS-PV Howlett+22 (clipped)")) pd.concat([dh, resh], axis=1).to_csv(os.path.join(F.DATA, "resid_howlett.csv"), index=False, float_format="%.8g") # ------------------------------------------------------------------ SDSS own (unclipped) ds, cfgs = S["sdss"] info, ress, css, cjs = fit_relations(ds) info["snr_par"] = cfgs["snr_par"] info["mu_lim"] = cfgs["mu_lim"] all_info["sdss"] = info for v in ("default", "snr_linear", "mlim1777", "ew_boost"): all_lim.append(sweep(ds, ress, css, cjs, cfgs, "SDSS-MGS own (unclipped)", variant=v)) qm = sdss_quality_mask(ds) rs = pd.concat([ds, ress], axis=1) rs["quality"] = qm.astype(int) rs.to_csv(F.resid_path("sdss"), index=False, float_format="%.8g") slim = rs[["id", "plate", "ra", "dec", "zdist", "logs", "e_logs", "logs_port", "theta_e", "logr", "logi", "msel", "snr", "fibermag_r", "gr", "ba", "d4000", "ew_ha_port", "clean", "in_howlett", "quality", "dI", "dI_sd", "dFJ"]].copy() for c in slim.columns: if c in ("ra", "dec", "zdist"): slim[c] = slim[c].round(5) elif c in ("msel", "mu50", "fibermag_r", "Mr", "gr", "d4000", "theta_e", "ba", "ew_ha_port"): slim[c] = slim[c].round(3) elif c == "snr": slim[c] = slim[c].round(1) elif slim[c].dtype.kind == "f": slim[c] = slim[c].round(4) slim.to_csv(os.path.join(F.DATA, "resid_sdss.csv"), index=False) all_lim.append(sweep(ds[qm].reset_index(drop=True), ress[qm].reset_index(drop=True), css, cjs, cfgs, "SDSS-MGS own, quality cuts", quality=True)) # sub-samples: passive (Portsmouth EW(Ha)<1A, as in Howlett+22) and galaxies present in Howlett+22 with np.errstate(invalid="ignore"): pm = (ds.ew_ha_port.values < 1.0) all_lim.append(sweep(ds[pm].reset_index(drop=True), ress[pm].reset_index(drop=True), css, cjs, cfgs, "SDSS-MGS own, EW(Ha)<1A", variant="ew_boost")) lim = pd.concat(all_lim, ignore_index=True) lim.to_csv(os.path.join(F.DATA, "limits_individual.csv"), index=False, float_format="%.6g") with open(os.path.join(F.DATA, "fp_fits.json"), "w") as fh: json.dump(all_info, fh, indent=1) # clipping diagnostics: what the published catalogues removed diag = [] for thr in (-0.30, -0.45, -0.60): diag.append(dict(threshold=thr, sdss_own_all=int((rs.dI <= thr).sum()), sdss_own_N=len(rs), sdss_own_in_howlett=int(((rs.dI <= thr) & (rs.in_howlett == 1)).sum()), sdss_own_in_howlett_N=int((rs.in_howlett == 1).sum()), sixdf_all=int((r6.dI <= thr).sum()), sixdf_N=len(r6), sixdf_published=int(((r6.dI <= thr) & (r6.in_6dFGSv_J == 1)).sum()), sixdf_published_N=int((r6.in_6dFGSv_J == 1).sum()))) pd.DataFrame(diag).to_csv(os.path.join(F.DATA, "clipping_diagnostics.csv"), index=False) best = best_of_thresholds(lim) best.to_csv(os.path.join(F.DATA, "limits_individual_best.csv"), index=False, float_format="%.6g") pd.set_option("display.width", 250) show = best[best.f.isin([0.5, 0.75, 0.9]) & (best.relation == "FP")] print(show.to_string(index=False)) show = best[best.f.isin([0.5, 0.75, 0.9]) & (best.relation == "FJ") & (best.variant == "default")] print(show.to_string(index=False)) print(pd.DataFrame(diag).to_string(index=False)) return all_info, lim if __name__ == "__main__": main()