Merge pull request #290 from cta-observatory/ang_res_multiple_quantiles

Multiple quantiles for angular_resolution

docs/changes/290.api.rst

@@ -0,0 +1,3 @@
+Make it possible to pass multiple quantiles to ``pyirf.benchmarks.angular_resolution``, calculating all of them.
+
+The column name(s) in the output now include(s) the percentage value of the calculated quantile, e.g. ``angular_resolution_68``.

docs/conf.py

@@ -84,8 +84,7 @@
 html_theme = "sphinx_rtd_theme"
 
 html_theme_options = {
-    'canonical_url': 'https://cta-observatory.github.io/pyirf',
-    'display_version': True,
+    'canonical_url': 'https://pyirf.readthedocs.io/',
 }
 
 intersphinx_mapping = {

examples/comparison_with_EventDisplay.ipynb

@@ -528,7 +528,7 @@
     "\n",
     "plt.errorbar(\n",
     "    0.5 * (ang_res['reco_energy_low'] + ang_res['reco_energy_high']).to_value(u.TeV),\n",
-    "    ang_res['angular_resolution'].to_value(u.deg),\n",
+    "    ang_res['angular_resolution_68'].to_value(u.deg),\n",
     "    xerr=0.5 * (ang_res['reco_energy_high'] - ang_res['reco_energy_low']).to_value(u.TeV),\n",
     "    ls='',\n",
     "    label='pyirf'\n",

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,13 @@ 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]
+
+    keys = [f"angular_resolution_{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 +72,9 @@ 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)
+        quantile_values = np.nanquantile(group["theta"], quantile)
+        for key, value in zip(keys, quantile_values):
+            result[key][bin_idx] = value
+
     return result

pyirf/benchmarks/tests/test_angular_resolution.py

@@ -8,17 +8,17 @@
 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,
+        }
     )
 
-    assert np.all(np.isnan(table["angular_resolution"]))
+    table = angular_resolution(events, [1, 10, 100] * u.TeV)
+
+    assert np.all(np.isnan(table["angular_resolution_68"]))
+
 
 @pytest.mark.parametrize("unit", (u.deg, u.rad))
 def test_angular_resolution(unit):
@@ -32,28 +32,32 @@ 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
     events["true_energy"].value[(2 * N) // 2] = np.nan
 
     bins = [1, 10, 100] * u.TeV
     table = angular_resolution(events, bins)
-    ang_res = table['angular_resolution'].to(u.deg)
+    ang_res = table["angular_resolution_68"].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_95"].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["angular_resolution_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["angular_resolution_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["angular_resolution_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)