diff --git a/pyirf/interpolation/base_extrapolators.py b/pyirf/interpolation/base_extrapolators.py index 79ff7a79c..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, bin_contents): + def __init__(self, grid_points, bin_edges, binned_pdf, normalization=PDFNormalization.AREA): """DiscretePDFExtrapolator Parameters @@ -92,7 +93,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 @@ -105,6 +106,7 @@ def __init__(self, grid_points, bin_edges, bin_contents): """ super().__init__(grid_points) + self.normalization = normalization 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..99ecb44af 100644 --- a/pyirf/interpolation/base_interpolators.py +++ b/pyirf/interpolation/base_interpolators.py @@ -1,10 +1,40 @@ """Base classes for interpolators""" from abc import ABCMeta, abstractmethod +import enum +import astropy.units as u 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): + """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 radian. + 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): @@ -84,21 +114,24 @@ 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, 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 - 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 the quantity that should be interpolated (e.g. the Migra axis for EDisp) - + normalization : PDFNormalization + How the PDF is normalized Note ---- @@ -108,4 +141,5 @@ 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 + self.normalization = normalization diff --git a/pyirf/interpolation/component_estimators.py b/pyirf/interpolation/component_estimators.py index c34f0a1f7..dc859c868 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, ParametrizedInterpolator +from .base_interpolators import ( + DiscretePDFInterpolator, + PDFNormalization, + ParametrizedInterpolator, +) from .griddata_interpolator import GridDataInterpolator from .quantile_interpolator import QuantileInterpolator @@ -152,7 +156,7 @@ def __init__( self, grid_points, bin_edges, - bin_contents, + binned_pdf, interpolator_cls=QuantileInterpolator, interpolator_kwargs=None, extrapolator_cls=None, @@ -168,7 +172,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 +195,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,28 +216,28 @@ 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 if self.grid_dim == 1: 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 +251,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 +262,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 ) @@ -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 @@ -596,7 +602,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, @@ -681,14 +687,23 @@ 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 = {} + + 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, - bin_contents=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, @@ -707,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) @@ -716,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 0d69bc802..e3f8520fe 100644 --- a/pyirf/interpolation/moment_morph_interpolator.py +++ b/pyirf/interpolation/moment_morph_interpolator.py @@ -2,23 +2,23 @@ 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, get_bin_width __all__ = [ "MomentMorphInterpolator", ] -def _estimate_mean_std(bin_edges, bin_contents): +def _estimate_mean_std(bin_edges, binned_pdf, normalization): """ 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,40 @@ 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 = get_bin_width(bin_edges, normalization) + + # 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.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 -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 +84,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 +177,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, normalization): """ 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 +204,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, normalization=normalization) 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 +226,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 +234,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 @@ -240,15 +245,18 @@ def moment_morph_estimation(bin_edges, bin_contents, coefficients): # 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, *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, normalization=PDFNormalization.AREA + ): """ Interpolator class using moment morphing. @@ -258,7 +266,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 +278,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, normalization) if self.grid_dim == 2: self.triangulation = Delaunay(self.grid_points) @@ -293,8 +301,9 @@ 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, + normalization=self.normalization, ) def _interpolate2D(self, target_point): @@ -312,8 +321,9 @@ 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, + normalization=self.normalization, ) def interpolate(self, target_point): diff --git a/pyirf/interpolation/nearest_neighbor_searcher.py b/pyirf/interpolation/nearest_neighbor_searcher.py index 8bcaaa108..6ddcf5df7 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 @@ -131,9 +129,7 @@ def __init__(self, grid_points, bin_edges, bin_contents, norm_ord=2): Also calls pyirf.interpolation.BaseNearestNeighborSearcher.__init__ """ - super().__init__( - grid_points=grid_points, contents=bin_contents, norm_ord=norm_ord - ) + super().__init__(grid_points=grid_points, values=binned_pdf, norm_ord=norm_ord) DiscretePDFInterpolator.register(DiscretePDFNearestNeighborSearcher) @@ -168,7 +164,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..c108766db 100644 --- a/pyirf/interpolation/nearest_simplex_extrapolator.py +++ b/pyirf/interpolation/nearest_simplex_extrapolator.py @@ -7,7 +7,9 @@ import numpy as np from scipy.spatial import Delaunay -from .base_extrapolators import DiscretePDFExtrapolator, ParametrizedExtrapolator +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, linesegment_1D_interpolation_coefficients, @@ -122,7 +124,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, normalization=PDFNormalization.AREA): """ Extrapolator class extending/reusing parts of Moment Morphing by allowing for negative extrapolation coefficients computed @@ -136,7 +138,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 +150,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, normalization) if self.grid_dim == 2: self.triangulation = Delaunay(self.grid_points) @@ -169,8 +171,9 @@ 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, + normalization=self.normalization, ) def _extrapolate2D(self, simplex_inds, target_point): @@ -185,8 +188,9 @@ 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, + normalization=self.normalization, ) def extrapolate(self, target_point): @@ -228,9 +232,10 @@ 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.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..bdc6c6d42 100644 --- a/pyirf/interpolation/quantile_interpolator.py +++ b/pyirf/interpolation/quantile_interpolator.py @@ -1,16 +1,17 @@ import numpy as np from scipy.interpolate import griddata, interp1d -from .base_interpolators import DiscretePDFInterpolator +from .base_interpolators import DiscretePDFInterpolator, PDFNormalization, get_bin_width __all__ = ["QuantileInterpolator"] -def cdf_values(bin_contents): +def cdf_values(binned_pdf, bin_edges, normalization): """ compute cdf values and assure they are normed to unity """ - cdfs = np.cumsum(bin_contents, axis=-1) + 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 with np.errstate(invalid="ignore"): @@ -95,7 +96,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. @@ -122,7 +123,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 @@ -131,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, ) @@ -149,28 +151,37 @@ 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, + 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 - 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 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 ------ 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 +189,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 +198,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 +208,14 @@ 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, + normalization=normalization, ) # Compute CDF values - self.cdfs = cdf_values(self.bin_contents) + 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) @@ -242,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/__init__.py b/pyirf/interpolation/tests/__init__.py new file mode 100644 index 000000000..e69de29bb 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..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(): @@ -239,9 +242,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 +252,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 +260,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 +277,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 +291,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 +308,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,12 +316,13 @@ 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, ) 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(): @@ -331,11 +335,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 +347,9 @@ 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, :]) + 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, @@ -352,7 +357,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}, @@ -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, @@ -406,7 +412,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 +421,24 @@ 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 + 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 + 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(): @@ -460,7 +468,9 @@ def extrapolate(self, target_point): 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])), params[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])), params[0, :]) + # Nearest neighbor is grid_point 1 at the index 1 of the original binned_pdf + 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 + 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 605f4310a..77915ccb5 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) @@ -11,13 +13,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 @@ -29,14 +31,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) + truth = binned_normal_pdf(target, 0, bins) return { "grid": grid, "target": target, - "bin_contents": bin_contents, + "binned_pdf": binned_pdf, "truth": truth, } @@ -45,14 +47,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) + truth = binned_normal_pdf(*target, bins) return { "grid": grid, "target": target, - "bin_contents": bin_contents, + "binned_pdf": binned_pdf, "truth": truth, } @@ -61,13 +63,13 @@ 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)], + [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 @@ -94,18 +96,18 @@ def test_estimate_mean_std(bins): ] ).squeeze() - mean, std = _estimate_mean_std(bins, bin_contents) + mean, std = _estimate_mean_std(bins, binned_pdf, PDFNormalization.AREA) # 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(): 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 +128,7 @@ def test_lookup(): ] ) - assert np.allclose(_lookup(bins, bin_contents, x), truth) + np.testing.assert_allclose(_lookup(bins, binned_pdf, x), truth) def test_linesegment_1D_interpolation_coefficients(): @@ -137,16 +139,16 @@ 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]), + 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]), + np.array([[0, 1]]), ) target_point = np.array([[25]]) @@ -156,11 +158,11 @@ 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) - 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 +175,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 +205,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): @@ -212,16 +214,16 @@ 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, 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 - 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): @@ -230,34 +232,34 @@ 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, 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 - 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): 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) - 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 - 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(): @@ -265,7 +267,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,41 +275,42 @@ 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]])) + 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): 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) - assert np.allclose(res, 0) + np.testing.assert_allclose(res, 0) 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) - assert np.allclose(res, 0) + np.testing.assert_allclose(res, 0) def test_MomentMorphInterpolator1D_mixed_data(bins): @@ -317,34 +320,31 @@ 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)], - [ - 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 ] ) # 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( [ [ - 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), ], ] ) @@ -354,17 +354,17 @@ 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) - expected_norms = np.array([[0, 1], [1, 0]]) - assert np.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, *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) + np.testing.assert_allclose(res.squeeze(), truth, atol=1e-3, rtol=1e-5) def test_MomentMorphInterpolator1D_extended_grid_extradims(bins): @@ -372,13 +372,13 @@ 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)], + [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 @@ -386,78 +386,78 @@ 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( [ [ - 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), ], ] ) res = interp(target) - assert np.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, *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) + np.testing.assert_allclose(res.squeeze(), truth, atol=1e-3, rtol=1e-5) 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) - 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 - 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): 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) - assert np.allclose(res, 0) + np.testing.assert_allclose(res, 0) 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) - assert np.allclose(res, 0) + np.testing.assert_allclose(res, 0) def test_MomentMorphInterpolator2D_mixed(bins): @@ -467,34 +467,31 @@ 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)], - [ - 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 ] ) # 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( [ [ - 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), ], ] ) @@ -504,17 +501,17 @@ 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) - expected_norms = np.array([[0, 1], [1, 0]]) - assert np.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, *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) + np.testing.assert_allclose(res.squeeze(), truth, atol=1e-3, rtol=1e-5) def test_MomentMorphInterpolator1D_extended_grid(bins): @@ -522,22 +519,22 @@ 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) - truth = binned_normal_pdf([target, 0], bins) + 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 - 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): @@ -545,68 +542,70 @@ 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(a, b, bins) for a, b in grid]) interp = MomentMorphInterpolator( grid_points=grid, bin_edges=bins, - bin_contents=bin_contents, + binned_pdf=binned_pdf, ) res = interp(target) - truth = binned_normal_pdf(target, bins) + 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 - 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): from pyirf.interpolation import MomentMorphInterpolator grid = np.array([[20, 20], [40, 20], [30, 40], [50, 20], [45, 40]]) - target = np.array([25, 25]) - bin_contents = np.array( + 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 ] ) interp = MomentMorphInterpolator( grid_points=grid, bin_edges=bins, - bin_contents=bin_contents, + binned_pdf=binned_pdf, ) 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), ], ] ) res = interp(target) - assert np.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, *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) + np.testing.assert_allclose(res.squeeze(), truth, atol=1e-3, rtol=1e-5) def test_MomentMorphInterpolator3D(): @@ -616,7 +615,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 +624,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..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 contents(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)]] @@ -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..313cdc727 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 @@ -45,17 +45,17 @@ 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, :] - 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]]) interpolant2 = interpolator(target_point2) - dummy_data_target2 = -2.5 * slope + 1 + dummy_data_target2 = (-2.5 * slope + 1)[np.newaxis, :] - 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]]) @@ -104,12 +104,14 @@ 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() + ) - assert np.allclose(interpolant2, dummy_data_target2) + np.testing.assert_allclose(interpolant2, dummy_data_target2) assert interpolant2.shape == (1, *dummy_data.shape[1:]) @@ -137,32 +139,34 @@ 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 - bin_contents = np.array( + 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), - np.zeros(len(bins) - 1), + binned_normal_pdf(x + 10, 0, bins), + 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] - bin_contents[0, 2, 0, 20] = 1 - bin_contents[1, 2, 0, 30] = 1 - bin_contents[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. - bin_contents[0, 0, 0, :] = np.zeros(len(bins) - 1) + binned_pdf[0, 0, 0, :] = 0.0 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]) @@ -171,84 +175,86 @@ def test_MomentMorphNearestSimplexExtrapolator_1DGrid(bins): truth1 = np.array( [ [ - np.zeros(len(bins) - 1), - binned_normal_pdf([target1 + 1, 0], bins), + np.zeros(n_bins), + binned_normal_pdf(target1 + 1, 0, bins), ], [ - binned_normal_pdf([target1 + 10, 0], bins), - np.zeros(len(bins) - 1), + binned_normal_pdf(target1 + 10, 0, bins), + 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]]) - assert np.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, *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) + np.testing.assert_allclose(res1.squeeze(), truth1, atol=4e-3, rtol=1e-5) target2 = np.array([45]) 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), - np.zeros(len(bins) - 1), + binned_normal_pdf(target2 + 10, 0, bins), + 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]]) - assert np.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, *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) + np.testing.assert_allclose(res2.squeeze(), truth2, atol=1e-3, rtol=1e-5) 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 - bin_contents = np.array( + 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(n_bins)], ] - for x in grid_points + for a, b in grid_points ] ) # 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(n_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 ) target1 = np.array([25, 10]) @@ -257,48 +263,48 @@ def test_MomentMorphNearestSimplexExtrapolator_2DGrid(bins): truth1 = np.array( [ [ - np.zeros(len(bins) - 1), - binned_normal_pdf([target1[0] + 1, target1[1]], bins), + 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), + binned_normal_pdf(target1[0] + 10, target1[1], bins), + np.zeros(n_bins), ], ] ) res1 = extrap(target1) - expected_norms1 = np.array([[0, 1], [1, 0]]) - assert np.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, *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) + np.testing.assert_allclose(res1.squeeze(), truth1, atol=1e-3, rtol=1e-5) target2 = np.array([45, 20]) 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), - np.zeros(len(bins) - 1), + binned_normal_pdf(target2[0] + 10, target2[1], bins), + np.zeros(n_bins), ], ] ) res2 = extrap(target2) - expected_norms2 = np.array([[1, 1], [1, 0]]) - assert np.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, *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) + np.testing.assert_allclose(res2.squeeze(), truth2, atol=1e-3, rtol=1e-5) def test_MomentMorphNearestSimplexExtrapolator_3DGrid(): @@ -315,7 +321,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,15 +336,15 @@ 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)] + [binned_normal_pdf(x, 0, bins), binned_normal_pdf(x + 10, 0, bins)] for x in grid_points ] ) 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( @@ -347,5 +353,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) diff --git a/pyirf/interpolation/tests/test_quantile_interpolator.py b/pyirf/interpolation/tests/test_quantile_interpolator.py index 195d95fdc..2a1892cd0 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,19 @@ 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 +52,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 +67,13 @@ 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 +88,38 @@ 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) + 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( + result, + data["binned_pdfs"][0], + ) + + result = norm_pdf( + np.zeros(len(data["bin_edges"]) - 1), + data["bin_edges"], + PDFNormalization.AREA, + ) + np.testing.assert_allclose(result, 0) def test_interpolate_binned_pdf(data): @@ -96,7 +128,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, ) 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", +]