Repository: gammapy/gammapy

Issue:
IRF.evaluate fails for IRFs with single energy bin
**Gammapy version**

Current master

**Bug description**

`IRF.evaluate` fails in the interpolation code when there is only a single energy bin, as for example is used with the global RAD_MAX:



**Expected behavior**

Returning the single global value.

**To Reproduce**

```
In [1]: from gammapy.irf import RadMax2D
In [2]: from gammapy.maps import MapAxis
In [3]: energy_axis = MapAxis.from_energy_bounds(0.01, 100, 1, unit="TeV")
In [4]: offset_axis = MapAxis.from_bounds(0., 5, 1, unit="deg", name="offset", )
In [5]: import astropy.units as u
In [6]: rad_max = RadMax2D(data=[[0.1]] * u.deg, axes=[energy_axis, offset_axis])
In [7]: rad_max.evaluate(1 * u.TeV)
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Input In [7], in <module>
----> 1 rad_max.evaluate(1 * u.TeV)

File ~/Projects/gammapy/gammapy/irf/core.py:239, in IRF.evaluate(self, method, **kwargs)
    237     if coord is not None:
    238         coords_default[key] = u.Quantity(coord, copy=False)
--> 239 data = self._interpolate(coords_default.values(), method=method)
    241 if self.interp_kwargs["fill_value"] is not None:
    242     idxs = self.axes.coord_to_idx(coords_default, clip=False)

File ~/Projects/gammapy/gammapy/utils/interpolation.py:106, in ScaledRegularGridInterpolator.__call__(self, points, method, clip, **kwargs)
    104 if self.axis is None:
    105     points = np.broadcast_arrays(*points)
--> 106     points_interp = np.stack([_.flat for _ in points]).T
    107     values = self._interpolate(points_interp, method, **kwargs)
    108     values = self.scale.inverse(values.reshape(points[0].shape))

File <__array_function__ internals>:180, in stack(*args, **kwargs)

File ~/.local/anaconda/envs/gammapy-dev/lib/python3.8/site-packages/numpy/core/shape_base.py:422, in stack(arrays, axis, out)
    420 arrays = [asanyarray(arr) for arr in arrays]
    421 if not arrays:
--> 422     raise ValueError('need at least one array to stack')
    424 shapes = {arr.shape for arr in arrays}
    425 if len(shapes) != 1:

ValueError: need at least one array to stack
```

Discussion:
codecov[bot]: # [Codecov](https://codecov.io/gh/gammapy/gammapy/pull/3837?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gammapy) Report
> Merging [#3837](https://codecov.io/gh/gammapy/gammapy/pull/3837?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gammapy) (f0fc779) into [master](https://codecov.io/gh/gammapy/gammapy/commit/dd35f7a27c2f939fa6af3080d9c0ec0d1029260e?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gammapy) (dd35f7a) will **increase** coverage by `0.00%`.
> The diff coverage is `100.00%`.

[![Impacted file tree graph](https://codecov.io/gh/gammapy/gammapy/pull/3837/graphs/tree.svg?width=650&height=150&src=pr&token=OUjhqZUZzq&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gammapy)](https://codecov.io/gh/gammapy/gammapy/pull/3837?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gammapy)

```diff
@@           Coverage Diff           @@
##           master    #3837   +/-   ##
=======================================
  Coverage   93.82%   93.83%           
=======================================
  Files         162      162           
  Lines       19834    19843    +9     
=======================================
+ Hits        18610    18619    +9     
  Misses       1224     1224           
```


| [Impacted Files](https://codecov.io/gh/gammapy/gammapy/pull/3837?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gammapy) | Coverage Δ | |
|---|---|---|
| [gammapy/irf/core.py](https://codecov.io/gh/gammapy/gammapy/pull/3837/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gammapy#diff-Z2FtbWFweS9pcmYvY29yZS5weQ==) | `89.38% <ø> (ø)` | |
| [gammapy/utils/interpolation.py](https://codecov.io/gh/gammapy/gammapy/pull/3837/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gammapy#diff-Z2FtbWFweS91dGlscy9pbnRlcnBvbGF0aW9uLnB5) | `94.94% <100.00%> (+0.50%)` | :arrow_up: |

------

[Continue to review full report at Codecov](https://codecov.io/gh/gammapy/gammapy/pull/3837?src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gammapy).
> **Legend** - [Click here to learn more](https://docs.codecov.io/docs/codecov-delta?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gammapy)
> `Δ = absolute <relative> (impact)`, `ø = not affected`, `? = missing data`
> Powered by [Codecov](https://codecov.io/gh/gammapy/gammapy/pull/3837?src=pr&el=footer&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gammapy). Last update [dd35f7a...f0fc779](https://codecov.io/gh/gammapy/gammapy/pull/3837?src=pr&el=lastupdated&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gammapy). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gammapy).

maxnoe: > Should we raise an error when scalar data is interpolated with method="linear" or just silently use "nearest"?

I'd say raise an error. Falling back to nearest is probably not the right choice every time, so it shouldn't be done silently

Patch:
diff --git a/gammapy/data/event_list.py b/gammapy/data/event_list.py
index 385a378ea..33c01a0ef 100644
--- a/gammapy/data/event_list.py
+++ b/gammapy/data/event_list.py
@@ -769,6 +769,34 @@ class EventList:
         mask &= offset < offset_band[1]
         return self.select_row_subset(mask)
 
+    def select_rad_max(self, rad_max, position=None):
+        """Select energy dependent offset
+
+        Parameters
+        ----------
+        rad_max : `~gamapy.irf.RadMax2D`
+            Rad max definition
+        position : `~astropy.coordinates.SkyCoord`
+            Center position. By default the pointing position is used.
+
+        Returns
+        -------
+        event_list : `EventList`
+            Copy of event list with selection applied.
+        """
+        if position is None:
+            position = self.pointing_radec
+
+        offset = position.separation(self.pointing_radec)
+        separation = position.separation(self.radec)
+
+        rad_max_for_events = rad_max.evaluate(
+            method="nearest", energy=self.energy, offset=offset
+        )
+
+        selected = separation <= rad_max_for_events
+        return self.select_row_subset(selected)
+
     @property
     def is_pointed_observation(self):
         """Whether observation is pointed"""
diff --git a/gammapy/irf/core.py b/gammapy/irf/core.py
index 5af7d6db2..0ec3a92ad 100644
--- a/gammapy/irf/core.py
+++ b/gammapy/irf/core.py
@@ -258,6 +258,7 @@ class IRF(metaclass=abc.ABCMeta):
         """
         # TODO: change to coord dict?
         non_valid_axis = set(kwargs).difference(self.axes.names)
+
         if non_valid_axis:
             raise ValueError(
                 f"Not a valid coordinate axis {non_valid_axis}"
@@ -270,6 +271,7 @@ class IRF(metaclass=abc.ABCMeta):
             coord = kwargs.get(key, value)
             if coord is not None:
                 coords_default[key] = u.Quantity(coord, copy=False)
+
         data = self._interpolate(coords_default.values(), method=method)
 
         if self.interp_kwargs["fill_value"] is not None:
diff --git a/gammapy/irf/rad_max.py b/gammapy/irf/rad_max.py
index cbe1a7344..32943ff87 100644
--- a/gammapy/irf/rad_max.py
+++ b/gammapy/irf/rad_max.py
@@ -32,11 +32,9 @@ class RadMax2D(IRF):
     required_axes = ["energy", "offset"]
     default_unit = u.deg
 
-
     @classmethod
     def from_irf(cls, irf):
-        '''
-        Create a RadMax2D instance from another IRF component.
+        """Create a RadMax2D instance from another IRF component.
 
         This reads the RAD_MAX metadata keyword from the irf and creates
         a RadMax2D with a single bin in energy and offset using the
@@ -57,21 +55,25 @@ class RadMax2D(IRF):
         -----
         This assumes the true energy axis limits are also valid for the
         reco energy limits.
-        '''
+        """
         if not irf.is_pointlike:
-            raise ValueError('RadMax2D.from_irf is only valid for point-like irfs')
+            raise ValueError("RadMax2D.from_irf requires a point-like irf")
 
-        if 'RAD_MAX' not in irf.meta:
-            raise ValueError('irf does not contain RAD_MAX keyword')
+        if "RAD_MAX" not in irf.meta:
+            raise ValueError("Irf does not contain RAD_MAX keyword")
 
         rad_max_value = irf.meta["RAD_MAX"]
         if not isinstance(rad_max_value, float):
-            raise ValueError('RAD_MAX must be a float')
+            raise ValueError(
+                f"RAD_MAX must be a float, got '{type(rad_max_value)}' instead"
+            )
 
         energy_axis = irf.axes["energy_true"].copy(name="energy").squash()
         offset_axis = irf.axes["offset"].squash()
 
         return cls(
-            data=u.Quantity([[rad_max_value]], u.deg),
+            data=rad_max_value,
             axes=[energy_axis, offset_axis],
+            unit="deg",
+            interp_kwargs={"method": "nearest", "fill_value": None},
         )
diff --git a/gammapy/makers/utils.py b/gammapy/makers/utils.py
index 3c2d8ed1c..975738506 100644
--- a/gammapy/makers/utils.py
+++ b/gammapy/makers/utils.py
@@ -428,34 +428,6 @@ def make_theta_squared_table(
     return table
 
 
-def get_rad_max_vs_energy(rad_max, pointing, geom):
-    """Obtain the values of `RAD_MAX` at a given offset and for an array of
-    estimated energy values (in the geom energy axis).
-
-    Parameters
-    ----------
-    rad_max : `~gammapy.irf.RadMax2D`
-        the RAD_MAX_2D table IRF
-    geom : `~gammapy.maps.Geom`
-        the map geom to be used
-    pointing : `~astropy.coordinates.SkyCoord`
-        pointing direction
-
-    Returns
-    -------
-    array : `~astropy.units.Quantity`
-        Values of the `RAD_MAX` corresponding to each estimated energy bin center.
-    """
-    on_center = geom.center_skydir
-    offset = on_center.separation(pointing)
-
-    rad_max_vals = rad_max.evaluate(
-        offset=offset, energy=geom.axes["energy"].center
-    )
-
-    return rad_max_vals
-
-
 def make_counts_rad_max(geom, rad_max, events):
     """Extract the counts using for the ON region size the values in the
     `RAD_MAX_2D` table.
@@ -474,31 +446,15 @@ def make_counts_rad_max(geom, rad_max, events):
     counts : `~gammapy.maps.RegionNDMap`
         Counts vs estimated energy extracted from the ON region.
     """
-    selected_events = apply_rad_max(events, rad_max, geom.region.center)
+    selected_events = events.select_rad_max(
+        rad_max=rad_max, position=geom.region.center
+    )
 
-    counts_geom = RegionGeom(geom.region, axes=[geom.axes['energy']])
-    counts = Map.from_geom(counts_geom)
+    counts = Map.from_geom(geom=geom)
     counts.fill_events(selected_events)
-
     return counts
 
 
-def apply_rad_max(events, rad_max, position):
-    '''Apply the RAD_MAX cut to the event list for given region'''
-    offset = position.separation(events.pointing_radec)
-    separation = position.separation(events.radec)
-
-    if rad_max.data.shape == (1, 1):
-        rad_max_for_events = rad_max.quantity[0, 0]
-    else:
-        rad_max_for_events = rad_max.evaluate(
-            method="nearest", energy=events.energy, offset=offset
-        )
-
-    selected = separation <= rad_max_for_events
-    return events.select_row_subset(selected)
-
-
 def are_regions_overlapping_rad_max(regions, rad_max, offset, e_min, e_max):
     """
     Calculate pair-wise separations between all regions and compare with rad_max
@@ -509,20 +465,15 @@ def are_regions_overlapping_rad_max(regions, rad_max, offset, e_min, e_max):
         for a, b in combinations(regions, 2)
     ])
 
-
-    # evaluate fails with a single bin somewhere trying to interpolate
-    if rad_max.data.shape == (1, 1):
-        rad_max_at_offset = rad_max.quantity[0, 0]
-    else:
-        rad_max_at_offset = rad_max.evaluate(offset=offset)
-        # do not check bins outside of energy range
-        edges_min = rad_max.axes['energy'].edges_min
-        edges_max = rad_max.axes['energy'].edges_max
-        # to be sure all possible values are included, we check
-        # for the *upper* energy bin to be larger than e_min and the *lower* edge
-        # to be larger than e_max
-        mask = (edges_max >= e_min) & (edges_min <= e_max)
-        rad_max_at_offset = rad_max_at_offset[mask]
+    rad_max_at_offset = rad_max.evaluate(offset=offset)
+    # do not check bins outside of energy range
+    edges_min = rad_max.axes['energy'].edges_min
+    edges_max = rad_max.axes['energy'].edges_max
+    # to be sure all possible values are included, we check
+    # for the *upper* energy bin to be larger than e_min and the *lower* edge
+    # to be larger than e_max
+    mask = (edges_max >= e_min) & (edges_min <= e_max)
+    rad_max_at_offset = rad_max_at_offset[mask]
 
     return np.any(separations[np.newaxis, :] < (2 * rad_max_at_offset))
 
@@ -564,7 +515,7 @@ def make_counts_off_rad_max(
     )
 
     if len(off_regions) == 0:
-        log.warn("RegionsFinder returned no regions")
+        log.warning("RegionsFinder returned no regions")
         # counts_off=None, acceptance_off=0
         return None, RegionNDMap.from_geom(on_geom, data=0)
 
@@ -592,7 +543,9 @@ def make_counts_off_rad_max(
     )
 
     for off_region in off_regions:
-        selected_events = apply_rad_max(events, rad_max, off_region.center)
+        selected_events = events.select_rad_max(
+            rad_max=rad_max, position=off_region.center
+        )
         counts_off.fill_events(selected_events)
 
     return counts_off, acceptance_off
diff --git a/gammapy/utils/interpolation.py b/gammapy/utils/interpolation.py
index 0475d139e..9c0450e88 100644
--- a/gammapy/utils/interpolation.py
+++ b/gammapy/utils/interpolation.py
@@ -3,6 +3,7 @@
 import numpy as np
 import scipy.interpolate
 from astropy import units as u
+from itertools import compress
 
 __all__ = [
     "interpolate_profile",
@@ -54,22 +55,29 @@ class ScaledRegularGridInterpolator:
 
         self.scale_points = [interpolation_scale(scale) for scale in points_scale]
         self.scale = interpolation_scale(values_scale)
-        self._include_dim = [len(p) > 1 for p in points]
-
-        points_scaled = tuple(
-            [
-                scale(p)
-                for p, scale, _ in zip(points, self.scale_points, self._include_dim)
-                if _
-            ]
-        )
-        values_scaled = self.scale(values).squeeze()
         self.axis = axis
 
+        self._include_dimensions = [len(p) > 1 for p in points]
+
+        values_scaled = self.scale(values)
+        points_scaled = self._scale_points(points=points)
+
         if extrapolate:
             kwargs.setdefault("bounds_error", False)
             kwargs.setdefault("fill_value", None)
 
+        method = kwargs.get("method", None)
+
+        if not np.any(self._include_dimensions):
+            if method != "nearest":
+                raise ValueError(
+                    "Interpolating scalar values requires using "
+                    "method='nearest' explicitely."
+                )
+
+        if np.any(self._include_dimensions):
+            values_scaled = np.squeeze(values_scaled)
+
         if axis is None:
             self._interpolate = scipy.interpolate.RegularGridInterpolator(
                 points=points_scaled, values=values_scaled, **kwargs
@@ -79,6 +87,14 @@ class ScaledRegularGridInterpolator:
                 points_scaled[0], values_scaled, axis=axis
             )
 
+    def _scale_points(self, points):
+        points_scaled = [scale(p) for p, scale in zip(points, self.scale_points)]
+
+        if np.any(self._include_dimensions):
+            points_scaled = compress(points_scaled, self._include_dimensions)
+
+        return tuple(points_scaled)
+
     def __call__(self, points, method="linear", clip=True, **kwargs):
         """Interpolate data points.
 
@@ -92,14 +108,7 @@ class ScaledRegularGridInterpolator:
         clip : bool
             Clip values at zero after interpolation.
         """
-
-        points = tuple(
-            [
-                scale(p)
-                for scale, p, _ in zip(self.scale_points, points, self._include_dim)
-                if _
-            ]
-        )
+        points = self._scale_points(points=points)
 
         if self.axis is None:
             points = np.broadcast_arrays(*points)


Test patch:
diff --git a/gammapy/irf/tests/test_rad_max.py b/gammapy/irf/tests/test_rad_max.py
index 0fdec6cd0..345d5ceee 100644
--- a/gammapy/irf/tests/test_rad_max.py
+++ b/gammapy/irf/tests/test_rad_max.py
@@ -1,12 +1,13 @@
 import numpy as np
+from numpy.testing import assert_allclose
 import astropy.units as u
 from gammapy.maps import MapAxis
+from gammapy.irf import RadMax2D, EffectiveAreaTable2D
+
 import pytest
 
 
 def test_rad_max_roundtrip(tmp_path):
-    from gammapy.irf import RadMax2D
-
     n_energy = 10
     energy_axis = MapAxis.from_energy_bounds(
         50 * u.GeV, 100 * u.TeV, n_energy, name="energy"
@@ -30,13 +31,11 @@ def test_rad_max_roundtrip(tmp_path):
     rad_max_2d.write(tmp_path / "rad_max.fits")
     rad_max_read = RadMax2D.read(tmp_path / "rad_max.fits")
 
-    assert np.all(rad_max_read.data.data == rad_max)
-    assert np.all(rad_max_read.data.data == rad_max_read.data.data)
+    assert np.all(rad_max_read.data == rad_max)
+    assert np.all(rad_max_read.data == rad_max_read.data)
 
 
 def test_rad_max_from_irf():
-    from gammapy.irf import RadMax2D, EffectiveAreaTable2D
-
     e_bins = 3
     o_bins = 2
     energy_axis = MapAxis.from_energy_bounds(1 * u.TeV, 10 * u.TeV, nbin=e_bins, name='energy_true')
@@ -71,3 +70,22 @@ def test_rad_max_from_irf():
     assert rad_max.axes['offset'].edges[1] == aeff.axes['offset'].edges[-1]
     assert rad_max.quantity.shape == (1, 1)
     assert rad_max.quantity[0, 0] == 0.2 * u.deg
+
+
+def test_rad_max_single_bin():
+    energy_axis = MapAxis.from_energy_bounds(0.01, 100, 1, unit="TeV")
+    offset_axis = MapAxis.from_bounds(0., 5, 1, unit="deg", name="offset", )
+
+    rad_max = RadMax2D(
+        data=[[0.1]] * u.deg,
+        axes=[energy_axis, offset_axis],
+        interp_kwargs={"method": "nearest", "fill_value": None}
+    )
+
+    value = rad_max.evaluate(energy=1 * u.TeV, offset=1 * u.deg)
+    assert_allclose(value, 0.1 * u.deg)
+
+    value = rad_max.evaluate(energy=[1, 2, 3] * u.TeV, offset=[[1]] * u.deg)
+    assert value.shape == (1, 3)
+    assert_allclose(value, 0.1 * u.deg)
+