#!/usr/bin/env python3 """Checks added in the revision (each answers a point of the referee report; see report.md, last section). 1. red_side(): galaxies excluded by the RED edge of the colour window (non-grey dimming of an old population reddens it, Lacki 2019), and limits for a sample defined without that edge. 2. mixture_bounds(): upper limit on the TOTAL fraction of dimmed galaxies when the population is spread over several values of f (linear programme over all pre-specified thresholds). 3. howlett_eligibility(): how many far-tail galaxies were eligible for the Howlett+22 catalogue at all, and how many would be expected in it without outlier rejection. 4. neff_population(): which galaxies make up N_eff(f) (the population each limit actually describes). 5. tail_asymmetry(): numbers of galaxies beyond +-0.30 / 0.45 / 0.60 dex on the dim and on the bright side. 6. cluster_decomposition(): the Coma sphere: finger-of-god, photometry, remainder versus the like-for-like (group-redshift) environmental trend; spheres centred on catalogue cluster positions. patch_identifications(): catalogued clusters inside the most significant dim spheres (SIMBAD). 7. wise_apertures(): AllWISE aperture photometry (IRSA) of the far-tail galaxies, because the profile-fit magnitudes used in vet_candidates.py underestimate the flux of resolved galaxies. 8. region_final(): the region limits with the tested share of the survey volume, the number of sphere volumes tested, Monte Carlo errors, redshift reach and galaxy-weighted limits. """ import io import json import os import time import urllib.parse import urllib.request import numpy as np import pandas as pd from scipy import optimize, stats from scipy.spatial import cKDTree import analysis_individual as A import build_samples as B import fpk_common as F import vet_candidates as V SIMBAD = "https://simbad.cds.unistra.fr/simbad/sim-tap/sync" IRSA = "https://irsa.ipac.caltech.edu/TAP/sync" CLUSTERS = ("ACO 1656", "ACO 2593", "ACO 2142") # Coma and the two other clusters named in the report F_MIX = np.array([0.5, 0.6, 0.7, 0.75, 0.8, 0.85, 0.9]) TAIL_THR = (0.30, 0.45, 0.60) LABEL_RED = "SDSS-MGS own, no red colour edge" def _get(url, params, dest, tries=5, timeout=300): if os.path.exists(dest) and os.path.getsize(dest) > 50: return open(dest).read() last = None for k in range(tries): try: with urllib.request.urlopen(url + "?" + urllib.parse.urlencode(params), timeout=timeout) as r: txt = r.read().decode() if txt.lstrip().startswith("<"): raise RuntimeError("service error: %s" % txt[:300]) with open(dest, "w") as fh: fh.write(txt) return txt except Exception as e: # noqa last = e time.sleep(10 * (k + 1)) raise RuntimeError("failed to fetch %s: %s" % (url, last)) def _plane_residuals(df, fit): """I-direction FP and Faber-Jackson residuals of any galaxy table with respect to the BASELINE fits.""" c, s0, r0, z0 = fit["fpI_coef"], fit["s0"], fit["r0"], fit["z0"] dI = df.logi.values - (c[0] + c[1] * (df.logs.values - s0) + c[2] * (df.logr.values - r0) + c[3] * (df.zdist.values - z0)) j = fit["fj_coef"] logL = df.logi.values + 2 * df.logr.values + np.log10(2 * np.pi) + 6.0 dFJ = logL - (j[0] + j[1] * (df.logs.values - s0) + j[2] * (df.zdist.values - z0)) return pd.DataFrame({"dI": dI, "dFJ": dFJ}, index=df.index) # ---------------------------------------------------------------------------------------------------------- # 1. the red side of the colour window # ---------------------------------------------------------------------------------------------------------- def red_side(S, fits): P = pd.read_csv(os.path.join(F.CACHE, "sdss_own_parent_derived.csv")) cut = B.sdss_cut_masks(P) lo, hi = B.sdss_colour_edges(P.Mr) res = _plane_residuals(P, fits["sdss"]) P["dI"], P["dFJ"] = res.dI.values, res.dFJ.values P["excess_red"] = P.gr - hi base = np.ones(len(P), bool) for k in B.SDSS_CUT_ORDER: if k not in ("colour", "d4000"): base &= cut[k].values d4 = cut["d4000"].values red = P.gr.values > hi.values blue = P.gr.values < lo.values inwin = cut["colour"].values rows = [] for lab, m in (("in colour window, D4000>=1.6 (= the sample)", base & inwin & d4), ("redder than window, D4000 cut not applied", base & red), ("redder than window, D4000>=1.6", base & red & d4), ("bluer than window, D4000 cut not applied", base & blue), ("in window, D4000<1.6 or missing", base & inwin & ~d4)): r = dict(subset=lab, N=int(m.sum()), median_dI=float(np.nanmedian(P.dI.values[m]))) for t in TAIL_THR: below = m & (P.dI.values <= -t) r["n_dI_le_-%.2f" % t] = int(below.sum()) r["n_dI_le_-%.2f_and_excess_gt_0.2mag" % t] = int((below & (P.excess_red.values > 0.2)).sum()) r["n_dI_le_-%.2f_and_excess_gt_0.5mag" % t] = int((below & (P.excess_red.values > 0.5)).sum()) rows.append(r) tab = pd.DataFrame(rows) tab.to_csv(os.path.join(F.DATA, "redside_check.csv"), index=False, float_format="%.5g") obj = P[base & red].sort_values("dI") cols = ["id", "ra", "dec", "z", "gr", "excess_red", "Mr", "dI", "veldisp", "theta_e", "msel", "snr", "d4000", "ba", "clean", "dmag_dev_petro", "ew_ha_port"] obj = obj[cols].copy() # stored visual notes for the objects that are both very red (> 0.2 mag beyond the edge) and >= 50% deficient man = os.path.join(F.HERE, "redside_manual.csv") if os.path.exists(man): obj = obj.merge(pd.read_csv(man)[["id", "visual_class", "visual_note"]], on="id", how="left") obj.to_csv(os.path.join(F.DATA, "redside_objects.csv"), index=False, float_format="%.6g") exc = P.excess_red.values[base & red] summ = dict(n_pass_noncolour_cuts_noD4000=int(base.sum()), n_red=int((base & red).sum()), n_red_d4000=int((base & red & d4).sum()), red_excess_p50=float(np.percentile(exc, 50)), red_excess_p90=float(np.percentile(exc, 90)), red_excess_max=float(exc.max()), sample_dist_to_red_edge_median=float(np.median(-P.excess_red.values[base & inwin & d4])), sample_dist_to_red_edge_p05=float(np.percentile(-P.excess_red.values[base & inwin & d4], 5))) # ---- limits for the "red-extended" sample: every cut of the sample except the red edge of the colour window. # f is then the fraction of r-band light removed, whatever happens to the colour (as long as the galaxy does # not become bluer than the blue edge). Residuals are taken with respect to the baseline plane. ext = base & d4 & cut["colour_blue_edge_only"].values dr = P[ext].reset_index(drop=True) dr["mfib"] = dr.fibermag_r rr = dr[["dI", "dFJ"]].copy() cfg = S["sdss"][1] cs, cj = fits["sdss"]["fpI_coef"][1], fits["sdss"]["fj_coef"][1] lim = A.sweep(dr, rr, cs, cj, cfg, LABEL_RED, variant="no_red_edge") best = A.best_of_thresholds(lim) best.to_csv(os.path.join(F.DATA, "limits_redext.csv"), index=False, float_format="%.6g") summ["n_redext_sample"] = int(ext.sum()) b = best[best.relation == "FP"].set_index("f") for f in (0.3, 0.4, 0.5, 0.6, 0.75, 0.9): summ["redext_q95_f%.2f" % f] = float(b.loc[f, "q95_best"]) summ["redext_Neff_f%.2f" % f] = int(b.loc[f, "N_eff"]) summ["redext_nobs_f%.2f" % f] = int(b.loc[f, "n_obs"]) print(tab.to_string(index=False)) print(best[(best.relation == "FP")].to_string(index=False)) return summ # ---------------------------------------------------------------------------------------------------------- # 2. mixtures of f # ---------------------------------------------------------------------------------------------------------- def _m_matrix(df, d, cs, cfg, fg, thr): """M[k, j] = expected number of sample galaxies that, dimmed by fg[j], survive and fall below thr[k].""" M = np.zeros((len(thr), len(fg))) for j, f in enumerate(fg): ok, g = A.survive(df, f, cfg) shift = np.log10(1 - f) extra = np.abs(cs) * df.e_logs.values * np.sqrt(np.clip(1.0 / np.maximum(g, 1e-3) ** 2 - 1.0, 0, None)) for k, t in enumerate(thr): 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[k, j] = p[ok].sum() return M def mixture_bounds(S, fits): """If a fraction q_j of galaxies is dimmed by f_j, the expected count below threshold t_k is at least sum_j q_j M[k, j]. With an upper limit U_k on each expected count, the largest total fraction sum_j q_j compatible with all thresholds is a linear programme. Two confidence conventions: 'per-threshold' U_k at 1 - 0.05/4 (the Bonferroni level used for the single-f limits); 'simultaneous' U_k at 1 - 0.05/K for all K thresholds at once (a rigorous 95% bound for any mixture).""" rows = [] for name, label in (("sdss", "SDSS-MGS own (unclipped)"), ("6df", "6dFGSv-J (unclipped)")): df, cfg = S[name] d = pd.read_csv(F.resid_path(name)).dI.values cs = fits[name]["fpI_coef"][1] mad = F.mad_sigma(d) thr = np.array(sorted(set(np.round([np.log10(1 - f) + k * mad for f in F_MIX for k in A.KAPPAS], 6)))) M = _m_matrix(df, d, cs, cfg, F_MIX, thr) nobs = np.array([(d <= t).sum() for t in thr]) for conv, cl in (("per-threshold", 1 - 0.05 / len(A.KAPPAS)), ("simultaneous", 1 - 0.05 / len(thr))): ul = F.poisson_ul(nobs, cl) for lo_f, hi_f in ((0.5, 0.9), (0.75, 0.9)): idx = (F_MIX >= lo_f - 1e-9) & (F_MIX <= hi_f + 1e-9) r = optimize.linprog(-np.ones(idx.sum()), A_ub=M[:, idx], b_ub=ul, bounds=[(0, None)] * int(idx.sum()), method="highs") single = [np.min(ul[M[:, j] > 0] / M[M[:, j] > 0, j]) if (M[:, j] > 0).any() else np.nan for j in np.where(idx)[0]] rows.append(dict(sample=label, convention=conv, confidence_per_threshold=cl, n_thresholds=len(thr), f_lo=lo_f, f_hi=hi_f, max_total_fraction=float(-r.fun), worst_single_f_limit=float(np.nanmax(single)), composition=";".join("f=%.2f:%.5f" % (f, x) for f, x in zip(F_MIX[idx], r.x) if x > 1e-7))) out = pd.DataFrame(rows) out.to_csv(os.path.join(F.DATA, "limits_mixture.csv"), index=False, float_format="%.6g") print(out.to_string(index=False)) return out # ---------------------------------------------------------------------------------------------------------- # 3. what the Howlett+22 catalogue could have contained # ---------------------------------------------------------------------------------------------------------- def howlett_eligibility(): """Howlett+22 require 10 <= m_r(deV) <= 17.0, z_CMB <= 0.1 and EW(Halpha) < 1 A (plus cuts similar to ours, a visual inspection that removed 20 per cent, and the FP-outlier rejection). 'Eligible' = passes those three cuts; the inclusion rate of eligible galaxies of our sample in their catalogue then gives the number of far-tail galaxies expected in it if the tail were treated like the rest of the sample.""" rs = pd.read_csv(F.resid_path("sdss")) with np.errstate(invalid="ignore"): elig = ((rs.devmag_r.values - rs.ext_r.values) <= 17.0) & (rs.z.values <= 0.1) & \ (rs.ew_ha_port.values < 1.0) rate = float(rs.in_howlett.values[elig].mean()) rows = [] for t in TAIL_THR: m = rs.dI.values <= -t n_el = int((m & elig).sum()) obs = int((m & (rs.in_howlett.values == 1)).sum()) rows.append(dict(threshold_dex=-t, n_tail=int(m.sum()), n_tail_eligible=n_el, inclusion_rate_of_eligible=rate, n_eligible_all=int(elig.sum()), expected_in_howlett=n_el * rate, observed_in_howlett=obs, p_binomial_le_observed=float(stats.binom.cdf(obs, n_el, rate)))) out = pd.DataFrame(rows) out.to_csv(os.path.join(F.DATA, "howlett_eligibility.csv"), index=False, float_format="%.5g") print(out.to_string(index=False)) return out # ---------------------------------------------------------------------------------------------------------- # 4. the population behind N_eff(f) # ---------------------------------------------------------------------------------------------------------- def neff_population(S): rows = [] for name, label in (("sdss", "SDSS"), ("6df", "6dFGSv")): df, cfg = S[name] d = pd.read_csv(F.resid_path(name)).dI.values for f in (0.3, 0.5, 0.75, 0.9): ok, _ = A.survive(df, f, cfg) okl, _ = A.survive(df, f, cfg, variant="snr_linear") if ok.sum() == 0: continue r = dict(sample=label, f=f, N=len(df), N_eff=int(ok.sum()), share_of_sample=float(ok.mean()), N_eff_snr_linear=int(okl.sum()), z_median=float(np.median(df.zdist.values[ok])), z_p95=float(np.percentile(df.zdist.values[ok], 95)), z_median_sample=float(np.median(df.zdist.values)), mag_faintest=float(df.msel.values[ok].max()), mag_median=float(np.median(df.msel.values[ok])), mag_median_sample=float(np.median(df.msel.values)), mean_dI=float(d[ok].mean()), mean_dI_sample=float(d.mean()), logsigma_median=float(np.median(df.logs.values[ok])), logsigma_median_sample=float(np.median(df.logs.values))) if "Mr" in df: r.update(Mr_median=float(np.median(df.Mr.values[ok])), Mr_median_sample=float(np.median(df.Mr.values))) rows.append(r) out = pd.DataFrame(rows) out.to_csv(os.path.join(F.DATA, "neff_population.csv"), index=False, float_format="%.5g") print(out.to_string(index=False)) return out # ---------------------------------------------------------------------------------------------------------- # 5. dim versus bright tail # ---------------------------------------------------------------------------------------------------------- def tail_asymmetry(): rows = [] for name, label in (("sdss", "SDSS"), ("6df", "6dFGSv")): d = pd.read_csv(F.resid_path(name)).dI.values for t in TAIL_THR: nd, nb = int((d <= -t).sum()), int((d >= t).sum()) rows.append(dict(sample=label, N=len(d), beyond_dex=t, n_dim=nd, n_bright=nb, rate_dim=nd / len(d), rate_bright=nb / len(d), p_binomial_dim_excess=float(stats.binomtest(nd, nd + nb, 0.5, alternative="greater").pvalue) if nd + nb > 0 else np.nan)) out = pd.DataFrame(rows) both = out.groupby("beyond_dex")[["n_dim", "n_bright"]].sum().reset_index() for _, r in both.iterrows(): rows.append(dict(sample="both", N=int(out.N.unique().sum()), beyond_dex=r.beyond_dex, n_dim=int(r.n_dim), n_bright=int(r.n_bright), rate_dim=np.nan, rate_bright=np.nan, p_binomial_dim_excess=float(stats.binomtest(int(r.n_dim), int(r.n_dim + r.n_bright), 0.5, alternative="greater").pvalue))) out = pd.DataFrame(rows) out.to_csv(os.path.join(F.DATA, "tail_asymmetry.csv"), index=False, float_format="%.5g") print(out.to_string(index=False)) # the objects at an apparent deficit of 75 per cent or more (dI <= -0.60) and their vetting status far = [] for name in ("sdss", "6df"): c = pd.read_csv(os.path.join(F.DATA, "candidates_%s_vetted.csv" % name), dtype={"id": str}) c = c[c.dI <= -0.60] for _, r in c.iterrows(): far.append(dict(label=r.label, id=r.id, ra=r.ra, dec=r["dec"], z=r.zdist, dI=r.dI, apparent_f=1 - 10 ** r.dI, explained=int(r.explained), vetting_class=r["class"], note=r.note)) far = pd.DataFrame(far) far.to_csv(os.path.join(F.DATA, "far_tail_objects.csv"), index=False, float_format="%.7g") print(far[["label", "dI", "apparent_f", "explained", "vetting_class"]].to_string(index=False)) return out, far # ---------------------------------------------------------------------------------------------------------- # 6. clusters # ---------------------------------------------------------------------------------------------------------- def simbad_clusters(): q = ("SELECT id, main_id, ra, dec, rvz_redshift FROM basic JOIN ident ON oidref=oid WHERE id IN (%s)" % ",".join("'%s'" % c for c in CLUSTERS)) txt = _get(SIMBAD, dict(REQUEST="doQuery", LANG="ADQL", FORMAT="csv", QUERY=q), os.path.join(F.RAW, "simbad_clusters.csv")) t = pd.read_csv(io.StringIO(txt)) t["id"] = t["id"].str.replace(r"\s+", " ", regex=True).str.strip() return t def cluster_decomposition(): s = pd.read_csv(F.resid_path("sdss")) s = s[(s.zdist >= 0.02) & (s.zdist <= 0.12)].reset_index(drop=True) xyz = F.radec_to_xyz(s.ra.values, s["dec"].values, F.comoving_distance(s.zdist.values)) tree = cKDTree(xyz) pt = pd.read_csv(os.path.join(F.DATA, "region_extreme_patches.csv")) fog = pd.read_csv(os.path.join(F.DATA, "region_patches_groupz.csv")) env = pd.read_csv(os.path.join(F.DATA, "env_trend.csv")) top = pt[(pt["sample"] == "sdss") & (pt.variant == "raw") & (pt.R_mpc == 10) & pt.kind.str.contains("negative")].iloc[0] f10 = fog[fog.R_mpc == 10].iloc[0] c = F.radec_to_xyz(np.array([top.ra]), np.array([top["dec"]]), F.comoving_distance(np.array([top.z])))[0] nb = np.array(tree.query_ball_point(c, 10.0), dtype=int) sem_sphere = float(s.dI.values[nb].std() / np.sqrt(len(nb))) fog_part = float(f10.mean_dI_ownz_subset - f10.mean_dI_groupz_subset) phot_part = float(-0.4 * top.mean_r_minus_J_resid_mag) phot_err = float(0.4 * top.sem_r_minus_J_resid_mag) rem = float(top.mean_dI_dex - fog_part - phot_part) hi_env = env[env.overdensity_lo == 8] tr_ind = hi_env[hi_env["sample"].str.contains("individual")].iloc[0] tr_grp = hi_env[hi_env["sample"].str.contains("Howlett")].iloc[0] rem_err = float(np.hypot(np.hypot(sem_sphere, phot_err), f10.sem_subset)) dec = dict(sphere_ra=float(top.ra), sphere_dec=float(top["dec"]), sphere_z=float(top.z), N=int(top.N), mean_dI=float(top.mean_dI_dex), sem_dI=sem_sphere, fog_part=fog_part, n_with_group_z=int(f10.n_with_group_z), phot_part=phot_part, phot_err=phot_err, identified_fraction=float((fog_part + phot_part) / top.mean_dI_dex), remainder=rem, remainder_err=rem_err, env_trend_individual_z=float(tr_ind.mean_dI_dex), env_trend_individual_z_err=float(tr_ind.sem_dex), env_trend_group_z=float(tr_grp.mean_dI_dex), env_trend_group_z_err=float(tr_grp.sem_dex), remainder_minus_group_trend=float(rem - tr_grp.mean_dI_dex), remainder_minus_group_trend_err=float(np.hypot(rem_err, tr_grp.sem_dex))) # ---- spheres centred on the CATALOGUE position and redshift of each cluster (not on the scan maximum): # symmetric about the cluster in redshift space, so the finger-of-god term largely cancels, and the centre # is not tuned to the data (no "most extreme of 1e5 spheres" selection of the centre). h = pd.read_csv(os.path.join(F.DATA, "sample_howlett.csv"))[["id", "zdist"]].rename(columns={"zdist": "z_group"}) s2 = s.merge(h, on="id", how="left") fits = json.load(open(os.path.join(F.DATA, "fp_fits.json"))) c_r, c_z = fits["sdss"]["fpI_coef"][2], fits["sdss"]["fpI_coef"][3] zg = s2.z_group.values with np.errstate(invalid="ignore"): d_logr = np.log10(F.angular_diameter_distance(zg) / F.angular_diameter_distance(s2.zdist.values)) d_logi = 4 * np.log10((1 + zg) / (1 + s2.zdist.values)) dI_g = s2.dI.values + d_logi - c_r * d_logr - c_z * (zg - s2.zdist.values) rows = [] for _, cl in simbad_clusters().iterrows(): zc = float(F.z_cmb_from_helio(cl.rvz_redshift, cl.ra, cl["dec"])) cc = F.radec_to_xyz(np.array([cl.ra]), np.array([cl["dec"]]), F.comoving_distance(np.array([zc])))[0] for R in (5.0, 10.0, 20.0): nb = np.array(tree.query_ball_point(cc, R), dtype=int) if len(nb) < 3: continue g = np.isfinite(dI_g[nb]) rows.append(dict(cluster=cl["id"], ra=cl.ra, dec=cl["dec"], z_helio_simbad=cl.rvz_redshift, z_cmb=zc, R_mpc=R, N=len(nb), mean_dI=float(s.dI.values[nb].mean()), sem_dI=float(s.dI.values[nb].std() / np.sqrt(len(nb))), n_group_z=int(g.sum()), mean_dI_ownz_subset=float(s.dI.values[nb][g].mean()) if g.any() else np.nan, mean_dI_groupz_subset=float(dI_g[nb][g].mean()) if g.any() else np.nan, sem_groupz_subset=float(dI_g[nb][g].std() / np.sqrt(g.sum())) if g.sum() > 1 else np.nan)) cl = pd.DataFrame(rows) cl.to_csv(os.path.join(F.DATA, "cluster_spheres.csv"), index=False, float_format="%.5g") pd.DataFrame([dec]).to_csv(os.path.join(F.DATA, "coma_decomposition.csv"), index=False, float_format="%.5g") print(json.dumps(dec, indent=1)) print(cl.to_string(index=False)) return dec, cl def patch_identifications(): """Catalogued clusters (SIMBAD object type ClG) inside each of the most significant dim SDSS spheres: within the sphere's angular radius and +-0.004 in redshift (+-1200 km/s) of its centre.""" pt = pd.read_csv(os.path.join(F.DATA, "region_extreme_patches.csv")) pt = pt[(pt["sample"] == "sdss") & (pt.variant == "raw") & pt.kind.str.contains("negative")].reset_index(drop=True) dest = os.path.join(F.RAW, "simbad_patch_clusters.csv") if os.path.exists(dest) and os.path.getsize(dest) > 100: allc = pd.read_csv(dest) else: parts = [] for k, p in pt.iterrows(): rad = float(np.degrees(p.R_mpc / F.comoving_distance(p.z))) q = ("SELECT TOP 300 main_id, otype, ra, dec, rvz_redshift, DISTANCE(POINT('ICRS',ra,dec)," "POINT('ICRS',%.5f,%.5f)) AS sep_deg FROM basic WHERE otype = 'ClG' AND " "CONTAINS(POINT('ICRS',ra,dec),CIRCLE('ICRS',%.5f,%.5f,%.4f))=1 AND rvz_redshift BETWEEN %.5f AND %.5f " "ORDER BY sep_deg" % (p.ra, p["dec"], p.ra, p["dec"], rad, p.z - 0.004, p.z + 0.004)) txt = _get(SIMBAD, dict(REQUEST="doQuery", LANG="ADQL", FORMAT="csv", QUERY=q), os.path.join(F.RAW, "simbad_patch_%02d.csv" % k)) t = pd.read_csv(io.StringIO(txt)) t["patch"] = k parts.append(t) allc = pd.concat(parts, ignore_index=True) allc.to_csv(dest, index=False) allc["main_id"] = allc.main_id.astype(str).str.replace(r"\s+", " ", regex=True).str.strip() rows = [] for k, p in pt.iterrows(): c = allc[allc.patch == k].sort_values("sep_deg") aco = c[c.main_id.str.startswith("ACO")] near = c.iloc[0] if len(c) else None rows.append(dict(R_mpc=p.R_mpc, S=p.S, N=p.N, ra=p.ra, dec=p["dec"], z=p.z, mean_dI_dex=p.mean_dI_dex, n_clusters_in_sphere=len(c), nearest_cluster=near.main_id if near is not None else "", nearest_sep_deg=float(near.sep_deg) if near is not None else np.nan, nearest_z=float(near.rvz_redshift) if near is not None else np.nan, abell_clusters_in_sphere=";".join(aco.main_id.head(6)))) out = pd.DataFrame(rows) out.to_csv(os.path.join(F.DATA, "region_patch_identifications.csv"), index=False, float_format="%.5g") print(out.to_string(index=False)) return out # ---------------------------------------------------------------------------------------------------------- # 7. AllWISE aperture photometry of the far-tail galaxies # ---------------------------------------------------------------------------------------------------------- def _irsa_row(ra, dec): aps = ",".join("w%dmag_%d,w%dsigm_%d" % (b, k, b, k) for b in (3, 4) for k in range(1, 9)) q = ("SELECT designation,ra,dec,w1mpro,w3mpro,w3sigmpro,w4mpro,w4sigmpro,w3gmag,w3gerr,w4gmag,w4gerr,%s," "ext_flg,cc_flags,ph_qual FROM allwise_p3as_psd WHERE " "CONTAINS(POINT('ICRS',ra,dec),CIRCLE('ICRS',%.6f,%.6f,%.6f))=1" % (aps, ra, dec, 5.0 / 3600)) last = None for k in range(4): try: with urllib.request.urlopen(IRSA + "?" + urllib.parse.urlencode(dict(QUERY=q, FORMAT="CSV")), timeout=120) as r: txt = r.read().decode() if txt.lstrip().startswith("<"): raise RuntimeError(txt[:300]) t = pd.read_csv(io.StringIO(txt)) if len(t) == 0: return None sep = np.hypot((t.ra - ra) * np.cos(np.radians(dec)), t["dec"] - dec) return t.iloc[int(np.argmin(sep.values))] except Exception as e: # noqa last = e time.sleep(5 * (k + 1)) raise RuntimeError("IRSA query failed: %s" % last) def wise_apertures(): """gamma_T (see vet_candidates.waste_heat) recomputed with aperture photometry. 'psf' profile-fit magnitudes (first pass; correct for point sources, too faint for resolved galaxies) 'ap8' largest AllWISE circular aperture, radius 24.75 arcsec (W3) / 49.5 arcsec (W4), no curve-of-growth correction needed at that size; includes neighbours and more sky noise -> generous upper bound 'max' the brightest of the profile-fit, the elliptical (gmag) and all eight aperture magnitudes, counting 95% upper limits as if they were detections -> the most pessimistic bound available.""" dest = os.path.join(F.RAW, "irsa_allwise_candidates.csv") cands = [] for name in ("sdss", "6df"): c = pd.read_csv(os.path.join(F.DATA, "candidates_%s_vetted.csv" % name), dtype={"id": str}) c["sample"] = name cands.append(c) if os.path.exists(dest) and os.path.getsize(dest) > 200: w = pd.read_csv(dest, dtype={"id": str}) else: # one small cone query per galaxy (about 15 s each at IRSA); four at a time, network-bound from concurrent.futures import ThreadPoolExecutor todo = [r for c in cands for _, r in c.iterrows()] def one(r): t = _irsa_row(r.ra, r["dec"]) d = dict(sample=r["sample"], label=r.label, id=r.id) if t is not None: d.update(t.to_dict()) return d with ThreadPoolExecutor(4) as ex: rows = list(ex.map(one, todo)) w = pd.DataFrame(rows) w.to_csv(dest, index=False) out = [] for c in cands: name = c["sample"].iloc[0] full = pd.read_csv(F.resid_path(name), dtype={"id": str}) keep = ["id", "msel"] + (["devmag_r", "ext_r"] if name == "sdss" else []) c = c.drop(columns=[k for k in keep if k in c.columns and k != "id"]).merge(full[keep], on="id", how="left") m = c.merge(w[w["sample"] == name].drop(columns=["sample", "id", "ra", "dec"]), on="label", how="left") for kind in ("psf", "ap8", "max"): t = m.copy() for b in (3, 4): mags = {"psf": [("w%dmpro" % b, "w%dsigmpro" % b)], "ap8": [("w%dmag_8" % b, "w%dsigm_8" % b)], "max": [("w%dmpro" % b, "w%dsigmpro" % b), ("w%dgmag" % b, "w%dgerr" % b)] + [("w%dmag_%d" % (b, k), "w%dsigm_%d" % (b, k)) for k in range(1, 9)]}[kind] mag = np.full(len(t), np.nan) det = np.zeros(len(t), int) for mc, sc in mags: v = t[mc].values.astype(float) better = np.isfinite(v) & ~(v >= mag) # brighter than the current choice (or first) mag = np.where(better, v, mag) det = np.where(better, np.isfinite(t[sc].values.astype(float)).astype(int), det) t["W%dmag" % b] = mag t["w%d_det" % b] = det g = V.waste_heat(t, name) m["gamma300K_%s" % kind] = g.gamma300K_W3.values m["gamma300K_%s_is_limit" % kind] = g.gamma300K_is_limit.values m["gamma150K_%s" % kind] = g.gamma150K_W4.values m["gamma150K_%s_is_limit" % kind] = g.gamma150K_is_limit.values m["sample"] = name out.append(m[["sample", "label", "id", "dI", "theta_e", "explained", "class", "ext_flg", "ph_qual"] + [c_ for c_ in m.columns if c_.startswith("gamma")]]) out = pd.concat(out, ignore_index=True) out.to_csv(os.path.join(F.DATA, "candidates_wise_apertures.csv"), index=False, float_format="%.5g") un = out[out.explained == 0] summ = dict(n_candidates=int(len(out)), n_unexplained=int(len(un)), n_unexplained_with_allwise=int(un.gamma300K_psf.notna().sum())) for k in ("psf", "ap8", "max"): summ["gamma300K_%s_max_unexplained" % k] = float(un["gamma300K_%s" % k].max()) summ["gamma150K_%s_max_unexplained" % k] = float(un["gamma150K_%s" % k].max()) summ["gamma300K_%s_median_unexplained" % k] = float(un["gamma300K_%s" % k].median()) summ["gamma150K_%s_median_unexplained" % k] = float(un["gamma150K_%s" % k].median()) print(un[["label", "dI", "theta_e", "gamma300K_psf", "gamma300K_ap8", "gamma300K_max", "gamma150K_psf", "gamma150K_ap8", "gamma150K_max"]].to_string(index=False)) print(json.dumps(summ, indent=1)) return summ # ---------------------------------------------------------------------------------------------------------- # 8. region limits with their scope # ---------------------------------------------------------------------------------------------------------- def region_final(): rows = pd.read_csv(os.path.join(F.DATA, "limits_regions.csv")) vol = pd.read_csv(os.path.join(F.DATA, "survey_volume.csv")) vsurv = {"sdss": float(vol.volume_mpc3[vol["sample"] == "SDSS own"].iloc[0]), "6df": float(vol.volume_mpc3[vol["sample"] == "6dFGSv"].iloc[0])} out = [] raw = rows[rows.variant == "raw"] for (samp, R, f), g in raw.groupby(["sample", "R_mpc", "f"]): v = g[g["mode"] == "volume"].iloc[0] q = g[g["mode"] == "galaxy"].iloc[0] vR = 4 * np.pi / 3 * R ** 3 a_w = not (v.fill95 < v.fill95_fix) # True if rule A gives the weaker (larger) limit fill = np.nanmax([v.fill95, v.fill95_fix]) vs = v.v_sens_mpc3 if a_w else v.v_sens_fix_mpc3 rel = v.mc_relerr if a_w else v.mc_relerr_fix zp90 = v.z_det_p90 if a_w else v.z_det_fix_p90 out.append(dict( sample=samp, R_mpc=R, f=f, n_inject_volume=int(v.n_inject), fill95_ruleA=v.fill95, mc_relerr_ruleA=v.mc_relerr, n_det_ruleA=int(v.n_detected), fill95_ruleB=v.fill95_fix, mc_relerr_ruleB=v.mc_relerr_fix, n_det_ruleB=int(v.n_detected_fix), k_obs_ruleB=int(v.k_obs_fix), fill95_weaker=fill, weaker_rule="A" if a_w else "B", mc_relerr_weaker=rel, n95_per_gpc3_weaker=fill / vR * 1e9, v_sens_mpc3_weaker=vs, tested_share_of_survey_volume=vs / vsurv[samp], n_sphere_volumes_tested=vs / vR, tested_share_ruleA=v.v_sens_mpc3 / vsurv[samp], n_sphere_volumes_ruleA=v.v_sens_mpc3 / vR, tested_share_ruleB=v.v_sens_fix_mpc3 / vsurv[samp], n_sphere_volumes_ruleB=v.v_sens_fix_mpc3 / vR, z_det_median_ruleA=v.z_det_median, z_det_p90_ruleA=v.z_det_p90, z_det_median_ruleB=v.z_det_fix_median, z_det_p90_ruleB=v.z_det_fix_p90, z_det_p90_weaker=zp90, # galaxy-weighted: fraction of sample galaxies that can lie inside such regions gal_eff_ruleA=q.frac_detected, gal_eff_ruleB=q.frac_detected_fix, mean_members=q.mean_n_in, frac_gal_inside95_ruleA=q.frac_gal_inside95, frac_gal_inside95_ruleB=q.frac_gal_inside95_fix, frac_gal_inside95_weaker=np.nanmax([q.frac_gal_inside95, q.frac_gal_inside95_fix]), gal_mc_relerr_ruleB=q.mc_relerr_fix, gal_z_det_p90_ruleB=q.z_det_fix_p90)) out = pd.DataFrame(out) out.to_csv(os.path.join(F.DATA, "limits_regions_final.csv"), index=False, float_format="%.5g") pd.set_option("display.width", 260) print(out[["sample", "R_mpc", "f", "fill95_ruleA", "fill95_ruleB", "weaker_rule", "mc_relerr_weaker", "tested_share_of_survey_volume", "n_sphere_volumes_tested", "z_det_p90_weaker", "frac_gal_inside95_ruleA", "frac_gal_inside95_ruleB"]].to_string(index=False)) return out def main(): pd.set_option("display.width", 260) S = A.load_samples() fits = json.load(open(os.path.join(F.DATA, "fp_fits.json"))) summary = {} print("-- red side of the colour window") summary["red_side"] = red_side(S, fits) print("-- mixtures of f") mix = mixture_bounds(S, fits) summary["mixture"] = mix.to_dict(orient="records") print("-- Howlett+22 eligibility of the far tail") summary["howlett_eligibility"] = howlett_eligibility().to_dict(orient="records") print("-- population behind N_eff") summary["neff_population"] = neff_population(S).to_dict(orient="records") print("-- tail asymmetry") ta, far = tail_asymmetry() summary["tail_asymmetry"] = ta.to_dict(orient="records") print("-- clusters") dec, cl = cluster_decomposition() summary["coma_decomposition"] = dec patch_identifications() print("-- AllWISE aperture photometry of the far tail") summary["wise_apertures"] = wise_apertures() print("-- region limits with scope") region_final() with open(os.path.join(F.DATA, "revision_summary.json"), "w") as fh: json.dump(summary, fh, indent=1) return True if __name__ == "__main__": main()