1.2 — ECDFS photo-z evaluation (4 bands, FlexZBoost)¶

Authors: Xiangchong Li, Tianqing Zhang, Eric Charles, Sam Schmidt LSST Science Pipelines: Weekly 2025_49 Kernel: Python-dp1

Purpose¶

Same evaluation pipeline as 1_1_ecdfs_photoz_evaluate.ipynb, but for the 4-band FlexZBoost variant trained on griz only (no u, no y). Used to quantify the photo-z degradation when the bluest / reddest bands are removed.

Steps¶

  1. Load the precomputed 4-band FlexZBoost output (fzb_4bands/.../output_estimate_fzboost.hdf5).
  2. Load the train/test/data tables.
  3. Color-histogram sanity checks for g-r, r-i, i-z, z-y.
  4. Plot i_mag_gauss2 vs FPFS trace to motivate the resolution cut.
  5. Apply the 4-band basic golden selection (griz mag cuts + trace + spec-z confidence > 0.82) and show 2-D (z_ref, z_phot).
  6. Compute and plot photo-z metrics vs i-mag and vs z_ref.

Inputs¶

  • photoz/fzb_4bands/projects/com_cam/data/gold_rubin/output_estimate_fzboost.hdf5
  • photoz/fzb_4bands/projects/com_cam/data/gold_rubin/model_inform_fzboost.pkl
  • rail_data/{data,train,test}_with_mag.hdf5

Outputs¶

Inline figures only (color histograms, mag–trace, 2-D z_ref-vs-z_phot, metric panels).

In [1]:
import os

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import colors

import xlens

import tables_io
from astropy.stats import biweight_location, biweight_scale
from scipy.stats import sigmaclip

# =============================
# Matplotlib defaults
# =============================
SMALL_SIZE = 8
MEDIUM_SIZE = 10
BIGGER_SIZE = 20
plt.rc("font", size=BIGGER_SIZE)
plt.rc("axes", titlesize=BIGGER_SIZE)
plt.rc("axes", labelsize=BIGGER_SIZE)
plt.rc("xtick", labelsize=BIGGER_SIZE)
plt.rc("ytick", labelsize=BIGGER_SIZE)
plt.rc("legend", fontsize=BIGGER_SIZE)
plt.rc("figure", titlesize=BIGGER_SIZE)
In [2]:
pscratch= "/gpfs02/astro/workarea/xli6/"
Dir = os.path.join(
    pscratch,
    "work/DP1/"
)
rail_pz_h5 = os.path.join(
    Dir, "photoz",
    "fzb_4bands/projects/com_cam/data/gold_rubin/output_estimate_fzboost.hdf5",
)

MODEL_PATH = os.path.join(
    Dir, "photoz",
    "fzb_4bands/projects/com_cam/data/gold_rubin/model_inform_fzboost.pkl",
)
In [3]:
data_fname = os.path.join(Dir, "rail_data", "data_with_mag.hdf5")
df = tables_io.read(
    data_fname,
    tType="astropyTable",
)
train_fname = os.path.join(Dir, "rail_data", "train_with_mag.hdf5")
df_train = tables_io.read(
    train_fname,
    tType="astropyTable",
)

test_fname = os.path.join(Dir, "rail_data", "test_with_mag.hdf5")
df_test = tables_io.read(
    test_fname,
    tType="pandasDataFrame",
)
In [4]:
fig, axs = plt.subplots(1, 4, figsize=(24, 5))
xbins = np.linspace(0, 4, 14)
magbins = np.linspace(14, 26.0, 52)
colorbins = np.linspace(-1, 2.5, 50)
bands = "grizy"
for i in range(4):
    bandi = bands[i]
    bandj = bands[i + 1]
    xi = f"{bandi}_mag_gauss2"
    xj = f"{bandj}_mag_gauss2"
    axs[i].hist(df_train[xi] - df_train[xj], bins=colorbins, alpha=0.35, color="b", density=True)
    axs[i].hist(df[xi] - df[xj], bins=colorbins, alpha=0.8, color="k", density=True, histtype="step", linewidth=2.0)
    axs[i].set_xlabel(f"{bandi}-{bandj}", fontsize=25)
plt.tight_layout()
No description has been provided for this image
In [5]:
pdfs = tables_io.read(rail_pz_h5)["data"]["yvals"]
ancil = tables_io.read(rail_pz_h5)["ancil"]
In [6]:
tables_io.read(rail_pz_h5)["meta"]["xvals"]
Out[6]:
array([[0.  , 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09, 0.1 ,
        0.11, 0.12, 0.13, 0.14, 0.15, 0.16, 0.17, 0.18, 0.19, 0.2 , 0.21,
        0.22, 0.23, 0.24, 0.25, 0.26, 0.27, 0.28, 0.29, 0.3 , 0.31, 0.32,
        0.33, 0.34, 0.35, 0.36, 0.37, 0.38, 0.39, 0.4 , 0.41, 0.42, 0.43,
        0.44, 0.45, 0.46, 0.47, 0.48, 0.49, 0.5 , 0.51, 0.52, 0.53, 0.54,
        0.55, 0.56, 0.57, 0.58, 0.59, 0.6 , 0.61, 0.62, 0.63, 0.64, 0.65,
        0.66, 0.67, 0.68, 0.69, 0.7 , 0.71, 0.72, 0.73, 0.74, 0.75, 0.76,
        0.77, 0.78, 0.79, 0.8 , 0.81, 0.82, 0.83, 0.84, 0.85, 0.86, 0.87,
        0.88, 0.89, 0.9 , 0.91, 0.92, 0.93, 0.94, 0.95, 0.96, 0.97, 0.98,
        0.99, 1.  , 1.01, 1.02, 1.03, 1.04, 1.05, 1.06, 1.07, 1.08, 1.09,
        1.1 , 1.11, 1.12, 1.13, 1.14, 1.15, 1.16, 1.17, 1.18, 1.19, 1.2 ,
        1.21, 1.22, 1.23, 1.24, 1.25, 1.26, 1.27, 1.28, 1.29, 1.3 , 1.31,
        1.32, 1.33, 1.34, 1.35, 1.36, 1.37, 1.38, 1.39, 1.4 , 1.41, 1.42,
        1.43, 1.44, 1.45, 1.46, 1.47, 1.48, 1.49, 1.5 , 1.51, 1.52, 1.53,
        1.54, 1.55, 1.56, 1.57, 1.58, 1.59, 1.6 , 1.61, 1.62, 1.63, 1.64,
        1.65, 1.66, 1.67, 1.68, 1.69, 1.7 , 1.71, 1.72, 1.73, 1.74, 1.75,
        1.76, 1.77, 1.78, 1.79, 1.8 , 1.81, 1.82, 1.83, 1.84, 1.85, 1.86,
        1.87, 1.88, 1.89, 1.9 , 1.91, 1.92, 1.93, 1.94, 1.95, 1.96, 1.97,
        1.98, 1.99, 2.  , 2.01, 2.02, 2.03, 2.04, 2.05, 2.06, 2.07, 2.08,
        2.09, 2.1 , 2.11, 2.12, 2.13, 2.14, 2.15, 2.16, 2.17, 2.18, 2.19,
        2.2 , 2.21, 2.22, 2.23, 2.24, 2.25, 2.26, 2.27, 2.28, 2.29, 2.3 ,
        2.31, 2.32, 2.33, 2.34, 2.35, 2.36, 2.37, 2.38, 2.39, 2.4 , 2.41,
        2.42, 2.43, 2.44, 2.45, 2.46, 2.47, 2.48, 2.49, 2.5 , 2.51, 2.52,
        2.53, 2.54, 2.55, 2.56, 2.57, 2.58, 2.59, 2.6 , 2.61, 2.62, 2.63,
        2.64, 2.65, 2.66, 2.67, 2.68, 2.69, 2.7 , 2.71, 2.72, 2.73, 2.74,
        2.75, 2.76, 2.77, 2.78, 2.79, 2.8 , 2.81, 2.82, 2.83, 2.84, 2.85,
        2.86, 2.87, 2.88, 2.89, 2.9 , 2.91, 2.92, 2.93, 2.94, 2.95, 2.96,
        2.97, 2.98, 2.99, 3.  , 3.01, 3.02, 3.03, 3.04, 3.05, 3.06, 3.07,
        3.08, 3.09, 3.1 , 3.11, 3.12, 3.13, 3.14, 3.15, 3.16, 3.17, 3.18,
        3.19, 3.2 , 3.21, 3.22, 3.23, 3.24, 3.25, 3.26, 3.27, 3.28, 3.29,
        3.3 , 3.31, 3.32, 3.33, 3.34, 3.35, 3.36, 3.37, 3.38, 3.39, 3.4 ,
        3.41, 3.42, 3.43, 3.44, 3.45, 3.46, 3.47, 3.48, 3.49, 3.5 , 3.51,
        3.52, 3.53, 3.54, 3.55, 3.56, 3.57, 3.58, 3.59, 3.6 , 3.61, 3.62,
        3.63, 3.64, 3.65, 3.66, 3.67, 3.68, 3.69, 3.7 , 3.71, 3.72, 3.73,
        3.74, 3.75, 3.76, 3.77, 3.78, 3.79, 3.8 , 3.81, 3.82, 3.83, 3.84,
        3.85, 3.86, 3.87, 3.88, 3.89, 3.9 , 3.91, 3.92, 3.93, 3.94, 3.95,
        3.96, 3.97, 3.98, 3.99, 4.  , 4.01, 4.02, 4.03, 4.04, 4.05, 4.06,
        4.07, 4.08, 4.09, 4.1 , 4.11, 4.12, 4.13, 4.14, 4.15, 4.16, 4.17,
        4.18, 4.19, 4.2 , 4.21, 4.22, 4.23, 4.24, 4.25, 4.26, 4.27, 4.28,
        4.29, 4.3 , 4.31, 4.32, 4.33, 4.34, 4.35, 4.36, 4.37, 4.38, 4.39,
        4.4 , 4.41, 4.42, 4.43, 4.44, 4.45, 4.46, 4.47, 4.48, 4.49, 4.5 ,
        4.51, 4.52, 4.53, 4.54, 4.55, 4.56, 4.57, 4.58, 4.59, 4.6 , 4.61,
        4.62, 4.63, 4.64, 4.65, 4.66, 4.67, 4.68, 4.69, 4.7 , 4.71, 4.72,
        4.73, 4.74, 4.75, 4.76, 4.77, 4.78, 4.79, 4.8 , 4.81, 4.82, 4.83,
        4.84, 4.85, 4.86, 4.87, 4.88, 4.89, 4.9 , 4.91, 4.92, 4.93, 4.94,
        4.95, 4.96, 4.97, 4.98, 4.99, 5.  ]])
In [7]:
# remove stars in the catalog with resolution cut
plt.scatter(df_train["i_mag_gauss2"], (df_train["i_fpfs1_m00"] + df_train["i_fpfs1_m20"]) / df_train["i_fpfs1_m00"], s=0.1)
Out[7]:
<matplotlib.collections.PathCollection at 0x7f63067c7500>
No description has been provided for this image
In [8]:
import xlens
out = xlens.catalog.redshift.get_point_estimates_from_pdfs(pdfs)
fzb_mode = np.ravel(ancil["zmode"])
fzb_mode = np.ravel(out["z500"])
fzb_mode = np.ravel(ancil["z_median"])
z_truth = np.ravel(ancil["redshift"])
msk = (
    (df_test["g_mag_gauss2"] < 27.8) & 
    (df_test["r_mag_gauss2"] < 27.1) & 
    (df_test["i_mag_gauss2"] < 24.5) & 
    (df_test["z_mag_gauss2"] < 25.8) & 
    ((df_test["i_fpfs1_m00"] + df_test["i_fpfs1_m20"]) / df_test["i_fpfs1_m00"] > 0.05)
).to_numpy()
fig, axes = plt.subplots(1,1, figsize=(5, 5))
plt.hist2d(z_truth[msk], fzb_mode[msk], range = ((0,4), (0,4)), bins = (100,100),
           cmap = 'jet', norm=colors.LogNorm())
plt.xlim(0,3.01)
plt.ylim(0,3.01)
plt.plot(np.linspace(0, 4.1), np.linspace(0, 4.1), ls="--", color="black")
plt.colorbar()

plt.xlabel(r'$z_\text{ref}$')
plt.ylabel(r'$z_\text{est}$')
Out[8]:
Text(0, 0.5, '$z_\\text{est}$')
No description has been provided for this image
In [9]:
def compute_nmad(dz: np.ndarray):
    med_dz = np.median(dz)
    return 1.48 * np.median(np.abs(dz - med_dz))


def get_biweight_mean_sigma_outlier(
    subset: np.ndarray, nclip: int = 0, abs_out_thresh: float = 0.15
):
    if sigmaclip is None or biweight_location is None or biweight_scale is None:
        raise RuntimeError("scipy and astropy.stats are required for this function")
    subset_clip, _, _ = sigmaclip(subset, low=99, high=99)
    for _ in range(nclip):
        subset_clip, _, _ = sigmaclip(subset_clip, low=5, high=5)

    mean = biweight_location(subset_clip)
    std = biweight_scale(subset_clip)
    nmad = compute_nmad(subset_clip)
    outlier_rate = np.sum(np.abs(subset) > np.min([3 * biweight_scale(subset_clip), 0.06])) / len(subset)
    abs_outlier_rate = np.sum(np.abs(subset) > abs_out_thresh) / len(subset)

    return (
        float(mean),
        float(std / np.sqrt(len(subset_clip))),
        float(std),
        float(nmad),
        float(outlier_rate),
        float(abs_outlier_rate),
    )


def process_data(
    zphot: np.ndarray,
    specz: np.ndarray,
    low: float = 0.01,
    high: float = 2.0,
    nclip: int = 0,
    nbin: int = 10,
    zbin_type: str = "spec",
    mag = None,
):
    if sigmaclip is None:
        raise RuntimeError("scipy is required for process_data")

    dz = (zphot - specz) / (1 + specz)
    z_bins = np.linspace(low, high, nbin)

    if zbin_type == "spec":
        zx = specz
    elif zbin_type == "mag":
        zx = mag if mag is not None else specz
    else:
        zx = zphot

    bin_indices = np.digitize(zx, bins=z_bins) - 1

    biweight_mean: list[float] = []
    biweight_std: list[float] = []
    biweight_sigma: list[float] = []
    nmad_list: list[float] = []
    biweight_outlier: list[float] = []
    z_mean: list[float] = []
    qt_95_low: list[float] = []
    qt_68_low: list[float] = []
    median: list[float] = []
    qt_68_high: list[float] = []
    qt_95_high: list[float] = []

    for i in range(len(z_bins) - 1):
        subset = dz[bin_indices == i]
        if len(subset) < 1:
            continue
        subset_clip, _, _ = sigmaclip(subset, low=5, high=5)

        biweight_mean.append(float(np.median(subset_clip)))
        biweight_std.append(float(np.std(subset_clip) / np.sqrt(len(subset_clip))))
        biweight_sigma.append(float(np.std(subset_clip)))
        nmad_list.append(compute_nmad(subset_clip))
        outlier_rate = np.sum(np.abs(subset) > 0.15) / len(subset)
        biweight_outlier.append(float(outlier_rate))

        qt_95_low.append(float(np.percentile(subset, 2.5)))
        qt_68_low.append(float(np.percentile(subset, 16)))
        median.append(float(np.percentile(subset, 50)))
        qt_68_high.append(float(np.percentile(subset, 84)))
        qt_95_high.append(float(np.percentile(subset, 97.5)))

        z_mean.append(float(np.mean(zx[bin_indices == i])))

    return {
        "z_mean": z_mean,
        "biweight_mean": biweight_mean,
        "biweight_std": biweight_std,
        "biweight_sigma": biweight_sigma,
        "biweight_outlier": biweight_outlier,
        "qt_95_low": qt_95_low,
        "qt_68_low": qt_68_low,
        "median": median,
        "qt_68_high": qt_68_high,
        "qt_95_high": qt_95_high,
        "nmad": nmad_list,
    }
In [10]:
algorithms = ['fzboost']
method_name_list = ['FlexZBoost', 'BPz', 'CMNN' ,'DNF']
color_list = [
    "#1f77b4", "#ff7f0e", "#2ca02c", "#d62728",
    "#9467bd", "#8c564b", "#e377c2", "#17becf"
]
fig, axes = plt.subplots(1, 3, figsize=(16, 5), sharex=True)
metric_names = [r'$< \Delta z >$', r'$\sigma_{NMAD}$', r'$P(\Delta z > 0.15)$']
title = ['Bias', 'Scatter', 'Outlier Rate']
y_keys = ['biweight_mean', 'nmad', 'biweight_outlier']
for algo, color, algo_name in zip(algorithms, color_list, method_name_list):
    result = process_data(
        zphot=fzb_mode[msk], 
        specz=z_truth[msk], 
        low=18.5, high=24.5, 
        nbin=10, 
        nclip = 0,  
        zbin_type = "mag", 
        mag = df_test['i_mag_gauss2'][msk]
    )
    for i in range(3):
        axes[i].plot(result['z_mean'], result[y_keys[i]], label=algo_name, color=color)
        axes[i].set_title(title[i])
        axes[i].set_xlabel(r'$i$-mag')

axes[0].set_ylim(-0.1, 0.1)
axes[0].set_xlim(18, 24.5)
axes[0].plot([0,2], [0, 0], '--', color = 'black', alpha = 0.3)
axes[1].axhspan(0.0, 0.03, alpha=0.8, color="gray")
axes[1].axhspan(0.0, 0.10, alpha=0.3, color="gray")

# --- Format plots ---
for i, ax in enumerate(axes):
    ax.set_ylabel(metric_names[i], fontsize=16)
axes[2].legend(loc='upper left', ncol=2, fontsize=14)
plt.tight_layout()
No description has been provided for this image
In [11]:
fig, axes = plt.subplots(1, 3, figsize=(16, 5), sharex=True)
metric_names = [r'$< \Delta z >$', r'$\sigma_{NMAD}$', r'$P(\Delta z > 0.15)$']
title = ['Bias', 'Scatter', 'Outlier Rate']
y_keys = ['biweight_mean', 'nmad', 'biweight_outlier']
for algo, color, algo_name in zip(algorithms, color_list, method_name_list):
    result = process_data(
        zphot=fzb_mode[msk], 
        specz=z_truth[msk], 
        low=0.2, high=2.0, 
        nbin=10, 
        nclip=0,  
        zbin_type="spec", 
    )
    for i in range(3):
        axes[i].plot(result['z_mean'], result[y_keys[i]], label=algo_name, color=color)
        axes[i].set_title(title[i])
        axes[i].set_xlabel(r'$z_\text{ref}$')

axes[0].set_ylim(-0.1, 0.1)
axes[0].set_xlim(-0.02, 2.02)
axes[0].plot([0,2], [0, 0], '--', color = 'black', alpha = 0.3)

# --- Format plots ---
for i, ax in enumerate(axes):
    ax.set_ylabel(metric_names[i], fontsize=16)
axes[2].legend(loc='upper left', ncol=2, fontsize=14)
axes[1].axhspan(0.0, 0.03, alpha=0.8, color="gray")
axes[1].axhspan(0.0, 0.10, alpha=0.3, color="gray")
plt.tight_layout()
No description has been provided for this image
In [ ]:
In [ ]: