display IRFs¶

In [1]:
import ctapipe 
print(ctapipe.__version__)

0.19.0
In [2]:
 from astropy.table import Table
In [3]:
from stereograph.plots import merge_stats, compute_auc, roc_auc_per_energy
from ctapipe.io import TableLoader
import pandas as pd

import numpy as np
import matplotlib.pyplot as plt
# from gammapy.irf import load_irf_dict_from_file
import astropy.units as u
from stereograph.plots import plot_auc_per_energy, combine_auc, compute_auc, roc_auc_per_energy
import matplotlib.pyplot as plt

Benchmarks¶

In [4]:
from astropy.table import Table
In [5]:
def plot_angular_resolution(ang_res_table, fraction=68, ax=None, **kwargs):
    ax = ax or plt.gca()
    
    energy_lo = ang_res_table['ENERG_LO'][0] * ang_res_table['ENERG_LO'].unit
    energy_hi = ang_res_table['ENERG_HI'][0] * ang_res_table['ENERG_HI'].unit
    energy_cen = (energy_lo + energy_hi) / 2
    energy_half_width = (energy_hi - energy_lo) / 2

    angular_resolution_68 = ang_res_table['ANGULAR_RESOLUTION_68'][0][0] * ang_res_table['ANGULAR_RESOLUTION_68'].unit


    kwargs.setdefault('label', f'Angular Resolution ({fraction}%)')
    kwargs.setdefault('fmt', 'o')
    kwargs.setdefault('linestyle', '')
    kwargs.setdefault('capsize', 3)
    ax.errorbar(
        energy_cen.to_value(u.TeV),
        angular_resolution_68.to_value(u.deg),
        xerr=energy_half_width.to_value(u.TeV),
        yerr=None,
        **kwargs
    )

    ax.set_xscale('log')
    ax.set_xlabel(f"True energy [{energy_cen.unit}]")
    ax.set_ylabel(f"Angular Resolution ({fraction}%) / {angular_resolution_68.unit}")

    ax.set_title("Angular Resolution vs Energy")
    ax.grid(True, which='both', axis='x', linestyle='--', linewidth=0.5)
    ax.grid(True, which='major', axis='y', linestyle='--', linewidth=0.5)
    ax.legend()
    return ax
In [6]:
def plot_energy_resolution(energy_res_table, ax=None, **kwargs):
    ax = ax or plt.gca()
    energy_lo = energy_res_table['ENERG_LO'][0] * energy_res_table['ENERG_LO'].unit
    energy_hi = energy_res_table['ENERG_HI'][0] * energy_res_table['ENERG_HI'].unit
    energy_cen = (energy_lo + energy_hi) / 2
    energy_half_width = (energy_hi - energy_lo) / 2
    resolution = energy_res_table['RESOLUTION'][0][0]


    kwargs.setdefault('label', 'Energy Resolution')
    kwargs.setdefault('fmt', 'o')
    kwargs.setdefault('linestyle', '')
    kwargs.setdefault('capsize', 3)
    ax.errorbar(
        energy_cen.to_value(u.TeV),
        resolution,
        xerr=energy_half_width.to_value(u.TeV),
        yerr=None,
        **kwargs
    )
    ax.set_xscale('log')
    ax.set_xlabel(f"True Energy / {energy_cen.unit}")
    ax.set_ylabel(f"Energy Resolution")
    ax.set_title("Energy Resolution vs Energy")
    ax.grid(True, which='both', axis='x', linestyle='--', linewidth=0.5)
    ax.grid(True, which='major', axis='y', linestyle='--', linewidth=0.5)
    ax.legend()
    return ax


def plot_energy_bias(energy_res_table, ax=None, **kwargs):
    ax = ax or plt.gca()
    energy_lo = energy_res_table['ENERG_LO'][0] * energy_res_table['ENERG_LO'].unit
    energy_hi = energy_res_table['ENERG_HI'][0] * energy_res_table['ENERG_HI'].unit
    energy_cen = (energy_lo + energy_hi) / 2
    energy_half_width = (energy_hi - energy_lo) / 2
    bias = energy_res_table['BIAS'][0][0] 

    kwargs.setdefault('label', 'Energy Bias')
    kwargs.setdefault('fmt', 'o')
    kwargs.setdefault('linestyle', '')
    kwargs.setdefault('capsize', 3)
    ax.errorbar(
        energy_cen.to_value(u.TeV),
        bias,
        xerr=energy_half_width.to_value(u.TeV),
        yerr=None,
        **kwargs
    )
    ax.set_xscale('log')
    ax.set_xlabel(f"True energy / {energy_cen.unit}")
    ax.set_ylabel(f"Energy Bias")
    ax.set_title("Energy Bias vs Energy")
    ax.grid(True, which='both', axis='x', linestyle='--', linewidth=0.5)
    ax.grid(True, which='major', axis='y', linestyle='--', linewidth=0.5)
    ax.legend()
    return ax
In [7]:
def plot_sensitivity(sensitivity_table, ax=None, **kwargs):
    ax = ax or plt.gca()
    energy_lo = sensitivity_table['ENERG_LO'][0] * sensitivity_table['ENERG_LO'].unit
    energy_hi = sensitivity_table['ENERG_HI'][0] * sensitivity_table['ENERG_HI'].unit
    energy_cen = (energy_lo + energy_hi) / 2
    energy_half_width = (energy_hi - energy_lo) / 2
    flux_sensitivity = sensitivity_table['ENERGY_FLUX_SENSITIVITY'][0][0] * sensitivity_table['ENERGY_FLUX_SENSITIVITY'].unit

    kwargs.setdefault('label', 'Energy Flux Sensitivity')
    kwargs.setdefault('fmt', 'o')
    kwargs.setdefault('linestyle', '-')
    kwargs.setdefault('capsize', 3)
    ax.errorbar(
        energy_cen.to_value(u.TeV),
        flux_sensitivity.to_value(sensitivity_table['ENERGY_FLUX_SENSITIVITY'].unit),
        xerr=energy_half_width.to_value(u.TeV),
        yerr=None,
        **kwargs
    )
    ax.set_xscale('log')
    ax.set_yscale('log')
    ax.set_xlabel(f" True energy / {energy_cen.unit}")
    ax.set_ylabel(f"Energy Flux Sensitivity / ({sensitivity_table['ENERGY_FLUX_SENSITIVITY'].unit})")
    ax.set_title("Energy Flux Sensitivity vs Energy")
    ax.grid(True, which='both', axis='x', linestyle='--', linewidth=0.5)
    ax.grid(True, which='major', axis='y', linestyle='--', linewidth=0.5)
    ax.legend()
    return ax
In [8]:
from astropy.table import Table
import numpy as np
import matplotlib.pyplot as plt

# Benchmark files
benchmark_files = {
    'Random Forests': "../results/RF/benchmarks_sensitivity.fits.gz",
    'FCN': "../results/FC/benchmarks_sensitivity.fits.gz",
    'GNNs': "../results/GNN/benchmarks_sensitivity.fits.gz",
}

model_colors = {
    'Random Forests': 'tab:blue',
    'FCN': 'tab:orange',
    'GNNs': 'tab:green',
}

energy_bins = np.logspace(-2, 2, 10) * u.TeV
fig, axes = plt.subplots(3, 2, figsize=(14, 12),
                         gridspec_kw={'height_ratios': [2, 2, 1]})

angular_res_plotted = False

for name, file in benchmark_files.items():
    energy_res_table = Table.read(file, hdu=1)
    ang_res_table = Table.read(file, hdu=2)
    sensitivity_table = Table.read(file, hdu=3)

    plot_angular_resolution(
        ang_res_table, ax=axes[0, 0],
        label="Hillas Reconstructor", 
        color=model_colors[name]
    )
    angular_res_plotted = True

    plot_energy_bias(energy_res_table, ax=axes[0, 1], label=name, color=model_colors[name])
    plot_energy_resolution(energy_res_table, ax=axes[1, 0], label=name, color=model_colors[name])
    plot_sensitivity(sensitivity_table, ax=axes[1, 1], label=name, color=model_colors[name])

energy_res_rf = Table.read(benchmark_files['Random Forests'], hdu=1)
energy_res_fc = Table.read(benchmark_files['FCN'], hdu=1)
energy_res_gnn = Table.read(benchmark_files['GNNs'], hdu=1)

sensitivity_rf = Table.read(benchmark_files['Random Forests'], hdu=3)
sensitivity_fc = Table.read(benchmark_files['FCN'], hdu=3)
sensitivity_gnn = Table.read(benchmark_files['GNNs'], hdu=3)

res_rf = np.array(list(energy_res_rf['RESOLUTION'][0][0]))
res_fc = np.array(list(energy_res_fc['RESOLUTION'][0][0]))
res_gnn = np.array(list(energy_res_gnn['RESOLUTION'][0][0]))

sens_rf = np.array(list(sensitivity_rf['FLUX_SENSITIVITY'][0][0]))
sens_fc = np.array(list(sensitivity_fc['FLUX_SENSITIVITY'][0][0]))
sens_gnn = np.array(list(sensitivity_gnn['FLUX_SENSITIVITY'][0][0]))
energy_cen_res = (energy_res_rf['ENERG_LO'][0] + energy_res_rf['ENERG_HI'][0]) / 2
energy_cen_sens = (sensitivity_rf['ENERG_LO'][0] + sensitivity_rf['ENERG_HI'][0]) / 2

# FC over RF and GNN over RF
rel_energy_fc_rf = (res_rf - res_fc) / res_rf * 100
rel_sens_fc_rf = (sens_rf - sens_fc) / sens_rf * 100

rel_energy_gnn_rf = (res_rf - res_gnn) / res_rf * 100
rel_sens_gnn_rf = (sens_rf - sens_gnn) / sens_rf * 100

# averages
avg_energy_fc_rf = np.mean(rel_energy_fc_rf)
avg_sens_fc_rf = np.mean(rel_sens_fc_rf)
avg_energy_gnn_rf = np.mean(rel_energy_gnn_rf)
avg_sens_gnn_rf = np.mean(rel_sens_gnn_rf)

props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)

# relative plots energy
axes[2, 0].axhline(0, color=model_colors['Random Forests'], linestyle=':', linewidth=1.5,)
axes[2, 0].semilogx(energy_cen_res, rel_energy_fc_rf, 'o-', label="FC over RF", color=model_colors['FCN'])
axes[2, 0].semilogx(energy_cen_res, rel_energy_gnn_rf, 's--', label="GNN over RF", color=model_colors['GNNs'])

axes[2, 0].set_title("Relative Energy Resolution")
axes[2, 0].set_xlabel("True energy [TeV]")
axes[2, 0].set_ylabel("Improvement [%]")
axes[2, 0].grid(True, which='both', ls='--', lw=0.5)
axes[2, 0].legend()

axes[2, 0].text(
    0.05, 0.05,
    f'FC over RF avg: {avg_energy_fc_rf:.2f}%\nGNN over RF avg: {avg_energy_gnn_rf:.2f}%',
    transform=axes[2, 0].transAxes,
    fontsize=9,
    verticalalignment='bottom',
    bbox=props
)

# relative plots sensitivity
axes[2, 1].axhline(0, color=model_colors['Random Forests'], linestyle=':', linewidth=1.5,)
axes[2, 1].semilogx(energy_cen_sens, rel_sens_fc_rf, 'o-', label="FC over RF", color=model_colors['FCN'])
axes[2, 1].semilogx(energy_cen_sens, rel_sens_gnn_rf, 's--', label="GNN over RF", color=model_colors['GNNs'])
axes[2, 1].set_title("Relative Sensitivity")
axes[2, 1].set_xlabel("True energy [TeV]")
axes[2, 1].set_ylabel("Improvement [%]")
axes[2, 1].grid(True, which='both', ls='--', lw=0.5)

axes[2, 1].legend(
    loc='upper center',
    bbox_to_anchor=(0.5, -0.25),
    ncol=3,
    frameon=True
)

axes[2, 1].text(
    0.5, 0.08,
    f'FC over RF avg: {avg_sens_fc_rf:.2f}%\nGNN over RF avg: {avg_sens_gnn_rf:.2f}%',
    transform=axes[2, 1].transAxes,
    fontsize=9,
    verticalalignment='bottom',
    horizontalalignment='center',
    bbox=props
)

plt.subplots_adjust(bottom=0.12)
plt.tight_layout()
plt.savefig("benchmark_RF_FCN_GNN.png", dpi=300)
plt.show()
No description has been provided for this image
In [9]:
gamma_gnn="../results/GNN/dl2_class_test.h5"
fc_gamma="../results/FC/dl2_class_test.h5"
In [10]:
table_fc=pd.read_hdf(fc_gamma,key="/dl2/event/subarray/FC")
table_gnn = pd.read_hdf(gamma_gnn, key="/dl2/event/subarray/Graph")
In [11]:
def compute_auc(dataframes, energy_bins, score_column=None):
    """
    Compute the AUC for each dataframe and return the results.

    Parameters
    ----------
    dataframes : list of pandas.DataFrame
        List of dataframes, each containing the data for a single experiment.
        Each dataframe must have the columns:
        - "true_energy": The true energy of the event.
        - "true_shower_primary_id": The true type of the shower.
        - score_column: Column containing the predicted score
    energy_bins : array-like
        Energy bins to use for computing the AUC.
    score_column : str
        Name of the column with the predicted score (example: "fcClassifier_prediction")

    Returns
    -------
    aucs : list of np.ndarray
        List of AUC arrays (one per dataframe)
    """
    aucs = []

    for df in dataframes:
        true_energy = df["true_energy"].values * u.TeV
        true_type = df["true_shower_id"].values
        gammaness = df[score_column].values  # <-- flexible column

        _, auc = roc_auc_per_energy(true_type, gammaness, true_energy, energy_bins=energy_bins)
        aucs.append(auc)

    return aucs
In [12]:
import matplotlib.pyplot as plt

energy_low = energy_bins[:-1]
energy_high = energy_bins[1:]
energy_mean = np.sqrt(energy_low * energy_high)
fig, (ax1, ax2) = plt.subplots(
    2, 1, figsize=(9, 6), sharex=True,
    gridspec_kw={'height_ratios': [3, 1], 'hspace': 0.05} 
)

# Fully Connected Network (FCN)
aucs_fc = compute_auc([table_fc], energy_bins, score_column="fc_gammaness")
avg_auc_fc = merge_stats(aucs_fc)
auc_fc, auc_min_fc, auc_max_fc = avg_auc_fc[0], avg_auc_fc[2], avg_auc_fc[3]
auc_err_fc = (auc_fc - auc_min_fc, auc_max_fc - auc_fc)
ax = plot_auc_per_energy(
    energy_low, energy_high, auc_fc, auc_err=auc_err_fc, ax=ax1, label="Fully Connected Network"
)

# Graph Neural Network (GNN)
aucs_gnn = compute_auc([table_gnn], energy_bins, score_column="Graph_reco_gammaness")
avg_auc_gnn = merge_stats(aucs_gnn)
auc_gnn, auc_min_gnn, auc_max_gnn = avg_auc_gnn[0], avg_auc_gnn[2], avg_auc_gnn[3]
auc_err_gnn = (auc_gnn - auc_min_gnn, auc_max_gnn - auc_gnn)
ax = plot_auc_per_energy(
    energy_low, energy_high, auc_gnn, auc_err=auc_err_gnn, ax=ax1, label="Graph Neural Network"
)

# Random Forest
auc_rf = None
if 'RandomForest_class_reco' in table_gnn.columns and table_gnn['RandomForest_class_reco'].notnull().all():
    _, auc_rf = roc_auc_per_energy(
        table_gnn["true_shower_id"].values,
        table_gnn['RandomForest_class_reco'].values,
        table_gnn['true_energy'].values * u.TeV,
        energy_bins=energy_bins
    )
    ax = plot_auc_per_energy(
        energy_low, energy_high, auc_rf, ax=ax1, label="Random Forest"
    )


handles, labels = ax.get_legend_handles_labels()
by_label = dict(zip(labels, handles))
ax1.legend(by_label.values(), by_label.keys(), loc='upper left', bbox_to_anchor=(1.02, 1))
ax1.set_ylabel("AUC")
ax1.grid(True, which="both")


if auc_rf is not None:
    rel_improvement_fcn = (auc_fc - auc_rf) / auc_rf * 100
    rel_improvement_gnn = (auc_gnn - auc_rf) / auc_rf * 100

    ax2.semilogx(energy_mean.value, rel_improvement_fcn, '-', color='red', label='FCN ')
    ax2.semilogx(energy_mean.value, rel_improvement_gnn, '--', color='blue', label='GNN ')

    ax2.set_xlabel(r"$E_\mathrm{True}$ / TeV")
    ax2.set_ylabel("Relative Performance [%]")
    ax2.grid(True, which="both")

    margin = 5
    y_min, y_max = ax2.get_ylim()
    ax2.set_ylim(
        min(y_min, np.min([rel_improvement_fcn.min(), rel_improvement_gnn.min()]) - margin),
        max(y_max, np.max([rel_improvement_fcn.max(), rel_improvement_gnn.max()]) + margin)
    )
    ax2.legend(loc='lower right', fontsize=8, framealpha=0.9)

    ax1.legend(by_label.values(), by_label.keys(), loc='lower right', fontsize=10)

    avg_improvement_fcn = np.mean(rel_improvement_fcn)
    avg_improvement_gnn = np.mean(rel_improvement_gnn)
    textstr = f'Relative Performance :\nFCN over RF: {avg_improvement_fcn:.2f}%\nGNN over RF: {avg_improvement_gnn:.2f}%'
    ax2.text(
        0.98, 0.95, textstr,
        transform=ax2.transAxes, fontsize=7,
        verticalalignment='top', horizontalalignment='right',
        bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.5)
    )


plt.subplots_adjust(left=0.08, right=0.78, top=0.95, bottom=0.12, hspace=0.12)  
plt.savefig('auc.png', dpi=200, bbox_inches='tight')
plt.show()
No description has been provided for this image
In [ ]:
In [ ]: