diff --git a/docs/changes/290.feature.rst b/docs/changes/290.feature.rst new file mode 100644 index 000000000..42dacd1b5 --- /dev/null +++ b/docs/changes/290.feature.rst @@ -0,0 +1,2 @@ +Make it possible to pass multiple quantiles to ``pyirf.benchmarks.angular_resolution``, calculating all of them. +If only one quantile is passed, the output is the same as before. diff --git a/pyirf/benchmarks/angular_resolution.py b/pyirf/benchmarks/angular_resolution.py index 2f3445cfe..577b9e588 100644 --- a/pyirf/benchmarks/angular_resolution.py +++ b/pyirf/benchmarks/angular_resolution.py @@ -1,3 +1,4 @@ +from collections.abc import Sequence import numpy as np from astropy.table import QTable from scipy.stats import norm @@ -10,7 +11,8 @@ def angular_resolution( - events, energy_bins, + events, + energy_bins, energy_type="true", quantile=ONE_SIGMA_QUANTILE, ): @@ -20,6 +22,8 @@ def angular_resolution( This implementation corresponds to the 68% containment of the angular distance distribution. + Passing a list of quantiles results in all the quantiles being calculated. + Parameters ---------- events : astropy.table.QTable @@ -29,8 +33,8 @@ def angular_resolution( energy_type: str Either "true" or "reco" energy. Default is "true". - quantile : float - Which quantile to use for the angular resolution, + quantile : list(float) + Which quantile(s) to use for the angular resolution, by default, the containment of the 1-sigma region of the normal distribution (~68%) is used. @@ -52,8 +56,16 @@ def angular_resolution( result[f"{energy_key}_center"] = 0.5 * (energy_bins[:-1] + energy_bins[1:]) result["n_events"] = 0 - key = "angular_resolution" - result[key] = u.Quantity(np.nan, table["theta"].unit) + if not isinstance(quantile, Sequence): + quantile = [quantile] + + if len(quantile) < 2: + keys = ["angular_resolution"] + else: + keys = [f"containment_quantile_{value * 100:.0f}" for value in quantile] + + for key in keys: + result[key] = u.Quantity(np.nan, table["theta"].unit) # if we get an empty input (no selected events available) # we return the table filled with NaNs @@ -63,6 +75,8 @@ def angular_resolution( # use groupby operations to calculate the percentile in each bin by_bin = table[valid].group_by(bin_index[valid]) for bin_idx, group in zip(by_bin.groups.keys, by_bin.groups): - result[key][bin_idx] = np.nanquantile(group["theta"], quantile) - result["n_events"][bin_idx] = len(group) + for key, value in zip(keys, quantile): + result[key][bin_idx] = np.nanquantile(group["theta"], value) + result["n_events"][bin_idx] = len(group) + return result diff --git a/pyirf/benchmarks/tests/test_angular_resolution.py b/pyirf/benchmarks/tests/test_angular_resolution.py index 269c6c46b..fc4941ae8 100644 --- a/pyirf/benchmarks/tests/test_angular_resolution.py +++ b/pyirf/benchmarks/tests/test_angular_resolution.py @@ -8,18 +8,18 @@ def test_empty_angular_resolution(): from pyirf.benchmarks import angular_resolution - events = QTable({ - 'true_energy': [] * u.TeV, - 'theta': [] * u.deg, - }) - - table = angular_resolution( - events, - [1, 10, 100] * u.TeV + events = QTable( + { + "true_energy": [] * u.TeV, + "theta": [] * u.deg, + } ) + table = angular_resolution(events, [1, 10, 100] * u.TeV) + assert np.all(np.isnan(table["angular_resolution"])) + @pytest.mark.parametrize("unit", (u.deg, u.rad)) def test_angular_resolution(unit): from pyirf.benchmarks import angular_resolution @@ -32,20 +32,24 @@ def test_angular_resolution(unit): true_resolution = np.append(np.full(N, TRUE_RES_1), np.full(N, TRUE_RES_2)) - rng = np.random.default_rng(0) - events = QTable({ - 'true_energy': np.concatenate([ - [0.5], # below bin 1 to test with underflow - np.full(N - 1, 5.0), - np.full(N - 1, 50.0), - [500], # above bin 2 to test with overflow - ]) * u.TeV, - 'theta': np.abs(rng.normal(0, true_resolution)) * u.deg - }) + events = QTable( + { + "true_energy": np.concatenate( + [ + [0.5], # below bin 1 to test with underflow + np.full(N - 1, 5.0), + np.full(N - 1, 50.0), + [500], # above bin 2 to test with overflow + ] + ) + * u.TeV, + "theta": np.abs(rng.normal(0, true_resolution)) * u.deg, + } + ) - events['theta'] = events['theta'].to(unit) + events["theta"] = events["theta"].to(unit) # add nans to test if nans are ignored events["true_energy"].value[N // 2] = np.nan @@ -53,7 +57,7 @@ def test_angular_resolution(unit): bins = [1, 10, 100] * u.TeV table = angular_resolution(events, bins) - ang_res = table['angular_resolution'].to(u.deg) + ang_res = table["angular_resolution"].to(u.deg) assert len(ang_res) == 2 assert u.isclose(ang_res[0], TRUE_RES_1 * u.deg, rtol=0.05) assert u.isclose(ang_res[1], TRUE_RES_2 * u.deg, rtol=0.05) @@ -64,9 +68,26 @@ def test_angular_resolution(unit): # 2 sigma coverage interval quantile = norm(0, 1).cdf(2) - norm(0, 1).cdf(-2) table = angular_resolution(events, bins, quantile=quantile) - ang_res = table['angular_resolution'].to(u.deg) - + ang_res = table["angular_resolution"].to(u.deg) assert len(ang_res) == 2 - assert u.isclose(ang_res[0], 2 * TRUE_RES_1 * u.deg, rtol=0.05) assert u.isclose(ang_res[1], 2 * TRUE_RES_2 * u.deg, rtol=0.05) + + # 25%, 50%, 90% coverage interval + table = angular_resolution(events, bins, quantile=[0.25, 0.5, 0.9]) + cov_25 = table["containment_quantile_25"].to(u.deg) + assert len(cov_25) == 2 + assert u.isclose( + cov_25[0], norm(0, TRUE_RES_1).interval(0.25)[1] * u.deg, rtol=0.05 + ) + assert u.isclose( + cov_25[1], norm(0, TRUE_RES_2).interval(0.25)[1] * u.deg, rtol=0.05 + ) + + cov_50 = table["containment_quantile_50"].to(u.deg) + assert u.isclose(cov_50[0], norm(0, TRUE_RES_1).interval(0.5)[1] * u.deg, rtol=0.05) + assert u.isclose(cov_50[1], norm(0, TRUE_RES_2).interval(0.5)[1] * u.deg, rtol=0.05) + + cov_90 = table["containment_quantile_90"].to(u.deg) + assert u.isclose(cov_90[0], norm(0, TRUE_RES_1).interval(0.9)[1] * u.deg, rtol=0.05) + assert u.isclose(cov_90[1], norm(0, TRUE_RES_2).interval(0.9)[1] * u.deg, rtol=0.05)