Extend and rework IRF interpolation docs

docs/interpolation.rst

@@ -1,22 +1,193 @@
 .. _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:
 
-Reference/API
+=============================================    ==================  ============    ==================================================================================================
+**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.
+:any:`ParametrizedNearestNeighborSearcher`       Nearest Neighbor    Arbitrary       Nearest neighbor finder usable instead of inter- and/or extrapolation. 
+=============================================    ==================  ============    ==================================================================================================
+
+For components represented by discretized PDFs (PSF and EDISP tables) these classes are:
+
+=============================================    ==================  ============    ==============================================================================
+**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.
+: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 pyirf.interpolation import MomentMorphNearestSimplexExtrapolator
+   from scipy.stats import norm 
+
+   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
+
+

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

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
+        )