#!/usr/bin/env python """Extended cross-tool comparison. Adds to compare_tools.py: - normalised flux ratios (renormalised by the median ratio so spread is what matters) - recovery fractions vs reference brightness - per-tool detection counts per filter - matched-source overlap matrix - mag-binned scatter Output: writeup/figures/* and writeup/tables/extended_summary.ecsv """ from __future__ import annotations import sys from pathlib import Path import numpy as np import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt from astropy.coordinates import SkyCoord from astropy.table import Table from astropy import units as u sys.path.insert(0, str(Path(__file__).resolve().parent)) import config # noqa: E402 FIGDIR = config.BENCHMARK_ROOT / "writeup" / "figures" TABDIR = config.BENCHMARK_ROOT / "writeup" / "tables" FIGDIR.mkdir(parents=True, exist_ok=True) TABDIR.mkdir(parents=True, exist_ok=True) REF_TOOL = "photutils_iter3" TOOLS = ["crowdsource", "starbug2", "dolphot", REF_TOOL] def load(tool): p = config.tool_bandmerged(tool) if not p.exists(): return None return Table.read(p) def safe_array(t, col): if col not in t.colnames: return np.full(len(t), np.nan) a = np.array(t[col]).astype(float) return a def main(): tabs = {tool: load(tool) for tool in TOOLS} ref = tabs[REF_TOOL] # Per-tool detection-count table det_rows = [] for tool, t in tabs.items(): if t is None: continue row = {"tool": tool, "n_band_merged": len(t)} for FILT in config.FILTERS: col = f"{FILT}_detected" row[f"{FILT}_n_detected"] = int(np.sum(t[col])) if col in t.colnames else 0 det_rows.append(row) det_tab = Table(det_rows) det_tab.write(TABDIR / "detection_counts.ecsv", overwrite=True) print("=== Detection counts per tool ===") det_tab.pprint_all() # Pairwise positional matches between every (tool, ref) pair c_ref = SkyCoord(ref["ra"] * u.deg, ref["dec"] * u.deg) summary_rows = [] for tool in ["crowdsource", "starbug2", "dolphot"]: t = tabs[tool] if t is None: continue c_oth = SkyCoord(t["ra"] * u.deg, t["dec"] * u.deg) idx, sep, _ = c_ref.match_to_catalog_sky(c_oth) match = sep.arcsec < config.MATCH_RADIUS_ARCSEC # Reverse match (tool -> ref) to count tool-only sources idx_b, sep_b, _ = c_oth.match_to_catalog_sky(c_ref) match_b = sep_b.arcsec < config.MATCH_RADIUS_ARCSEC for FILT in config.FILTERS: ref_flux = safe_array(ref, f"{FILT}_flux") ref_err = safe_array(ref, f"{FILT}_flux_err") tool_flux_all = safe_array(t, f"{FILT}_flux") tool_err_all = safe_array(t, f"{FILT}_flux_err") tool_flux = tool_flux_all[idx] tool_err = tool_err_all[idx] both = match & np.isfinite(ref_flux) & np.isfinite(tool_flux) & (ref_flux > 0) & (tool_flux > 0) ratio = tool_flux[both] / ref_flux[both] med = float(np.median(ratio)) if both.sum() else np.nan # normalised log-ratio log_ratio = np.log10(ratio / med) if both.sum() else np.array([]) mad_log = float(np.median(np.abs(log_ratio))) if both.sum() else np.nan sig_log = 1.4826 * mad_log if both.sum() else np.nan # rel uncertainty median rel_ref = ref_err[both] / ref_flux[both] rel_tool = tool_err[both] / tool_flux[both] row = { "tool": tool, "filter": FILT, "n_ref_total": int(np.sum(np.isfinite(ref_flux))), "n_tool_total": int(np.sum(np.isfinite(tool_flux_all))), "n_matched_both": int(both.sum()), "median_flux_ratio_tool_over_ref": med, "log10_scatter_dex": sig_log, # 1-sigma in dex after re-normalising "rel_err_median_ref": float(np.nanmedian(rel_ref)) if both.sum() else np.nan, "rel_err_median_tool": float(np.nanmedian(rel_tool)) if both.sum() else np.nan, "n_tool_only": int(np.sum(~match_b & np.isfinite(tool_flux_all))), } summary_rows.append(row) if both.sum() == 0: continue # Normalised flux scatter (log10 ratio, mean-subtracted) fig, ax = plt.subplots(figsize=(5, 5)) ax.hist(log_ratio, bins=80, range=(-1.0, 1.0), histtype="step") ax.axvline(0.0, color="r", lw=0.8) ax.set_xlabel(r"log10(tool flux / ref flux) - log10(median)") ax.set_ylabel("N") ax.set_title(f"{tool} vs {REF_TOOL}, {FILT}\nN={both.sum()}, sigma={sig_log:.3f} dex") fig.tight_layout() fig.savefig(FIGDIR / f"flux_scatter_norm_{tool}_{FILT}.png", dpi=120) plt.close(fig) # Mag-binned scatter (mag = -2.5*log10(ref_flux)) mag_ref = -2.5 * np.log10(ref_flux[both]) order = np.argsort(mag_ref) mag_sorted = mag_ref[order] lr_sorted = log_ratio[order] # bin in 25 quantile bins n_per_bin = max(50, len(mag_sorted) // 25) bins = [] for i in range(0, len(mag_sorted) - n_per_bin + 1, n_per_bin): chunk = lr_sorted[i:i + n_per_bin] bins.append((np.mean(mag_sorted[i:i + n_per_bin]), np.median(chunk), 1.4826 * np.median(np.abs(chunk - np.median(chunk))))) if bins: bins = np.array(bins) fig, ax = plt.subplots(figsize=(6, 4)) ax.errorbar(bins[:, 0], bins[:, 1], yerr=bins[:, 2], fmt=".", capsize=2) ax.axhline(0, color="r", lw=0.5) ax.set_xlabel(f"-2.5*log10(ref_flux), {FILT}") ax.set_ylabel("median(log_ratio) +/- MAD") ax.set_title(f"{tool}: mag-binned bias and scatter, {FILT}") fig.tight_layout() fig.savefig(FIGDIR / f"flux_magbin_{tool}_{FILT}.png", dpi=120) plt.close(fig) summary = Table(summary_rows) summary.write(TABDIR / "extended_summary.ecsv", overwrite=True) print() print("=== Per-(tool, filter) extended summary ===") summary.pprint_all() # Recovery fraction plot (matched / ref total) per filter and tool fig, ax = plt.subplots(figsize=(7, 4)) width = 0.25 x = np.arange(len(config.FILTERS)) for i, tool in enumerate(["crowdsource", "starbug2", "dolphot"]): rec = [] for FILT in config.FILTERS: row = next((r for r in summary_rows if r["tool"] == tool and r["filter"] == FILT), None) if row is None or row["n_ref_total"] == 0: rec.append(0.0) else: rec.append(row["n_matched_both"] / row["n_ref_total"]) ax.bar(x + (i - 1) * width, rec, width, label=tool) ax.set_xticks(x) ax.set_xticklabels(config.FILTERS) ax.set_ylabel("matched fraction (vs photutils_iter3 detections)") ax.set_title("Recovery fraction (within %.2f arcsec)" % config.MATCH_RADIUS_ARCSEC) ax.legend() fig.tight_layout() fig.savefig(FIGDIR / "recovery_fraction.png", dpi=120) plt.close(fig) print(f"Wrote recovery_fraction.png") if __name__ == "__main__": main()