From b252ad56692e5f920191318d89ae2078025d57bd Mon Sep 17 00:00:00 2001 From: Rune Michael Dominik Date: Mon, 4 Sep 2023 16:47:19 +0200 Subject: [PATCH 1/4] Start reworking interpolation main page --- docs/interpolation.rst | 196 ++++++++++++++++++++++++++++++++++++++--- 1 file changed, 184 insertions(+), 12 deletions(-) diff --git a/docs/interpolation.rst b/docs/interpolation.rst index 2a356c538..4d30f99c8 100644 --- a/docs/interpolation.rst +++ b/docs/interpolation.rst @@ -1,22 +1,194 @@ .. _interpolation: -Interpolation of IRFs -===================== +Interpolation and Extrapolation of IRFs +======================================= -This module contains functions to interpolate from a set of IRFs for different -conditions to a new IRF. Implementations of interpolation algorithms exist as interpolator -classes and are applied by top-level estimator classes to IRF components. -Direct usage of the interpolator classes is discuraged, as only the estimator classes -check the data for consistency. +.. currentmodule:: pyirf.interpolation -This can e.g. be used to interpolate IRFs for zenith angles of 20° and 40° -to 30°. +This module contains functions to inter- or extrapolate from a set of IRFs for different +conditions to a new IRF. Implementations of interpolation and extrapolation algorithms +exist as interpolator and extrapolator classes and are applied by top-level estimator +classes to IRF components. +Direct usage of the inter- and extrapolator classes is discouraged, as only the estimator classes +check the data for consistency. Most methods support an arbitrary number of interpolation dimensions although it -is strongly advised to limit those for resonable results. +is strongly advised to limit those for reasonable results. +The herein provided functionalities can e.g. be used to interpolate IRFs for zenith +angles of 20° and 40° to 30°. + + +IRF Component Estimator Classes +------------------------------- + +.. autosummary:: + :nosignatures: + + EffectiveAreaEstimator Estimate AEFF tables. + RadMaxEstimator Estimate RadMax tables. + EnergyDispersionEstimator Estimate 2D EDISPs. + PSFTableEstimator Estimate PSF tables. + +Each of them is either tailored to parametrized or discrete PDF components though +inheritance from the respective base class + +.. autosummary:: + :nosignatures: + + ParametrizedComponentEstimator Parametrized components + DiscretePDFComponentEstimator Discrete PDF components + + +Inter- and Extrapolation Classes +-------------------------------- + +This module provides inter- and extrapolation classes that can be +plugged into the estimator classes. +Not all of these classes support arbitrary grid-dimensions where the grid +in this context is the grid of e.g. observation parameters like zenith angle and +magnetic field inclination (this would be a 2D grid) on which template IRFs exist +and are meant to be inter- or extrapolated. + +For parametrized components (Effective Areas and Rad-Max tables) these classes are: +============================================= ================== ============ ================================================================================================== +**Name** **Component Type** **Grid-Dim** **Note** +============================================= ================== ============ ================================================================================================== +:any:`GridDataInterpolator` Interpolation Arbitrary See also :any:`scipy.interpolate.griddata`. +:any:`ParametrizedNearestSimplexExtrapolator` Extrapolation 1D or 2D Linear (1D) or baryzentric (2D) extension outside the grid's convex hull from the nearest simplex. +:any:`ParametrizedNearestNeighborSearcher` Nearest Neighbor Arbitrary Nearest neighbor finder usable instead of inter- and/or extrapolation. +============================================= ================== ============ ================================================================================================== -Reference/API +For components represented by discretized PDFs (PSF and EDISP tables) these classes are: + +============================================= ================== ============ ============================================================================== +**Name** **Component Type** **Grid-Dim** **Note** +============================================= ================== ============ ============================================================================== +:any:`QuantileInterpolator` Interpolation Arbitrary Adaption of [Hol+13]_ and [Rea99]_ to discretized PDFs. +:any:`MomentMorphInterpolator` Interpolation 1D or 2D Adaption of [Baa+15]_ to discretized PDFs. +:any:`MomentMorphNearestSimplexExtrapolator` Extrapolation 1D or 2D Extension of [Baa+15]_ beyond the grid's convex hull from the nearest simplex. +:any:`DiscretePDFNearestNeighborSearcher` Nearest Neighbor Arbitrary Nearest neighbor finder usable instead of inter- and/or extrapolation. +============================================= ================== ============ ============================================================================== + +.. [Hol+13] B. E. Hollister and A. T. Pang (2013). Interpolation of Non-Gaussian Probability Distributions for Ensemble Visualization. + https://engineering.ucsc.edu/sites/default/files/technical-reports/UCSC-SOE-13-13.pdf +.. [Rea99] A. L. Read (1999). Linear Interpolation of Histograms. + Nucl. Instrum. Methods Phys. Res. A 425, 357-360. https://doi.org/10.1016/S0168-9002(98)01347-3 +.. [Baa+15] M. Baak, S. Gadatsch, R. Harrington and W. Verkerke (2015). Interpolation between + multi-dimensional histograms using a new non-linear moment morphing method + Nucl. Instrum. Methods Phys. Res. A 771, 39-48. https://doi.org/10.1016/j.nima.2014.10.033 + + +Usage Example ------------- -.. automodapi:: pyirf.interpolation +To create a estimator class for an IRF component not yet implemented, one can simply +inherit from respective base class. +Consider an example, where one is interested in an estimator for simple Gaussians. +As this is already the scope of the ``DiscretePDFComponentEstimator`` base class and +for the sake of this demonstration, let the Gaussians come with some +units attached that need handling: + +.. code-block:: python + + import astropy.units as u + + from pyirf.interpolation import DiscretePDFComponentEstimator, MomentMorphInterpolator + + class GaussianEstimatior(DiscretePDFComponentEstimator): + @u.quantity_input(gaussians=u.m) + def __init__( + self, + grid_points, + bin_edges, + gaussians, + interpolator_cls=MomentMorphInterpolator, + interpolator_kwargs=None, + extrapolator_cls=None, + extrapolator_kwargs=None, + ): + if interpolator_kwargs is None: + interpolator_kwargs = {} + + if extrapolator_kwargs is None: + extrapolator_kwargs = {} + + self.unit = gaussians.unit + + super().__init__( + grid_points=grid_points, + bin_edges=bin_edges, + binned_pdf=gaussians.to_value(u.m), + interpolator_cls=interpolator_cls, + interpolator_kwargs=interpolator_kwargs, + extrapolator_cls=extrapolator_cls, + extrapolator_kwargs=extrapolator_kwargs, + ) + + def __call__(self, target_point): + res = super().__call__(target_point) + + # Return result with correct unit + return u.Quantity(res, u.m, copy=False).to(self.unit) + +This new estimator class can now be used just like any other estimator class already +implemented in ``pyirf.interpolation``. +While the ``extrapolator_cls`` argument can be empty when creating an instance of +``GaussianEstimator``, effectively disabling extrapolation and raising an error in +case it would be needed regardless, assume the desired extrapolation method to be +``MomentMorphNearestSimplexExtrapolator``: + +.. code-block:: python + + import numpy as np + + from scipy.stats import norm + from pyirf.interpolation import MomentMorphNearestSimplexExtrapolator + + bins = np.linspace(-10, 10, 51) + grid = np.array([[1], [2], [3]]) + + gaussians = np.array([np.diff(norm(loc=x, scale=1/x).cdf(bins))/np.diff(bins) for x in grid]) + + estimator = GaussianEstimatior( + grid_points = grid, + bin_edges = bins, + gaussians = gaussians * u.m, + interpolator_cls = MomentMorphInterpolator, + extrapolator_cls = MomentMorphNearestSimplexExtrapolator + ) + +This estimator object can now easily be used to estimate Gaussians at arbitrary target points: + +.. code-block:: python + + targets = np.array([[0.9], [1.5]]) + + results = u.Quantity([estimator(target).squeeze() for target in targets]) + + +Helper Classes +-------------- + +.. autosummary:: + :nosignatures: + + PDFNormalization + + +Base Classes +------------ + +.. autosummary:: + :nosignatures: + + BaseComponentEstimator + BaseInterpolator + ParametrizedInterpolator + DiscretePDFInterpolator + BaseExtrapolator + ParametrizedExtrapolator + DiscretePDFExtrapolator + BaseNearestNeighborSearcher + + From e8fb1e4b295849c5b4c866de6e7c87d34df150cb Mon Sep 17 00:00:00 2001 From: Rune Michael Dominik Date: Wed, 6 Sep 2023 11:48:04 +0200 Subject: [PATCH 2/4] Streamline shape indications in docstrings --- docs/interpolation.rst | 13 +++--- pyirf/interpolation/base_extrapolators.py | 16 ++++--- pyirf/interpolation/component_estimators.py | 28 +++++------ pyirf/interpolation/griddata_interpolator.py | 11 +++-- .../moment_morph_interpolator.py | 46 ++++++++++--------- .../nearest_neighbor_searcher.py | 8 ++-- .../nearest_simplex_extrapolator.py | 33 ++++++++----- pyirf/interpolation/quantile_interpolator.py | 35 +++++++------- pyirf/interpolation/utils.py | 12 ++--- 9 files changed, 112 insertions(+), 90 deletions(-) diff --git a/docs/interpolation.rst b/docs/interpolation.rst index 4d30f99c8..76407530b 100644 --- a/docs/interpolation.rst +++ b/docs/interpolation.rst @@ -52,7 +52,7 @@ and are meant to be inter- or extrapolated. For parametrized components (Effective Areas and Rad-Max tables) these classes are: ============================================= ================== ============ ================================================================================================== -**Name** **Component Type** **Grid-Dim** **Note** +**Name** **Type** **Grid-Dim** **Note** ============================================= ================== ============ ================================================================================================== :any:`GridDataInterpolator` Interpolation Arbitrary See also :any:`scipy.interpolate.griddata`. :any:`ParametrizedNearestSimplexExtrapolator` Extrapolation 1D or 2D Linear (1D) or baryzentric (2D) extension outside the grid's convex hull from the nearest simplex. @@ -62,7 +62,7 @@ For parametrized components (Effective Areas and Rad-Max tables) these classes a For components represented by discretized PDFs (PSF and EDISP tables) these classes are: ============================================= ================== ============ ============================================================================== -**Name** **Component Type** **Grid-Dim** **Note** +**Name** **Type** **Grid-Dim** **Note** ============================================= ================== ============ ============================================================================== :any:`QuantileInterpolator` Interpolation Arbitrary Adaption of [Hol+13]_ and [Rea99]_ to discretized PDFs. :any:`MomentMorphInterpolator` Interpolation 1D or 2D Adaption of [Baa+15]_ to discretized PDFs. @@ -92,8 +92,8 @@ units attached that need handling: .. code-block:: python import astropy.units as u - - from pyirf.interpolation import DiscretePDFComponentEstimator, MomentMorphInterpolator + from pyirf.interpolation import (DiscretePDFComponentEstimator, + MomentMorphInterpolator) class GaussianEstimatior(DiscretePDFComponentEstimator): @u.quantity_input(gaussians=u.m) @@ -141,9 +141,8 @@ case it would be needed regardless, assume the desired extrapolation method to b .. code-block:: python import numpy as np - - from scipy.stats import norm - from pyirf.interpolation import MomentMorphNearestSimplexExtrapolator + from pyirf.interpolation import MomentMorphNearestSimplexExtrapolator + from scipy.stats import norm bins = np.linspace(-10, 10, 51) grid = np.array([[1], [2], [3]]) diff --git a/pyirf/interpolation/base_extrapolators.py b/pyirf/interpolation/base_extrapolators.py index 28b3b38a8..9edd9933b 100644 --- a/pyirf/interpolation/base_extrapolators.py +++ b/pyirf/interpolation/base_extrapolators.py @@ -26,7 +26,7 @@ def __init__(self, grid_points): self.grid_points = grid_points if self.grid_points.ndim == 1: self.grid_points = self.grid_points.reshape(*self.grid_points.shape, 1) - self.n_points = self.grid_points.shape[0] + self.N = self.grid_points.shape[0] self.grid_dim = self.grid_points.shape[1] @abstractmethod @@ -60,7 +60,7 @@ def __init__(self, grid_points, params): 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 params: np.ndarray, shape=(n_points, ..., n_params) Corresponding parameter values at each point in grid_points. @@ -84,21 +84,23 @@ class DiscretePDFExtrapolator(BaseExtrapolator): Derived from pyirf.interpolation.BaseExtrapolator """ - def __init__(self, grid_points, bin_edges, binned_pdf, normalization=PDFNormalization.AREA): + def __init__( + self, grid_points, bin_edges, binned_pdf, normalization=PDFNormalization.AREA + ): """DiscretePDFExtrapolator 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 + normalization: PDFNormalization How the PDF is normalized Note diff --git a/pyirf/interpolation/component_estimators.py b/pyirf/interpolation/component_estimators.py index dc859c868..1d102d933 100644 --- a/pyirf/interpolation/component_estimators.py +++ b/pyirf/interpolation/component_estimators.py @@ -9,8 +9,8 @@ from .base_extrapolators import DiscretePDFExtrapolator, ParametrizedExtrapolator from .base_interpolators import ( DiscretePDFInterpolator, - PDFNormalization, ParametrizedInterpolator, + PDFNormalization, ) from .griddata_interpolator import GridDataInterpolator from .quantile_interpolator import QuantileInterpolator @@ -37,7 +37,7 @@ class BaseComponentEstimator: def __init__(self, grid_points): """ Base __init__, doing sanity checks on the grid, building a - triangulated version of the grid and intantiating inter- and extrapolator. + triangulated version of the grid and instantiating inter- and extrapolator. Parameters ---------- @@ -46,7 +46,7 @@ def __init__(self, grid_points): Raises ------ TypeError: - When grid_points is not a np.ndarray + When grid_points is not a np.ndarray of float compatible values TypeError: When grid_point has dtype object ValueError: @@ -370,6 +370,8 @@ def __init__( class EffectiveAreaEstimator(ParametrizedComponentEstimator): + """Estimator class for effective area tables (AEFF_2D).""" + @u.quantity_input(effective_area=u.m**2, min_effective_area=u.m**2) def __init__( self, @@ -382,8 +384,6 @@ def __init__( min_effective_area=1 * u.m**2, ): """ - Estimator class for effective areas. - Takes a grid of effective areas for a bunch of different parameters and inter-/extrapolates (log) effective areas to given value of those parameters. @@ -471,6 +471,8 @@ def __call__(self, target_point): class RadMaxEstimator(ParametrizedComponentEstimator): + """Estimator class for rad-max tables (RAD_MAX, RAD_MAX_2D).""" + def __init__( self, grid_points, @@ -481,8 +483,6 @@ def __init__( extrapolator_kwargs=None, ): """ - Estimator class for RAD_MAX tables. - Takes a grid of rad max values for a bunch of different parameters and inter-/extrapolates rad max values to given value of those parameters. @@ -545,6 +545,8 @@ def __call__(self, target_point): class EnergyDispersionEstimator(DiscretePDFComponentEstimator): + """Estimator class for energy dispersions (EDISP_2D).""" + def __init__( self, grid_points, @@ -557,8 +559,6 @@ def __init__( axis=-2, ): """ - Estimator class for energy dispersions. - Takes a grid of energy dispersions for a bunch of different parameters and inter-/extrapolates energy dispersions to given value of those parameters. @@ -631,6 +631,8 @@ def __call__(self, target_point): class PSFTableEstimator(DiscretePDFComponentEstimator): + """Estimator class for point spread function tables (PSF_TABLE).""" + @u.quantity_input(psf=u.sr**-1, source_offset_bins=u.deg) def __init__( self, @@ -644,8 +646,6 @@ def __init__( axis=-1, ): """ - Estimator class for point spread functions. - Takes a grid of psfs or a bunch of different parameters and inter-/extrapolates psfs to given value of those parameters. @@ -712,7 +712,7 @@ def __init__( def __call__(self, target_point): """ - Estimating energy dispersions at target_point, inter-/extrapolates as needed and + Estimating psf tables at target_point, inter-/extrapolates as needed and specified in __init__. Parameters @@ -731,4 +731,6 @@ def __call__(self, target_point): interpolated_psf_normed = super().__call__(target_point) # Undo normalisation to get a proper PSF and return - return u.Quantity(np.swapaxes(interpolated_psf_normed, -1, self.axis), u.sr**-1, copy=False) + return u.Quantity( + np.swapaxes(interpolated_psf_normed, -1, self.axis), u.sr**-1, copy=False + ) diff --git a/pyirf/interpolation/griddata_interpolator.py b/pyirf/interpolation/griddata_interpolator.py index 9130a5701..7d50359a9 100644 --- a/pyirf/interpolation/griddata_interpolator.py +++ b/pyirf/interpolation/griddata_interpolator.py @@ -9,14 +9,16 @@ class GridDataInterpolator(ParametrizedInterpolator): + """ "Wrapper arounf scipy.interpolate.griddata.""" + def __init__(self, grid_points, params, **griddata_kwargs): """Parametrized Interpolator using scipy.interpolate.griddata Parameters ---------- - grid_points: np.ndarray, shape=(N, O) + grid_points: np.ndarray, shape=(n_points, n_dims) Grid points at which interpolation templates exist - params: np.ndarray, shape=(N, ...) + params: np.ndarray, shape=(n_points, ..., n_params) Structured array of corresponding parameter values at each point in grid_points. First dimesion has to correspond to number of grid_points @@ -48,12 +50,13 @@ def interpolate(self, target_point): Parameters ---------- - target_point: np.ndarray, shape=(1, O) + target_point: np.ndarray, shape=(1, n_dims) Target point for interpolation Returns ------- - interpolant: np.ndarray + interpolant: np.ndarray, shape=(1, ..., n_params) + Interpolated parameter values References ---------- diff --git a/pyirf/interpolation/moment_morph_interpolator.py b/pyirf/interpolation/moment_morph_interpolator.py index f561e050d..77c5000b9 100644 --- a/pyirf/interpolation/moment_morph_interpolator.py +++ b/pyirf/interpolation/moment_morph_interpolator.py @@ -16,18 +16,18 @@ def _estimate_mean_std(bin_edges, binned_pdf, normalization): Parameters ---------- - bin_edges: np.ndarray, shape=(M+1) + bin_edges: np.ndarray, shape=(n_bins+1) Array of common bin-edges for binned_pdf - binned_pdf: np.ndarray, shape=(N, ..., M) + binned_pdf: np.ndarray, shape=(n_points, ..., n_bins) PDF values from which to compute mean and std normalization : PDFNormalization How the PDF is normalized Returns ------- - mean: np.ndarray, shape=(N, ...) + mean: np.ndarray, shape=(n_points, ...) Estimated mean for each input template - std: np.ndarray, shape=(N, ...) + std: np.ndarray, shape=(n_points, ...) Estimated standard deviation for each input template. Set to width/2 if only one bin in the input template is =/= 0 """ @@ -51,7 +51,7 @@ def _estimate_mean_std(bin_edges, binned_pdf, normalization): if np.any(mask): 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) + 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 @@ -63,17 +63,17 @@ def _lookup(bin_edges, binned_pdf, x): Parameters ---------- - bin_edges: np.ndarray, shape=(M+1) + bin_edges: np.ndarray, shape=(n_bins+1) Array of common bin-edges for binned_pdf - binned_pdf: np.ndarray, shape=(N, ..., M) + binned_pdf: np.ndarray, shape=(n_points, ..., n_bins) 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 + x: numpy.ndarray, shape=(n_points, ..., n_bins) + Array of n_bins points for each input template, where the histogram-value (bin-height) should be found Returns ------- - y: numpy.ndarray, shape=(N, ..., M) - Array of the bin-heights at the M points x, set to 0 at each point outside the histogram + y: numpy.ndarray, shape=(n_points, ..., n_bins) + Array of the bin-heights at the n_bins points x, set to 0 at each point outside the histogram """ # Find the bin where each point x is located in @@ -186,18 +186,18 @@ def moment_morph_estimation(bin_edges, binned_pdf, coefficients, normalization): Parameters ---------- - bin_edges: np.ndarray, shape=(M+1) + bin_edges: np.ndarray, shape=(n_bins+1) Array of common bin-edges for binned_pdf - binned_pdf: np.ndarray, shape=(N, ..., M) + binned_pdf: np.ndarray, shape=(n_points, ..., n_bins) Array of bin-entries, actual - coefficients: np.ndarray, shape=(N) + coefficients: np.ndarray, shape=(n_points) Estimation coefficients for each entry in binned_pdf normalization : PDFNormalization How the PDF is normalized Returns ------- - f_new: numpy.ndarray, shape=(1, M) + f_new: numpy.ndarray, shape=(1, n_bins) Interpolated histogram References @@ -216,7 +216,9 @@ def moment_morph_estimation(bin_edges, binned_pdf, coefficients, normalization): # 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, normalization=normalization) + 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") ) @@ -259,6 +261,8 @@ def moment_morph_estimation(bin_edges, binned_pdf, coefficients, normalization): class MomentMorphInterpolator(DiscretePDFInterpolator): + """Interpolator class providing Moment Morphing to interpolate discretized PDFs.""" + def __init__( self, grid_points, bin_edges, binned_pdf, normalization=PDFNormalization.AREA ): @@ -267,11 +271,11 @@ def __init__( Parameters ---------- - grid_points: np.ndarray, shape=(N, ...) + grid_points: np.ndarray, shape=(n_points, n_dims) Grid points at which interpolation templates exist. May be one ot two dimensional. - bin_edges: np.ndarray, shape=(M+1) + bin_edges: np.ndarray, shape=(n_bins+1) Edges of the data binning - binned_pdf: np.ndarray, shape=(N, ..., M) + 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. Interpolation dimension, meaning the @@ -342,12 +346,12 @@ def interpolate(self, target_point): Parameters ---------- - target_point: numpy.ndarray + target_point: numpy.ndarray, shape=(1, n_dims) Value for which the interpolation is performed (target point) Returns ------- - f_new: numpy.ndarray, shape=(1,...,M,...) + f_new: numpy.ndarray, shape=(1,...,n_bins) Interpolated and binned pdf References diff --git a/pyirf/interpolation/nearest_neighbor_searcher.py b/pyirf/interpolation/nearest_neighbor_searcher.py index 6ddcf5df7..4543d4c07 100644 --- a/pyirf/interpolation/nearest_neighbor_searcher.py +++ b/pyirf/interpolation/nearest_neighbor_searcher.py @@ -71,12 +71,12 @@ def interpolate(self, target_point): Parameters ---------- - target_point: numpy.ndarray + target_point: numpy.ndarray, shape=(1, n_dims) Value for which the nearest neighbor should be found (target point) Returns ------- - content_new: numpy.ndarray, shape=(1,...,M,...) + content_new: numpy.ndarray, shape=(1, ...) values at nearest neighbor Note @@ -100,7 +100,7 @@ def interpolate(self, target_point): class DiscretePDFNearestNeighborSearcher(BaseNearestNeighborSearcher): """ Dummy NearestNeighbor approach usable instead of - actual Interpolation/Extrapolation. + actual interpolation/extrapolation. Compatible with discretized PDF IRF component API. """ @@ -139,7 +139,7 @@ def __init__(self, grid_points, bin_edges, binned_pdf, norm_ord=2): class ParametrizedNearestNeighborSearcher(BaseNearestNeighborSearcher): """ Dummy NearestNeighbor approach usable instead of - actual Interpolation/Extrapolation + actual interpolation/extrapolation Compatible with parametrized IRF component API. """ diff --git a/pyirf/interpolation/nearest_simplex_extrapolator.py b/pyirf/interpolation/nearest_simplex_extrapolator.py index 3dc8c7aa8..a04481be0 100644 --- a/pyirf/interpolation/nearest_simplex_extrapolator.py +++ b/pyirf/interpolation/nearest_simplex_extrapolator.py @@ -7,8 +7,11 @@ import numpy as np from scipy.spatial import Delaunay - -from .base_extrapolators import DiscretePDFExtrapolator, ParametrizedExtrapolator, PDFNormalization +from .base_extrapolators import ( + DiscretePDFExtrapolator, + ParametrizedExtrapolator, + PDFNormalization, +) from .moment_morph_interpolator import ( barycentric_2D_interpolation_coefficients, linesegment_1D_interpolation_coefficients, @@ -23,17 +26,19 @@ class ParametrizedNearestSimplexExtrapolator(ParametrizedExtrapolator): + """Extrapolator class extending linear or baryzentric interpolation outside a grid's convex hull.""" + def __init__(self, grid_points, params): """ - Extrapolator class using linear extrapolation in one ore two + Extrapolator class using linear or baryzentric extrapolation in one ore two grid-dimensions. Parameters ---------- - grid_points: np.ndarray, shape=(N, ...) + grid_points: np.ndarray, shape=(n_points, n_dims) Grid points at which templates exist. May be one ot two dimensional. Have to be sorted in accending order for 1D. - params: np.ndarray, shape=(N, ...) + params: np.ndarray, shape=(n_points, ...) Array of corresponding parameter values at each point in grid_points. First dimesion has to correspond to number of grid_points @@ -95,7 +100,7 @@ def extrapolate(self, target_point): Parameters ---------- - target_point: numpy.ndarray + target_point: numpy.ndarray, shape=(1, n_dims) Value for which the extrapolation is performed (target point) Returns @@ -123,7 +128,11 @@ def extrapolate(self, target_point): class MomentMorphNearestSimplexExtrapolator(DiscretePDFExtrapolator): - def __init__(self, grid_points, bin_edges, binned_pdf, normalization=PDFNormalization.AREA): + """Extrapolator class extending moment morphing interpolation outside a grid's convex hull.""" + + 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 @@ -132,12 +141,12 @@ def __init__(self, grid_points, bin_edges, binned_pdf, normalization=PDFNormaliz Parameters ---------- - grid_points: np.ndarray, shape=(N, ...) + grid_points: np.ndarray, shape=(n_points, n_dims) Grid points at which templates exist. May be one ot two dimensional. Have to be sorted in accending order for 1D. - bin_edges: np.ndarray, shape=(M+1) + bin_edges: np.ndarray, shape=(n_bins+1) Edges of the data binning - binned_pdf: np.ndarray, shape=(N, ..., M) + 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. Extrapolation dimension, meaning the @@ -201,12 +210,12 @@ def extrapolate(self, target_point): Parameters ---------- - target_point: numpy.ndarray + target_point: numpy.ndarray, shape=(1, n_dims) Value for which the extrapolation is performed (target point) Returns ------- - values: numpy.ndarray, shape=(1, ..., M) + values: numpy.ndarray, shape=(1, ..., n_bins) Extrapolated discretized PDFs """ diff --git a/pyirf/interpolation/quantile_interpolator.py b/pyirf/interpolation/quantile_interpolator.py index 890624191..deaeae184 100644 --- a/pyirf/interpolation/quantile_interpolator.py +++ b/pyirf/interpolation/quantile_interpolator.py @@ -27,18 +27,18 @@ def ppf_values(bin_mids, cdfs, quantiles): Parameters ---------- - bin_mids: numpy.ndarray, shape=(M) + bin_mids: numpy.ndarray, shape=(n_bins) Bin-mids for each bin along interpolation axis - cdfs: numpy.ndarray, shape=(N,...,M) + cdfs: numpy.ndarray, shape=(n_points,...,n_bins) Corresponding cdf-values for all quantiles - quantiles: numpy.ndarray, shape=(L) + quantiles: numpy.ndarray, shape=(n_quantiles) Quantiles for which the ppf-values should be estimated Returns ------- - ppfs: numpy.ndarray, shape=(1,...,L) + ppfs: numpy.ndarray, shape=(1,...,n_quantiles) Corresponding ppf-values for all quantiles at the target interpolation point """ @@ -93,20 +93,20 @@ def pdf_from_ppf(bin_edges, interp_ppfs, quantiles): Parameters ---------- - bin_edges: numpy.ndarray, shape=(M+1) + bin_edges: numpy.ndarray, shape=(n_bins+1) Edges of the bins in which the final pdf should be binned - interp_ppfs: numpy.ndarray, shape=(1,...,L) + interp_ppfs: numpy.ndarray, shape=(1,...,n_quantiles) 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. - quantiles: numpy.ndarray, shape=(L) + quantiles: numpy.ndarray, shape=(n_quantiles) Quantiles corresponding to the ppf-values in interp_ppfs Returns ------- - pdf_values: numpy.ndarray, shape=(1,...,M) + pdf_values: numpy.ndarray, shape=(1,...,n_bins) Recomputed, binned pdf at target_point """ # recalculate pdf values through numerical differentiation @@ -152,6 +152,8 @@ def norm_pdf(pdf_values, bin_edges, normalization): class QuantileInterpolator(DiscretePDFInterpolator): + """Interpolator class providing quantile interpoalation.""" + def __init__( self, grid_points, @@ -160,15 +162,14 @@ def __init__( quantile_resolution=1e-3, normalization=PDFNormalization.AREA, ): - """BinnedInterpolator constructor - + """ Parameters ---------- - grid_points : np.ndarray + grid_points : np.ndarray, shape=(n_points, n_dims) Grid points at which interpolation templates exist - bin_edges : np.ndarray + bin_edges : np.ndarray, shape=(n_bins+1) Edges of the data binning - binned_pdf : np.ndarray + 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, the last axis has to correspond to number @@ -232,12 +233,12 @@ def interpolate(self, target_point): Parameters ---------- - target_point: numpy.ndarray, shape=(O) + target_point: numpy.ndarray, shape=(1, n_dims) Value for which the interpolation is performed (target point) Returns ------- - f_new: numpy.ndarray, shape=(1,...,M,...) + f_new: numpy.ndarray, shape=(1, ..., n_bins) Interpolated and binned pdf References @@ -257,7 +258,9 @@ def interpolate(self, target_point): )[..., 1:] # Renormalize pdf to sum of 1 - normed_interpolated_pdfs = norm_pdf(interpolated_pdfs, self.bin_edges[1:], self.normalization) + 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/utils.py b/pyirf/interpolation/utils.py index 5936ec63a..5792a60c8 100644 --- a/pyirf/interpolation/utils.py +++ b/pyirf/interpolation/utils.py @@ -22,9 +22,9 @@ def plumb_point_dist(line, target): Parameters ---------- - line: np.ndarray, shape=(2, M) - Array of two points spanning a line segment. Might be in two or three dims M. - target: np.ndarray, shape=(M) + line: np.ndarray, shape=(2, n_dims) + Array of two points spanning a line segment. Might be in two or three dims n_dims. + target: np.ndarray, shape=(n_dims) Target point, of which the minimal distance to line segement is needed Returns @@ -75,9 +75,9 @@ def point_facet_angle(line, target): Parameters ---------- - line: np.ndarray, shape=(2, M) - Array of two points spanning a line segment. Might be in two or three dims M. - target: np.ndarray, shape=(M) + line: np.ndarray, shape=(2, n_dims) + Array of two points spanning a line segment. Might be in two or three dims n_dims. + target: np.ndarray, shape=(n_dims) Target point, of which the angle is needed. Returns From 669d3a50dd399bc5bd6c54942a147e53596bb828 Mon Sep 17 00:00:00 2001 From: Rune Michael Dominik Date: Wed, 6 Sep 2023 14:23:47 +0200 Subject: [PATCH 3/4] Have full pyirf.interpolation API be created --- docs/interpolation.rst | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/docs/interpolation.rst b/docs/interpolation.rst index 76407530b..f869391b6 100644 --- a/docs/interpolation.rst +++ b/docs/interpolation.rst @@ -191,3 +191,13 @@ Base Classes BaseNearestNeighborSearcher +Full API +-------- + +.. automodapi:: pyirf.interpolation + :no-heading: + :no-main-docstr: + :inherited-members: + :no-inheritance-diagram: + + From 59b94087f1fd324bceabec5e840ee3600936f182 Mon Sep 17 00:00:00 2001 From: Rune Michael Dominik Date: Thu, 7 Sep 2023 14:59:49 +0200 Subject: [PATCH 4/4] Add basic usage code example --- docs/interpolation.rst | 84 +++++++++++++++++++++++++++++++++++------- 1 file changed, 70 insertions(+), 14 deletions(-) diff --git a/docs/interpolation.rst b/docs/interpolation.rst index f869391b6..6bed34075 100644 --- a/docs/interpolation.rst +++ b/docs/interpolation.rst @@ -14,8 +14,8 @@ check the data for consistency. Most methods support an arbitrary number of interpolation dimensions although it is strongly advised to limit those for reasonable results. -The herein provided functionalities can e.g. be used to interpolate IRFs for zenith -angles of 20° and 40° to 30°. +The herein provided functionalities can e.g. be used to interpolate the IRF +for a zenith angle of 30° from available IRFs at 20° and 40°. IRF Component Estimator Classes @@ -29,14 +29,6 @@ IRF Component Estimator Classes EnergyDispersionEstimator Estimate 2D EDISPs. PSFTableEstimator Estimate PSF tables. -Each of them is either tailored to parametrized or discrete PDF components though -inheritance from the respective base class - -.. autosummary:: - :nosignatures: - - ParametrizedComponentEstimator Parametrized components - DiscretePDFComponentEstimator Discrete PDF components Inter- and Extrapolation Classes @@ -79,11 +71,75 @@ For components represented by discretized PDFs (PSF and EDISP tables) these clas Nucl. Instrum. Methods Phys. Res. A 771, 39-48. https://doi.org/10.1016/j.nima.2014.10.033 -Usage Example -------------- +Using Estimator Classes +----------------------- + +Usage of the estimator classes is simple. +As an example, consider CTA's Prod5 IRFs [CTA+21]_, they can be downloaded manually or by executing +``download_irfs.py`` in ``pyirf's`` root directory, which downloads them to ``.../pyirf/irfs/``. +The estimator classes can simply be used by first creating an instance of the respective class with all +relevant information and then using the object's ``__call__`` interface the obtain results for a specific +target point. +As the energy dispersion represents one of the discretized PDF IRF components, one can use the +``MomentMorphInterpolator`` for interpolation and the ``DiscretePDFNearestNeighborSearcher`` +for extrapolation. + +.. code-block:: python + + import numpy as np + + from gammapy.irf import load_irf_dict_from_file + from glob import glob + from pyirf.interpolation import ( + EnergyDispersionEstimator, + MomentMorphInterpolator, + DiscretePDFNearestNeighborSearcher + ) + + # Load IRF data, replace path with actual path + PROD5_IRF_PATH = "pyirf/irfs/*.fits.gz" + + irfs = [load_irf_dict_from_file(path) for path in sorted(glob(PROD5_IRF_PATH))] + + edisps = np.array([irf["edisp"].quantity for irf in irfs]) + bin_edges = irfs[0]["edisp"].axes["migra"].edges + # IRFs are for zenith distances of 20, 40 and 60 deg + zen_pnt = np.array([[20], [40], [60]]) + + # Create estimator instance + edisp_estimator = EnergyDispersionEstimator( + grid_points=zen_pnt, + migra_bins=bin_edges, + energy_dispersion=edisps, + interpolator_cls=MomentMorphInterpolator, + interpolator_kwargs=None, + extrapolator_cls=DiscretePDFNearestNeighborSearcher, + extrapolator_kwargs=None, + ) + + # Estimate energy dispersions + interpolated_edisp = edisp_estimator(np.array([[30]])) + extrapolated_edisp = edisp_estimator(np.array([[10]])) + + +.. [CTA+21] Cherenkov Telescope Array Observatory & Cherenkov Telescope Array Consortium. (2021). + CTAO Instrument Response Functions - prod5 version v0.1 (v0.1) [Data set]. Zenodo. + https://doi.org/10.5281/zenodo.5499840 + + +Creating new Estimator Classes +------------------------------ To create a estimator class for an IRF component not yet implemented, one can simply inherit from respective base class. +There are two, tailored to either parametrized or discrete PDF components. + +.. autosummary:: + :nosignatures: + + ParametrizedComponentEstimator Parametrized components + DiscretePDFComponentEstimator Discrete PDF components + Consider an example, where one is interested in an estimator for simple Gaussians. As this is already the scope of the ``DiscretePDFComponentEstimator`` base class and for the sake of this demonstration, let the Gaussians come with some @@ -182,6 +238,8 @@ Base Classes :nosignatures: BaseComponentEstimator + ParametrizedComponentEstimator + DiscretePDFComponentEstimator BaseInterpolator ParametrizedInterpolator DiscretePDFInterpolator @@ -199,5 +257,3 @@ Full API :no-main-docstr: :inherited-members: :no-inheritance-diagram: - -