From 8275a433cb7270481ffb2a237bda48f3af113166 Mon Sep 17 00:00:00 2001 From: Maximilian Linhoff Date: Wed, 2 Aug 2023 13:51:10 +0900 Subject: [PATCH 01/14] Normalize edisp to integral of 1, not sum of 1 --- pyirf/irf/energy_dispersion.py | 10 ++++---- pyirf/irf/tests/test_energy_dispersion.py | 30 ++++++++++------------- 2 files changed, 18 insertions(+), 22 deletions(-) diff --git a/pyirf/irf/energy_dispersion.py b/pyirf/irf/energy_dispersion.py index 9bb203999..de631e426 100644 --- a/pyirf/irf/energy_dispersion.py +++ b/pyirf/irf/energy_dispersion.py @@ -10,9 +10,10 @@ ] -def _normalize_hist(hist): +def _normalize_hist(hist, migration_bins): # make sure we do not mutate the input array hist = hist.copy() + bin_width = np.diff(migration_bins) # calculate number of events along the N_MIGRA axis to get events # per energy per fov @@ -21,8 +22,8 @@ def _normalize_hist(hist): with np.errstate(invalid="ignore"): # hist shape is (N_E, N_MIGRA, N_FOV), norm shape is (N_E, N_FOV) # so we need to add a new axis in the middle to get (N_E, 1, N_FOV) - # for broadcasting - hist = hist / norm[:, np.newaxis, :] + # bin_width is 1d, so we need newaxis, use the values, newaxis + hist = hist / norm[:, np.newaxis, :] / bin_width[np.newaxis, :, np.newaxis] return np.nan_to_num(hist) @@ -74,8 +75,7 @@ def energy_dispersion( ], ) - n_events_per_energy = energy_dispersion.sum(axis=1) - energy_dispersion = _normalize_hist(energy_dispersion) + energy_dispersion = _normalize_hist(energy_dispersion, migration_bins) return energy_dispersion diff --git a/pyirf/irf/tests/test_energy_dispersion.py b/pyirf/irf/tests/test_energy_dispersion.py index b7dba84a5..036d106de 100644 --- a/pyirf/irf/tests/test_energy_dispersion.py +++ b/pyirf/irf/tests/test_energy_dispersion.py @@ -50,35 +50,31 @@ def test_energy_dispersion(): ) assert result.shape == (3, 1000, 2) - assert np.isclose(result.sum(), 6.0) - cumulative_sum = np.cumsum(result, axis=1) + bin_width = np.diff(migration_bins) bin_centers = 0.5 * (migration_bins[1:] + migration_bins[:-1]) + # edisp shape is (N_E, N_MIGRA, N_FOV), we need to integrate over migration axis + integral = (result * bin_width[np.newaxis, :, np.newaxis]).sum(axis=1) + np.testing.assert_allclose(integral, 1.0) + + cdf = np.cumsum(result * bin_width[np.newaxis, :, np.newaxis], axis=1) + + def ppf(cdf, bins, value): + return np.interp(value, cdf, bins[1:]) + assert np.isclose( TRUE_SIGMA_1, - ( - bin_centers[np.where(cumulative_sum[0, :, :] >= 0.84)[0][0]] - - bin_centers[np.where(cumulative_sum[0, :, :] >= 0.16)[0][0]] - ) - / 2, + 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, - ( - bin_centers[np.where(cumulative_sum[1, :, :] >= 0.84)[0][0]] - - bin_centers[np.where(cumulative_sum[1, :, :] >= 0.16)[0][0]] - ) - / 2, + 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, - ( - bin_centers[np.where(cumulative_sum[2, :, :] >= 0.84)[0][0]] - - bin_centers[np.where(cumulative_sum[2, :, :] >= 0.16)[0][0]] - ) - / 2, + 0.5 * (ppf(cdf[2, :, 0], migration_bins, 0.84) - ppf(cdf[2, :, 0], migration_bins, 0.16)), rtol=0.1, ) From a3477c9d8d78ee8de45e6925138a1788aa59b114 Mon Sep 17 00:00:00 2001 From: Maximilian Linhoff Date: Thu, 3 Aug 2023 22:59:48 +0900 Subject: [PATCH 02/14] Start working on fixing interpolation for fixed edisp --- pyirf/interpolation/base_extrapolators.py | 6 +- pyirf/interpolation/base_interpolators.py | 6 +- pyirf/interpolation/component_estimators.py | 32 +++--- .../moment_morph_interpolator.py | 66 ++++++------ .../nearest_neighbor_searcher.py | 28 +++-- .../nearest_simplex_extrapolator.py | 12 +-- pyirf/interpolation/quantile_interpolator.py | 20 ++-- .../tests/test_base_extrapolators.py | 8 +- .../tests/test_base_interpolators.py | 8 +- .../test_component_estimator_base_classes.py | 56 +++++----- .../tests/test_moment_morph_interpolator.py | 100 +++++++++--------- .../tests/test_nearest_neighbor_searcher.py | 54 +++++----- .../test_nearest_simplex_extrapolator.py | 32 +++--- .../tests/test_quantile_interpolator.py | 2 +- 14 files changed, 216 insertions(+), 214 deletions(-) diff --git a/pyirf/interpolation/base_extrapolators.py b/pyirf/interpolation/base_extrapolators.py index 79ff7a79c..f63800166 100644 --- a/pyirf/interpolation/base_extrapolators.py +++ b/pyirf/interpolation/base_extrapolators.py @@ -83,7 +83,7 @@ class DiscretePDFExtrapolator(BaseExtrapolator): Derived from pyirf.interpolation.BaseExtrapolator """ - def __init__(self, grid_points, bin_edges, bin_contents): + def __init__(self, grid_points, bin_edges, binned_pdf): """DiscretePDFExtrapolator Parameters @@ -92,7 +92,7 @@ def __init__(self, grid_points, bin_edges, bin_contents): Grid points at which templates exist bin_edges: np.ndarray, shape=(n_bins+1) Edges of the data binning - bin_content: np.ndarray, shape=(n_points, ..., n_bins) + binned_pdf: np.ndarray, shape=(n_points, ..., n_bins) Content of each bin in bin_edges for each point in grid_points. First dimesion has to correspond to number of grid_points, last dimension has to correspond to number of bins for @@ -107,4 +107,4 @@ def __init__(self, grid_points, bin_edges, bin_contents): self.bin_edges = bin_edges self.bin_mids = bin_center(self.bin_edges) - self.bin_contents = bin_contents + self.binned_pdf = binned_pdf diff --git a/pyirf/interpolation/base_interpolators.py b/pyirf/interpolation/base_interpolators.py index 0c53ed2c2..269f528e7 100644 --- a/pyirf/interpolation/base_interpolators.py +++ b/pyirf/interpolation/base_interpolators.py @@ -84,7 +84,7 @@ class DiscretePDFInterpolator(BaseInterpolator): Derived from pyirf.interpolation.BaseInterpolator """ - def __init__(self, grid_points, bin_edges, bin_contents): + def __init__(self, grid_points, bin_edges, binned_pdf): """DiscretePDFInterpolator Parameters @@ -93,7 +93,7 @@ def __init__(self, grid_points, bin_edges, bin_contents): Grid points at which interpolation templates exist bin_edges: np.ndarray, shape=(n_bins+1) Edges of the data binning - bin_content: np.ndarray, shape=(n_points, ..., n_bins) + binned_pdf: np.ndarray, shape=(n_points, ..., n_bins) Content of each bin in bin_edges for each point in grid_points. First dimesion has to correspond to number of grid_points, last dimension has to correspond to number of bins for @@ -108,4 +108,4 @@ def __init__(self, grid_points, bin_edges, bin_contents): self.bin_edges = bin_edges self.bin_mids = bin_center(self.bin_edges) - self.bin_contents = bin_contents + self.binned_pdf = binned_pdf diff --git a/pyirf/interpolation/component_estimators.py b/pyirf/interpolation/component_estimators.py index c34f0a1f7..e2caf7ef3 100644 --- a/pyirf/interpolation/component_estimators.py +++ b/pyirf/interpolation/component_estimators.py @@ -152,7 +152,7 @@ def __init__( self, grid_points, bin_edges, - bin_contents, + binned_pdf, interpolator_cls=QuantileInterpolator, interpolator_kwargs=None, extrapolator_cls=None, @@ -168,7 +168,7 @@ def __init__( Grid points at which interpolation templates exist bin_edges: np.ndarray, shape=(n_bins+1) Common set of bin-edges for all discretized PDFs. - bin_contents: np.ndarray, shape=(n_points, ..., n_bins) + binned_pdf: np.ndarray, shape=(n_points, ..., n_bins) Discretized PDFs for all grid points and arbitrary further dimensions (in IRF term e.g. field-of-view offset bins). Actual interpolation dimension, meaning the dimensions that contains actual histograms, has to be along @@ -191,16 +191,16 @@ def __init__( TypeError: When bin_edges is not a np.ndarray. TypeError: - When bin_content is not a np.ndarray.. + When binned_pdf is not a np.ndarray.. TypeError: When interpolator_cls is not a DiscretePDFInterpolator subclass. TypeError: When extrapolator_cls is not a DiscretePDFExtrapolator subclass. ValueError: - When number of bins in bin_edges and contents in bin_contents is + When number of bins in bin_edges and contents in binned_pdf is not matching. ValueError: - When number of histograms in bin_contents and points in grid_points + When number of histograms in binned_pdf and points in grid_points is not matching. Note @@ -212,20 +212,20 @@ def __init__( grid_points, ) - if not isinstance(bin_contents, np.ndarray): - raise TypeError("Input bin_contents is not a numpy array.") - elif self.n_points != bin_contents.shape[0]: + if not isinstance(binned_pdf, np.ndarray): + raise TypeError("Input binned_pdf is not a numpy array.") + elif self.n_points != binned_pdf.shape[0]: raise ValueError( f"Shape missmatch, number of grid_points ({self.n_points}) and " - f"number of histograms in bin_contents ({bin_contents.shape[0]}) " + f"number of histograms in binned_pdf ({binned_pdf.shape[0]}) " "not matching." ) elif not isinstance(bin_edges, np.ndarray): raise TypeError("Input bin_edges is not a numpy array.") - elif bin_contents.shape[-1] != (bin_edges.shape[0] - 1): + elif binned_pdf.shape[-1] != (bin_edges.shape[0] - 1): raise ValueError( f"Shape missmatch, bin_edges ({bin_edges.shape[0] - 1} bins) " - f"and bin_contents ({bin_contents.shape[-1]} bins) not matching." + f"and binned_pdf ({binned_pdf.shape[-1]} bins) not matching." ) # Make sure that 1D input is sorted in increasing order @@ -233,7 +233,7 @@ def __init__( sorting_inds = np.argsort(self.grid_points.squeeze()) self.grid_points = self.grid_points[sorting_inds] - bin_contents = bin_contents[sorting_inds] + binned_pdf = binned_pdf[sorting_inds] if interpolator_kwargs is None: interpolator_kwargs = {} @@ -247,7 +247,7 @@ def __init__( ) self.interpolator = interpolator_cls( - self.grid_points, bin_edges, bin_contents, **interpolator_kwargs + self.grid_points, bin_edges, binned_pdf, **interpolator_kwargs ) if extrapolator_cls is None: @@ -258,7 +258,7 @@ def __init__( ) else: self.extrapolator = extrapolator_cls( - self.grid_points, bin_edges, bin_contents, **extrapolator_kwargs + self.grid_points, bin_edges, binned_pdf, **extrapolator_kwargs ) @@ -596,7 +596,7 @@ def __init__( super().__init__( grid_points=grid_points, bin_edges=migra_bins, - bin_contents=np.swapaxes(energy_dispersion, axis, -1), + binned_pdf=np.swapaxes(energy_dispersion, axis, -1), interpolator_cls=interpolator_cls, interpolator_kwargs=interpolator_kwargs, extrapolator_cls=extrapolator_cls, @@ -688,7 +688,7 @@ def __init__( super().__init__( grid_points=grid_points, bin_edges=source_offset_bins, - bin_contents=psf_normed, + binned_pdf=psf_normed, interpolator_cls=interpolator_cls, interpolator_kwargs=interpolator_kwargs, extrapolator_cls=extrapolator_cls, diff --git a/pyirf/interpolation/moment_morph_interpolator.py b/pyirf/interpolation/moment_morph_interpolator.py index 0d69bc802..e94a81f60 100644 --- a/pyirf/interpolation/moment_morph_interpolator.py +++ b/pyirf/interpolation/moment_morph_interpolator.py @@ -9,16 +9,16 @@ ] -def _estimate_mean_std(bin_edges, bin_contents): +def _estimate_mean_std(bin_edges, binned_pdf): """ Function to roughly estimate mean and standard deviation from a histogram. Parameters ---------- bin_edges: np.ndarray, shape=(M+1) - Array of common bin-edges for bin_contents - bin_contents: np.ndarray, shape=(N, ..., M) - Array of bin-entries, actual + Array of common bin-edges for binned_pdf + binned_pdf: np.ndarray, shape=(N, ..., M) + PDF values from which to compute mean and std Returns ------- @@ -29,35 +29,39 @@ def _estimate_mean_std(bin_edges, bin_contents): the input template is =/= 0 """ # Create an 2darray where the 1darray mids is repeated n_template times - mids = np.broadcast_to(bin_center(bin_edges), bin_contents.shape) + mids = np.broadcast_to(bin_center(bin_edges), binned_pdf.shape) + width = np.diff(bin_edges) + + # integrate pdf to get probability in each bin + probability = binned_pdf * width # Weighted averages to compute mean and std - mean = np.average(mids, weights=bin_contents, axis=-1) + mean = np.average(mids, weights=probability, axis=-1) std = np.sqrt( - np.average((mids - mean[..., np.newaxis]) ** 2, weights=bin_contents, axis=-1) + np.average((mids - mean[..., np.newaxis])**2, weights=probability, axis=-1) ) # Set std to 0.5*width for all those templates that have only one bin =/= 0. In those # cases mids-mean = 0 and therefore std = 0. Uses the width of the one bin with - # bin_content!=0 + # binned_pdf!=0 mask = std == 0 if np.any(mask): # Create array of bin_widths and clone for each template width = np.diff(bin_edges) / 2 - width = np.repeat(width[np.newaxis, :], bin_contents.shape[0], axis=0) - std[mask] = width[bin_contents[mask, :] != 0] + width = np.repeat(width[np.newaxis, :], binned_pdf.shape[0], axis=0) + std[mask] = width[binned_pdf[mask, :] != 0] return mean, std -def _lookup(bin_edges, bin_contents, x): +def _lookup(bin_edges, binned_pdf, x): """ Function to return the bin-height at a desired point. Parameters ---------- bin_edges: np.ndarray, shape=(M+1) - Array of common bin-edges for bin_contents - bin_contents: np.ndarray, shape=(N, ..., M) + Array of common bin-edges for binned_pdf + binned_pdf: np.ndarray, shape=(N, ..., M) Array of bin-entries, actual x: numpy.ndarray, shape=(N, ..., M) Array of M points for each input template, where the histogram-value (bin-height) should be found @@ -79,11 +83,11 @@ def _lookup(bin_edges, bin_contents, x): [ cont[binnr_row] for cont, binnr_row in zip( - bin_contents.reshape(-1, bin_contents.shape[-1]), + binned_pdf.reshape(-1, binned_pdf.shape[-1]), binnr.reshape(-1, binnr.shape[-1]), ) ] - ).reshape(bin_contents.shape) + ).reshape(binned_pdf.shape) # Set all invalid bins to 0 lu[~valid] = 0 @@ -172,18 +176,18 @@ def barycentric_2D_interpolation_coefficients(grid_points, target_point): return coefficients -def moment_morph_estimation(bin_edges, bin_contents, coefficients): +def moment_morph_estimation(bin_edges, binned_pdf, coefficients): """ Function that wraps up the moment morph procedure [1] adopted for histograms. Parameters ---------- bin_edges: np.ndarray, shape=(M+1) - Array of common bin-edges for bin_contents - bin_contents: np.ndarray, shape=(N, ..., M) + Array of common bin-edges for binned_pdf + binned_pdf: np.ndarray, shape=(N, ..., M) Array of bin-entries, actual coefficients: np.ndarray, shape=(N) - Estimation coefficients for each entry in bin_contents + Estimation coefficients for each entry in binned_pdf Returns ------- @@ -199,16 +203,16 @@ def moment_morph_estimation(bin_edges, bin_contents, coefficients): bin_mids = bin_center(bin_edges) # Catch all those templates, where at least one template histogram is all zeros. - zero_templates = ~np.all(~np.isclose(np.sum(bin_contents, axis=-1), 0), 0) + zero_templates = ~np.all(~np.isclose(np.sum(binned_pdf, axis=-1), 0), 0) # Manipulate those templates so that computations pass without error - bin_contents[:, zero_templates] = np.full(len(bin_mids), 1 / len(bin_mids)) + binned_pdf[:, zero_templates] = np.full(len(bin_mids), 1 / len(bin_mids)) # Estimate mean and std for each input template histogram. First adaption needed to extend # the moment morph procedure to histograms - mus, sigs = _estimate_mean_std(bin_edges=bin_edges, bin_contents=bin_contents) + mus, sigs = _estimate_mean_std(bin_edges=bin_edges, binned_pdf=binned_pdf) coefficients = coefficients.reshape( - bin_contents.shape[0], *np.ones(mus.ndim - 1, "int") + binned_pdf.shape[0], *np.ones(mus.ndim - 1, "int") ) # Transform mean and std as in eq. (11) and (12) in [1] @@ -221,7 +225,7 @@ def moment_morph_estimation(bin_edges, bin_contents, coefficients): bij = mus - mu_prime * aij # Transformation as in eq. (13) in [1] - mids = np.broadcast_to(bin_mids, bin_contents.shape) + mids = np.broadcast_to(bin_mids, binned_pdf.shape) transf_mids = aij[..., np.newaxis] * mids + bij[..., np.newaxis] # Compute the morphed historgram according to eq. (18) in [1]. The function "lookup" "resamples" @@ -229,7 +233,7 @@ def moment_morph_estimation(bin_edges, bin_contents, coefficients): # bin-mid as new value for a whole transformed bin. Second adaption needed to extend # the moment morph procedure to histograms, adaptes the behaviour of eq. (16) - transf_hist = _lookup(bin_edges=bin_edges, bin_contents=bin_contents, x=transf_mids) + transf_hist = _lookup(bin_edges=bin_edges, binned_pdf=binned_pdf, x=transf_mids) f_new = np.sum( np.expand_dims(coefficients, -1) * transf_hist * np.expand_dims(aij, -1), axis=0 @@ -243,12 +247,12 @@ def moment_morph_estimation(bin_edges, bin_contents, coefficients): norm = np.expand_dims(np.sum(f_new, axis=-1), -1) return np.divide(f_new, norm, out=np.zeros_like(f_new), where=norm != 0).reshape( - 1, *bin_contents.shape[1:] + 1, *binned_pdf.shape[1:] ) class MomentMorphInterpolator(DiscretePDFInterpolator): - def __init__(self, grid_points, bin_edges, bin_contents): + def __init__(self, grid_points, bin_edges, binned_pdf): """ Interpolator class using moment morphing. @@ -258,7 +262,7 @@ def __init__(self, grid_points, bin_edges, bin_contents): Grid points at which interpolation templates exist. May be one ot two dimensional. bin_edges: np.ndarray, shape=(M+1) Edges of the data binning - bin_content: np.ndarray, shape=(N, ..., M) + binned_pdf: np.ndarray, shape=(N, ..., M) Content of each bin in bin_edges for each point in grid_points. First dimesion has to correspond to number of grid_points. Interpolation dimension, meaning the @@ -270,7 +274,7 @@ def __init__(self, grid_points, bin_edges, bin_contents): ---- Also calls pyirf.interpolation.DiscretePDFInterpolator.__init__. """ - super().__init__(grid_points, bin_edges, bin_contents) + super().__init__(grid_points, bin_edges, binned_pdf) if self.grid_dim == 2: self.triangulation = Delaunay(self.grid_points) @@ -293,7 +297,7 @@ def _interpolate1D(self, target_point): return moment_morph_estimation( bin_edges=self.bin_edges, - bin_contents=self.bin_contents[segment_inds], + binned_pdf=self.binned_pdf[segment_inds], coefficients=coefficients, ) @@ -312,7 +316,7 @@ def _interpolate2D(self, target_point): return moment_morph_estimation( bin_edges=self.bin_edges, - bin_contents=self.bin_contents[simplex_inds], + binned_pdf=self.binned_pdf[simplex_inds], coefficients=coefficients, ) diff --git a/pyirf/interpolation/nearest_neighbor_searcher.py b/pyirf/interpolation/nearest_neighbor_searcher.py index 8bcaaa108..4db01d14f 100644 --- a/pyirf/interpolation/nearest_neighbor_searcher.py +++ b/pyirf/interpolation/nearest_neighbor_searcher.py @@ -20,7 +20,7 @@ class BaseNearestNeighborSearcher(BaseInterpolator): actual Interpolation/Extrapolation """ - def __init__(self, grid_points, contents, norm_ord=2): + def __init__(self, grid_points, values, norm_ord=2): """ BaseNearestNeighborSearcher @@ -28,8 +28,8 @@ def __init__(self, grid_points, contents, norm_ord=2): ---------- grid_points: np.ndarray, shape=(n_points, n_dims) Grid points at which templates exist - contents: np.ndarray, shape=(n_points, ...) - Corresponding IRF contents at grid_points + values: np.ndarray, shape=(n_points, ...) + Corresponding IRF values at grid_points norm_ord: non-zero int Order of the norm which is used to compute the distances, passed to numpy.linalg.norm [1]. Defaults to 2, @@ -44,10 +44,9 @@ def __init__(self, grid_points, contents, norm_ord=2): ---- Also calls pyirf.interpolation.BaseInterpolators.__init__ """ - super().__init__(grid_points) - self.contents = contents + self.values = values # Test wether norm_ord is a number try: @@ -67,9 +66,8 @@ def __init__(self, grid_points, contents, norm_ord=2): def interpolate(self, target_point): """ - Takes a grid of IRF contents for a bunch of different parameters - and returns the contents at the nearest grid point - as seen from the target point. + Takes a grid of IRF values for a bunch of different parameters and returns + the values at the nearest grid point as seen from the target point. Parameters ---------- @@ -79,11 +77,11 @@ def interpolate(self, target_point): Returns ------- content_new: numpy.ndarray, shape=(1,...,M,...) - Contents at nearest neighbor + values at nearest neighbor Note ---- - In case of multiple nearest neighbors, the contents corresponding + In case of multiple nearest neighbors, the values corresponding to the first one are returned. """ @@ -96,7 +94,7 @@ def interpolate(self, target_point): index = np.argmin(distances) - return self.contents[index, :] + return self.values[index, :] class DiscretePDFNearestNeighborSearcher(BaseNearestNeighborSearcher): @@ -106,7 +104,7 @@ class DiscretePDFNearestNeighborSearcher(BaseNearestNeighborSearcher): Compatible with discretized PDF IRF component API. """ - def __init__(self, grid_points, bin_edges, bin_contents, norm_ord=2): + def __init__(self, grid_points, bin_edges, binned_pdf, norm_ord=2): """ NearestNeighborSearcher compatible with discretized PDF IRF components API @@ -116,7 +114,7 @@ def __init__(self, grid_points, bin_edges, bin_contents, norm_ord=2): Grid points at which templates exist bin_edges: np.ndarray, shape=(n_bins+1) Edges of the data binning. Ignored for nearest neighbor searching. - bin_content: np.ndarray, shape=(n_points, ..., n_bins) + binned_pdf: np.ndarray, shape=(n_points, ..., n_bins) Content of each bin in bin_edges for each point in grid_points. First dimesion has to correspond to number of grid_points, last dimension has to correspond to number of bins for @@ -132,7 +130,7 @@ def __init__(self, grid_points, bin_edges, bin_contents, norm_ord=2): """ super().__init__( - grid_points=grid_points, contents=bin_contents, norm_ord=norm_ord + grid_points=grid_points, values=binned_pdf, norm_ord=norm_ord ) @@ -168,7 +166,7 @@ def __init__(self, grid_points, params, norm_ord=2): Also calls pyirf.interpolation.BaseNearestNeighborSearcher.__init__ """ - super().__init__(grid_points=grid_points, contents=params, norm_ord=norm_ord) + super().__init__(grid_points=grid_points, values=params, norm_ord=norm_ord) ParametrizedInterpolator.register(ParametrizedNearestNeighborSearcher) diff --git a/pyirf/interpolation/nearest_simplex_extrapolator.py b/pyirf/interpolation/nearest_simplex_extrapolator.py index 321ac6b11..e29db3176 100644 --- a/pyirf/interpolation/nearest_simplex_extrapolator.py +++ b/pyirf/interpolation/nearest_simplex_extrapolator.py @@ -122,7 +122,7 @@ def extrapolate(self, target_point): class MomentMorphNearestSimplexExtrapolator(DiscretePDFExtrapolator): - def __init__(self, grid_points, bin_edges, bin_contents): + def __init__(self, grid_points, bin_edges, binned_pdf): """ Extrapolator class extending/reusing parts of Moment Morphing by allowing for negative extrapolation coefficients computed @@ -136,7 +136,7 @@ def __init__(self, grid_points, bin_edges, bin_contents): Have to be sorted in accending order for 1D. bin_edges: np.ndarray, shape=(M+1) Edges of the data binning - bin_content: np.ndarray, shape=(N, ..., M) + binned_pdf: np.ndarray, shape=(N, ..., M) Content of each bin in bin_edges for each point in grid_points. First dimesion has to correspond to number of grid_points. Extrapolation dimension, meaning the @@ -148,7 +148,7 @@ def __init__(self, grid_points, bin_edges, bin_contents): ---- Also calls pyirf.interpolation.DiscretePDFExtrapolator.__init__. """ - super().__init__(grid_points, bin_edges, bin_contents) + super().__init__(grid_points, bin_edges, binned_pdf) if self.grid_dim == 2: self.triangulation = Delaunay(self.grid_points) @@ -169,7 +169,7 @@ def _extrapolate1D(self, segment_inds, target_point): return moment_morph_estimation( bin_edges=self.bin_edges, - bin_contents=self.bin_contents[segment_inds], + binned_pdf=self.binned_pdf[segment_inds], coefficients=coefficients, ) @@ -185,7 +185,7 @@ def _extrapolate2D(self, simplex_inds, target_point): return moment_morph_estimation( bin_edges=self.bin_edges, - bin_contents=self.bin_contents[simplex_inds], + binned_pdf=self.binned_pdf[simplex_inds], coefficients=coefficients, ) @@ -231,6 +231,6 @@ def extrapolate(self, target_point): norm = np.expand_dims(np.sum(extrapolant, axis=-1), -1) return np.divide( extrapolant, norm, out=np.zeros_like(extrapolant), where=norm != 0 - ).reshape(1, *self.bin_contents.shape[1:]) + ).reshape(1, *self.binned_pdf.shape[1:]) else: return extrapolant diff --git a/pyirf/interpolation/quantile_interpolator.py b/pyirf/interpolation/quantile_interpolator.py index 7faeb5cb5..11e0c10e3 100644 --- a/pyirf/interpolation/quantile_interpolator.py +++ b/pyirf/interpolation/quantile_interpolator.py @@ -6,11 +6,11 @@ __all__ = ["QuantileInterpolator"] -def cdf_values(bin_contents): +def cdf_values(binned_pdf, bin_edges): """ compute cdf values and assure they are normed to unity """ - cdfs = np.cumsum(bin_contents, axis=-1) + cdfs = np.cumsum(binned_pdf * bin_edges, axis=-1) # assure the last cdf value is 1, ignore errors for empty pdfs as they are reset to 0 by nan_to_num with np.errstate(invalid="ignore"): @@ -149,7 +149,7 @@ def norm_pdf(pdf_values): class QuantileInterpolator(DiscretePDFInterpolator): - def __init__(self, grid_points, bin_edges, bin_contents, quantile_resolution=1e-3): + def __init__(self, grid_points, bin_edges, binned_pdf, quantile_resolution=1e-3): """BinnedInterpolator constructor Parameters @@ -158,7 +158,7 @@ def __init__(self, grid_points, bin_edges, bin_contents, quantile_resolution=1e- Grid points at which interpolation templates exist bin_edges: np.ndarray Edges of the data binning - bin_contents: np.ndarray + binned_pdf: np.ndarray Content of each bin in bin_edges for each point in grid_points. First dimesion has to correspond to number of grid_points, the last axis has to correspond to number @@ -170,7 +170,7 @@ def __init__(self, grid_points, bin_edges, bin_contents, quantile_resolution=1e- Raises ------ ValueError: - When last axis in bin_contents and number of bins are not equal. + When last axis in binned_pdf and number of bins are not equal. Note ---- @@ -178,7 +178,7 @@ def __init__(self, grid_points, bin_edges, bin_contents, quantile_resolution=1e- DiscretePDFInterpolator """ # Remember input shape - self.input_shape = bin_contents.shape + self.input_shape = binned_pdf.shape if self.input_shape[-1] != len(bin_edges) - 1: raise ValueError( @@ -187,8 +187,8 @@ def __init__(self, grid_points, bin_edges, bin_contents, quantile_resolution=1e- # include 0-bin at first position in each pdf to avoid edge-effects where the CDF would otherwise # start at a value != 0, also extend edges with one bin to the left - fill_zeros = np.zeros(shape=bin_contents.shape[:-1])[..., np.newaxis] - bin_contents = np.concatenate((fill_zeros, bin_contents), axis=-1) + fill_zeros = np.zeros(shape=binned_pdf.shape[:-1])[..., np.newaxis] + binned_pdf = np.concatenate((fill_zeros, binned_pdf), axis=-1) bin_edges = np.append(bin_edges[0] - np.diff(bin_edges)[0], bin_edges) # compute quantiles from quantile_resolution @@ -197,11 +197,11 @@ def __init__(self, grid_points, bin_edges, bin_contents, quantile_resolution=1e- ) super().__init__( - grid_points=grid_points, bin_edges=bin_edges, bin_contents=bin_contents + grid_points=grid_points, bin_edges=bin_edges, binned_pdf=binned_pdf ) # Compute CDF values - self.cdfs = cdf_values(self.bin_contents) + self.cdfs = cdf_values(self.binned_pdf) # compute ppf values at quantiles, determine quantile step of [1] self.ppfs = ppf_values(self.bin_mids, self.cdfs, self.quantiles) diff --git a/pyirf/interpolation/tests/test_base_extrapolators.py b/pyirf/interpolation/tests/test_base_extrapolators.py index b7997ea4a..9331a26cf 100644 --- a/pyirf/interpolation/tests/test_base_extrapolators.py +++ b/pyirf/interpolation/tests/test_base_extrapolators.py @@ -67,19 +67,19 @@ def test_DiscretePDFExtrapolator_instantiation(): target2D = np.array([[0.25, 0.25]]) bin_edges = np.linspace(-1, 1, 11) - bin_content = np.ones(shape=(len(grid_points1D), len(bin_edges) - 1)) + binned_pdf = np.ones(shape=(len(grid_points1D), len(bin_edges) - 1)) with pytest.raises(TypeError): # Abstract class, cannot be instantiated - DiscretePDFExtrapolator(grid_points1D, bin_edges, bin_content) + DiscretePDFExtrapolator(grid_points1D, bin_edges, binned_pdf) class DummyBinnedExtrapolator(DiscretePDFExtrapolator): def extrapolate(self, target_point): return 42 - interp1D = DummyBinnedExtrapolator(grid_points1D, bin_edges, bin_content) + interp1D = DummyBinnedExtrapolator(grid_points1D, bin_edges, binned_pdf) assert interp1D(target1D) == 42 - interp2D = DummyBinnedExtrapolator(grid_points2D, bin_edges, bin_content) + interp2D = DummyBinnedExtrapolator(grid_points2D, bin_edges, binned_pdf) assert interp2D(target2D) == 42 diff --git a/pyirf/interpolation/tests/test_base_interpolators.py b/pyirf/interpolation/tests/test_base_interpolators.py index c9bb63203..b2b319c3a 100644 --- a/pyirf/interpolation/tests/test_base_interpolators.py +++ b/pyirf/interpolation/tests/test_base_interpolators.py @@ -62,19 +62,19 @@ def test_DiscretePDFInterpolator_instantiation(): target2D = np.array([[1.25, 1.25]]) bin_edges = np.linspace(-1, 1, 11) - bin_content = np.ones(shape=(len(grid_points1D), len(bin_edges) - 1)) + binned_pdf = np.ones(shape=(len(grid_points1D), len(bin_edges) - 1)) with pytest.raises(TypeError): # Abstract class, cannot be instantiated - DiscretePDFInterpolator(grid_points1D, bin_edges, bin_content) + DiscretePDFInterpolator(grid_points1D, bin_edges, binned_pdf) class DummyBinnedInterpolator(DiscretePDFInterpolator): def interpolate(self, target_point, **kwargs): return 42 - interp1D = DummyBinnedInterpolator(grid_points1D, bin_edges, bin_content) + interp1D = DummyBinnedInterpolator(grid_points1D, bin_edges, binned_pdf) assert interp1D(target1D) == 42 - interp2D = DummyBinnedInterpolator(grid_points2D, bin_edges, bin_content) + interp2D = DummyBinnedInterpolator(grid_points2D, bin_edges, binned_pdf) assert interp2D(target2D) == 42 diff --git a/pyirf/interpolation/tests/test_component_estimator_base_classes.py b/pyirf/interpolation/tests/test_component_estimator_base_classes.py index 083a67fdb..534fdd209 100644 --- a/pyirf/interpolation/tests/test_component_estimator_base_classes.py +++ b/pyirf/interpolation/tests/test_component_estimator_base_classes.py @@ -239,9 +239,9 @@ def extrapolate(self, target_point): grid_points = np.array([1, 2, 3]) bin_edges = np.linspace(-1, 1, 11) - bin_content_good_shape = np.ones(shape=(len(grid_points), len(bin_edges) - 1)) - bin_content_bad_nhist = np.ones(shape=(len(grid_points) - 1, len(bin_edges) - 1)) - bin_content_bad_nbins = np.ones(shape=(len(grid_points), len(bin_edges))) + binned_pdf_good_shape = np.ones(shape=(len(grid_points), len(bin_edges) - 1)) + binned_pdf_bad_nhist = np.ones(shape=(len(grid_points) - 1, len(bin_edges) - 1)) + binned_pdf_bad_nbins = np.ones(shape=(len(grid_points), len(bin_edges))) with pytest.raises( TypeError, @@ -249,7 +249,7 @@ def extrapolate(self, target_point): ): DiscretePDFComponentEstimator( grid_points=grid_points, - bin_contents=bin_content_good_shape, + binned_pdf=binned_pdf_good_shape, bin_edges=bin_edges, interpolator_cls=WrongInterpolator, ) @@ -257,15 +257,15 @@ def extrapolate(self, target_point): with pytest.raises(TypeError, match="Input bin_edges is not a numpy array."): DiscretePDFComponentEstimator( grid_points=grid_points, - bin_contents=bin_content_good_shape, + binned_pdf=binned_pdf_good_shape, bin_edges=bin_edges.tolist(), interpolator_cls=DummyInterpolator, ) - with pytest.raises(TypeError, match="Input bin_contents is not a numpy array."): + with pytest.raises(TypeError, match="Input binned_pdf is not a numpy array."): DiscretePDFComponentEstimator( grid_points=grid_points, - bin_contents=bin_content_good_shape.tolist(), + binned_pdf=binned_pdf_good_shape.tolist(), bin_edges=bin_edges, interpolator_cls=DummyInterpolator, ) @@ -274,12 +274,12 @@ def extrapolate(self, target_point): ValueError, match=re.escape( "Shape missmatch, number of grid_points (3) and " - "number of histograms in bin_contents (2) not matching." + "number of histograms in binned_pdf (2) not matching." ), ): DiscretePDFComponentEstimator( grid_points=grid_points, - bin_contents=bin_content_bad_nhist, + binned_pdf=binned_pdf_bad_nhist, bin_edges=bin_edges, interpolator_cls=DummyInterpolator, ) @@ -288,12 +288,12 @@ def extrapolate(self, target_point): ValueError, match=re.escape( "Shape missmatch, bin_edges (10 bins) " - "and bin_contents (11 bins) not matching." + "and binned_pdf (11 bins) not matching." ), ): DiscretePDFComponentEstimator( grid_points=grid_points, - bin_contents=bin_content_bad_nbins, + binned_pdf=binned_pdf_bad_nbins, bin_edges=bin_edges, interpolator_cls=DummyInterpolator, ) @@ -305,7 +305,7 @@ def extrapolate(self, target_point): DiscretePDFComponentEstimator( grid_points=grid_points, bin_edges=bin_edges, - bin_contents=bin_content_good_shape, + binned_pdf=binned_pdf_good_shape, interpolator_cls=DummyInterpolator, extrapolator_cls=WrongExtrapolator, ) @@ -313,7 +313,7 @@ def extrapolate(self, target_point): estim = DiscretePDFComponentEstimator( grid_points=grid_points, bin_edges=bin_edges, - bin_contents=bin_content_good_shape, + binned_pdf=binned_pdf_good_shape, interpolator_cls=DummyInterpolator, extrapolator_cls=DummyExtrapolator, ) @@ -331,11 +331,11 @@ def test_DiscretePDFComponentEstimator_NearestNeighbors(): grid_points = np.array([1, 2, 3]) bin_edges = np.linspace(0, 1, 11) - bin_contents = np.array([np.full(10, x) for x in grid_points]) + binned_pdf = np.array([np.full(10, x) for x in grid_points]) estim = DiscretePDFComponentEstimator( grid_points=grid_points, - bin_contents=bin_contents, + binned_pdf=binned_pdf, bin_edges=bin_edges, interpolator_cls=DiscretePDFNearestNeighborSearcher, interpolator_kwargs={"norm_ord": 2}, @@ -343,8 +343,8 @@ def test_DiscretePDFComponentEstimator_NearestNeighbors(): extrapolator_kwargs={"norm_ord": 1}, ) - assert np.allclose(estim(target_point=np.array([1.1])), bin_contents[0, :]) - assert np.allclose(estim(target_point=np.array([4.1])), bin_contents[2, :]) + assert np.allclose(estim(target_point=np.array([1.1])), binned_pdf[0, :]) + assert np.allclose(estim(target_point=np.array([4.1])), binned_pdf[2, :]) with pytest.raises( TypeError, @@ -352,7 +352,7 @@ def test_DiscretePDFComponentEstimator_NearestNeighbors(): ): DiscretePDFComponentEstimator( grid_points=grid_points, - bin_contents=bin_contents, + binned_pdf=binned_pdf, bin_edges=bin_edges, interpolator_cls=ParametrizedNearestNeighborSearcher, interpolator_kwargs={"norm_ord": 3}, @@ -406,7 +406,7 @@ def test_DiscretePDFComponentEstimator_1Dsorting(): grid_points = np.array([[3], [1], [2]]) bin_edges = np.linspace(0, 1, 11) - bin_contents = np.array([np.full(10, x) for x in grid_points]) + binned_pdf = np.array([np.full(10, x) for x in grid_points]) class DummyInterpolator(DiscretePDFInterpolator): def interpolate(self, target_point): @@ -415,22 +415,22 @@ def interpolate(self, target_point): class DummyExtrapolator(DiscretePDFExtrapolator): def extrapolate(self, target_point): if target_point < self.grid_points.min(): - return self.bin_contents[0] + return self.binned_pdf[0] elif target_point > self.grid_points.max(): - return self.bin_contents[-1] + return self.binned_pdf[-1] estim = DiscretePDFComponentEstimator( grid_points=grid_points, - bin_contents=bin_contents, + binned_pdf=binned_pdf, bin_edges=bin_edges, interpolator_cls=DummyInterpolator, extrapolator_cls=DummyExtrapolator, ) - # Nearest neighbor is grid_point 1 at the index 1 of the original bin_contents - assert np.allclose(estim(target_point=np.array([0])), bin_contents[1, :]) - # Nearest neighbor is grid_point 3 at the index 0 of the original bin_contents - assert np.allclose(estim(target_point=np.array([4])), bin_contents[0, :]) + # Nearest neighbor is grid_point 1 at the index 1 of the original binned_pdf + assert np.allclose(estim(target_point=np.array([0])), binned_pdf[1, :]) + # Nearest neighbor is grid_point 3 at the index 0 of the original binned_pdf + assert np.allclose(estim(target_point=np.array([4])), binned_pdf[0, :]) def test_ParametrizedComponentEstimator_1Dsorting(): @@ -460,7 +460,7 @@ def extrapolate(self, target_point): extrapolator_cls=DummyExtrapolator, ) - # Nearest neighbor is grid_point 1 at the index 1 of the original bin_contents + # Nearest neighbor is grid_point 1 at the index 1 of the original binned_pdf assert np.allclose(estim(target_point=np.array([0])), params[1, :]) - # Nearest neighbor is grid_point 3 at the index 0 of the original bin_contents + # Nearest neighbor is grid_point 3 at the index 0 of the original binned_pdf assert np.allclose(estim(target_point=np.array([4])), params[0, :]) diff --git a/pyirf/interpolation/tests/test_moment_morph_interpolator.py b/pyirf/interpolation/tests/test_moment_morph_interpolator.py index 605f4310a..a452ed29c 100644 --- a/pyirf/interpolation/tests/test_moment_morph_interpolator.py +++ b/pyirf/interpolation/tests/test_moment_morph_interpolator.py @@ -29,14 +29,14 @@ def bins(): def simple_1D_data(bins): grid = np.array([[20], [40]]) target = np.array([30]) - bin_contents = np.array([binned_normal_pdf([x, 0], bins) for x in grid]) + binned_pdf = np.array([binned_normal_pdf([x, 0], bins) for x in grid]) truth = binned_normal_pdf([target, 0], bins) return { "grid": grid, "target": target, - "bin_contents": bin_contents, + "binned_pdf": binned_pdf, "truth": truth, } @@ -45,14 +45,14 @@ def simple_1D_data(bins): def simple_2D_data(bins): grid = np.array([[20, 20], [60, 20], [40, 60]]) target = np.array([25, 25]) - bin_contents = np.array([binned_normal_pdf(x, bins) for x in grid]) + binned_pdf = np.array([binned_normal_pdf(x, bins) for x in grid]) truth = binned_normal_pdf(target, bins) return { "grid": grid, "target": target, - "bin_contents": bin_contents, + "binned_pdf": binned_pdf, "truth": truth, } @@ -61,7 +61,7 @@ def test_estimate_mean_std(bins): from pyirf.interpolation.moment_morph_interpolator import _estimate_mean_std grid = np.array([[20], [40]]) - bin_contents = np.array( + binned_pdf = np.array( [ [ [binned_normal_pdf([x, 0], bins), binned_normal_pdf([x + 1, 0], bins)], @@ -94,7 +94,7 @@ def test_estimate_mean_std(bins): ] ).squeeze() - mean, std = _estimate_mean_std(bins, bin_contents) + mean, std = _estimate_mean_std(bins, binned_pdf) # Assert estimation and truth within one bin assert np.allclose(mean, true_mean, atol=np.diff(bins)[0] / 2) @@ -105,7 +105,7 @@ def test_lookup(): from pyirf.interpolation.moment_morph_interpolator import _lookup bins = np.array([0, 0.1, 0.2, 0.3, 0.4, 0.5]) - bin_contents = np.array( + binned_pdf = np.array( [ [[1, 2, 3, 4, 5], [6, 7, 8, 9, 10]], [[6, 7, 8, 9, 10], [1, 2, 3, 4, 5]], @@ -126,7 +126,7 @@ def test_lookup(): ] ) - assert np.allclose(_lookup(bins, bin_contents, x), truth) + assert np.allclose(_lookup(bins, binned_pdf, x), truth) def test_linesegment_1D_interpolation_coefficients(): @@ -212,10 +212,10 @@ def test_moment_morph_estimation1D(bins, simple_1D_data): moment_morph_estimation, ) - grid, target, bin_contents, truth = simple_1D_data.values() + grid, target, binned_pdf, truth = simple_1D_data.values() coeffs = linesegment_1D_interpolation_coefficients(grid, target) - res = moment_morph_estimation(bins, bin_contents, coeffs) + res = moment_morph_estimation(bins, binned_pdf, coeffs) assert np.isclose(np.sum(res), 1) assert np.all(np.isfinite(res)) @@ -230,10 +230,10 @@ def test_moment_morph_estimation2D(bins, simple_2D_data): moment_morph_estimation, ) - grid, target, bin_contents, truth = simple_2D_data.values() + grid, target, binned_pdf, truth = simple_2D_data.values() coeffs = barycentric_2D_interpolation_coefficients(grid, target) - res = moment_morph_estimation(bins, bin_contents, coeffs) + res = moment_morph_estimation(bins, binned_pdf, coeffs) assert np.isclose(np.sum(res), 1) assert np.all(np.isfinite(res)) @@ -245,10 +245,10 @@ def test_moment_morph_estimation2D(bins, simple_2D_data): def test_MomentMorphInterpolator1D(bins, simple_1D_data): from pyirf.interpolation import MomentMorphInterpolator - grid, target, bin_contents, truth = simple_1D_data.values() + grid, target, binned_pdf, truth = simple_1D_data.values() interp = MomentMorphInterpolator( - grid_points=grid, bin_edges=bins, bin_contents=bin_contents + grid_points=grid, bin_edges=bins, binned_pdf=binned_pdf ) res = interp(target) @@ -265,7 +265,7 @@ def test_MomentMorphInterpolator1D_dirac_delta_input(): grid = np.array([[1], [3]]) bin_edges = np.array([0, 1, 2, 3, 4]) - bin_contents = np.array( + binned_pdf = np.array( [ [[0, 1, 0, 0], [0.25, 0.25, 0.25, 0.25]], [[0, 0, 0, 1], [0.25, 0.25, 0.25, 0.25]], @@ -273,7 +273,7 @@ def test_MomentMorphInterpolator1D_dirac_delta_input(): ) target = np.array([2]) - interp = MomentMorphInterpolator(grid, bin_edges, bin_contents) + interp = MomentMorphInterpolator(grid, bin_edges, binned_pdf) res = interp(target) assert np.allclose(res, np.array([[0, 0, 1, 0], [0.25, 0.25, 0.25, 0.25]])) @@ -283,10 +283,10 @@ def test_MomentMorphInterpolator1D_all_empty(bins, simple_1D_data): from pyirf.interpolation import MomentMorphInterpolator grid, target, _, _ = simple_1D_data.values() - bin_contents = np.array([np.zeros(len(bins) - 1) for _ in grid]) + binned_pdf = np.array([np.zeros(len(bins) - 1) for _ in grid]) interp = MomentMorphInterpolator( - grid_points=grid, bin_edges=bins, bin_contents=bin_contents + grid_points=grid, bin_edges=bins, binned_pdf=binned_pdf ) res = interp(target) @@ -297,12 +297,12 @@ def test_MomentMorphInterpolator1D_all_empty(bins, simple_1D_data): def test_MomentMorphInterpolator1D_partially_empty(bins, simple_1D_data): from pyirf.interpolation import MomentMorphInterpolator - grid, target, bin_contents, _ = simple_1D_data.values() + grid, target, binned_pdf, _ = simple_1D_data.values() - bin_contents[0, :] = np.zeros(len(bins) - 1) + binned_pdf[0, :] = np.zeros(len(bins) - 1) interp = MomentMorphInterpolator( - grid_points=grid, bin_edges=bins, bin_contents=bin_contents + grid_points=grid, bin_edges=bins, binned_pdf=binned_pdf ) res = interp(target) @@ -317,7 +317,7 @@ def test_MomentMorphInterpolator1D_mixed_data(bins): target = np.array([30]) # Create template histograms - bin_contents = np.array( + binned_pdf = np.array( [ [ [binned_normal_pdf([x, 0], bins), binned_normal_pdf([x + 1, 0], bins)], @@ -331,10 +331,10 @@ def test_MomentMorphInterpolator1D_mixed_data(bins): ) # Make template histograms at indizes [:, 1, 1, :] all zeroed - bin_contents[:, 1, 1, :] = np.zeros(len(bins) - 1) + binned_pdf[:, 1, 1, :] = np.zeros(len(bins) - 1) # Zero template histogram at index [1, 0, 0, :] - bin_contents[1, 0, 0, :] = np.zeros(len(bins) - 1) + binned_pdf[1, 0, 0, :] = np.zeros(len(bins) - 1) truth = np.array( [ @@ -354,7 +354,7 @@ def test_MomentMorphInterpolator1D_mixed_data(bins): truth[1, 1, :] = np.zeros(len(bins) - 1) interp = MomentMorphInterpolator( - grid_points=grid, bin_edges=bins, bin_contents=bin_contents + grid_points=grid, bin_edges=bins, binned_pdf=binned_pdf ) res = interp(target) @@ -362,7 +362,7 @@ def test_MomentMorphInterpolator1D_mixed_data(bins): expected_norms = np.array([[0, 1], [1, 0]]) assert np.allclose(np.sum(res, axis=-1), expected_norms) assert np.all(np.isfinite(res)) - assert res.shape == (1, *bin_contents.shape[1:]) + assert res.shape == (1, *binned_pdf.shape[1:]) # Assert truth and result matching within +- 0.1%, atol dominates comparison assert np.allclose(res.squeeze(), truth, atol=1e-3, rtol=1e-5) @@ -372,7 +372,7 @@ def test_MomentMorphInterpolator1D_extended_grid_extradims(bins): grid = np.array([[20], [40], [60], [80]]) target = np.array([25]) - bin_contents = np.array( + binned_pdf = np.array( [ [ [binned_normal_pdf([x, 0], bins), binned_normal_pdf([x + 1, 0], bins)], @@ -386,7 +386,7 @@ def test_MomentMorphInterpolator1D_extended_grid_extradims(bins): ) interp = MomentMorphInterpolator( - grid_points=grid, bin_edges=bins, bin_contents=bin_contents + grid_points=grid, bin_edges=bins, binned_pdf=binned_pdf ) truth = np.array( @@ -406,7 +406,7 @@ def test_MomentMorphInterpolator1D_extended_grid_extradims(bins): assert np.allclose(np.sum(res, axis=-1), 1) assert np.all(np.isfinite(res)) - assert res.shape == (1, *bin_contents.shape[1:]) + assert res.shape == (1, *binned_pdf.shape[1:]) # Assert truth and result matching within +- 0.1%, atol dominates comparison assert np.allclose(res.squeeze(), truth, atol=1e-4, rtol=1e-5) @@ -414,10 +414,10 @@ def test_MomentMorphInterpolator1D_extended_grid_extradims(bins): def test_MomentMorphInterpolator2D(bins, simple_2D_data): from pyirf.interpolation import MomentMorphInterpolator - grid, target, bin_contents, truth = simple_2D_data.values() + grid, target, binned_pdf, truth = simple_2D_data.values() interp = MomentMorphInterpolator( - grid_points=grid, bin_edges=bins, bin_contents=bin_contents + grid_points=grid, bin_edges=bins, binned_pdf=binned_pdf ) res = interp(target) @@ -432,12 +432,12 @@ def test_MomentMorphInterpolator2D(bins, simple_2D_data): def test_MomentMorphInterpolator2D_partially_empty(bins, simple_2D_data): from pyirf.interpolation import MomentMorphInterpolator - grid, target, bin_contents, _ = simple_2D_data.values() + grid, target, binned_pdf, _ = simple_2D_data.values() - bin_contents[0, :] = np.zeros(len(bins) - 1) + binned_pdf[0, :] = np.zeros(len(bins) - 1) interp = MomentMorphInterpolator( - grid_points=grid, bin_edges=bins, bin_contents=bin_contents + grid_points=grid, bin_edges=bins, binned_pdf=binned_pdf ) res = interp(target) @@ -449,10 +449,10 @@ def test_MomentMorphInterpolator2D_all_empty(bins, simple_2D_data): from pyirf.interpolation import MomentMorphInterpolator grid, target, _, _ = simple_2D_data.values() - bin_contents = np.array([np.zeros(len(bins) - 1) for _ in grid]) + binned_pdf = np.array([np.zeros(len(bins) - 1) for _ in grid]) interp = MomentMorphInterpolator( - grid_points=grid, bin_edges=bins, bin_contents=bin_contents + grid_points=grid, bin_edges=bins, binned_pdf=binned_pdf ) res = interp(target) @@ -467,7 +467,7 @@ def test_MomentMorphInterpolator2D_mixed(bins): target = np.array([25, 25]) # Create template histograms - bin_contents = np.array( + binned_pdf = np.array( [ [ [binned_normal_pdf(x, bins), binned_normal_pdf([x[0] + 1, x[1]], bins)], @@ -481,10 +481,10 @@ def test_MomentMorphInterpolator2D_mixed(bins): ) # Make template histograms at indizes [:, 1, 1, :] all zeroed - bin_contents[:, 1, 1, :] = np.zeros(len(bins) - 1) + binned_pdf[:, 1, 1, :] = np.zeros(len(bins) - 1) # Zero template histogram at index [1, 0, 0, :] - bin_contents[1, 0, 0, :] = np.zeros(len(bins) - 1) + binned_pdf[1, 0, 0, :] = np.zeros(len(bins) - 1) truth = np.array( [ @@ -504,7 +504,7 @@ def test_MomentMorphInterpolator2D_mixed(bins): truth[1, 1, :] = np.zeros(len(bins) - 1) interp = MomentMorphInterpolator( - grid_points=grid, bin_edges=bins, bin_contents=bin_contents + grid_points=grid, bin_edges=bins, binned_pdf=binned_pdf ) res = interp(target) @@ -512,7 +512,7 @@ def test_MomentMorphInterpolator2D_mixed(bins): expected_norms = np.array([[0, 1], [1, 0]]) assert np.allclose(np.sum(res, axis=-1), expected_norms) assert np.all(np.isfinite(res)) - assert res.shape == (1, *bin_contents.shape[1:]) + assert res.shape == (1, *binned_pdf.shape[1:]) # Assert truth and result matching within +- 0.1%, atol dominates comparison assert np.allclose(res.squeeze(), truth, atol=1e-3, rtol=1e-5) @@ -522,12 +522,12 @@ def test_MomentMorphInterpolator1D_extended_grid(bins): grid = np.array([[20], [40], [60], [80]]) target = np.array([25]) - bin_contents = np.array([binned_normal_pdf([x, 0], bins) for x in grid]) + binned_pdf = np.array([binned_normal_pdf([x, 0], bins) for x in grid]) interp = MomentMorphInterpolator( grid_points=grid, bin_edges=bins, - bin_contents=bin_contents, + binned_pdf=binned_pdf, ) res = interp(target) @@ -545,12 +545,12 @@ def test_MomentMorphInterpolator2D_extended_grid(bins): grid = np.array([[20, 20], [40, 20], [30, 40], [50, 20], [45, 40]]) target = np.array([25, 25]) - bin_contents = np.array([binned_normal_pdf(x, bins) for x in grid]) + binned_pdf = np.array([binned_normal_pdf(x, bins) for x in grid]) interp = MomentMorphInterpolator( grid_points=grid, bin_edges=bins, - bin_contents=bin_contents, + binned_pdf=binned_pdf, ) res = interp(target) @@ -568,7 +568,7 @@ def test_MomentMorphInterpolator2D_extended_grid_extradims(bins): grid = np.array([[20, 20], [40, 20], [30, 40], [50, 20], [45, 40]]) target = np.array([25, 25]) - bin_contents = np.array( + binned_pdf = np.array( [ [ [binned_normal_pdf(x, bins), binned_normal_pdf([x[0] + 1, x[1]], bins)], @@ -584,7 +584,7 @@ def test_MomentMorphInterpolator2D_extended_grid_extradims(bins): interp = MomentMorphInterpolator( grid_points=grid, bin_edges=bins, - bin_contents=bin_contents, + binned_pdf=binned_pdf, ) truth = np.array( @@ -604,7 +604,7 @@ def test_MomentMorphInterpolator2D_extended_grid_extradims(bins): assert np.allclose(np.sum(res, axis=-1), 1) assert np.all(np.isfinite(res)) - assert res.shape == (1, *bin_contents.shape[1:]) + assert res.shape == (1, *binned_pdf.shape[1:]) # Assert truth and result matching within +- 0.1%, atol dominates comparison assert np.allclose(res.squeeze(), truth, atol=1e-3, rtol=1e-5) @@ -616,7 +616,7 @@ def test_MomentMorphInterpolator3D(): grid = np.array([[0, 0, 0], [0, 20, 0], [20, 0, 0], [20, 20, 0], [10, 10, 10]]) - bin_contents = np.array([np.ones(len(bins) - 1) / (len(bins) - 1) for _ in grid]) + binned_pdf = np.array([np.ones(len(bins) - 1) / (len(bins) - 1) for _ in grid]) with pytest.raises( NotImplementedError, @@ -625,5 +625,5 @@ def test_MomentMorphInterpolator3D(): MomentMorphInterpolator( grid_points=grid, bin_edges=bins, - bin_contents=bin_contents, + binned_pdf=binned_pdf, ) diff --git a/pyirf/interpolation/tests/test_nearest_neighbor_searcher.py b/pyirf/interpolation/tests/test_nearest_neighbor_searcher.py index 66c25f45d..6394b8ae4 100644 --- a/pyirf/interpolation/tests/test_nearest_neighbor_searcher.py +++ b/pyirf/interpolation/tests/test_nearest_neighbor_searcher.py @@ -13,7 +13,7 @@ def grid_2d(): @pytest.fixture -def contents(grid_1d): +def values(grid_1d): return np.array( [ [[np.full(10, x), np.full(10, x / 2)], [np.full(10, 2 * x), np.zeros(10)]] @@ -22,83 +22,83 @@ def contents(grid_1d): ) -def test_BaseNearestNeighborSearcher_1DGrid(grid_1d, contents): +def test_BaseNearestNeighborSearcher_1DGrid(grid_1d, binned_pdf): from pyirf.interpolation import BaseNearestNeighborSearcher - searcher = BaseNearestNeighborSearcher(grid_1d, contents, norm_ord=2) + searcher = BaseNearestNeighborSearcher(grid_1d, binned_pdf, norm_ord=2) target = np.array([0]) - assert np.array_equal(searcher(target), contents[0, :]) + assert np.array_equal(searcher(target), binned_pdf[0, :]) target = np.array([[1.9]]) - assert np.array_equal(searcher(target), contents[1, :]) + assert np.array_equal(searcher(target), binned_pdf[1, :]) -def test_BaseNearestNeighborSearcher_2DGrid(grid_2d, contents): +def test_BaseNearestNeighborSearcher_2DGrid(grid_2d, binned_pdf): from pyirf.interpolation import BaseNearestNeighborSearcher - searcher = BaseNearestNeighborSearcher(grid_2d, contents, norm_ord=2) + searcher = BaseNearestNeighborSearcher(grid_2d, binned_pdf, norm_ord=2) target = np.array([[0, 1]]) - assert np.array_equal(searcher(target), contents[0, :]) + assert np.array_equal(searcher(target), binned_pdf[0, :]) target = np.array([3, 3]) - assert np.array_equal(searcher(target), contents[-1, :]) + assert np.array_equal(searcher(target), binned_pdf[-1, :]) -def test_BaseNearestNeighborSearcher_manhatten_norm(grid_2d, contents): +def test_BaseNearestNeighborSearcher_manhatten_norm(grid_2d, binned_pdf): from pyirf.interpolation import BaseNearestNeighborSearcher - searcher = BaseNearestNeighborSearcher(grid_2d, contents, norm_ord=1) + searcher = BaseNearestNeighborSearcher(grid_2d, binned_pdf, norm_ord=1) target = np.array([[0, 1]]) - assert np.array_equal(searcher(target), contents[0, :]) + assert np.array_equal(searcher(target), binned_pdf[0, :]) target = np.array([[3, 3]]) - assert np.array_equal(searcher(target), contents[-1, :]) + assert np.array_equal(searcher(target), binned_pdf[-1, :]) -def test_BaseNearestNeighborSearcher_wrong_norm(grid_1d, contents): +def test_BaseNearestNeighborSearcher_wrong_norm(grid_1d, binned_pdf): from pyirf.interpolation import BaseNearestNeighborSearcher with pytest.raises(ValueError, match="Only positiv integers allowed for norm_ord"): - BaseNearestNeighborSearcher(grid_1d, contents, norm_ord=-2) + BaseNearestNeighborSearcher(grid_1d, binned_pdf, norm_ord=-2) with pytest.raises(ValueError, match="Only positiv integers allowed for norm_ord"): - BaseNearestNeighborSearcher(grid_1d, contents, norm_ord=1.5) + BaseNearestNeighborSearcher(grid_1d, binned_pdf, norm_ord=1.5) with pytest.raises(ValueError, match="Only positiv integers allowed for norm_ord"): - BaseNearestNeighborSearcher(grid_1d, contents, norm_ord=np.inf) + BaseNearestNeighborSearcher(grid_1d, binned_pdf, norm_ord=np.inf) with pytest.raises(ValueError, match="Only positiv integers allowed for norm_ord"): - BaseNearestNeighborSearcher(grid_1d, contents, norm_ord="nuc") + BaseNearestNeighborSearcher(grid_1d, binned_pdf, norm_ord="nuc") -def test_DiscretePDFNearestNeighborSearcher(grid_2d, contents): +def test_DiscretePDFNearestNeighborSearcher(grid_2d, binned_pdf): from pyirf.interpolation import DiscretePDFNearestNeighborSearcher - bin_edges = np.linspace(0, 1, contents.shape[-1] + 1) + bin_edges = np.linspace(0, 1, binned_pdf.shape[-1] + 1) searcher = DiscretePDFNearestNeighborSearcher( - grid_points=grid_2d, bin_edges=bin_edges, bin_contents=contents, norm_ord=1 + grid_points=grid_2d, bin_edges=bin_edges, binned_pdf=binned_pdf, norm_ord=1 ) target = np.array([[0, 1]]) - assert np.array_equal(searcher(target), contents[0, :]) + assert np.array_equal(searcher(target), binned_pdf[0, :]) target = np.array([[3, 3]]) - assert np.array_equal(searcher(target), contents[-1, :]) + assert np.array_equal(searcher(target), binned_pdf[-1, :]) -def test_ParametrizedNearestNeighborSearcher(grid_2d, contents): +def test_ParametrizedNearestNeighborSearcher(grid_2d, binned_pdf): from pyirf.interpolation import ParametrizedNearestNeighborSearcher searcher = ParametrizedNearestNeighborSearcher( - grid_points=grid_2d, params=contents, norm_ord=1 + grid_points=grid_2d, params=binned_pdf, norm_ord=1 ) target = np.array([[0, 1]]) - assert np.array_equal(searcher(target), contents[0, :]) + assert np.array_equal(searcher(target), binned_pdf[0, :]) target = np.array([[3, 3]]) - assert np.array_equal(searcher(target), contents[-1, :]) + assert np.array_equal(searcher(target), binned_pdf[-1, :]) diff --git a/pyirf/interpolation/tests/test_nearest_simplex_extrapolator.py b/pyirf/interpolation/tests/test_nearest_simplex_extrapolator.py index 2ffc5bd8c..389cfa2f7 100644 --- a/pyirf/interpolation/tests/test_nearest_simplex_extrapolator.py +++ b/pyirf/interpolation/tests/test_nearest_simplex_extrapolator.py @@ -138,7 +138,7 @@ def test_MomentMorphNearestSimplexExtrapolator_1DGrid(bins): grid_points = np.array([[20], [30], [40]]) # Create template histograms - bin_contents = np.array( + binned_pdf = np.array( [ [ [binned_normal_pdf([x, 0], bins), binned_normal_pdf([x + 1, 0], bins)], @@ -153,16 +153,16 @@ def test_MomentMorphNearestSimplexExtrapolator_1DGrid(bins): ) # Include dirac-delta moving with grid_point x at [:, 2, 0, x] - bin_contents[0, 2, 0, 20] = 1 - bin_contents[1, 2, 0, 30] = 1 - bin_contents[2, 2, 0, 40] = 1 + binned_pdf[0, 2, 0, 20] = 1 + binned_pdf[1, 2, 0, 30] = 1 + binned_pdf[2, 2, 0, 40] = 1 # Zero template histogram at index [0, 0, 0, :], extrapolations "lower" # then this bin has consequently to be zero in these bins. - bin_contents[0, 0, 0, :] = np.zeros(len(bins) - 1) + binned_pdf[0, 0, 0, :] = np.zeros(len(bins) - 1) extrap = MomentMorphNearestSimplexExtrapolator( - grid_points=grid_points, bin_contents=bin_contents, bin_edges=bins + grid_points=grid_points, binned_pdf=binned_pdf, bin_edges=bins ) target1 = np.array([10]) @@ -190,7 +190,7 @@ def test_MomentMorphNearestSimplexExtrapolator_1DGrid(bins): expected_norms1 = np.array([[0, 1], [1, 0], [1, 1]]) assert np.allclose(np.sum(res1, axis=-1), expected_norms1) assert np.all(np.isfinite(res1)) - assert res1.shape == (1, *bin_contents.shape[1:]) + assert res1.shape == (1, *binned_pdf.shape[1:]) # Assert truth and result matching within +- 0.1%, atol dominates comparison assert np.allclose(res1.squeeze(), truth1, atol=1e-3, rtol=1e-5) @@ -218,7 +218,7 @@ def test_MomentMorphNearestSimplexExtrapolator_1DGrid(bins): expected_norms2 = np.array([[1, 1], [1, 0], [1, 1]]) assert np.allclose(np.sum(res2, axis=-1), expected_norms2) assert np.all(np.isfinite(res2)) - assert res2.shape == (1, *bin_contents.shape[1:]) + assert res2.shape == (1, *binned_pdf.shape[1:]) # Assert truth and result matching within +- 0.1%, atol dominates comparison assert np.allclose(res2.squeeze(), truth2, atol=1e-3, rtol=1e-5) @@ -230,7 +230,7 @@ def test_MomentMorphNearestSimplexExtrapolator_2DGrid(bins): grid_points = np.array([[20, 20], [30, 10], [40, 20], [50, 10]]) # Create template histograms - bin_contents = np.array( + binned_pdf = np.array( [ [ [binned_normal_pdf(x, bins), binned_normal_pdf([x[0] + 1, x[1]], bins)], @@ -245,10 +245,10 @@ def test_MomentMorphNearestSimplexExtrapolator_2DGrid(bins): # Zero template histogram at index [0, 0, 0, :] (the [20, 20] point), # so that targets interpolated from left simplex should also be zero at [0, 0, :] - bin_contents[0, 0, 0, :] = np.zeros(len(bins) - 1) + binned_pdf[0, 0, 0, :] = np.zeros(len(bins) - 1) extrap = MomentMorphNearestSimplexExtrapolator( - grid_points=grid_points, bin_contents=bin_contents, bin_edges=bins + grid_points=grid_points, binned_pdf=binned_pdf, bin_edges=bins ) target1 = np.array([25, 10]) @@ -272,7 +272,7 @@ def test_MomentMorphNearestSimplexExtrapolator_2DGrid(bins): expected_norms1 = np.array([[0, 1], [1, 0]]) assert np.allclose(np.sum(res1, axis=-1), expected_norms1) assert np.all(np.isfinite(res1)) - assert res1.shape == (1, *bin_contents.shape[1:]) + assert res1.shape == (1, *binned_pdf.shape[1:]) # Assert truth and result matching within +- 0.1%, atol dominates comparison assert np.allclose(res1.squeeze(), truth1, atol=1e-3, rtol=1e-5) @@ -296,7 +296,7 @@ def test_MomentMorphNearestSimplexExtrapolator_2DGrid(bins): expected_norms2 = np.array([[1, 1], [1, 0]]) assert np.allclose(np.sum(res2, axis=-1), expected_norms2) assert np.all(np.isfinite(res2)) - assert res2.shape == (1, *bin_contents.shape[1:]) + assert res2.shape == (1, *binned_pdf.shape[1:]) # Assert truth and result matching within +- 0.1%, atol dominates comparison assert np.allclose(res2.squeeze(), truth2, atol=1e-3, rtol=1e-5) @@ -315,7 +315,7 @@ def test_MomentMorphNearestSimplexExtrapolator_3DGrid(): match="Extrapolation in more then two dimension not impemented.", ): MomentMorphNearestSimplexExtrapolator( - grid_points=grid_points, bin_contents=dummy_data, bin_edges=bin_edges + grid_points=grid_points, binned_pdf=dummy_data, bin_edges=bin_edges ) @@ -330,7 +330,7 @@ def test_MomentMorphNearestSimplexExtrapolator_below_zero_warning(bins): grid_points = np.array([[30], [40]]) bins = np.linspace(-10, 40, 51) - bin_contents = np.array( + binned_pdf = np.array( [ [binned_normal_pdf([x, 0], bins), binned_normal_pdf([x + 10, 0], bins)] for x in grid_points @@ -338,7 +338,7 @@ def test_MomentMorphNearestSimplexExtrapolator_below_zero_warning(bins): ) extrap = MomentMorphNearestSimplexExtrapolator( - grid_points=grid_points, bin_contents=bin_contents, bin_edges=bins + grid_points=grid_points, binned_pdf=binned_pdf, bin_edges=bins ) with pytest.warns( diff --git a/pyirf/interpolation/tests/test_quantile_interpolator.py b/pyirf/interpolation/tests/test_quantile_interpolator.py index 195d95fdc..a5a9cc3fd 100644 --- a/pyirf/interpolation/tests/test_quantile_interpolator.py +++ b/pyirf/interpolation/tests/test_quantile_interpolator.py @@ -96,7 +96,7 @@ def test_interpolate_binned_pdf(data): interpolator = QuantileInterpolator( grid_points=data["grid_points"][[0, 2]], bin_edges=data["bin_edges"], - bin_contents=data["binned_pdfs"][[0, 2], :], + binned_pdf=data["binned_pdfs"][[0, 2], :], quantile_resolution=1e-3, ) From b9b2081e4270dbf7fc0bbc326a4766bebba48e6a Mon Sep 17 00:00:00 2001 From: Maximilian Linhoff Date: Mon, 21 Aug 2023 12:52:56 +0200 Subject: [PATCH 03/14] Add normalization arg to PDF interpolator class --- pyirf/interpolation/base_interpolators.py | 22 +++++++++--- .../moment_morph_interpolator.py | 6 ++-- pyirf/interpolation/quantile_interpolator.py | 35 +++++++++++++------ 3 files changed, 45 insertions(+), 18 deletions(-) diff --git a/pyirf/interpolation/base_interpolators.py b/pyirf/interpolation/base_interpolators.py index 269f528e7..9b4227f38 100644 --- a/pyirf/interpolation/base_interpolators.py +++ b/pyirf/interpolation/base_interpolators.py @@ -1,5 +1,6 @@ """Base classes for interpolators""" from abc import ABCMeta, abstractmethod +import enum import numpy as np from pyirf.binning import bin_center @@ -7,6 +8,15 @@ __all__ = ["BaseInterpolator", "ParametrizedInterpolator", "DiscretePDFInterpolator"] +class PDFNormalization(enum.Enum): + """How a discrete PDF is normalized""" + #: PDF is normalized to a "normal" area integral of 1 + AREA = enum.auto() + #: PDF is normalized to 1 over the solid angle integral where the bin + #: edges represent the opening angles of cones in radean. + CONE_SOLID_ANGLE = enum.auto() + + class BaseInterpolator(metaclass=ABCMeta): """ Base class for all interpolators, only knowing grid-points, @@ -84,21 +94,22 @@ class DiscretePDFInterpolator(BaseInterpolator): Derived from pyirf.interpolation.BaseInterpolator """ - def __init__(self, grid_points, bin_edges, binned_pdf): + def __init__(self, grid_points, bin_edges, binned_pdf, normalization=PDFNormalization.AREA): """DiscretePDFInterpolator Parameters ---------- - grid_points: np.ndarray, shape=(n_points, n_dims) + grid_points : np.ndarray, shape=(n_points, n_dims) Grid points at which interpolation templates exist - bin_edges: np.ndarray, shape=(n_bins+1) + bin_edges : np.ndarray, shape=(n_bins+1) Edges of the data binning - binned_pdf: np.ndarray, shape=(n_points, ..., n_bins) + binned_pdf : np.ndarray, shape=(n_points, ..., n_bins) Content of each bin in bin_edges for each point in grid_points. First dimesion has to correspond to number of grid_points, last dimension has to correspond to number of bins for the quantity that should be interpolated (e.g. the Migra axis for EDisp) - + normalization : PDFNormalization + How the PDF is normalized Note ---- @@ -109,3 +120,4 @@ def __init__(self, grid_points, bin_edges, binned_pdf): self.bin_edges = bin_edges self.bin_mids = bin_center(self.bin_edges) self.binned_pdf = binned_pdf + self.normalization = normalization diff --git a/pyirf/interpolation/moment_morph_interpolator.py b/pyirf/interpolation/moment_morph_interpolator.py index e94a81f60..fafcfccb4 100644 --- a/pyirf/interpolation/moment_morph_interpolator.py +++ b/pyirf/interpolation/moment_morph_interpolator.py @@ -2,7 +2,7 @@ from pyirf.binning import bin_center, calculate_bin_indices from scipy.spatial import Delaunay -from .base_interpolators import DiscretePDFInterpolator +from .base_interpolators import DiscretePDFInterpolator, PDFNormalization __all__ = [ "MomentMorphInterpolator", @@ -252,7 +252,7 @@ def moment_morph_estimation(bin_edges, binned_pdf, coefficients): class MomentMorphInterpolator(DiscretePDFInterpolator): - def __init__(self, grid_points, bin_edges, binned_pdf): + def __init__(self, grid_points, bin_edges, binned_pdf, normalization=PDFNormalization.AREA): """ Interpolator class using moment morphing. @@ -274,7 +274,7 @@ def __init__(self, grid_points, bin_edges, binned_pdf): ---- Also calls pyirf.interpolation.DiscretePDFInterpolator.__init__. """ - super().__init__(grid_points, bin_edges, binned_pdf) + super().__init__(grid_points, bin_edges, binned_pdf, normalization) if self.grid_dim == 2: self.triangulation = Delaunay(self.grid_points) diff --git a/pyirf/interpolation/quantile_interpolator.py b/pyirf/interpolation/quantile_interpolator.py index 11e0c10e3..8170107ba 100644 --- a/pyirf/interpolation/quantile_interpolator.py +++ b/pyirf/interpolation/quantile_interpolator.py @@ -1,16 +1,26 @@ import numpy as np from scipy.interpolate import griddata, interp1d +import enum -from .base_interpolators import DiscretePDFInterpolator +from ..utils import cone_solid_angle +from .base_interpolators import DiscretePDFInterpolator, PDFNormalization __all__ = ["QuantileInterpolator"] -def cdf_values(binned_pdf, bin_edges): +def cdf_values(binned_pdf, bin_edges, normalization): """ compute cdf values and assure they are normed to unity """ - cdfs = np.cumsum(binned_pdf * bin_edges, axis=-1) + + if normalization is PDFNormalization.AREA: + bin_widths = np.diff(bin_edges) + elif normalization is PDFNormalization.CONE_SOLID_ANGLE: + bin_widths = np.diff(cone_solid_angle(bin_edges)) + else: + raise ValueError(f"Invalid PDF normalization: {normalization}") + + cdfs = np.cumsum(binned_pdf * bin_widths, axis=-1) # assure the last cdf value is 1, ignore errors for empty pdfs as they are reset to 0 by nan_to_num with np.errstate(invalid="ignore"): @@ -149,23 +159,25 @@ def norm_pdf(pdf_values): class QuantileInterpolator(DiscretePDFInterpolator): - def __init__(self, grid_points, bin_edges, binned_pdf, quantile_resolution=1e-3): + def __init__(self, grid_points, bin_edges, binned_pdf, quantile_resolution=1e-3, normalization=PDFNormalization.AREA): """BinnedInterpolator constructor Parameters ---------- - grid_points: np.ndarray + grid_points : np.ndarray Grid points at which interpolation templates exist - bin_edges: np.ndarray + bin_edges : np.ndarray Edges of the data binning - binned_pdf: np.ndarray + binned_pdf : np.ndarray Content of each bin in bin_edges for each point in grid_points. First dimesion has to correspond to number of grid_points, the last axis has to correspond to number of bins for the quantity that should be interpolated (e.g. the Migra axis for EDisp) - quantile_resolution: float + quantile_resolution : float Spacing between quantiles + normalization : PDFNormalization + How the discrete PDF is normalized Raises ------ @@ -197,11 +209,14 @@ def __init__(self, grid_points, bin_edges, binned_pdf, quantile_resolution=1e-3) ) super().__init__( - grid_points=grid_points, bin_edges=bin_edges, binned_pdf=binned_pdf + grid_points=grid_points, + bin_edges=bin_edges, + binned_pdf=binned_pdf, + normalization=normalization, ) # Compute CDF values - self.cdfs = cdf_values(self.binned_pdf) + self.cdfs = cdf_values(self.binned_pdf, self.bin_edges, self.normalization) # compute ppf values at quantiles, determine quantile step of [1] self.ppfs = ppf_values(self.bin_mids, self.cdfs, self.quantiles) From bb66247d1d43ea3c5cae201ff7f36fed634b2e27 Mon Sep 17 00:00:00 2001 From: Maximilian Linhoff Date: Mon, 21 Aug 2023 16:19:13 +0200 Subject: [PATCH 04/14] Change binned_normal_pdf to have integral 1 --- pyirf/interpolation/tests/__init__.py | 0 .../tests/test_moment_morph_interpolator.py | 98 +++++++++---------- .../test_nearest_simplex_extrapolator.py | 47 +++++---- 3 files changed, 69 insertions(+), 76 deletions(-) create mode 100644 pyirf/interpolation/tests/__init__.py diff --git a/pyirf/interpolation/tests/__init__.py b/pyirf/interpolation/tests/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/pyirf/interpolation/tests/test_moment_morph_interpolator.py b/pyirf/interpolation/tests/test_moment_morph_interpolator.py index a452ed29c..715597aab 100644 --- a/pyirf/interpolation/tests/test_moment_morph_interpolator.py +++ b/pyirf/interpolation/tests/test_moment_morph_interpolator.py @@ -11,13 +11,13 @@ def expected_std(a, b): return 1 + 0.05 * (a + b) -def binned_normal_pdf(mean_std_args, bins): - pdf = np.diff( - norm(loc=expected_mean(*mean_std_args), scale=expected_std(*mean_std_args)).cdf( - bins - ) +def binned_normal_pdf(a, b, bins): + dist = norm( + loc=expected_mean(a, b), + scale=expected_std(a, b), ) - return pdf / pdf.sum() + pdf = np.diff(dist.cdf(bins)) + return pdf / np.diff(bins) @pytest.fixture @@ -31,7 +31,7 @@ def simple_1D_data(bins): target = np.array([30]) binned_pdf = np.array([binned_normal_pdf([x, 0], bins) for x in grid]) - truth = binned_normal_pdf([target, 0], bins) + truth = binned_normal_pdf(target, 0, bins) return { "grid": grid, @@ -45,9 +45,9 @@ def simple_1D_data(bins): def simple_2D_data(bins): grid = np.array([[20, 20], [60, 20], [40, 60]]) target = np.array([25, 25]) - binned_pdf = np.array([binned_normal_pdf(x, bins) for x in grid]) + binned_pdf = np.array([binned_normal_pdf(*x, bins) for x in grid]) - truth = binned_normal_pdf(target, bins) + truth = binned_normal_pdf(*target, bins) return { "grid": grid, @@ -64,10 +64,10 @@ def test_estimate_mean_std(bins): binned_pdf = np.array( [ [ - [binned_normal_pdf([x, 0], bins), binned_normal_pdf([x + 1, 0], bins)], + [binned_normal_pdf(x, 0, bins), binned_normal_pdf(x + 1, 0, bins)], [ - binned_normal_pdf([x + 1.5, 0], bins), - binned_normal_pdf([x + 10, 0], bins), + binned_normal_pdf(x + 1.5, 0, bins), + binned_normal_pdf(x + 10, 0, bins), ], ] for x in grid @@ -320,11 +320,8 @@ def test_MomentMorphInterpolator1D_mixed_data(bins): binned_pdf = np.array( [ [ - [binned_normal_pdf([x, 0], bins), binned_normal_pdf([x + 1, 0], bins)], - [ - binned_normal_pdf([x + 10, 0], bins), - binned_normal_pdf([x + 2, 0], bins), - ], + [binned_normal_pdf(x, 0, bins), binned_normal_pdf(x + 1, 0, bins)], + [binned_normal_pdf(x + 10, 0, bins), binned_normal_pdf(x + 2, 0, bins)], ] for x in grid ] @@ -339,12 +336,12 @@ def test_MomentMorphInterpolator1D_mixed_data(bins): truth = np.array( [ [ - binned_normal_pdf([target, 0], bins), - binned_normal_pdf([target + 1, 0], bins), + binned_normal_pdf(target, 0, bins), + binned_normal_pdf(target + 1, 0, bins), ], [ - binned_normal_pdf([target + 10, 0], bins), - binned_normal_pdf([target + 2, 0], bins), + binned_normal_pdf(target + 10, 0, bins), + binned_normal_pdf(target + 2, 0, bins), ], ] ) @@ -375,10 +372,10 @@ def test_MomentMorphInterpolator1D_extended_grid_extradims(bins): binned_pdf = np.array( [ [ - [binned_normal_pdf([x, 0], bins), binned_normal_pdf([x + 1, 0], bins)], + [binned_normal_pdf(x, 0, bins), binned_normal_pdf(x + 1, 0, bins)], [ - binned_normal_pdf([x + 10, 0], bins), - binned_normal_pdf([x + 2, 0], bins), + binned_normal_pdf(x + 10, 0, bins), + binned_normal_pdf(x + 2, 0, bins), ], ] for x in grid @@ -392,12 +389,12 @@ def test_MomentMorphInterpolator1D_extended_grid_extradims(bins): truth = np.array( [ [ - binned_normal_pdf([target, 0], bins), - binned_normal_pdf([target + 1, 0], bins), + binned_normal_pdf(target, 0, bins), + binned_normal_pdf(target + 1, 0, bins), ], [ - binned_normal_pdf([target + 10, 0], bins), - binned_normal_pdf([target + 2, 0], bins), + binned_normal_pdf(target + 10, 0, bins), + binned_normal_pdf(target + 2, 0, bins), ], ] ) @@ -470,13 +467,10 @@ def test_MomentMorphInterpolator2D_mixed(bins): binned_pdf = np.array( [ [ - [binned_normal_pdf(x, bins), binned_normal_pdf([x[0] + 1, x[1]], bins)], - [ - binned_normal_pdf([x[0] + 10, x[1]], bins), - binned_normal_pdf([x[0] + 2, x[1]], bins), - ], + [binned_normal_pdf(a, b, bins), binned_normal_pdf(a + 1, b, bins)], + [binned_normal_pdf(a + 10, b, bins), binned_normal_pdf(a + 2, b, bins)], ] - for x in grid + for a, b in grid ] ) @@ -489,12 +483,12 @@ def test_MomentMorphInterpolator2D_mixed(bins): truth = np.array( [ [ - binned_normal_pdf(target, bins), - binned_normal_pdf([target[0] + 1, target[1]], bins), + binned_normal_pdf(*target, bins), + binned_normal_pdf(target[0] + 1, target[1], bins), ], [ - binned_normal_pdf([target[0] + 10, target[1]], bins), - binned_normal_pdf([target[0] + 2, target[1]], bins), + binned_normal_pdf(target[0] + 10, target[1], bins), + binned_normal_pdf(target[0] + 2, target[1], bins), ], ] ) @@ -522,7 +516,7 @@ def test_MomentMorphInterpolator1D_extended_grid(bins): grid = np.array([[20], [40], [60], [80]]) target = np.array([25]) - binned_pdf = np.array([binned_normal_pdf([x, 0], bins) for x in grid]) + binned_pdf = np.array([binned_normal_pdf(x, 0, bins) for x in grid]) interp = MomentMorphInterpolator( grid_points=grid, @@ -531,7 +525,7 @@ def test_MomentMorphInterpolator1D_extended_grid(bins): ) res = interp(target) - truth = binned_normal_pdf([target, 0], bins) + truth = binned_normal_pdf(target, 0, bins) assert np.isclose(np.sum(res), 1) assert np.all(np.isfinite(res)) @@ -545,7 +539,7 @@ def test_MomentMorphInterpolator2D_extended_grid(bins): grid = np.array([[20, 20], [40, 20], [30, 40], [50, 20], [45, 40]]) target = np.array([25, 25]) - binned_pdf = np.array([binned_normal_pdf(x, bins) for x in grid]) + binned_pdf = np.array([binned_normal_pdf(a, b, bins) for a, b in grid]) interp = MomentMorphInterpolator( grid_points=grid, @@ -554,7 +548,7 @@ def test_MomentMorphInterpolator2D_extended_grid(bins): ) res = interp(target) - truth = binned_normal_pdf(target, bins) + truth = binned_normal_pdf(*target, bins) assert np.isclose(np.sum(res), 1) assert np.all(np.isfinite(res)) @@ -567,17 +561,19 @@ def test_MomentMorphInterpolator2D_extended_grid_extradims(bins): from pyirf.interpolation import MomentMorphInterpolator grid = np.array([[20, 20], [40, 20], [30, 40], [50, 20], [45, 40]]) - target = np.array([25, 25]) + a = 25 + b = 25 + target = np.array([a, b]) binned_pdf = np.array( [ [ - [binned_normal_pdf(x, bins), binned_normal_pdf([x[0] + 1, x[1]], bins)], + [binned_normal_pdf(a, b, bins), binned_normal_pdf(a + 1, b, bins)], [ - binned_normal_pdf([x[0] + 10, x[1]], bins), - binned_normal_pdf([x[0] + 2, x[1]+5], bins), + binned_normal_pdf(a + 10, b, bins), + binned_normal_pdf(a + 2, b + 5, bins), ], ] - for x in grid + for a, b in grid ] ) @@ -590,12 +586,12 @@ def test_MomentMorphInterpolator2D_extended_grid_extradims(bins): truth = np.array( [ [ - binned_normal_pdf(target, bins), - binned_normal_pdf([target[0] + 1, target[1]], bins), + binned_normal_pdf(a, b, bins), + binned_normal_pdf(a + 1, b, bins), ], [ - binned_normal_pdf([target[0] + 10, target[1]], bins), - binned_normal_pdf([target[0] + 2, target[1]+5], bins), + binned_normal_pdf(a + 10, b, bins), + binned_normal_pdf(a + 2, b + 5], bins), ], ] ) diff --git a/pyirf/interpolation/tests/test_nearest_simplex_extrapolator.py b/pyirf/interpolation/tests/test_nearest_simplex_extrapolator.py index 389cfa2f7..ca4aa3ca4 100644 --- a/pyirf/interpolation/tests/test_nearest_simplex_extrapolator.py +++ b/pyirf/interpolation/tests/test_nearest_simplex_extrapolator.py @@ -11,13 +11,13 @@ def expected_std(a, b): return 1 + 0.05 * (a + b) -def binned_normal_pdf(mean_std_args, bins): - pdf = np.diff( - norm(loc=expected_mean(*mean_std_args), scale=expected_std(*mean_std_args)).cdf( - bins - ) +def binned_normal_pdf(a, b, bins): + dist = norm( + loc=expected_mean(a, b), + scale=expected_std(a, b), ) - return pdf / pdf.sum() + pdf = np.diff(dist.cdf(bins)) + return pdf / np.diff(bins) @pytest.fixture @@ -141,9 +141,9 @@ def test_MomentMorphNearestSimplexExtrapolator_1DGrid(bins): binned_pdf = np.array( [ [ - [binned_normal_pdf([x, 0], bins), binned_normal_pdf([x + 1, 0], bins)], + [binned_normal_pdf(x, 0, bins), binned_normal_pdf(x + 1, 0, bins)], [ - binned_normal_pdf([x + 10, 0], bins), + binned_normal_pdf(x + 10, 0, bins), np.zeros(len(bins) - 1), ], [np.zeros(len(bins) - 1), np.ones(len(bins) - 1) / (len(bins) - 1)], @@ -172,10 +172,10 @@ def test_MomentMorphNearestSimplexExtrapolator_1DGrid(bins): [ [ np.zeros(len(bins) - 1), - binned_normal_pdf([target1 + 1, 0], bins), + binned_normal_pdf(target1 + 1, 0, bins), ], [ - binned_normal_pdf([target1 + 10, 0], bins), + binned_normal_pdf(target1 + 10, 0, bins), np.zeros(len(bins) - 1), ], [np.zeros(len(bins) - 1), np.ones(len(bins) - 1) / (len(bins) - 1)], @@ -199,11 +199,11 @@ def test_MomentMorphNearestSimplexExtrapolator_1DGrid(bins): truth2 = np.array( [ [ - binned_normal_pdf([target2, 0], bins), - binned_normal_pdf([target2 + 1, 0], bins), + binned_normal_pdf(target2, 0, bins), + binned_normal_pdf(target2 + 1, 0, bins), ], [ - binned_normal_pdf([target2 + 10, 0], bins), + binned_normal_pdf(target2 + 10, 0, bins), np.zeros(len(bins) - 1), ], [np.zeros(len(bins) - 1), np.ones(len(bins) - 1) / (len(bins) - 1)], @@ -233,13 +233,10 @@ def test_MomentMorphNearestSimplexExtrapolator_2DGrid(bins): binned_pdf = np.array( [ [ - [binned_normal_pdf(x, bins), binned_normal_pdf([x[0] + 1, x[1]], bins)], - [ - binned_normal_pdf([x[0] + 10, x[1]], bins), - np.zeros(len(bins) - 1), - ], + [binned_normal_pdf(a, b, bins), binned_normal_pdf(a + 1, b, bins)], + [binned_normal_pdf(a + 10, b, bins), np.zeros(len(bins) - 1)], ] - for x in grid_points + for a, b in grid_points ] ) @@ -258,10 +255,10 @@ def test_MomentMorphNearestSimplexExtrapolator_2DGrid(bins): [ [ np.zeros(len(bins) - 1), - binned_normal_pdf([target1[0] + 1, target1[1]], bins), + binned_normal_pdf(target1[0] + 1, target1[1], bins), ], [ - binned_normal_pdf([target1[0] + 10, target1[1]], bins), + binned_normal_pdf(target1[0] + 10, target1[1], bins), np.zeros(len(bins) - 1), ], ] @@ -281,11 +278,11 @@ def test_MomentMorphNearestSimplexExtrapolator_2DGrid(bins): truth2 = np.array( [ [ - binned_normal_pdf(target2, bins), - binned_normal_pdf([target2[0] + 1, target2[1]], bins), + binned_normal_pdf(*target2, bins), + binned_normal_pdf(target2[0] + 1, target2[1], bins), ], [ - binned_normal_pdf([target2[0] + 10, target2[1]], bins), + binned_normal_pdf(target2[0] + 10, target2[1], bins), np.zeros(len(bins) - 1), ], ] @@ -332,7 +329,7 @@ def test_MomentMorphNearestSimplexExtrapolator_below_zero_warning(bins): binned_pdf = np.array( [ - [binned_normal_pdf([x, 0], bins), binned_normal_pdf([x + 10, 0], bins)] + [binned_normal_pdf(x, 0, bins), binned_normal_pdf(x + 10, 0, bins)] for x in grid_points ] ) From 6d51dc61123078fdac54a13792e09b9bc5ce8034 Mon Sep 17 00:00:00 2001 From: Maximilian Linhoff Date: Mon, 21 Aug 2023 16:40:16 +0200 Subject: [PATCH 05/14] Start working on fixing interpolation... --- pyirf/interpolation/component_estimators.py | 11 +++- pyirf/interpolation/quantile_interpolator.py | 24 +++---- .../tests/test_moment_morph_interpolator.py | 62 +++++++++---------- .../test_nearest_simplex_extrapolator.py | 2 +- .../tests/test_quantile_interpolator.py | 24 ++++--- 5 files changed, 69 insertions(+), 54 deletions(-) diff --git a/pyirf/interpolation/component_estimators.py b/pyirf/interpolation/component_estimators.py index e2caf7ef3..93b90fd5f 100644 --- a/pyirf/interpolation/component_estimators.py +++ b/pyirf/interpolation/component_estimators.py @@ -7,7 +7,7 @@ from scipy.spatial import Delaunay from .base_extrapolators import DiscretePDFExtrapolator, ParametrizedExtrapolator -from .base_interpolators import DiscretePDFInterpolator, ParametrizedInterpolator +from .base_interpolators import DiscretePDFInterpolator, PDFNormalization, ParametrizedInterpolator from .griddata_interpolator import GridDataInterpolator from .quantile_interpolator import QuantileInterpolator @@ -685,6 +685,15 @@ def __init__( self.omegas = np.diff(cone_solid_angle(source_offset_bins)) psf_normed = psf * self.omegas + if interpolator_kwargs is None: + interpolator_kwargs = {} + + if extrapolator_kwargs is None: + extrapolator_kwargs = {} + + interpolator_kwargs.setdefault("normalization", PDFNormalization.CONE_SOLID_ANGLE) + extrapolator_kwargs.setdefault("normalization", PDFNormalization.CONE_SOLID_ANGLE) + super().__init__( grid_points=grid_points, bin_edges=source_offset_bins, diff --git a/pyirf/interpolation/quantile_interpolator.py b/pyirf/interpolation/quantile_interpolator.py index 8170107ba..844028cb9 100644 --- a/pyirf/interpolation/quantile_interpolator.py +++ b/pyirf/interpolation/quantile_interpolator.py @@ -1,6 +1,5 @@ import numpy as np from scipy.interpolate import griddata, interp1d -import enum from ..utils import cone_solid_angle from .base_interpolators import DiscretePDFInterpolator, PDFNormalization @@ -8,18 +7,21 @@ __all__ = ["QuantileInterpolator"] +def get_bin_width(bin_edges, normalization): + if normalization is PDFNormalization.AREA: + return np.diff(bin_edges) + + if normalization is PDFNormalization.CONE_SOLID_ANGLE: + return np.diff(cone_solid_angle(bin_edges)) + + raise ValueError(f"Invalid PDF normalization: {normalization}") + + def cdf_values(binned_pdf, bin_edges, normalization): """ compute cdf values and assure they are normed to unity """ - - if normalization is PDFNormalization.AREA: - bin_widths = np.diff(bin_edges) - elif normalization is PDFNormalization.CONE_SOLID_ANGLE: - bin_widths = np.diff(cone_solid_angle(bin_edges)) - else: - raise ValueError(f"Invalid PDF normalization: {normalization}") - + bin_widths = get_bin_width(bin_edges, normalization) cdfs = np.cumsum(binned_pdf * bin_widths, axis=-1) # assure the last cdf value is 1, ignore errors for empty pdfs as they are reset to 0 by nan_to_num @@ -105,7 +107,7 @@ def pdf_from_ppf(bin_edges, interp_ppfs, quantiles): Edges of the bins in which the final pdf should be binned interp_ppfs: numpy.ndarray, shape=(1,...,L) - Corresponding ppf-values for all self.quantiles at the target_point, + Corresponding ppf-values for all quantiles at the target_point, not to be confused with QunatileInterpolators self.ppfs, the ppfs computed from the input distributions. @@ -132,7 +134,7 @@ def interpolate_ppf(xy): pdf = xy[len(xy) // 2 :] interpolate = interp1d(ppf, pdf, bounds_error=False, fill_value=(0, 0)) result = np.nan_to_num(interpolate(bin_edges[:-1])) - return np.diff(bin_edges) * result + return result # Interpolate pdf samples and evaluate at bin edges, weight with the bin_width to estimate # correct bin height via the midpoint rule formulation of the trapezoidal rule diff --git a/pyirf/interpolation/tests/test_moment_morph_interpolator.py b/pyirf/interpolation/tests/test_moment_morph_interpolator.py index 715597aab..f496dfc31 100644 --- a/pyirf/interpolation/tests/test_moment_morph_interpolator.py +++ b/pyirf/interpolation/tests/test_moment_morph_interpolator.py @@ -97,8 +97,8 @@ def test_estimate_mean_std(bins): mean, std = _estimate_mean_std(bins, binned_pdf) # Assert estimation and truth within one bin - assert np.allclose(mean, true_mean, atol=np.diff(bins)[0] / 2) - assert np.allclose(std, true_std, atol=np.diff(bins)[0] / 2) + np.testing.assert_allclose(mean, true_mean, atol=np.diff(bins)[0] / 2) + np.testing.assert_allclose(std, true_std, atol=np.diff(bins)[0] / 2) def test_lookup(): @@ -126,7 +126,7 @@ def test_lookup(): ] ) - assert np.allclose(_lookup(bins, binned_pdf, x), truth) + np.testing.assert_allclose(_lookup(bins, binned_pdf, x), truth) def test_linesegment_1D_interpolation_coefficients(): @@ -137,14 +137,14 @@ def test_linesegment_1D_interpolation_coefficients(): grid_points = np.array([[20], [40]]) target_point = np.array([[20]]) - assert np.allclose( + np.testing.assert_allclose( linesegment_1D_interpolation_coefficients(grid_points, target_point), np.array([1, 0]), ) target_point = np.array([[40]]) - assert np.allclose( + np.testing.assert_allclose( linesegment_1D_interpolation_coefficients(grid_points, target_point), np.array([0, 1]), ) @@ -160,7 +160,7 @@ def test_linesegment_1D_interpolation_coefficients(): res = linesegment_1D_interpolation_coefficients(grid_points, target_point) - assert np.allclose(res, cs) + np.testing.assert_allclose(res, cs) assert np.sum(res) == 1 assert np.all(np.logical_and(res < 1, res > 0)) @@ -173,28 +173,28 @@ def test_barycentric_2D_interpolation_coefficients(): grid_points = np.array([[20, 20], [60, 20], [40, 60]]) target_point = np.array([20, 20]) - assert np.allclose( + np.testing.assert_allclose( barycentric_2D_interpolation_coefficients(grid_points, target_point), np.array([1, 0, 0]), ) target_point = np.array([60, 20]) - assert np.allclose( + np.testing.assert_allclose( barycentric_2D_interpolation_coefficients(grid_points, target_point), np.array([0, 1, 0]), ) target_point = np.array([40, 60]) - assert np.allclose( + np.testing.assert_allclose( barycentric_2D_interpolation_coefficients(grid_points, target_point), np.array([0, 0, 1]), ) target_point = np.array([40, 20]) - assert np.allclose( + np.testing.assert_allclose( barycentric_2D_interpolation_coefficients(grid_points, target_point), np.array([0.5, 0.5, 0]), ) @@ -203,7 +203,7 @@ def test_barycentric_2D_interpolation_coefficients(): target_point = np.array([40, 100 / 3]) res = barycentric_2D_interpolation_coefficients(grid_points, target_point) - assert np.allclose(res, np.array([1, 1, 1]) / 3) + np.testing.assert_allclose(res, np.array([1, 1, 1]) / 3) def test_moment_morph_estimation1D(bins, simple_1D_data): @@ -221,7 +221,7 @@ def test_moment_morph_estimation1D(bins, simple_1D_data): assert np.all(np.isfinite(res)) assert res.shape == (1, len(bins) - 1) # Assert truth and result matching within +- 0.1%, atol dominates comparison - assert np.allclose(res.squeeze(), truth, atol=1e-3, rtol=1e-5) + np.testing.assert_allclose(res.squeeze(), truth, atol=1e-3, rtol=1e-5) def test_moment_morph_estimation2D(bins, simple_2D_data): @@ -239,7 +239,7 @@ def test_moment_morph_estimation2D(bins, simple_2D_data): assert np.all(np.isfinite(res)) assert res.shape == (1, len(bins) - 1) # Assert truth and result matching within +- 0.1%, atol dominates comparison - assert np.allclose(res.squeeze(), truth, atol=1e-3, rtol=1e-5) + np.testing.assert_allclose(res.squeeze(), truth, atol=1e-3, rtol=1e-5) def test_MomentMorphInterpolator1D(bins, simple_1D_data): @@ -257,7 +257,7 @@ def test_MomentMorphInterpolator1D(bins, simple_1D_data): assert np.all(np.isfinite(res)) assert res.shape == (1, len(bins) - 1) # Assert truth and result matching within +- 0.1%, atol dominates comparison - assert np.allclose(res.squeeze(), truth, atol=1e-3, rtol=1e-5) + np.testing.assert_allclose(res.squeeze(), truth, atol=1e-3, rtol=1e-5) def test_MomentMorphInterpolator1D_dirac_delta_input(): @@ -276,7 +276,7 @@ def test_MomentMorphInterpolator1D_dirac_delta_input(): interp = MomentMorphInterpolator(grid, bin_edges, binned_pdf) res = interp(target) - assert np.allclose(res, np.array([[0, 0, 1, 0], [0.25, 0.25, 0.25, 0.25]])) + np.testing.assert_allclose(res, np.array([[0, 0, 1, 0], [0.25, 0.25, 0.25, 0.25]])) def test_MomentMorphInterpolator1D_all_empty(bins, simple_1D_data): @@ -291,7 +291,7 @@ def test_MomentMorphInterpolator1D_all_empty(bins, simple_1D_data): res = interp(target) - assert np.allclose(res, 0) + np.testing.assert_allclose(res, 0) def test_MomentMorphInterpolator1D_partially_empty(bins, simple_1D_data): @@ -307,7 +307,7 @@ def test_MomentMorphInterpolator1D_partially_empty(bins, simple_1D_data): res = interp(target) - assert np.allclose(res, 0) + np.testing.assert_allclose(res, 0) def test_MomentMorphInterpolator1D_mixed_data(bins): @@ -357,11 +357,11 @@ def test_MomentMorphInterpolator1D_mixed_data(bins): res = interp(target) expected_norms = np.array([[0, 1], [1, 0]]) - assert np.allclose(np.sum(res, axis=-1), expected_norms) + np.testing.assert_allclose(np.sum(res, axis=-1), expected_norms) assert np.all(np.isfinite(res)) assert res.shape == (1, *binned_pdf.shape[1:]) # Assert truth and result matching within +- 0.1%, atol dominates comparison - assert np.allclose(res.squeeze(), truth, atol=1e-3, rtol=1e-5) + np.testing.assert_allclose(res.squeeze(), truth, atol=1e-3, rtol=1e-5) def test_MomentMorphInterpolator1D_extended_grid_extradims(bins): @@ -401,11 +401,11 @@ def test_MomentMorphInterpolator1D_extended_grid_extradims(bins): res = interp(target) - assert np.allclose(np.sum(res, axis=-1), 1) + np.testing.assert_allclose(np.sum(res, axis=-1), 1) assert np.all(np.isfinite(res)) assert res.shape == (1, *binned_pdf.shape[1:]) # Assert truth and result matching within +- 0.1%, atol dominates comparison - assert np.allclose(res.squeeze(), truth, atol=1e-4, rtol=1e-5) + np.testing.assert_allclose(res.squeeze(), truth, atol=1e-4, rtol=1e-5) def test_MomentMorphInterpolator2D(bins, simple_2D_data): @@ -423,7 +423,7 @@ def test_MomentMorphInterpolator2D(bins, simple_2D_data): assert np.all(np.isfinite(res)) assert res.shape == (1, len(bins) - 1) # Assert truth and result matching within +- 0.1%, atol dominates comparison - assert np.allclose(res.squeeze(), truth, atol=1e-3, rtol=1e-5) + np.testing.assert_allclose(res.squeeze(), truth, atol=1e-3, rtol=1e-5) def test_MomentMorphInterpolator2D_partially_empty(bins, simple_2D_data): @@ -439,7 +439,7 @@ def test_MomentMorphInterpolator2D_partially_empty(bins, simple_2D_data): res = interp(target) - assert np.allclose(res, 0) + np.testing.assert_allclose(res, 0) def test_MomentMorphInterpolator2D_all_empty(bins, simple_2D_data): @@ -454,7 +454,7 @@ def test_MomentMorphInterpolator2D_all_empty(bins, simple_2D_data): res = interp(target) - assert np.allclose(res, 0) + np.testing.assert_allclose(res, 0) def test_MomentMorphInterpolator2D_mixed(bins): @@ -504,11 +504,11 @@ def test_MomentMorphInterpolator2D_mixed(bins): res = interp(target) expected_norms = np.array([[0, 1], [1, 0]]) - assert np.allclose(np.sum(res, axis=-1), expected_norms) + np.testing.assert_allclose(np.sum(res, axis=-1), expected_norms) assert np.all(np.isfinite(res)) assert res.shape == (1, *binned_pdf.shape[1:]) # Assert truth and result matching within +- 0.1%, atol dominates comparison - assert np.allclose(res.squeeze(), truth, atol=1e-3, rtol=1e-5) + np.testing.assert_allclose(res.squeeze(), truth, atol=1e-3, rtol=1e-5) def test_MomentMorphInterpolator1D_extended_grid(bins): @@ -531,7 +531,7 @@ def test_MomentMorphInterpolator1D_extended_grid(bins): assert np.all(np.isfinite(res)) assert res.shape == (1, len(bins) - 1) # Assert truth and result matching within +- 0.1%, atol dominates comparison - assert np.allclose(res.squeeze(), truth, atol=1e-3, rtol=1e-5) + np.testing.assert_allclose(res.squeeze(), truth, atol=1e-3, rtol=1e-5) def test_MomentMorphInterpolator2D_extended_grid(bins): @@ -554,7 +554,7 @@ def test_MomentMorphInterpolator2D_extended_grid(bins): assert np.all(np.isfinite(res)) assert res.shape == (1, len(bins) - 1) # Assert truth and result matching within +- 0.1%, atol dominates comparison - assert np.allclose(res.squeeze(), truth, atol=1e-3, rtol=1e-5) + np.testing.assert_allclose(res.squeeze(), truth, atol=1e-3, rtol=1e-5) def test_MomentMorphInterpolator2D_extended_grid_extradims(bins): @@ -591,18 +591,18 @@ def test_MomentMorphInterpolator2D_extended_grid_extradims(bins): ], [ binned_normal_pdf(a + 10, b, bins), - binned_normal_pdf(a + 2, b + 5], bins), + binned_normal_pdf(a + 2, b + 5, bins), ], ] ) res = interp(target) - assert np.allclose(np.sum(res, axis=-1), 1) + np.testing.assert_allclose(np.sum(res, axis=-1), 1) assert np.all(np.isfinite(res)) assert res.shape == (1, *binned_pdf.shape[1:]) # Assert truth and result matching within +- 0.1%, atol dominates comparison - assert np.allclose(res.squeeze(), truth, atol=1e-3, rtol=1e-5) + np.testing.assert_allclose(res.squeeze(), truth, atol=1e-3, rtol=1e-5) def test_MomentMorphInterpolator3D(): diff --git a/pyirf/interpolation/tests/test_nearest_simplex_extrapolator.py b/pyirf/interpolation/tests/test_nearest_simplex_extrapolator.py index ca4aa3ca4..384c32570 100644 --- a/pyirf/interpolation/tests/test_nearest_simplex_extrapolator.py +++ b/pyirf/interpolation/tests/test_nearest_simplex_extrapolator.py @@ -47,7 +47,7 @@ def test_ParametrizedNearestSimplexExtrapolator_1DGrid(): dummy_data_target1 = 3 * slope + 1 - assert np.allclose(interpolant1, dummy_data_target1) + np.testing.assert_allclose(interpolant1, dummy_data_target1) assert interpolant1.shape == (1, *dummy_data.shape[1:]) target_point2 = np.array([[-2.5]]) diff --git a/pyirf/interpolation/tests/test_quantile_interpolator.py b/pyirf/interpolation/tests/test_quantile_interpolator.py index a5a9cc3fd..28a129bd5 100644 --- a/pyirf/interpolation/tests/test_quantile_interpolator.py +++ b/pyirf/interpolation/tests/test_quantile_interpolator.py @@ -3,6 +3,8 @@ from pyirf.binning import bin_center from scipy.stats import norm +from pyirf.interpolation.base_interpolators import PDFNormalization + @pytest.fixture def data(): @@ -12,6 +14,8 @@ def data(): # create binned pdfs by interpolation of bin content binned_pdfs = np.array([np.diff(dist.cdf(bin_edges)) for dist in distributions]) + bin_width = np.diff(bin_edges) + binned_pdfs /= bin_width[np.newaxis, :] dataset = { "bin_edges": bin_edges, @@ -28,11 +32,10 @@ def data(): def test_cdf_values(data): from pyirf.interpolation.quantile_interpolator import cdf_values - cdf_est = cdf_values(data["binned_pdfs"][0]) - # Assert empty histograms result in cdf containing zeros - assert np.all(cdf_values(np.zeros(shape=5)) == 0) + np.testing.assert_array_equal(cdf_values(np.zeros_like(data["binned_pdfs"])[0], data["bin_edges"], PDFNormalization.AREA), 0) + cdf_est = cdf_values(data["binned_pdfs"][0], data["bin_edges"], PDFNormalization.AREA) # Assert cdf is increasing or constant for actual pdfs assert np.all(np.diff(cdf_est) >= 0) @@ -40,7 +43,8 @@ def test_cdf_values(data): assert np.max(cdf_est) == 1 # Assert estimated and true cdf are matching - assert np.allclose(cdf_est, data["distributions"][0].cdf(data["bin_edges"][1:])) + true_cdf = data["distributions"][0].cdf(data["bin_edges"][1:]) + np.testing.assert_allclose(cdf_est, true_cdf, atol=1e-12) def test_ppf_values(data): @@ -54,11 +58,11 @@ def test_ppf_values(data): bin_mids = bin_center(data["bin_edges"]) # Estimated ppf-values - cdf_est = cdf_values(data["binned_pdfs"][0]) + cdf_est = cdf_values(data["binned_pdfs"][0], data["bin_edges"], PDFNormalization.AREA) ppf_est = ppf_values(bin_mids, cdf_est, quantiles) # Assert truth and estimation match allowing for +- bin_width deviation - assert np.allclose(ppf_true, ppf_est, atol=np.diff(data["bin_edges"])[0]) + np.testing.assert_allclose(ppf_true, ppf_est, atol=np.diff(data["bin_edges"])[0]) def test_pdf_from_ppf(data): @@ -73,21 +77,21 @@ def test_pdf_from_ppf(data): bin_mids = bin_center(data["bin_edges"]) # Estimate ppf-values - cdf_est = cdf_values(data["binned_pdfs"][0]) + cdf_est = cdf_values(data["binned_pdfs"][0], data["bin_edges"], PDFNormalization.AREA) ppf_est = ppf_values(bin_mids, cdf_est, quantiles) # Compute pdf_values pdf_est = pdf_from_ppf(data["bin_edges"], ppf_est, quantiles) # Assert pdf-values matching true pdf within +-1% - assert np.allclose(pdf_est, data["binned_pdfs"][0], atol=1e-2) + np.testing.assert_allclose(pdf_est, data["binned_pdfs"][0], atol=1e-2) def test_norm_pdf(data): from pyirf.interpolation.quantile_interpolator import norm_pdf - assert np.allclose(norm_pdf(2 * data["binned_pdfs"][0]), data["binned_pdfs"][0]) - assert np.allclose(norm_pdf(np.zeros(5)), 0) + np.testing.assert_allclose(norm_pdf(2 * data["binned_pdfs"][0]), data["binned_pdfs"][0]) + np.testing.assert_allclose(norm_pdf(np.zeros(5)), 0) def test_interpolate_binned_pdf(data): From 8fe8b875ce14f528bf7a65fdc8e2238440dbdadd Mon Sep 17 00:00:00 2001 From: Maximilian Linhoff Date: Mon, 21 Aug 2023 16:43:05 +0200 Subject: [PATCH 06/14] Fix mean / std estimation for psf --- pyirf/interpolation/base_interpolators.py | 16 ++++++++++++++-- pyirf/interpolation/moment_morph_interpolator.py | 9 +++++---- pyirf/interpolation/quantile_interpolator.py | 13 +------------ 3 files changed, 20 insertions(+), 18 deletions(-) diff --git a/pyirf/interpolation/base_interpolators.py b/pyirf/interpolation/base_interpolators.py index 9b4227f38..14467e593 100644 --- a/pyirf/interpolation/base_interpolators.py +++ b/pyirf/interpolation/base_interpolators.py @@ -3,9 +3,11 @@ import enum import numpy as np -from pyirf.binning import bin_center -__all__ = ["BaseInterpolator", "ParametrizedInterpolator", "DiscretePDFInterpolator"] +from ..binning import bin_center +from ..utils import cone_solid_angle + +__all__ = ["BaseInterpolator", "ParametrizedInterpolator", "DiscretePDFInterpolator", "PDFNormalization", "get_bin_width"] class PDFNormalization(enum.Enum): @@ -17,6 +19,16 @@ class PDFNormalization(enum.Enum): CONE_SOLID_ANGLE = enum.auto() +def get_bin_width(bin_edges, normalization): + if normalization is PDFNormalization.AREA: + return np.diff(bin_edges) + + if normalization is PDFNormalization.CONE_SOLID_ANGLE: + return np.diff(cone_solid_angle(bin_edges)) + + raise ValueError(f"Invalid PDF normalization: {normalization}") + + class BaseInterpolator(metaclass=ABCMeta): """ Base class for all interpolators, only knowing grid-points, diff --git a/pyirf/interpolation/moment_morph_interpolator.py b/pyirf/interpolation/moment_morph_interpolator.py index fafcfccb4..3480c37a8 100644 --- a/pyirf/interpolation/moment_morph_interpolator.py +++ b/pyirf/interpolation/moment_morph_interpolator.py @@ -2,14 +2,14 @@ from pyirf.binning import bin_center, calculate_bin_indices from scipy.spatial import Delaunay -from .base_interpolators import DiscretePDFInterpolator, PDFNormalization +from .base_interpolators import DiscretePDFInterpolator, PDFNormalization, get_bin_width __all__ = [ "MomentMorphInterpolator", ] -def _estimate_mean_std(bin_edges, binned_pdf): +def _estimate_mean_std(bin_edges, binned_pdf, normalization): """ Function to roughly estimate mean and standard deviation from a histogram. @@ -30,10 +30,11 @@ def _estimate_mean_std(bin_edges, binned_pdf): """ # Create an 2darray where the 1darray mids is repeated n_template times mids = np.broadcast_to(bin_center(bin_edges), binned_pdf.shape) - width = np.diff(bin_edges) + + width = get_bin_width(bin_edges, normalization) # integrate pdf to get probability in each bin - probability = binned_pdf * width + probability = binned_pdf * width # Weighted averages to compute mean and std mean = np.average(mids, weights=probability, axis=-1) std = np.sqrt( diff --git a/pyirf/interpolation/quantile_interpolator.py b/pyirf/interpolation/quantile_interpolator.py index 844028cb9..e36092ca7 100644 --- a/pyirf/interpolation/quantile_interpolator.py +++ b/pyirf/interpolation/quantile_interpolator.py @@ -1,22 +1,11 @@ import numpy as np from scipy.interpolate import griddata, interp1d -from ..utils import cone_solid_angle -from .base_interpolators import DiscretePDFInterpolator, PDFNormalization +from .base_interpolators import DiscretePDFInterpolator, PDFNormalization, get_bin_width __all__ = ["QuantileInterpolator"] -def get_bin_width(bin_edges, normalization): - if normalization is PDFNormalization.AREA: - return np.diff(bin_edges) - - if normalization is PDFNormalization.CONE_SOLID_ANGLE: - return np.diff(cone_solid_angle(bin_edges)) - - raise ValueError(f"Invalid PDF normalization: {normalization}") - - def cdf_values(binned_pdf, bin_edges, normalization): """ compute cdf values and assure they are normed to unity From 0705afa760e24e2d61d3be4d35b4121b8ac5e6eb Mon Sep 17 00:00:00 2001 From: Maximilian Linhoff Date: Mon, 21 Aug 2023 16:43:58 +0200 Subject: [PATCH 07/14] Run black --- pyirf/interpolation/base_interpolators.py | 15 ++++++++--- pyirf/interpolation/component_estimators.py | 22 +++++++++++----- .../moment_morph_interpolator.py | 6 +++-- .../nearest_neighbor_searcher.py | 4 +-- pyirf/interpolation/quantile_interpolator.py | 11 ++++++-- .../tests/test_quantile_interpolator.py | 25 +++++++++++++++---- 6 files changed, 62 insertions(+), 21 deletions(-) diff --git a/pyirf/interpolation/base_interpolators.py b/pyirf/interpolation/base_interpolators.py index 14467e593..cf4cf2ffa 100644 --- a/pyirf/interpolation/base_interpolators.py +++ b/pyirf/interpolation/base_interpolators.py @@ -7,11 +7,18 @@ from ..binning import bin_center from ..utils import cone_solid_angle -__all__ = ["BaseInterpolator", "ParametrizedInterpolator", "DiscretePDFInterpolator", "PDFNormalization", "get_bin_width"] +__all__ = [ + "BaseInterpolator", + "ParametrizedInterpolator", + "DiscretePDFInterpolator", + "PDFNormalization", + "get_bin_width", +] class PDFNormalization(enum.Enum): """How a discrete PDF is normalized""" + #: PDF is normalized to a "normal" area integral of 1 AREA = enum.auto() #: PDF is normalized to 1 over the solid angle integral where the bin @@ -25,7 +32,7 @@ def get_bin_width(bin_edges, normalization): if normalization is PDFNormalization.CONE_SOLID_ANGLE: return np.diff(cone_solid_angle(bin_edges)) - + raise ValueError(f"Invalid PDF normalization: {normalization}") @@ -106,7 +113,9 @@ class DiscretePDFInterpolator(BaseInterpolator): Derived from pyirf.interpolation.BaseInterpolator """ - def __init__(self, grid_points, bin_edges, binned_pdf, normalization=PDFNormalization.AREA): + def __init__( + self, grid_points, bin_edges, binned_pdf, normalization=PDFNormalization.AREA + ): """DiscretePDFInterpolator Parameters diff --git a/pyirf/interpolation/component_estimators.py b/pyirf/interpolation/component_estimators.py index 93b90fd5f..6ed0edc5d 100644 --- a/pyirf/interpolation/component_estimators.py +++ b/pyirf/interpolation/component_estimators.py @@ -7,7 +7,11 @@ from scipy.spatial import Delaunay from .base_extrapolators import DiscretePDFExtrapolator, ParametrizedExtrapolator -from .base_interpolators import DiscretePDFInterpolator, PDFNormalization, ParametrizedInterpolator +from .base_interpolators import ( + DiscretePDFInterpolator, + PDFNormalization, + ParametrizedInterpolator, +) from .griddata_interpolator import GridDataInterpolator from .quantile_interpolator import QuantileInterpolator @@ -231,7 +235,7 @@ def __init__( # Make sure that 1D input is sorted in increasing order if self.grid_dim == 1: sorting_inds = np.argsort(self.grid_points.squeeze()) - + self.grid_points = self.grid_points[sorting_inds] binned_pdf = binned_pdf[sorting_inds] @@ -334,7 +338,7 @@ def __init__( # Make sure that 1D input is sorted in increasing order if self.grid_dim == 1: sorting_inds = np.argsort(self.grid_points.squeeze()) - + self.grid_points = self.grid_points[sorting_inds] params = params[sorting_inds] @@ -349,7 +353,9 @@ def __init__( f"interpolator_cls must be a ParametrizedInterpolator subclass, got {interpolator_cls}" ) - self.interpolator = interpolator_cls(self.grid_points, params, **interpolator_kwargs) + self.interpolator = interpolator_cls( + self.grid_points, params, **interpolator_kwargs + ) if extrapolator_cls is None: self.extrapolator = None @@ -691,8 +697,12 @@ def __init__( if extrapolator_kwargs is None: extrapolator_kwargs = {} - interpolator_kwargs.setdefault("normalization", PDFNormalization.CONE_SOLID_ANGLE) - extrapolator_kwargs.setdefault("normalization", PDFNormalization.CONE_SOLID_ANGLE) + interpolator_kwargs.setdefault( + "normalization", PDFNormalization.CONE_SOLID_ANGLE + ) + extrapolator_kwargs.setdefault( + "normalization", PDFNormalization.CONE_SOLID_ANGLE + ) super().__init__( grid_points=grid_points, diff --git a/pyirf/interpolation/moment_morph_interpolator.py b/pyirf/interpolation/moment_morph_interpolator.py index 3480c37a8..c4c40827a 100644 --- a/pyirf/interpolation/moment_morph_interpolator.py +++ b/pyirf/interpolation/moment_morph_interpolator.py @@ -38,7 +38,7 @@ def _estimate_mean_std(bin_edges, binned_pdf, normalization): # Weighted averages to compute mean and std mean = np.average(mids, weights=probability, axis=-1) std = np.sqrt( - np.average((mids - mean[..., np.newaxis])**2, weights=probability, axis=-1) + np.average((mids - mean[..., np.newaxis]) ** 2, weights=probability, axis=-1) ) # Set std to 0.5*width for all those templates that have only one bin =/= 0. In those @@ -253,7 +253,9 @@ def moment_morph_estimation(bin_edges, binned_pdf, coefficients): class MomentMorphInterpolator(DiscretePDFInterpolator): - def __init__(self, grid_points, bin_edges, binned_pdf, normalization=PDFNormalization.AREA): + def __init__( + self, grid_points, bin_edges, binned_pdf, normalization=PDFNormalization.AREA + ): """ Interpolator class using moment morphing. diff --git a/pyirf/interpolation/nearest_neighbor_searcher.py b/pyirf/interpolation/nearest_neighbor_searcher.py index 4db01d14f..6ddcf5df7 100644 --- a/pyirf/interpolation/nearest_neighbor_searcher.py +++ b/pyirf/interpolation/nearest_neighbor_searcher.py @@ -129,9 +129,7 @@ def __init__(self, grid_points, bin_edges, binned_pdf, norm_ord=2): Also calls pyirf.interpolation.BaseNearestNeighborSearcher.__init__ """ - super().__init__( - grid_points=grid_points, values=binned_pdf, norm_ord=norm_ord - ) + super().__init__(grid_points=grid_points, values=binned_pdf, norm_ord=norm_ord) DiscretePDFInterpolator.register(DiscretePDFNearestNeighborSearcher) diff --git a/pyirf/interpolation/quantile_interpolator.py b/pyirf/interpolation/quantile_interpolator.py index e36092ca7..3bc9583fc 100644 --- a/pyirf/interpolation/quantile_interpolator.py +++ b/pyirf/interpolation/quantile_interpolator.py @@ -150,7 +150,14 @@ def norm_pdf(pdf_values): class QuantileInterpolator(DiscretePDFInterpolator): - def __init__(self, grid_points, bin_edges, binned_pdf, quantile_resolution=1e-3, normalization=PDFNormalization.AREA): + def __init__( + self, + grid_points, + bin_edges, + binned_pdf, + quantile_resolution=1e-3, + normalization=PDFNormalization.AREA, + ): """BinnedInterpolator constructor Parameters @@ -168,7 +175,7 @@ def __init__(self, grid_points, bin_edges, binned_pdf, quantile_resolution=1e-3, quantile_resolution : float Spacing between quantiles normalization : PDFNormalization - How the discrete PDF is normalized + How the discrete PDF is normalized Raises ------ diff --git a/pyirf/interpolation/tests/test_quantile_interpolator.py b/pyirf/interpolation/tests/test_quantile_interpolator.py index 28a129bd5..c7a19163c 100644 --- a/pyirf/interpolation/tests/test_quantile_interpolator.py +++ b/pyirf/interpolation/tests/test_quantile_interpolator.py @@ -33,9 +33,18 @@ def test_cdf_values(data): from pyirf.interpolation.quantile_interpolator import cdf_values # Assert empty histograms result in cdf containing zeros - np.testing.assert_array_equal(cdf_values(np.zeros_like(data["binned_pdfs"])[0], data["bin_edges"], PDFNormalization.AREA), 0) + np.testing.assert_array_equal( + cdf_values( + np.zeros_like(data["binned_pdfs"])[0], + data["bin_edges"], + PDFNormalization.AREA, + ), + 0, + ) - cdf_est = cdf_values(data["binned_pdfs"][0], data["bin_edges"], PDFNormalization.AREA) + cdf_est = cdf_values( + data["binned_pdfs"][0], data["bin_edges"], PDFNormalization.AREA + ) # Assert cdf is increasing or constant for actual pdfs assert np.all(np.diff(cdf_est) >= 0) @@ -58,7 +67,9 @@ def test_ppf_values(data): bin_mids = bin_center(data["bin_edges"]) # Estimated ppf-values - cdf_est = cdf_values(data["binned_pdfs"][0], data["bin_edges"], PDFNormalization.AREA) + cdf_est = cdf_values( + data["binned_pdfs"][0], data["bin_edges"], PDFNormalization.AREA + ) ppf_est = ppf_values(bin_mids, cdf_est, quantiles) # Assert truth and estimation match allowing for +- bin_width deviation @@ -77,7 +88,9 @@ def test_pdf_from_ppf(data): bin_mids = bin_center(data["bin_edges"]) # Estimate ppf-values - cdf_est = cdf_values(data["binned_pdfs"][0], data["bin_edges"], PDFNormalization.AREA) + cdf_est = cdf_values( + data["binned_pdfs"][0], data["bin_edges"], PDFNormalization.AREA + ) ppf_est = ppf_values(bin_mids, cdf_est, quantiles) # Compute pdf_values @@ -90,7 +103,9 @@ def test_pdf_from_ppf(data): def test_norm_pdf(data): from pyirf.interpolation.quantile_interpolator import norm_pdf - np.testing.assert_allclose(norm_pdf(2 * data["binned_pdfs"][0]), data["binned_pdfs"][0]) + np.testing.assert_allclose( + norm_pdf(2 * data["binned_pdfs"][0]), data["binned_pdfs"][0] + ) np.testing.assert_allclose(norm_pdf(np.zeros(5)), 0) From 3c81ad9df6a312a4605efa0dbde748ad6974786b Mon Sep 17 00:00:00 2001 From: Maximilian Linhoff Date: Mon, 21 Aug 2023 16:55:09 +0200 Subject: [PATCH 08/14] Fix some errors in tests --- pyirf/interpolation/base_extrapolators.py | 4 +++- .../moment_morph_interpolator.py | 6 +++-- .../nearest_simplex_extrapolator.py | 8 ++++--- .../tests/test_moment_morph_interpolator.py | 14 ++++++----- .../tests/test_nearest_neighbor_searcher.py | 2 +- .../test_nearest_simplex_extrapolator.py | 24 +++++++++---------- 6 files changed, 33 insertions(+), 25 deletions(-) diff --git a/pyirf/interpolation/base_extrapolators.py b/pyirf/interpolation/base_extrapolators.py index f63800166..d848009c4 100644 --- a/pyirf/interpolation/base_extrapolators.py +++ b/pyirf/interpolation/base_extrapolators.py @@ -3,6 +3,7 @@ import numpy as np from pyirf.binning import bin_center +from pyirf.interpolation.base_interpolators import PDFNormalization __all__ = ["BaseExtrapolator", "ParametrizedExtrapolator", "DiscretePDFExtrapolator"] @@ -83,7 +84,7 @@ class DiscretePDFExtrapolator(BaseExtrapolator): Derived from pyirf.interpolation.BaseExtrapolator """ - def __init__(self, grid_points, bin_edges, binned_pdf): + def __init__(self, grid_points, bin_edges, binned_pdf, normalization=PDFNormalization.AREA): """DiscretePDFExtrapolator Parameters @@ -105,6 +106,7 @@ def __init__(self, grid_points, bin_edges, binned_pdf): """ super().__init__(grid_points) + self.normalization = normalization self.bin_edges = bin_edges self.bin_mids = bin_center(self.bin_edges) self.binned_pdf = binned_pdf diff --git a/pyirf/interpolation/moment_morph_interpolator.py b/pyirf/interpolation/moment_morph_interpolator.py index c4c40827a..6c102d66b 100644 --- a/pyirf/interpolation/moment_morph_interpolator.py +++ b/pyirf/interpolation/moment_morph_interpolator.py @@ -177,7 +177,7 @@ def barycentric_2D_interpolation_coefficients(grid_points, target_point): return coefficients -def moment_morph_estimation(bin_edges, binned_pdf, coefficients): +def moment_morph_estimation(bin_edges, binned_pdf, coefficients, normalization): """ Function that wraps up the moment morph procedure [1] adopted for histograms. @@ -211,7 +211,7 @@ def moment_morph_estimation(bin_edges, binned_pdf, coefficients): # Estimate mean and std for each input template histogram. First adaption needed to extend # the moment morph procedure to histograms - mus, sigs = _estimate_mean_std(bin_edges=bin_edges, binned_pdf=binned_pdf) + mus, sigs = _estimate_mean_std(bin_edges=bin_edges, binned_pdf=binned_pdf, normalization=normalization) coefficients = coefficients.reshape( binned_pdf.shape[0], *np.ones(mus.ndim - 1, "int") ) @@ -302,6 +302,7 @@ def _interpolate1D(self, target_point): bin_edges=self.bin_edges, binned_pdf=self.binned_pdf[segment_inds], coefficients=coefficients, + normalization=self.normalization, ) def _interpolate2D(self, target_point): @@ -321,6 +322,7 @@ def _interpolate2D(self, target_point): bin_edges=self.bin_edges, binned_pdf=self.binned_pdf[simplex_inds], coefficients=coefficients, + normalization=self.normalization, ) def interpolate(self, target_point): diff --git a/pyirf/interpolation/nearest_simplex_extrapolator.py b/pyirf/interpolation/nearest_simplex_extrapolator.py index e29db3176..9d473d576 100644 --- a/pyirf/interpolation/nearest_simplex_extrapolator.py +++ b/pyirf/interpolation/nearest_simplex_extrapolator.py @@ -7,7 +7,7 @@ import numpy as np from scipy.spatial import Delaunay -from .base_extrapolators import DiscretePDFExtrapolator, ParametrizedExtrapolator +from .base_extrapolators import DiscretePDFExtrapolator, ParametrizedExtrapolator, PDFNormalization from .moment_morph_interpolator import ( barycentric_2D_interpolation_coefficients, linesegment_1D_interpolation_coefficients, @@ -122,7 +122,7 @@ def extrapolate(self, target_point): class MomentMorphNearestSimplexExtrapolator(DiscretePDFExtrapolator): - def __init__(self, grid_points, bin_edges, binned_pdf): + def __init__(self, grid_points, bin_edges, binned_pdf, normalization=PDFNormalization.AREA): """ Extrapolator class extending/reusing parts of Moment Morphing by allowing for negative extrapolation coefficients computed @@ -148,7 +148,7 @@ def __init__(self, grid_points, bin_edges, binned_pdf): ---- Also calls pyirf.interpolation.DiscretePDFExtrapolator.__init__. """ - super().__init__(grid_points, bin_edges, binned_pdf) + super().__init__(grid_points, bin_edges, binned_pdf, normalization) if self.grid_dim == 2: self.triangulation = Delaunay(self.grid_points) @@ -171,6 +171,7 @@ def _extrapolate1D(self, segment_inds, target_point): bin_edges=self.bin_edges, binned_pdf=self.binned_pdf[segment_inds], coefficients=coefficients, + normalization=self.normalization, ) def _extrapolate2D(self, simplex_inds, target_point): @@ -187,6 +188,7 @@ def _extrapolate2D(self, simplex_inds, target_point): bin_edges=self.bin_edges, binned_pdf=self.binned_pdf[simplex_inds], coefficients=coefficients, + normalization=self.normalization, ) def extrapolate(self, target_point): diff --git a/pyirf/interpolation/tests/test_moment_morph_interpolator.py b/pyirf/interpolation/tests/test_moment_morph_interpolator.py index f496dfc31..8112b6994 100644 --- a/pyirf/interpolation/tests/test_moment_morph_interpolator.py +++ b/pyirf/interpolation/tests/test_moment_morph_interpolator.py @@ -2,6 +2,8 @@ import pytest from scipy.stats import norm +from pyirf.interpolation.base_interpolators import PDFNormalization + def expected_mean(a, b): return 5 + (a / 5) + (b / 15) @@ -29,7 +31,7 @@ def bins(): def simple_1D_data(bins): grid = np.array([[20], [40]]) target = np.array([30]) - binned_pdf = np.array([binned_normal_pdf([x, 0], bins) for x in grid]) + binned_pdf = np.array([binned_normal_pdf(x, 0, bins) for x in grid]) truth = binned_normal_pdf(target, 0, bins) @@ -94,7 +96,7 @@ def test_estimate_mean_std(bins): ] ).squeeze() - mean, std = _estimate_mean_std(bins, binned_pdf) + mean, std = _estimate_mean_std(bins, binned_pdf, PDFNormalization.AREA) # Assert estimation and truth within one bin np.testing.assert_allclose(mean, true_mean, atol=np.diff(bins)[0] / 2) @@ -139,14 +141,14 @@ def test_linesegment_1D_interpolation_coefficients(): np.testing.assert_allclose( linesegment_1D_interpolation_coefficients(grid_points, target_point), - np.array([1, 0]), + np.array([[1, 0]]), ) target_point = np.array([[40]]) np.testing.assert_allclose( linesegment_1D_interpolation_coefficients(grid_points, target_point), - np.array([0, 1]), + np.array([[0, 1]]), ) target_point = np.array([[25]]) @@ -215,7 +217,7 @@ def test_moment_morph_estimation1D(bins, simple_1D_data): grid, target, binned_pdf, truth = simple_1D_data.values() coeffs = linesegment_1D_interpolation_coefficients(grid, target) - res = moment_morph_estimation(bins, binned_pdf, coeffs) + res = moment_morph_estimation(bins, binned_pdf, coeffs, PDFNormalization.AREA) assert np.isclose(np.sum(res), 1) assert np.all(np.isfinite(res)) @@ -233,7 +235,7 @@ def test_moment_morph_estimation2D(bins, simple_2D_data): grid, target, binned_pdf, truth = simple_2D_data.values() coeffs = barycentric_2D_interpolation_coefficients(grid, target) - res = moment_morph_estimation(bins, binned_pdf, coeffs) + res = moment_morph_estimation(bins, binned_pdf, coeffs, PDFNormalization.AREA) assert np.isclose(np.sum(res), 1) assert np.all(np.isfinite(res)) diff --git a/pyirf/interpolation/tests/test_nearest_neighbor_searcher.py b/pyirf/interpolation/tests/test_nearest_neighbor_searcher.py index 6394b8ae4..6a3fc57bf 100644 --- a/pyirf/interpolation/tests/test_nearest_neighbor_searcher.py +++ b/pyirf/interpolation/tests/test_nearest_neighbor_searcher.py @@ -13,7 +13,7 @@ def grid_2d(): @pytest.fixture -def values(grid_1d): +def binned_pdf(grid_1d): return np.array( [ [[np.full(10, x), np.full(10, x / 2)], [np.full(10, 2 * x), np.zeros(10)]] diff --git a/pyirf/interpolation/tests/test_nearest_simplex_extrapolator.py b/pyirf/interpolation/tests/test_nearest_simplex_extrapolator.py index 384c32570..b8ca7c4a2 100644 --- a/pyirf/interpolation/tests/test_nearest_simplex_extrapolator.py +++ b/pyirf/interpolation/tests/test_nearest_simplex_extrapolator.py @@ -55,7 +55,7 @@ def test_ParametrizedNearestSimplexExtrapolator_1DGrid(): dummy_data_target2 = -2.5 * slope + 1 - assert np.allclose(interpolant2, dummy_data_target2) + np.testing.assert_allclose(interpolant2, dummy_data_target2) assert interpolant2.shape == (1, *dummy_data.shape[1:]) @@ -96,7 +96,7 @@ def test_ParametrizedNearestSimplexExtrapolator_2DGrid(): ] ).squeeze() - assert np.allclose(interpolant1.squeeze(), dummy_data_target1) + np.testing.assert_allclose(interpolant1.squeeze(), dummy_data_target1) assert interpolant1.shape == (1, *dummy_data.shape[1:]) target_point2 = np.array([[-2.5, -5.5]]) @@ -109,7 +109,7 @@ def test_ParametrizedNearestSimplexExtrapolator_2DGrid(): ] ).squeeze() - assert np.allclose(interpolant2, dummy_data_target2) + np.testing.assert_allclose(interpolant2, dummy_data_target2) assert interpolant2.shape == (1, *dummy_data.shape[1:]) @@ -188,11 +188,11 @@ def test_MomentMorphNearestSimplexExtrapolator_1DGrid(bins): res1 = extrap(target1) expected_norms1 = np.array([[0, 1], [1, 0], [1, 1]]) - assert np.allclose(np.sum(res1, axis=-1), expected_norms1) + np.testing.assert_allclose(np.sum(res1, axis=-1), expected_norms1) assert np.all(np.isfinite(res1)) assert res1.shape == (1, *binned_pdf.shape[1:]) # Assert truth and result matching within +- 0.1%, atol dominates comparison - assert np.allclose(res1.squeeze(), truth1, atol=1e-3, rtol=1e-5) + np.testing.assert_allclose(res1.squeeze(), truth1, atol=1e-3, rtol=1e-5) target2 = np.array([45]) @@ -216,11 +216,11 @@ def test_MomentMorphNearestSimplexExtrapolator_1DGrid(bins): res2 = extrap(target2) expected_norms2 = np.array([[1, 1], [1, 0], [1, 1]]) - assert np.allclose(np.sum(res2, axis=-1), expected_norms2) + np.testing.assert_allclose(np.sum(res2, axis=-1), expected_norms2) assert np.all(np.isfinite(res2)) assert res2.shape == (1, *binned_pdf.shape[1:]) # Assert truth and result matching within +- 0.1%, atol dominates comparison - assert np.allclose(res2.squeeze(), truth2, atol=1e-3, rtol=1e-5) + np.testing.assert_allclose(res2.squeeze(), truth2, atol=1e-3, rtol=1e-5) def test_MomentMorphNearestSimplexExtrapolator_2DGrid(bins): @@ -267,11 +267,11 @@ def test_MomentMorphNearestSimplexExtrapolator_2DGrid(bins): res1 = extrap(target1) expected_norms1 = np.array([[0, 1], [1, 0]]) - assert np.allclose(np.sum(res1, axis=-1), expected_norms1) + np.testing.assert_allclose(np.sum(res1, axis=-1), expected_norms1) assert np.all(np.isfinite(res1)) assert res1.shape == (1, *binned_pdf.shape[1:]) # Assert truth and result matching within +- 0.1%, atol dominates comparison - assert np.allclose(res1.squeeze(), truth1, atol=1e-3, rtol=1e-5) + np.testing.assert_allclose(res1.squeeze(), truth1, atol=1e-3, rtol=1e-5) target2 = np.array([45, 20]) @@ -291,11 +291,11 @@ def test_MomentMorphNearestSimplexExtrapolator_2DGrid(bins): res2 = extrap(target2) expected_norms2 = np.array([[1, 1], [1, 0]]) - assert np.allclose(np.sum(res2, axis=-1), expected_norms2) + np.testing.assert_allclose(np.sum(res2, axis=-1), expected_norms2) assert np.all(np.isfinite(res2)) assert res2.shape == (1, *binned_pdf.shape[1:]) # Assert truth and result matching within +- 0.1%, atol dominates comparison - assert np.allclose(res2.squeeze(), truth2, atol=1e-3, rtol=1e-5) + np.testing.assert_allclose(res2.squeeze(), truth2, atol=1e-3, rtol=1e-5) def test_MomentMorphNearestSimplexExtrapolator_3DGrid(): @@ -344,5 +344,5 @@ def test_MomentMorphNearestSimplexExtrapolator_below_zero_warning(bins): ): res = extrap(np.array([0])) - assert np.allclose(np.sum(res, axis=-1), 1) + np.testing.assert_allclose(np.sum(res, axis=-1), 1) assert np.all(res >= 0) From ba7c0ddddf1a593d699074c6e28ba85daf9bdb1b Mon Sep 17 00:00:00 2001 From: Maximilian Linhoff Date: Tue, 22 Aug 2023 12:44:42 +0200 Subject: [PATCH 09/14] Adapt interpolation to new normalization of binned pdfs --- pyirf/interpolation/base_interpolators.py | 3 +- pyirf/interpolation/component_estimators.py | 12 +-- .../moment_morph_interpolator.py | 3 +- .../nearest_simplex_extrapolator.py | 5 +- pyirf/interpolation/quantile_interpolator.py | 7 +- .../test_component_estimator_base_classes.py | 30 +++++--- ...st_component_estimator_specific_classes.py | 16 ++-- .../tests/test_moment_morph_interpolator.py | 31 ++++---- .../test_nearest_simplex_extrapolator.py | 73 +++++++++++-------- .../tests/test_quantile_interpolator.py | 17 ++++- pyproject.toml | 9 ++- 11 files changed, 124 insertions(+), 82 deletions(-) diff --git a/pyirf/interpolation/base_interpolators.py b/pyirf/interpolation/base_interpolators.py index cf4cf2ffa..6b0cb8294 100644 --- a/pyirf/interpolation/base_interpolators.py +++ b/pyirf/interpolation/base_interpolators.py @@ -1,6 +1,7 @@ """Base classes for interpolators""" from abc import ABCMeta, abstractmethod import enum +import astropy.units as u import numpy as np @@ -31,7 +32,7 @@ def get_bin_width(bin_edges, normalization): return np.diff(bin_edges) if normalization is PDFNormalization.CONE_SOLID_ANGLE: - return np.diff(cone_solid_angle(bin_edges)) + return np.diff(cone_solid_angle(bin_edges).to_value(u.sr)) raise ValueError(f"Invalid PDF normalization: {normalization}") diff --git a/pyirf/interpolation/component_estimators.py b/pyirf/interpolation/component_estimators.py index 6ed0edc5d..dc859c868 100644 --- a/pyirf/interpolation/component_estimators.py +++ b/pyirf/interpolation/component_estimators.py @@ -687,10 +687,6 @@ def __init__( psf = np.swapaxes(psf, axis, -1) - # Renormalize along the source offset axis to have a proper PDF - self.omegas = np.diff(cone_solid_angle(source_offset_bins)) - psf_normed = psf * self.omegas - if interpolator_kwargs is None: interpolator_kwargs = {} @@ -706,8 +702,8 @@ def __init__( super().__init__( grid_points=grid_points, - bin_edges=source_offset_bins, - binned_pdf=psf_normed, + bin_edges=source_offset_bins.to_value(u.rad), + binned_pdf=psf, interpolator_cls=interpolator_cls, interpolator_kwargs=interpolator_kwargs, extrapolator_cls=extrapolator_cls, @@ -726,7 +722,7 @@ def __call__(self, target_point): Returns ------- - psf_interp: np.ndarray, shape=(n_points, ..., n_source_offset_bins) + psf_interp: u.Quantity[sr-1], shape=(n_points, ..., n_source_offset_bins) Interpolated psf table with same shape as input matrices. For PSF_TABLE of shape (n_points, n_energy_bins, n_fov_offset_bins, n_source_offset_bins) @@ -735,4 +731,4 @@ def __call__(self, target_point): interpolated_psf_normed = super().__call__(target_point) # Undo normalisation to get a proper PSF and return - return np.swapaxes(interpolated_psf_normed / self.omegas, -1, self.axis) + return u.Quantity(np.swapaxes(interpolated_psf_normed, -1, self.axis), u.sr**-1, copy=False) diff --git a/pyirf/interpolation/moment_morph_interpolator.py b/pyirf/interpolation/moment_morph_interpolator.py index 6c102d66b..cb0c748b8 100644 --- a/pyirf/interpolation/moment_morph_interpolator.py +++ b/pyirf/interpolation/moment_morph_interpolator.py @@ -245,7 +245,8 @@ def moment_morph_estimation(bin_edges, binned_pdf, coefficients, normalization): # Re-Normalize, needed, as the estimation of the std used above is not exact but the result is scaled with # the estimated std - norm = np.expand_dims(np.sum(f_new, axis=-1), -1) + bin_widths = get_bin_width(bin_edges, normalization) + norm = np.expand_dims(np.sum(f_new * bin_widths, axis=-1), -1) return np.divide(f_new, norm, out=np.zeros_like(f_new), where=norm != 0).reshape( 1, *binned_pdf.shape[1:] diff --git a/pyirf/interpolation/nearest_simplex_extrapolator.py b/pyirf/interpolation/nearest_simplex_extrapolator.py index 9d473d576..c108766db 100644 --- a/pyirf/interpolation/nearest_simplex_extrapolator.py +++ b/pyirf/interpolation/nearest_simplex_extrapolator.py @@ -7,6 +7,8 @@ import numpy as np from scipy.spatial import Delaunay +from pyirf.interpolation.base_interpolators import get_bin_width + from .base_extrapolators import DiscretePDFExtrapolator, ParametrizedExtrapolator, PDFNormalization from .moment_morph_interpolator import ( barycentric_2D_interpolation_coefficients, @@ -230,7 +232,8 @@ def extrapolate(self, target_point): # cut off values below 0 and re-normalize extrapolant[extrapolant < 0] = 0 - norm = np.expand_dims(np.sum(extrapolant, axis=-1), -1) + bin_width = get_bin_width(self.bin_edges, self.normalization) + norm = np.expand_dims(np.sum(extrapolant * bin_width, axis=-1), -1) return np.divide( extrapolant, norm, out=np.zeros_like(extrapolant), where=norm != 0 ).reshape(1, *self.binned_pdf.shape[1:]) diff --git a/pyirf/interpolation/quantile_interpolator.py b/pyirf/interpolation/quantile_interpolator.py index 3bc9583fc..bdc6c6d42 100644 --- a/pyirf/interpolation/quantile_interpolator.py +++ b/pyirf/interpolation/quantile_interpolator.py @@ -132,16 +132,17 @@ def interpolate_ppf(xy): return pdf_values -def norm_pdf(pdf_values): +def norm_pdf(pdf_values, bin_edges, normalization): """ Normalize binned_pdf to a sum of 1 """ norm = np.sum(pdf_values, axis=-1) + width = get_bin_width(bin_edges, normalization) # Norm all binned_pdfs to unity that are not empty normed_pdf_values = np.divide( pdf_values, - norm[..., np.newaxis], + norm[..., np.newaxis] * width, out=np.zeros_like(pdf_values), where=norm[..., np.newaxis] != 0, ) @@ -255,7 +256,7 @@ def interpolate(self, target_point): )[..., 1:] # Renormalize pdf to sum of 1 - normed_interpolated_pdfs = norm_pdf(interpolated_pdfs) + normed_interpolated_pdfs = norm_pdf(interpolated_pdfs, self.bin_edges[1:], self.normalization) # Re-swap axes and set all nans to zero return np.nan_to_num(normed_interpolated_pdfs).reshape(1, *self.input_shape[1:]) diff --git a/pyirf/interpolation/tests/test_component_estimator_base_classes.py b/pyirf/interpolation/tests/test_component_estimator_base_classes.py index 534fdd209..816c2ed35 100644 --- a/pyirf/interpolation/tests/test_component_estimator_base_classes.py +++ b/pyirf/interpolation/tests/test_component_estimator_base_classes.py @@ -127,10 +127,12 @@ def __init__(self, grid_points): assert estim2D(target2D_inGrid) == 42 estim1D = DummyEstimator(grid_points1D_good) - assert estim1D(target1D_outofGrid) == 43 + with pytest.warns(UserWarning, match="has to be extrapolated"): + assert estim1D(target1D_outofGrid) == 43 estim2D = DummyEstimator(grid_points2D_good) - assert estim2D(target2D_outofGrid) == 43 + with pytest.warns(UserWarning, match="has to be extrapolated"): + assert estim2D(target2D_outofGrid) == 43 def test_ParametrizedComponentEstimator_checks(): @@ -208,7 +210,8 @@ def extrapolate(self, target_point): extrapolator_cls=DummyExtrapolator, ) assert estim(np.array([[1.5]])) == 42 - assert estim(np.array([[0]])) == 43 + with pytest.warns(UserWarning, match="has to be extrapolated"): + assert estim(np.array([[0]])) == 43 def test_DiscretePDFComponentEstimator_checks(): @@ -318,7 +321,8 @@ def extrapolate(self, target_point): extrapolator_cls=DummyExtrapolator, ) assert estim(np.array([[1.5]])) == 42 - assert estim(np.array([[0]])) == 43 + with pytest.warns(UserWarning, match="has to be extrapolated"): + assert estim(np.array([[0]])) == 43 def test_DiscretePDFComponentEstimator_NearestNeighbors(): @@ -344,7 +348,8 @@ def test_DiscretePDFComponentEstimator_NearestNeighbors(): ) assert np.allclose(estim(target_point=np.array([1.1])), binned_pdf[0, :]) - assert np.allclose(estim(target_point=np.array([4.1])), binned_pdf[2, :]) + with pytest.warns(UserWarning, match="has to be extrapolated"): + assert np.allclose(estim(target_point=np.array([4.1])), binned_pdf[2, :]) with pytest.raises( TypeError, @@ -382,7 +387,8 @@ def test_ParametrizedComponentEstimator_NearestNeighbors(): ) assert np.allclose(estim(target_point=np.array([1.1])), params[0, :]) - assert np.allclose(estim(target_point=np.array([4.1])), params[2, :]) + with pytest.warns(UserWarning, match="has to be extrapolated"): + assert np.allclose(estim(target_point=np.array([4.1])), params[2, :]) with pytest.raises( TypeError, @@ -428,9 +434,11 @@ def extrapolate(self, target_point): ) # Nearest neighbor is grid_point 1 at the index 1 of the original binned_pdf - assert np.allclose(estim(target_point=np.array([0])), binned_pdf[1, :]) + with pytest.warns(UserWarning, match="has to be extrapolated"): + assert np.allclose(estim(target_point=np.array([0])), binned_pdf[1, :]) # Nearest neighbor is grid_point 3 at the index 0 of the original binned_pdf - assert np.allclose(estim(target_point=np.array([4])), binned_pdf[0, :]) + with pytest.warns(UserWarning, match="has to be extrapolated"): + assert np.allclose(estim(target_point=np.array([4])), binned_pdf[0, :]) def test_ParametrizedComponentEstimator_1Dsorting(): @@ -461,6 +469,8 @@ def extrapolate(self, target_point): ) # Nearest neighbor is grid_point 1 at the index 1 of the original binned_pdf - assert np.allclose(estim(target_point=np.array([0])), params[1, :]) + with pytest.warns(UserWarning, match="has to be extrapolated"): + assert np.allclose(estim(target_point=np.array([0])), params[1, :]) # Nearest neighbor is grid_point 3 at the index 0 of the original binned_pdf - assert np.allclose(estim(target_point=np.array([4])), params[0, :]) + with pytest.warns(UserWarning, match="has to be extrapolated"): + assert np.allclose(estim(target_point=np.array([4])), params[0, :]) diff --git a/pyirf/interpolation/tests/test_component_estimator_specific_classes.py b/pyirf/interpolation/tests/test_component_estimator_specific_classes.py index 49531e4ce..c1f77913b 100644 --- a/pyirf/interpolation/tests/test_component_estimator_specific_classes.py +++ b/pyirf/interpolation/tests/test_component_estimator_specific_classes.py @@ -10,6 +10,7 @@ def test_EnergyDispersionEstimator(prod5_irfs): zen_pnt = np.array([key.value for key in prod5_irfs.keys()]) edisps = np.array([irf["edisp"].data for irf in prod5_irfs.values()]) bin_edges = list(prod5_irfs.values())[0]["edisp"].axes["migra"].edges + bin_width = np.diff(bin_edges) estimator = EnergyDispersionEstimator( grid_points=zen_pnt[[0, 2]], @@ -24,13 +25,12 @@ def test_EnergyDispersionEstimator(prod5_irfs): interp = estimator(target_point=zen_pnt[[1]]) - assert np.max(interp) <= 1 assert np.min(interp) >= 0 assert np.all(np.isfinite(interp)) assert np.all( np.logical_or( - np.isclose(np.sum(interp, axis=-2), 1), - np.isclose(np.sum(interp, axis=-2), 0), + np.isclose(np.sum(interp * bin_width[:, np.newaxis], axis=-2), 1), + np.isclose(np.sum(interp * bin_width[:, np.newaxis], axis=-2), 0), ) ) assert interp.shape == edisps[[1]].shape @@ -74,15 +74,15 @@ def hist(pnt): interp = estimator(target_point=zen_pnt[[1]]) - interp *= omegas[np.newaxis, np.newaxis, np.newaxis, ...] + probability = (interp * omegas[np.newaxis, np.newaxis, np.newaxis, ...]).to_value(u.one) - assert np.max(interp) <= 1 - assert np.min(interp) >= 0 + assert np.max(probability) <= 1 + assert np.min(probability) >= 0 assert np.all(np.isfinite(interp)) assert np.all( np.logical_or( - np.isclose(np.sum(interp, axis=-1), 1), - np.isclose(np.sum(interp, axis=-1), 0), + np.isclose(np.sum(probability, axis=-1), 1), + np.isclose(np.sum(probability, axis=-1), 0), ) ) assert interp.shape == dummy_psfs[[1]].shape diff --git a/pyirf/interpolation/tests/test_moment_morph_interpolator.py b/pyirf/interpolation/tests/test_moment_morph_interpolator.py index 8112b6994..77915ccb5 100644 --- a/pyirf/interpolation/tests/test_moment_morph_interpolator.py +++ b/pyirf/interpolation/tests/test_moment_morph_interpolator.py @@ -158,7 +158,7 @@ def test_linesegment_1D_interpolation_coefficients(): mfrac = (target_point[0, 0] - grid_points[0, 0]) / ( grid_points[1, 0] - grid_points[0, 0] ) - cs = np.array([1 - mfrac, mfrac]) + cs = np.array([[1 - mfrac, mfrac]]) res = linesegment_1D_interpolation_coefficients(grid_points, target_point) @@ -219,7 +219,7 @@ def test_moment_morph_estimation1D(bins, simple_1D_data): coeffs = linesegment_1D_interpolation_coefficients(grid, target) res = moment_morph_estimation(bins, binned_pdf, coeffs, PDFNormalization.AREA) - assert np.isclose(np.sum(res), 1) + np.testing.assert_almost_equal(np.sum(res * np.diff(bins)), 1) assert np.all(np.isfinite(res)) assert res.shape == (1, len(bins) - 1) # Assert truth and result matching within +- 0.1%, atol dominates comparison @@ -237,7 +237,7 @@ def test_moment_morph_estimation2D(bins, simple_2D_data): coeffs = barycentric_2D_interpolation_coefficients(grid, target) res = moment_morph_estimation(bins, binned_pdf, coeffs, PDFNormalization.AREA) - assert np.isclose(np.sum(res), 1) + np.testing.assert_almost_equal(np.sum(res * np.diff(bins)), 1) assert np.all(np.isfinite(res)) assert res.shape == (1, len(bins) - 1) # Assert truth and result matching within +- 0.1%, atol dominates comparison @@ -255,7 +255,7 @@ def test_MomentMorphInterpolator1D(bins, simple_1D_data): res = interp(target) - assert np.isclose(np.sum(res), 1) + np.testing.assert_almost_equal(np.sum(res * np.diff(bins)), 1) assert np.all(np.isfinite(res)) assert res.shape == (1, len(bins) - 1) # Assert truth and result matching within +- 0.1%, atol dominates comparison @@ -278,7 +278,8 @@ def test_MomentMorphInterpolator1D_dirac_delta_input(): interp = MomentMorphInterpolator(grid, bin_edges, binned_pdf) res = interp(target) - np.testing.assert_allclose(res, np.array([[0, 0, 1, 0], [0.25, 0.25, 0.25, 0.25]])) + expected = np.array([[[0, 0, 1, 0], [0.25, 0.25, 0.25, 0.25]]]) + np.testing.assert_allclose(res, expected) def test_MomentMorphInterpolator1D_all_empty(bins, simple_1D_data): @@ -358,8 +359,8 @@ def test_MomentMorphInterpolator1D_mixed_data(bins): res = interp(target) - expected_norms = np.array([[0, 1], [1, 0]]) - np.testing.assert_allclose(np.sum(res, axis=-1), expected_norms) + expected_norms = np.array([[[0, 1], [1, 0]]]) + np.testing.assert_allclose(np.sum(res * np.diff(bins), axis=-1), expected_norms) assert np.all(np.isfinite(res)) assert res.shape == (1, *binned_pdf.shape[1:]) # Assert truth and result matching within +- 0.1%, atol dominates comparison @@ -403,11 +404,11 @@ def test_MomentMorphInterpolator1D_extended_grid_extradims(bins): res = interp(target) - np.testing.assert_allclose(np.sum(res, axis=-1), 1) + np.testing.assert_allclose(np.sum(res * np.diff(bins), axis=-1), 1) assert np.all(np.isfinite(res)) assert res.shape == (1, *binned_pdf.shape[1:]) # Assert truth and result matching within +- 0.1%, atol dominates comparison - np.testing.assert_allclose(res.squeeze(), truth, atol=1e-4, rtol=1e-5) + np.testing.assert_allclose(res.squeeze(), truth, atol=1e-3, rtol=1e-5) def test_MomentMorphInterpolator2D(bins, simple_2D_data): @@ -421,7 +422,7 @@ def test_MomentMorphInterpolator2D(bins, simple_2D_data): res = interp(target) - assert np.isclose(np.sum(res), 1) + np.testing.assert_almost_equal(np.sum(res * np.diff(bins)), 1) assert np.all(np.isfinite(res)) assert res.shape == (1, len(bins) - 1) # Assert truth and result matching within +- 0.1%, atol dominates comparison @@ -505,8 +506,8 @@ def test_MomentMorphInterpolator2D_mixed(bins): res = interp(target) - expected_norms = np.array([[0, 1], [1, 0]]) - np.testing.assert_allclose(np.sum(res, axis=-1), expected_norms) + expected_norms = np.array([[[0, 1], [1, 0]]]) + np.testing.assert_allclose(np.sum(res * np.diff(bins), axis=-1), expected_norms) assert np.all(np.isfinite(res)) assert res.shape == (1, *binned_pdf.shape[1:]) # Assert truth and result matching within +- 0.1%, atol dominates comparison @@ -529,7 +530,7 @@ def test_MomentMorphInterpolator1D_extended_grid(bins): res = interp(target) truth = binned_normal_pdf(target, 0, bins) - assert np.isclose(np.sum(res), 1) + np.testing.assert_almost_equal(np.sum(res * np.diff(bins)), 1) assert np.all(np.isfinite(res)) assert res.shape == (1, len(bins) - 1) # Assert truth and result matching within +- 0.1%, atol dominates comparison @@ -552,7 +553,7 @@ def test_MomentMorphInterpolator2D_extended_grid(bins): res = interp(target) truth = binned_normal_pdf(*target, bins) - assert np.isclose(np.sum(res), 1) + np.testing.assert_almost_equal(np.sum(res * np.diff(bins)), 1) assert np.all(np.isfinite(res)) assert res.shape == (1, len(bins) - 1) # Assert truth and result matching within +- 0.1%, atol dominates comparison @@ -600,7 +601,7 @@ def test_MomentMorphInterpolator2D_extended_grid_extradims(bins): res = interp(target) - np.testing.assert_allclose(np.sum(res, axis=-1), 1) + np.testing.assert_allclose(np.sum(res * np.diff(bins), axis=-1), 1) assert np.all(np.isfinite(res)) assert res.shape == (1, *binned_pdf.shape[1:]) # Assert truth and result matching within +- 0.1%, atol dominates comparison diff --git a/pyirf/interpolation/tests/test_nearest_simplex_extrapolator.py b/pyirf/interpolation/tests/test_nearest_simplex_extrapolator.py index b8ca7c4a2..313cdc727 100644 --- a/pyirf/interpolation/tests/test_nearest_simplex_extrapolator.py +++ b/pyirf/interpolation/tests/test_nearest_simplex_extrapolator.py @@ -45,7 +45,7 @@ def test_ParametrizedNearestSimplexExtrapolator_1DGrid(): target_point1 = np.array([3]) interpolant1 = interpolator(target_point1) - dummy_data_target1 = 3 * slope + 1 + dummy_data_target1 = (3 * slope + 1)[np.newaxis, :] np.testing.assert_allclose(interpolant1, dummy_data_target1) assert interpolant1.shape == (1, *dummy_data.shape[1:]) @@ -53,7 +53,7 @@ def test_ParametrizedNearestSimplexExtrapolator_1DGrid(): target_point2 = np.array([[-2.5]]) interpolant2 = interpolator(target_point2) - dummy_data_target2 = -2.5 * slope + 1 + dummy_data_target2 = (-2.5 * slope + 1)[np.newaxis, :] np.testing.assert_allclose(interpolant2, dummy_data_target2) assert interpolant2.shape == (1, *dummy_data.shape[1:]) @@ -104,10 +104,12 @@ def test_ParametrizedNearestSimplexExtrapolator_2DGrid(): dummy_data_target2 = np.array( [ - np.dot((m.T * target_point2 + n), np.array([1, 1])) - for m, n in zip(slope, intercept) + [ + np.dot((m.T * target_point2 + n), np.array([1, 1])) + for m, n in zip(slope, intercept) + ] ] - ).squeeze() + ) np.testing.assert_allclose(interpolant2, dummy_data_target2) assert interpolant2.shape == (1, *dummy_data.shape[1:]) @@ -137,6 +139,9 @@ def test_MomentMorphNearestSimplexExtrapolator_1DGrid(bins): grid_points = np.array([[20], [30], [40]]) + n_bins = len(bins) - 1 + bin_width = np.diff(bins) + # Create template histograms binned_pdf = np.array( [ @@ -144,22 +149,21 @@ def test_MomentMorphNearestSimplexExtrapolator_1DGrid(bins): [binned_normal_pdf(x, 0, bins), binned_normal_pdf(x + 1, 0, bins)], [ binned_normal_pdf(x + 10, 0, bins), - np.zeros(len(bins) - 1), + np.zeros(n_bins), ], - [np.zeros(len(bins) - 1), np.ones(len(bins) - 1) / (len(bins) - 1)], + [np.zeros(n_bins), np.ones(n_bins) / n_bins / bin_width], ] for x in grid_points ] ) # Include dirac-delta moving with grid_point x at [:, 2, 0, x] - binned_pdf[0, 2, 0, 20] = 1 - binned_pdf[1, 2, 0, 30] = 1 - binned_pdf[2, 2, 0, 40] = 1 + for i, bin_idx in enumerate(grid_points.ravel()): + binned_pdf[i, 2, 0, bin_idx] = 1 / bin_width[bin_idx] # Zero template histogram at index [0, 0, 0, :], extrapolations "lower" # then this bin has consequently to be zero in these bins. - binned_pdf[0, 0, 0, :] = np.zeros(len(bins) - 1) + binned_pdf[0, 0, 0, :] = 0.0 extrap = MomentMorphNearestSimplexExtrapolator( grid_points=grid_points, binned_pdf=binned_pdf, bin_edges=bins @@ -171,28 +175,31 @@ def test_MomentMorphNearestSimplexExtrapolator_1DGrid(bins): truth1 = np.array( [ [ - np.zeros(len(bins) - 1), + np.zeros(n_bins), binned_normal_pdf(target1 + 1, 0, bins), ], [ binned_normal_pdf(target1 + 10, 0, bins), - np.zeros(len(bins) - 1), + np.zeros_like(bin_width), ], - [np.zeros(len(bins) - 1), np.ones(len(bins) - 1) / (len(bins) - 1)], + [np.zeros(n_bins), np.ones(n_bins) / n_bins / bin_width], ] ) # Inlcude dirac-delta - truth1[2, 0, 10] = 1 + delta = np.zeros(len(bins) - 1) + delta[target1[0]] = 1 + delta /= bin_width + truth1[2, 0] = delta res1 = extrap(target1) - expected_norms1 = np.array([[0, 1], [1, 0], [1, 1]]) - np.testing.assert_allclose(np.sum(res1, axis=-1), expected_norms1) + expected_norms1 = np.array([[[0, 1], [1, 0], [1, 1]]]) + np.testing.assert_allclose(np.sum(res1 * np.diff(bins), axis=-1), expected_norms1) assert np.all(np.isfinite(res1)) assert res1.shape == (1, *binned_pdf.shape[1:]) # Assert truth and result matching within +- 0.1%, atol dominates comparison - np.testing.assert_allclose(res1.squeeze(), truth1, atol=1e-3, rtol=1e-5) + np.testing.assert_allclose(res1.squeeze(), truth1, atol=4e-3, rtol=1e-5) target2 = np.array([45]) @@ -204,19 +211,19 @@ def test_MomentMorphNearestSimplexExtrapolator_1DGrid(bins): ], [ binned_normal_pdf(target2 + 10, 0, bins), - np.zeros(len(bins) - 1), + np.zeros(n_bins), ], - [np.zeros(len(bins) - 1), np.ones(len(bins) - 1) / (len(bins) - 1)], + [np.zeros(n_bins), np.ones(n_bins) / n_bins / bin_width], ] ) # Inlcude dirac-delta - truth2[2, 0, 45] = 1 + truth2[2, 0, 45] = 1 / bin_width[45] res2 = extrap(target2) - expected_norms2 = np.array([[1, 1], [1, 0], [1, 1]]) - np.testing.assert_allclose(np.sum(res2, axis=-1), expected_norms2) + expected_norms2 = np.array([[[1, 1], [1, 0], [1, 1]]]) + np.testing.assert_allclose(np.sum(res2 * bin_width, axis=-1), expected_norms2) assert np.all(np.isfinite(res2)) assert res2.shape == (1, *binned_pdf.shape[1:]) # Assert truth and result matching within +- 0.1%, atol dominates comparison @@ -227,6 +234,8 @@ def test_MomentMorphNearestSimplexExtrapolator_2DGrid(bins): """Test ParametrizedNearestSimplexExtrapolator on a 2D Grid""" from pyirf.interpolation import MomentMorphNearestSimplexExtrapolator + n_bins = len(bins) - 1 + bin_width = np.diff(bins) grid_points = np.array([[20, 20], [30, 10], [40, 20], [50, 10]]) # Create template histograms @@ -234,7 +243,7 @@ def test_MomentMorphNearestSimplexExtrapolator_2DGrid(bins): [ [ [binned_normal_pdf(a, b, bins), binned_normal_pdf(a + 1, b, bins)], - [binned_normal_pdf(a + 10, b, bins), np.zeros(len(bins) - 1)], + [binned_normal_pdf(a + 10, b, bins), np.zeros(n_bins)], ] for a, b in grid_points ] @@ -242,7 +251,7 @@ def test_MomentMorphNearestSimplexExtrapolator_2DGrid(bins): # Zero template histogram at index [0, 0, 0, :] (the [20, 20] point), # so that targets interpolated from left simplex should also be zero at [0, 0, :] - binned_pdf[0, 0, 0, :] = np.zeros(len(bins) - 1) + binned_pdf[0, 0, 0, :] = np.zeros(n_bins) extrap = MomentMorphNearestSimplexExtrapolator( grid_points=grid_points, binned_pdf=binned_pdf, bin_edges=bins @@ -254,20 +263,20 @@ def test_MomentMorphNearestSimplexExtrapolator_2DGrid(bins): truth1 = np.array( [ [ - np.zeros(len(bins) - 1), + np.zeros(n_bins), binned_normal_pdf(target1[0] + 1, target1[1], bins), ], [ binned_normal_pdf(target1[0] + 10, target1[1], bins), - np.zeros(len(bins) - 1), + np.zeros(n_bins), ], ] ) res1 = extrap(target1) - expected_norms1 = np.array([[0, 1], [1, 0]]) - np.testing.assert_allclose(np.sum(res1, axis=-1), expected_norms1) + expected_norms1 = np.array([[[0, 1], [1, 0]]]) + np.testing.assert_allclose(np.sum(res1 * bin_width, axis=-1), expected_norms1) assert np.all(np.isfinite(res1)) assert res1.shape == (1, *binned_pdf.shape[1:]) # Assert truth and result matching within +- 0.1%, atol dominates comparison @@ -283,15 +292,15 @@ def test_MomentMorphNearestSimplexExtrapolator_2DGrid(bins): ], [ binned_normal_pdf(target2[0] + 10, target2[1], bins), - np.zeros(len(bins) - 1), + np.zeros(n_bins), ], ] ) res2 = extrap(target2) - expected_norms2 = np.array([[1, 1], [1, 0]]) - np.testing.assert_allclose(np.sum(res2, axis=-1), expected_norms2) + expected_norms2 = np.array([[[1, 1], [1, 0]]]) + np.testing.assert_allclose(np.sum(res2 * bin_width, axis=-1), expected_norms2) assert np.all(np.isfinite(res2)) assert res2.shape == (1, *binned_pdf.shape[1:]) # Assert truth and result matching within +- 0.1%, atol dominates comparison diff --git a/pyirf/interpolation/tests/test_quantile_interpolator.py b/pyirf/interpolation/tests/test_quantile_interpolator.py index c7a19163c..2a1892cd0 100644 --- a/pyirf/interpolation/tests/test_quantile_interpolator.py +++ b/pyirf/interpolation/tests/test_quantile_interpolator.py @@ -103,10 +103,23 @@ def test_pdf_from_ppf(data): def test_norm_pdf(data): from pyirf.interpolation.quantile_interpolator import norm_pdf + result = norm_pdf( + 2 * data["binned_pdfs"][0], + data["bin_edges"], + PDFNormalization.AREA, + ) + np.testing.assert_allclose(np.sum(result * np.diff(data["bin_edges"])), 1) np.testing.assert_allclose( - norm_pdf(2 * data["binned_pdfs"][0]), data["binned_pdfs"][0] + result, + data["binned_pdfs"][0], + ) + + result = norm_pdf( + np.zeros(len(data["bin_edges"]) - 1), + data["bin_edges"], + PDFNormalization.AREA, ) - np.testing.assert_allclose(norm_pdf(np.zeros(5)), 0) + np.testing.assert_allclose(result, 0) def test_interpolate_binned_pdf(data): diff --git a/pyproject.toml b/pyproject.toml index 51f3b9859..95c0fe5d2 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -32,4 +32,11 @@ build-backend = "setuptools.build_meta" [[tool.towncrier.section]] name = "" - path = "" \ No newline at end of file + path = "" + + +[tool.pytest.ini_options] +filterwarnings = [ + "ignore:pkg_resources is deprecated:DeprecationWarning", + "ignore:Deprecated call to `pkg_resources:DeprecationWarning", +] From 1cd2c58487d4d610146fe2fecf4d6e2cca127107 Mon Sep 17 00:00:00 2001 From: Maximilian Linhoff Date: Tue, 22 Aug 2023 12:54:18 +0200 Subject: [PATCH 10/14] Fix typo --- pyirf/interpolation/base_interpolators.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyirf/interpolation/base_interpolators.py b/pyirf/interpolation/base_interpolators.py index 6b0cb8294..99ecb44af 100644 --- a/pyirf/interpolation/base_interpolators.py +++ b/pyirf/interpolation/base_interpolators.py @@ -23,7 +23,7 @@ class PDFNormalization(enum.Enum): #: PDF is normalized to a "normal" area integral of 1 AREA = enum.auto() #: PDF is normalized to 1 over the solid angle integral where the bin - #: edges represent the opening angles of cones in radean. + #: edges represent the opening angles of cones in radian. CONE_SOLID_ANGLE = enum.auto() From f4a929ca5359c87cb938e0ee69d83a3f8a34a09b Mon Sep 17 00:00:00 2001 From: Maximilian Linhoff Date: Tue, 22 Aug 2023 16:02:58 +0200 Subject: [PATCH 11/14] Assume uniform distribution for pdfs where only a single bin has a value != 0 in moment morphing --- pyirf/interpolation/moment_morph_interpolator.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/pyirf/interpolation/moment_morph_interpolator.py b/pyirf/interpolation/moment_morph_interpolator.py index cb0c748b8..e3f8520fe 100644 --- a/pyirf/interpolation/moment_morph_interpolator.py +++ b/pyirf/interpolation/moment_morph_interpolator.py @@ -46,10 +46,10 @@ def _estimate_mean_std(bin_edges, binned_pdf, normalization): # binned_pdf!=0 mask = std == 0 if np.any(mask): - # Create array of bin_widths and clone for each template - width = np.diff(bin_edges) / 2 - width = np.repeat(width[np.newaxis, :], binned_pdf.shape[0], axis=0) - std[mask] = width[binned_pdf[mask, :] != 0] + width = np.diff(bin_edges) + # std of a uniform distribution inside the bin + uniform_std = np.broadcast_to(np.sqrt(1/12) * width, binned_pdf[mask].shape) + std[mask] = uniform_std[binned_pdf[mask, :] != 0] return mean, std From e02405cfb300b36c775ff5995744405a3eb65b07 Mon Sep 17 00:00:00 2001 From: Maximilian Linhoff Date: Tue, 22 Aug 2023 16:36:34 +0200 Subject: [PATCH 12/14] Add changelog entry --- docs/changes/250.api.rst | 6 ++++++ 1 file changed, 6 insertions(+) create mode 100644 docs/changes/250.api.rst diff --git a/docs/changes/250.api.rst b/docs/changes/250.api.rst new file mode 100644 index 000000000..97a5b1d44 --- /dev/null +++ b/docs/changes/250.api.rst @@ -0,0 +1,6 @@ +In prior versions of pyirf, the energy dispersion matrix was normalized to a +sum of 1 over the migration axis. +This is wrong, the correct normalization is to an integral of 1, which is fixed now. + +The internal API of the interpolation functions had to be adapted to take in additional +keywords, mainly the bin edges and the kind of normalization (standard or solid angle cone sections). From db105c1836bee06fae1ce8268038177779db7fb0 Mon Sep 17 00:00:00 2001 From: Maximilian Linhoff Date: Tue, 22 Aug 2023 17:02:23 +0200 Subject: [PATCH 13/14] Add missing tests/__init__.py --- pyirf/irf/tests/__init__.py | 0 1 file changed, 0 insertions(+), 0 deletions(-) create mode 100644 pyirf/irf/tests/__init__.py diff --git a/pyirf/irf/tests/__init__.py b/pyirf/irf/tests/__init__.py new file mode 100644 index 000000000..e69de29bb From 116464413c9696a33e86d45b6b10db4ef29e5ccf Mon Sep 17 00:00:00 2001 From: Maximilian Linhoff Date: Tue, 22 Aug 2023 17:26:10 +0200 Subject: [PATCH 14/14] Address review comments --- pyirf/interpolation/__init__.py | 2 ++ pyirf/interpolation/base_extrapolators.py | 9 +++++---- pyirf/interpolation/base_interpolators.py | 13 ------------- pyirf/interpolation/moment_morph_interpolator.py | 9 ++++++++- .../interpolation/nearest_simplex_extrapolator.py | 5 +++-- pyirf/interpolation/quantile_interpolator.py | 3 ++- pyirf/interpolation/tests/test_utils.py | 12 ++++++++++++ pyirf/interpolation/utils.py | 14 ++++++++++++++ 8 files changed, 46 insertions(+), 21 deletions(-) diff --git a/pyirf/interpolation/__init__.py b/pyirf/interpolation/__init__.py index 47d71cb84..2863650b5 100644 --- a/pyirf/interpolation/__init__.py +++ b/pyirf/interpolation/__init__.py @@ -11,6 +11,7 @@ BaseInterpolator, DiscretePDFInterpolator, ParametrizedInterpolator, + PDFNormalization, ) from .component_estimators import ( BaseComponentEstimator, @@ -39,6 +40,7 @@ "BaseInterpolator", "BaseNearestNeighborSearcher", "BaseExtrapolator", + "PDFNormalization", "DiscretePDFExtrapolator", "ParametrizedExtrapolator", "DiscretePDFComponentEstimator", diff --git a/pyirf/interpolation/base_extrapolators.py b/pyirf/interpolation/base_extrapolators.py index d848009c4..28b3b38a8 100644 --- a/pyirf/interpolation/base_extrapolators.py +++ b/pyirf/interpolation/base_extrapolators.py @@ -89,16 +89,17 @@ def __init__(self, grid_points, bin_edges, binned_pdf, normalization=PDFNormaliz Parameters ---------- - grid_points: np.ndarray, shape=(n_points, n_dims) + grid_points : np.ndarray, shape=(n_points, n_dims) Grid points at which templates exist - bin_edges: np.ndarray, shape=(n_bins+1) + bin_edges : np.ndarray, shape=(n_bins+1) Edges of the data binning - binned_pdf: np.ndarray, shape=(n_points, ..., n_bins) + binned_pdf : np.ndarray, shape=(n_points, ..., n_bins) Content of each bin in bin_edges for each point in grid_points. First dimesion has to correspond to number of grid_points, last dimension has to correspond to number of bins for the quantity that should be extrapolated (e.g. the Migra axis for EDisp) - + normalization : PDFNormalization + How the PDF is normalized Note ---- diff --git a/pyirf/interpolation/base_interpolators.py b/pyirf/interpolation/base_interpolators.py index 99ecb44af..9f47f0c31 100644 --- a/pyirf/interpolation/base_interpolators.py +++ b/pyirf/interpolation/base_interpolators.py @@ -1,19 +1,16 @@ """Base classes for interpolators""" from abc import ABCMeta, abstractmethod import enum -import astropy.units as u import numpy as np from ..binning import bin_center -from ..utils import cone_solid_angle __all__ = [ "BaseInterpolator", "ParametrizedInterpolator", "DiscretePDFInterpolator", "PDFNormalization", - "get_bin_width", ] @@ -27,16 +24,6 @@ class PDFNormalization(enum.Enum): CONE_SOLID_ANGLE = enum.auto() -def get_bin_width(bin_edges, normalization): - if normalization is PDFNormalization.AREA: - return np.diff(bin_edges) - - if normalization is PDFNormalization.CONE_SOLID_ANGLE: - return np.diff(cone_solid_angle(bin_edges).to_value(u.sr)) - - raise ValueError(f"Invalid PDF normalization: {normalization}") - - class BaseInterpolator(metaclass=ABCMeta): """ Base class for all interpolators, only knowing grid-points, diff --git a/pyirf/interpolation/moment_morph_interpolator.py b/pyirf/interpolation/moment_morph_interpolator.py index e3f8520fe..f561e050d 100644 --- a/pyirf/interpolation/moment_morph_interpolator.py +++ b/pyirf/interpolation/moment_morph_interpolator.py @@ -2,7 +2,8 @@ from pyirf.binning import bin_center, calculate_bin_indices from scipy.spatial import Delaunay -from .base_interpolators import DiscretePDFInterpolator, PDFNormalization, get_bin_width +from .base_interpolators import DiscretePDFInterpolator, PDFNormalization +from .utils import get_bin_width __all__ = [ "MomentMorphInterpolator", @@ -19,6 +20,8 @@ def _estimate_mean_std(bin_edges, binned_pdf, normalization): Array of common bin-edges for binned_pdf binned_pdf: np.ndarray, shape=(N, ..., M) PDF values from which to compute mean and std + normalization : PDFNormalization + How the PDF is normalized Returns ------- @@ -189,6 +192,8 @@ def moment_morph_estimation(bin_edges, binned_pdf, coefficients, normalization): Array of bin-entries, actual coefficients: np.ndarray, shape=(N) Estimation coefficients for each entry in binned_pdf + normalization : PDFNormalization + How the PDF is normalized Returns ------- @@ -273,6 +278,8 @@ def __init__( the quantity that should be interpolated (e.g. the Migra axis for EDisp) has to be at axis specified by axis-keyword as well as having entries corresponding to the number of bins given through bin_edges keyword. + normalization : PDFNormalization + How the PDF is normalized Note ---- diff --git a/pyirf/interpolation/nearest_simplex_extrapolator.py b/pyirf/interpolation/nearest_simplex_extrapolator.py index c108766db..3dc8c7aa8 100644 --- a/pyirf/interpolation/nearest_simplex_extrapolator.py +++ b/pyirf/interpolation/nearest_simplex_extrapolator.py @@ -7,7 +7,6 @@ import numpy as np from scipy.spatial import Delaunay -from pyirf.interpolation.base_interpolators import get_bin_width from .base_extrapolators import DiscretePDFExtrapolator, ParametrizedExtrapolator, PDFNormalization from .moment_morph_interpolator import ( @@ -15,7 +14,7 @@ linesegment_1D_interpolation_coefficients, moment_morph_estimation, ) -from .utils import find_nearest_simplex +from .utils import find_nearest_simplex, get_bin_width __all__ = [ "MomentMorphNearestSimplexExtrapolator", @@ -145,6 +144,8 @@ def __init__(self, grid_points, bin_edges, binned_pdf, normalization=PDFNormaliz the quantity that should be interpolated (e.g. the Migra axis for EDisp) has to be at the last axis as well as having entries corresponding to the number of bins given through bin_edges keyword. + normalization : PDFNormalization + How the PDF is normalized Note ---- diff --git a/pyirf/interpolation/quantile_interpolator.py b/pyirf/interpolation/quantile_interpolator.py index bdc6c6d42..890624191 100644 --- a/pyirf/interpolation/quantile_interpolator.py +++ b/pyirf/interpolation/quantile_interpolator.py @@ -1,7 +1,8 @@ import numpy as np from scipy.interpolate import griddata, interp1d -from .base_interpolators import DiscretePDFInterpolator, PDFNormalization, get_bin_width +from .base_interpolators import DiscretePDFInterpolator, PDFNormalization +from .utils import get_bin_width __all__ = ["QuantileInterpolator"] diff --git a/pyirf/interpolation/tests/test_utils.py b/pyirf/interpolation/tests/test_utils.py index 60fe5f7df..cba06d6cf 100644 --- a/pyirf/interpolation/tests/test_utils.py +++ b/pyirf/interpolation/tests/test_utils.py @@ -147,3 +147,15 @@ def test_find_nearest_simplex(non_rect_grid): assert find_nearest_simplex(non_rect_grid, np.array([-10, -10])) == 0 assert find_nearest_simplex(non_rect_grid, np.array([10, 30])) == 1 assert find_nearest_simplex(non_rect_grid, np.array([20.00000000001, -10])) == 2 + + +def test_get_bin_width(): + from pyirf.interpolation.utils import get_bin_width + from pyirf.interpolation import PDFNormalization + + bins = np.array([0, 1, 3]) + np.testing.assert_allclose(get_bin_width(bins, PDFNormalization.AREA), [1, 2]) + + bins = np.array([0, np.pi / 3, np.pi / 2]) + width = get_bin_width(bins, PDFNormalization.CONE_SOLID_ANGLE) + np.testing.assert_allclose(width, [np.pi, np.pi]) diff --git a/pyirf/interpolation/utils.py b/pyirf/interpolation/utils.py index aa54ad100..5936ec63a 100644 --- a/pyirf/interpolation/utils.py +++ b/pyirf/interpolation/utils.py @@ -1,5 +1,19 @@ +import astropy.units as u import numpy as np +from ..utils import cone_solid_angle +from .base_interpolators import PDFNormalization + + +def get_bin_width(bin_edges, normalization): + if normalization is PDFNormalization.AREA: + return np.diff(bin_edges) + + if normalization is PDFNormalization.CONE_SOLID_ANGLE: + return np.diff(cone_solid_angle(bin_edges).to_value(u.sr)) + + raise ValueError(f"Invalid PDF normalization: {normalization}") + def plumb_point_dist(line, target): """