From 11cd4f3f97aff95ec3969b7110fbcbde7f3a83d2 Mon Sep 17 00:00:00 2001 From: Michele Peresano Date: Tue, 12 Dec 2023 11:58:52 +0100 Subject: [PATCH 1/5] Update energy_dispersion_to_migration Remove stale debugging code Add news fragment Update energy_dispersion_to_migration - divide by the migration bin widths before resampling - minor style and format changes Fix pdf conversion to probability, remove final norm --- docs/changes/273.bugfix.rst | 4 ++++ pyirf/irf/energy_dispersion.py | 34 +++++++++++++++++++++------------- 2 files changed, 25 insertions(+), 13 deletions(-) create mode 100644 docs/changes/273.bugfix.rst diff --git a/docs/changes/273.bugfix.rst b/docs/changes/273.bugfix.rst new file mode 100644 index 000000000..607ffa71f --- /dev/null +++ b/docs/changes/273.bugfix.rst @@ -0,0 +1,4 @@ +Fix ``pyirf.irfs.energy_dispersion.energy_dispersion_to_migration``. +This function was not adapted to the now correct normalization of the +energy dispersion matrix, resulting in wrong results on the now correct +matrices. \ No newline at end of file diff --git a/pyirf/irf/energy_dispersion.py b/pyirf/irf/energy_dispersion.py index de631e426..f6100adda 100644 --- a/pyirf/irf/energy_dispersion.py +++ b/pyirf/irf/energy_dispersion.py @@ -29,7 +29,10 @@ def _normalize_hist(hist, migration_bins): def energy_dispersion( - selected_events, true_energy_bins, fov_offset_bins, migration_bins, + selected_events, + true_energy_bins, + fov_offset_bins, + migration_bins, ): """ Calculate energy dispersion for the given DL2 event list. @@ -88,8 +91,8 @@ def energy_dispersion_to_migration( new_reco_energy_edges, ): """ - Construct a energy migration matrix from a energy - dispersion matrix. + Construct a energy migration matrix from an energy dispersion matrix. + Depending on the new energy ranges, the sum over the first axis can be smaller than 1. The new true energy bins need to be a subset of the old range, @@ -110,21 +113,25 @@ def energy_dispersion_to_migration( new_reco_energy_edges: astropy.units.Quantity[energy] Reco energy edges matching the second dimension of the output - Returns: - -------- + Returns + ------- migration_matrix: numpy.ndarray Three-dimensional energy migration matrix. The third dimension equals the fov offset dimension of the energy dispersion matrix. """ + migration_matrix = np.zeros( + ( + len(new_true_energy_edges) - 1, + len(new_reco_energy_edges) - 1, + dispersion_matrix.shape[2], + ) + ) - migration_matrix = np.zeros(( - len(new_true_energy_edges) - 1, - len(new_reco_energy_edges) - 1, - dispersion_matrix.shape[2], - )) + migra_width = np.diff(disp_migration_edges) + probability = dispersion_matrix * migra_width[np.newaxis, :, np.newaxis] true_energy_interpolation = resample_histogram1d( - dispersion_matrix, + probability, disp_true_energy_edges, new_true_energy_edges, axis=0, @@ -137,13 +144,14 @@ def energy_dispersion_to_migration( for idx, e_true in enumerate( (new_true_energy_edges[1:] + new_true_energy_edges[:-1]) / 2 ): - # get migration for the new true energy bin e_true_dispersion = true_energy_interpolation[idx] with warnings.catch_warnings(): # silence inf/inf division warning - warnings.filterwarnings('ignore', 'invalid value encountered in true_divide') + warnings.filterwarnings( + "ignore", "invalid value encountered in true_divide" + ) interpolation_edges = new_reco_energy_edges / e_true y = resample_histogram1d( From 8507697cf34b0bd9829ccc8d718a3ee6ae6b9d43 Mon Sep 17 00:00:00 2001 From: Michele Peresano Date: Tue, 12 Dec 2023 14:53:06 +0100 Subject: [PATCH 2/5] Add function to compute energy migration matrix from events Simplify normalization code for migra matrix --- pyirf/irf/energy_dispersion.py | 53 +++++++++++++++++++++++++++++++++- 1 file changed, 52 insertions(+), 1 deletion(-) diff --git a/pyirf/irf/energy_dispersion.py b/pyirf/irf/energy_dispersion.py index f6100adda..c3d7bcac8 100644 --- a/pyirf/irf/energy_dispersion.py +++ b/pyirf/irf/energy_dispersion.py @@ -6,7 +6,8 @@ __all__ = [ "energy_dispersion", - "energy_dispersion_to_migration" + "energy_migration_matrix", + "energy_dispersion_to_migration", ] @@ -83,6 +84,56 @@ def energy_dispersion( return energy_dispersion +@u.quantity_input(true_energy_bins=u.TeV, reco_energy_bins=u.TeV, fov_offset_bins=u.deg) +def energy_migration_matrix( + events, true_energy_bins, reco_energy_bins, fov_offset_bins +): + """Compute the energy migration matrix directly from the events. + + Parameters + ---------- + events : `~astropy.table.QTable` + Table of the DL2 events. + Required columns: ``reco_energy``, ``true_energy``, ``true_source_fov_offset``. + true_energy_bins : `~astropy.units.Quantity` + Bin edges in true energy. + reco_energy_bins : `~astropy.units.Quantity` + Bin edges in reconstructed energy. + + Returns + ------- + matrix : array-like + Migration matrix as probabilities along the true + energy acis with shape + (n_true_energy_bins, n_reco_energy_bins, n_fov_offset_bins) + containing energies in TeV. + """ + + hist, _ = np.histogramdd( + np.column_stack( + [ + events["true_energy"].to_value(u.TeV), + events["reco_energy"].to_value(u.TeV), + events["true_source_fov_offset"].to_value(u.deg), + ] + ), + bins=[ + true_energy_bins.to_value(u.TeV), + reco_energy_bins.to_value(u.TeV), + fov_offset_bins.to_value(u.deg), + ], + ) + + with np.errstate(invalid="ignore"): + hist /= hist.sum(axis=1)[:, np.newaxis, :] + # the nans come from the fact that the sum along the reconstructed energy axis + # might sometimes be 0 when there are no events in that given true energy bin + # and fov offset bin + hist[np.isnan(hist)] = 0 + + return hist + + def energy_dispersion_to_migration( dispersion_matrix, disp_true_energy_edges, From 3be9e03d5aa2fd9e9b8bc1bc80f957e9405c25f0 Mon Sep 17 00:00:00 2001 From: Michele Peresano Date: Tue, 12 Dec 2023 15:03:55 +0100 Subject: [PATCH 3/5] Add unit test for energy migration matrix from events update unit test --- pyirf/irf/tests/test_energy_dispersion.py | 104 +++++++++++++++++----- 1 file changed, 80 insertions(+), 24 deletions(-) diff --git a/pyirf/irf/tests/test_energy_dispersion.py b/pyirf/irf/tests/test_energy_dispersion.py index 036d106de..ec0e444fb 100644 --- a/pyirf/irf/tests/test_energy_dispersion.py +++ b/pyirf/irf/tests/test_energy_dispersion.py @@ -64,17 +64,29 @@ def ppf(cdf, bins, value): assert np.isclose( TRUE_SIGMA_1, - 0.5 * (ppf(cdf[0, :, 0], migration_bins, 0.84) - ppf(cdf[0, :, 0], migration_bins, 0.16)), + 0.5 + * ( + ppf(cdf[0, :, 0], migration_bins, 0.84) + - ppf(cdf[0, :, 0], migration_bins, 0.16) + ), rtol=0.1, ) assert np.isclose( TRUE_SIGMA_2, - 0.5 * (ppf(cdf[1, :, 0], migration_bins, 0.84) - ppf(cdf[1, :, 0], migration_bins, 0.16)), + 0.5 + * ( + ppf(cdf[1, :, 0], migration_bins, 0.84) + - ppf(cdf[1, :, 0], migration_bins, 0.16) + ), rtol=0.1, ) assert np.isclose( TRUE_SIGMA_3, - 0.5 * (ppf(cdf[2, :, 0], migration_bins, 0.84) - ppf(cdf[2, :, 0], migration_bins, 0.16)), + 0.5 + * ( + ppf(cdf[2, :, 0], migration_bins, 0.84) + - ppf(cdf[2, :, 0], migration_bins, 0.16) + ), rtol=0.1, ) @@ -85,16 +97,15 @@ def test_energy_dispersion_to_migration(): np.random.seed(0) N = 10000 - true_energy_bins = 10**np.arange(np.log10(0.2), np.log10(200), 1/10) * u.TeV + true_energy_bins = 10 ** np.arange(np.log10(0.2), np.log10(200), 1 / 10) * u.TeV fov_offset_bins = np.array([0, 1, 2]) * u.deg migration_bins = np.linspace(0, 2, 101) - true_energy = np.random.uniform( - true_energy_bins[0].value, - true_energy_bins[-1].value, - size=N - ) * u.TeV + true_energy = ( + np.random.uniform(true_energy_bins[0].value, true_energy_bins[-1].value, size=N) + * u.TeV + ) reco_energy = true_energy * np.random.uniform(0.5, 1.5, size=N) selected_events = QTable( @@ -116,15 +127,14 @@ def test_energy_dispersion_to_migration(): ) # migration matrix selecting a limited energy band with different binsizes - new_true_energy_bins = 10**np.arange(np.log10(2), np.log10(20), 1/5) * u.TeV - new_reco_energy_bins = 10**np.arange(np.log10(2), np.log10(20), 1/5) * u.TeV + new_true_energy_bins = 10 ** np.arange(np.log10(2), np.log10(20), 1 / 5) * u.TeV + new_reco_energy_bins = 10 ** np.arange(np.log10(2), np.log10(20), 1 / 5) * u.TeV migration_matrix = energy_dispersion_to_migration( dispersion_matrix, true_energy_bins, migration_bins, new_true_energy_bins, new_reco_energy_bins, - ) # test dimension @@ -133,14 +143,57 @@ def test_energy_dispersion_to_migration(): assert migration_matrix.shape[2] == dispersion_matrix.shape[2] # test that all migrations are included for central energies - assert np.isclose(migration_matrix.sum(axis=1).max(), 1, rtol=0.01) + np.allclose(migration_matrix.sum(axis=1).max(), 1, rtol=0.1) # test that migrations dont always sum to 1 (since some energies are # not included in the matrix) - assert migration_matrix.sum() < (len(new_true_energy_bins) - 1) * (len(fov_offset_bins) - 1) + assert migration_matrix.sum() < (len(new_true_energy_bins) - 1) * ( + len(fov_offset_bins) - 1 + ) assert np.all(np.isfinite(migration_matrix)) +def test_energy_migration_matrix_from_events(): + from pyirf.irf.energy_dispersion import energy_migration_matrix + + np.random.seed(0) + N = 10000 + true_energy_bins = 10 ** np.arange(np.log10(0.2), np.log10(200), 1 / 10) * u.TeV + reco_energy_bins = 10 ** np.arange(np.log10(2), np.log10(20), 1 / 5) * u.TeV + fov_offset_bins = np.array([0, 1, 2]) * u.deg + + true_energy = ( + np.random.uniform(true_energy_bins[0].value, true_energy_bins[-1].value, size=N) + * u.TeV + ) + reco_energy = true_energy * np.random.uniform(0.5, 1.5, size=N) + + events = QTable( + { + "reco_energy": reco_energy, + "true_energy": true_energy, + "true_source_fov_offset": np.concatenate( + [ + np.full(N // 2, 0.2), + np.full(N // 2, 1.5), + ] + ) + * u.deg, + } + ) + + matrix = energy_migration_matrix( + events, true_energy_bins, reco_energy_bins, fov_offset_bins + ) + + assert matrix.shape == ( + len(true_energy_bins) - 1, + len(reco_energy_bins) - 1, + len(fov_offset_bins) - 1, + ) + assert np.allclose(matrix.sum(axis=1).max(), 1, rtol=0.1) + + def test_overflow_bins_migration_matrix(): from pyirf.irf import energy_dispersion from pyirf.irf.energy_dispersion import energy_dispersion_to_migration @@ -150,21 +203,24 @@ def test_overflow_bins_migration_matrix(): N = 10000 # add under/overflow bins - bins = 10 ** np.arange( - np.log10(0.2), - np.log10(200), - 1 / 10, - ) * u.TeV + bins = ( + 10 + ** np.arange( + np.log10(0.2), + np.log10(200), + 1 / 10, + ) + * u.TeV + ) true_energy_bins = add_overflow_bins(bins, positive=False) fov_offset_bins = np.array([0, 1, 2]) * u.deg migration_bins = np.linspace(0, 2, 101) - true_energy = np.random.uniform( - true_energy_bins[1].value, - true_energy_bins[-2].value, - size=N - ) * u.TeV + true_energy = ( + np.random.uniform(true_energy_bins[1].value, true_energy_bins[-2].value, size=N) + * u.TeV + ) reco_energy = true_energy * np.random.uniform(0.5, 1.5, size=N) selected_events = QTable( From 883274136941cb6a8e149b96874d650df0dc7976 Mon Sep 17 00:00:00 2001 From: Michele Peresano Date: Tue, 12 Dec 2023 15:10:23 +0100 Subject: [PATCH 4/5] Update test_energy_dispersion_to_migration Fix test_energy_dispersion_to_migration --- pyirf/irf/tests/test_energy_dispersion.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyirf/irf/tests/test_energy_dispersion.py b/pyirf/irf/tests/test_energy_dispersion.py index ec0e444fb..451511ffc 100644 --- a/pyirf/irf/tests/test_energy_dispersion.py +++ b/pyirf/irf/tests/test_energy_dispersion.py @@ -143,7 +143,7 @@ def test_energy_dispersion_to_migration(): assert migration_matrix.shape[2] == dispersion_matrix.shape[2] # test that all migrations are included for central energies - np.allclose(migration_matrix.sum(axis=1).max(), 1, rtol=0.1) + assert np.isclose(migration_matrix.sum(axis=1).max(), 1, rtol=0.01) # test that migrations dont always sum to 1 (since some energies are # not included in the matrix) From 4e59fb1c8138fedbcf29bb5d67bca2813b0c0ca9 Mon Sep 17 00:00:00 2001 From: Michele Peresano Date: Thu, 7 Nov 2024 11:28:39 +0100 Subject: [PATCH 5/5] fix: docstring irf.energy_dispersion.energy_migration_matrix --- pyirf/irf/energy_dispersion.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pyirf/irf/energy_dispersion.py b/pyirf/irf/energy_dispersion.py index c3d7bcac8..f34c3546d 100644 --- a/pyirf/irf/energy_dispersion.py +++ b/pyirf/irf/energy_dispersion.py @@ -103,8 +103,8 @@ def energy_migration_matrix( Returns ------- matrix : array-like - Migration matrix as probabilities along the true - energy acis with shape + Migration matrix as probabilities along the reconstructed energy axis. + energy axis with shape (n_true_energy_bins, n_reco_energy_bins, n_fov_offset_bins) containing energies in TeV. """